diff --git a/include/libgmshdg/gmsh_io.h b/include/libgmshdg/gmsh_io.h
index 81d959b50f647796988a042dad898b3a5a9f8453..c642d40f4dcd5dc6480d2f48407426ca516bcf8d 100644
--- a/include/libgmshdg/gmsh_io.h
+++ b/include/libgmshdg/gmsh_io.h
@@ -183,11 +183,8 @@ public:
model(const char *, int, int);
~model();
- void build(void);
- void build(double);
+ void build();
- void generate_mesh(void);
- void generate_mesh(double);
void partition_mesh(int);
void populate_from_gmsh(void);
diff --git a/include/libgmshdg/param_loader.h b/include/libgmshdg/param_loader.h
index 7ca629455d77b287ecebedd1cc1b031cdb6dc507..4390e96258208fdac076bb74cc2b9bc6ad4b2c49 100644
--- a/include/libgmshdg/param_loader.h
+++ b/include/libgmshdg/param_loader.h
@@ -66,6 +66,8 @@ public:
void call_timestep_callback(size_t);
void call_initialization_callback();
+ void call_finalization_callback();
+
sol::state& lua_state(void) { return lua; }
};
diff --git a/include/maxwell/maxwell_postpro.h b/include/maxwell/maxwell_postpro.h
index fe345ade591e8aedd4380882c93e0db6768a4020..943bf50aedc0dc1c22e84ec334858d54653e0fea 100644
--- a/include/maxwell/maxwell_postpro.h
+++ b/include/maxwell/maxwell_postpro.h
@@ -8,5 +8,7 @@ field_values eval_field(const field&, size_t, size_t, const vecxd&);
field_values compute_error(const model&, const solver_state&, const parameter_loader&);
field_values compute_energy(const model&, const solver_state&, const parameter_loader&);
+void compare_at_gauss_points(const model&, const solver_state&, const parameter_loader&);
+
} // namespace maxwell
diff --git a/share/examples/capacitor/capacitor.geo b/share/examples/capacitor/capacitor.geo
index 3a2306b212ef89647147b81b13ab669689d57781..b20a1d3c68b5d4940068f681f601063814b87ab2 100644
--- a/share/examples/capacitor/capacitor.geo
+++ b/share/examples/capacitor/capacitor.geo
@@ -11,4 +11,10 @@ Cylinder(3) = {0, 0, -d/2-t, 0, 0, t, R};
Sphere(4) = {0, 0, 0, 2*R, -Pi/2, Pi/2, 2*Pi};
Coherence;
+Physical Volume(1000) = {1, 3};
+Physical Volume(1001) = {2};
+Physical Volume(1002) = {4};
+
+Physical Surface(2000) = {1};
+
MeshSize{ PointsOf{ Volume{2}; } } = 0.0003;
diff --git a/share/examples/capacitor2/capacitor.geo b/share/examples/capacitor2/capacitor.geo
index 3a2306b212ef89647147b81b13ab669689d57781..b20a1d3c68b5d4940068f681f601063814b87ab2 100644
--- a/share/examples/capacitor2/capacitor.geo
+++ b/share/examples/capacitor2/capacitor.geo
@@ -11,4 +11,10 @@ Cylinder(3) = {0, 0, -d/2-t, 0, 0, t, R};
Sphere(4) = {0, 0, 0, 2*R, -Pi/2, Pi/2, 2*Pi};
Coherence;
+Physical Volume(1000) = {1, 3};
+Physical Volume(1001) = {2};
+Physical Volume(1002) = {4};
+
+Physical Surface(2000) = {1};
+
MeshSize{ PointsOf{ Volume{2}; } } = 0.0003;
diff --git a/share/examples/capacitor2/params_capacitor.lua b/share/examples/capacitor2/params_capacitor.lua
index c1bc00acb52754f87738f07d82ec3a2a78b86a51..6e7787a38ac3f97ebde149a1cc0421fa0cd150e8 100644
--- a/share/examples/capacitor2/params_capacitor.lua
+++ b/share/examples/capacitor2/params_capacitor.lua
@@ -1,11 +1,11 @@
sim.name = "capacitor" -- simulation name
sim.dt = 1e-14 -- timestep size
-sim.timesteps = 1001 -- num of iterations
-sim.gmsh_model = "capacitor.geo" -- gmsh model filename
+sim.timesteps = 101 -- num of iterations
+sim.gmsh_model = "capacitor.msh" -- gmsh model filename
sim.use_gpu = 0 -- 0: cpu, 1: gpu
-sim.approx_order = 1 -- approximation order
+sim.approx_order = 2 -- approximation order
sim.time_integrator = "leapfrog" -- time integration method
-postpro.silo_output_rate = 10 -- rate at which to write silo files
+postpro.silo_output_rate = 100 -- rate at which to write silo files
postpro.cycle_print_rate = 10 -- console print rate
local diel_epsr = 10 -- Permittivity of the capacitor dielectric
@@ -38,11 +38,11 @@ for i,v in ipairs(air) do
end
-- ** Boundary conditions **
-local absorbing_bcs = {1}
-for i,v in ipairs(absorbing_bcs) do
- bndconds[v] = {}
- bndconds[v].kind = "impedance"
-end
+--local absorbing_bcs = {1}
+--for i,v in ipairs(absorbing_bcs) do
+-- bndconds[v] = {}
+-- bndconds[v].kind = "impedance"
+--end
local R = 0.01 -- Capacitor disk radius (must match capacitor.geo)
local d = 0.001 -- Dielectric thickness (must match capacitor.geo)
@@ -103,6 +103,15 @@ end
-- Apply the current source to entity 2
sources[2] = current_source
+
+function on_timestepping_finished()
+ if ( do_IO ) then
+ print("\x1b[32m*** GMSH-FEM COMPARISON DATA ***\x1b[0m")
+ --compare_at_gauss_points() -- NOT MPI SAFE
+ end
+end
+
+
-- Dump charging current profile for debug
if ( do_IO ) then
local file = io.open("current_profile.dat", "w")
diff --git a/share/validation/params_resonator_gmshfem_comparison.lua b/share/validation/params_resonator_gmshfem_comparison.lua
new file mode 100644
index 0000000000000000000000000000000000000000..7d8a785d3003b3b55a0d3a1c10d153b2eaf8d0cc
--- /dev/null
+++ b/share/validation/params_resonator_gmshfem_comparison.lua
@@ -0,0 +1,97 @@
+-- Resonant cavity test. See Pozar, section 6.3 for theory.
+
+-----------------------
+--[[ Problem setup ]]--
+sim.name = "test_maxwell_resonator" -- simulation name
+sim.dt = 1e-11 -- timestep size
+sim.timesteps = 11 -- num of iterations
+sim.gmsh_model = "resonator.geo" -- gmsh model filename
+sim.use_gpu = 0 -- 0: cpu, 1: gpu
+sim.approx_order = 2 -- approximation order
+sim.time_integrator = "leapfrog"
+postpro.silo_output_rate = 100
+postpro.cycle_print_rate = 10 -- console print rate
+
+postpro["E"].silo_mode = "zonal"
+postpro["H"].silo_mode = "none"
+postpro["J"].silo_mode = "none"
+
+local epsr = 1
+local mur = 1
+
+materials[1] = {}
+materials[1].epsilon = epsr
+materials[1].mu = mur
+materials[1].sigma = 0
+
+function electric_initial_condition(x, y, z)
+ local Ex = 0
+ local Ey = math.sin(math.pi*x) * math.sin(math.pi*z)
+ local Ez = 0
+ return Ex, Ey, Ez
+end
+
+--------------------------
+--[[ Validation stuff ]]--
+debug = {}
+
+-- determine if we should do I/O
+local do_IO = (not parallel) or (parallel and parallel.comm_rank == 0)
+
+local c0 = 1/math.sqrt(const.eps0*const.mu0)
+-- Mode
+local m = 1 -- along x
+local n = 0 -- along y
+local l = 1 -- along z
+-- Cavity dimensions (must match sim.gmsh_model)
+local a = 1 -- along x
+local b = 0.1 -- along y
+local d = 1 -- along z
+
+local u = m*math.pi/a
+local v = n*math.pi/b
+local w = l*math.pi/d
+-- Compute resonant frequency
+local omega = c0*math.sqrt(u*u + v*v + w*w)/math.sqrt(epsr*mur)
+local resonance_f = omega/(2*math.pi)
+local resonance_MHz = resonance_f/1e6
+local cycle_timesteps = 1/(resonance_f*sim.dt)
+-- Compute impedance
+local eps = materials[1].epsilon * const.eps0
+local mu = materials[1].mu * const.mu0
+local k_sq = (omega*mu)*(omega*eps)
+local kc_sq = u*u + v*v
+local beta = math.sqrt(k_sq - kc_sq)
+local Zte = omega*mu/beta
+local We = eps*a*b*d/16
+
+if ( do_IO ) then
+ print("\x1b[1mANALYTICAL SOLUTION DATA:")
+ print(" Resonance frequency: " .. resonance_MHz .. " Mhz")
+ print(" Cavity impedance: " .. Zte .. " Ohm")
+ print(" Timesteps for 1 cycle: " .. cycle_timesteps)
+ print(" Expected energy " .. 2*We .. " J \x1b[0m")
+end
+
+function ansol(tag, x, y, z, t)
+ local Ex = 0.0
+ local Ey = math.cos(omega*t)*math.sin(math.pi*x)*math.sin(math.pi*z)
+ local Ez = 0.0
+ local Hx = math.sin(omega*t)*math.sin(math.pi*x)*math.cos(math.pi*z)/Zte
+ local Hy = 0.0
+ local Hz = -math.sin(omega*t)*math.cos(math.pi*x)*math.sin(math.pi*z)/Zte
+
+ return Ex, Ey, Ez, Hx, Hy, Hz
+end
+
+debug.analytical_solution = ansol
+--debug.dump_cell_ranks = true
+
+function on_timestepping_finished()
+ err = compute_error()
+ if ( do_IO ) then
+ print("\x1b[32m*** GMSH-FEM COMPARISON DATA ***\x1b[0m")
+ print("Error on Ey: " .. err.Ey)
+ compare_at_gauss_points() -- NOT MPI SAFE
+ end
+end
diff --git a/share/validation/resonator.geo b/share/validation/resonator.geo
index 3c3f82936d969c3b7a9a80eeb9c048c2b67cbfac..8c80fad0745ead4e1c2e54b635226b093822b3b7 100644
--- a/share/validation/resonator.geo
+++ b/share/validation/resonator.geo
@@ -1,4 +1,7 @@
SetFactory("OpenCASCADE");
-Box(1) = {0, 0, 0, 1, 0.1, 1};
+Box(1) = {0, 0, 0, 1, 1, 1};
+
+Physical Volume(1000) = {1};
+Physical Surface(2000) = {1,2,3,4,5,6};
MeshSize{:} = 0.05;
diff --git a/src/libgmshdg/gmsh_io.cpp b/src/libgmshdg/gmsh_io.cpp
index 23f5fb7225fb0b35689990da4999362e00b1eb97..056aee4b9ed20cbbf31e1c81a18338462ee6df88 100644
--- a/src/libgmshdg/gmsh_io.cpp
+++ b/src/libgmshdg/gmsh_io.cpp
@@ -10,6 +10,9 @@
#include "common/mpi_helpers.h"
#endif /* USE_MPI */
+#include "sgr.hpp"
+using namespace sgr;
+
std::string
quadrature_name(int order)
{
@@ -252,13 +255,14 @@ model::import_gmsh_entities_rank0(void)
flux_base_world += e.num_fluxes();
index_base_world += e.num_cells();
- /* Send entity to the appropriate process */
assert(tag_partition > 0);
int rank = tag_partition - 1;
- e.mpi_send(rank, MPI_COMM_WORLD);
#ifdef ENABLE_DEBUG_PRINT
- std::cout << "Sending entity " << tag << " to " << rank << " " << e.num_dofs() << std::endl;
+ std::cout << "Sending entity " << tag << " to " << rank << ", ";
+ std::cout << e.num_dofs() << " DoFs" << std::endl;
#endif /* ENABLE_DEBUG_PRINT */
+ /* Send entity to the appropriate process */
+ e.mpi_send(rank, MPI_COMM_WORLD);
/* Save locally the entity */
remote_entities[tag_partition].push_back( std::move(e) );
all_entities_tags.push_back(tag);
@@ -341,13 +345,14 @@ model::import_gmsh_entities_rankN(void)
for (auto& rtag : partition_map.at(mp))
{
entity e;
- e.mpi_recv(0, MPI_COMM_WORLD);
-
#ifdef ENABLE_DEBUG_PRINT
- std::cout << "Receiving " << rtag << ", rank " << mp-1 << " ";
- std::cout << e.num_dofs() << std::endl;
+ std::cout << "Receiving " << rtag << ", rank " << mp-1 << ", ";
#endif /* ENABLE_DEBUG_PRINT */
(void) rtag;
+ e.mpi_recv(0, MPI_COMM_WORLD);
+#ifdef ENABLE_DEBUG_PRINT
+ std::cout << e.num_dofs() << " DoFs" << std::endl;
+#endif /* ENABLE_DEBUG_PRINT */
for (size_t i = 0; i < e.num_cells(); i++)
{
@@ -405,17 +410,38 @@ model::get_facekey_tag_pairs(void)
void
model::map_boundaries(void)
{
-#ifdef USE_MPI
- //ASSERT_MPI_RANK_0;
-#endif /* USE_MPI */
std::vector<bfk_t> bfk = get_facekey_tag_pairs();
for (size_t i = 1; i < bfk.size(); i++)
{
if ( bfk[i-1].first == bfk[i].first )
{
- std::cout << bfk[i-1].second << " " << bfk[i].second << std::endl;
+ std::cout << redfg << Bon << "WARNING: " << nofg << __FILE__;
+ std::cout << "(" << __LINE__ << "): Face identified by ";
+ std::cout << bfk[i].first << " was found on interfaces ";
+ std::cout << bfk[i-1].second;
+#ifdef USE_MPI
+ if ( is_interprocess_boundary(bfk[i-1].second) )
+ std::cout << " (IP)";
+#endif
+ std::cout << " and " << bfk[i].second;
+#ifdef USE_MPI
+ if ( is_interprocess_boundary(bfk[i].second) )
+ std::cout << " (IP)";
+ int rank;
+ MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+ std::cout << " [MPI rank " << rank << "]";
+#endif
+ std::cout << ". " << reset << std::endl;
+ /*
+#ifdef USE_MPI
+ int rank;
+ MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+ if (rank == 0)
+#endif
+ gmsh::write("crash.msh");
throw 42;
+ */
}
}
@@ -668,7 +694,7 @@ void
model::populate_from_gmsh_rank0()
{
#ifdef USE_MPI
- ASSERT_MPI_RANK_0;
+ ASSERT_MPI_RANK_0;
#endif /* USE_MPI */
gmm::setOrder( geometric_order );
@@ -764,36 +790,10 @@ model::populate_from_gmsh()
void
model::build()
{
- generate_mesh();
- make_boundary_to_faces_map();
+ //make_boundary_to_faces_map();
populate_from_gmsh_rank0();
}
-void
-model::build(double sf)
-{
- generate_mesh(sf);
- make_boundary_to_faces_map();
- populate_from_gmsh_rank0();
-}
-
-void
-model::generate_mesh()
-{
- gmm::generate( DIMENSION(3) );
-}
-
-void
-model::generate_mesh(double sf)
-{
- gmm::generate( DIMENSION(3) );
-
- std::vector<double> tr = {sf, 0, 0, 0,
- 0, sf, 0, 0,
- 0, 0, sf, 0 };
- gmm::affineTransform(tr);
-}
-
#ifdef USE_MPI
void
model::partition_mesh(int parts)
diff --git a/src/libgmshdg/param_loader.cpp b/src/libgmshdg/param_loader.cpp
index 238da20ec851af53dad694de80ccb29394348196..831aec733bbb5265218b0c7b6df0ff5e7f231102 100644
--- a/src/libgmshdg/param_loader.cpp
+++ b/src/libgmshdg/param_loader.cpp
@@ -247,6 +247,21 @@ parameter_loader_base::call_initialization_callback()
}
}
+void
+parameter_loader_base::call_finalization_callback()
+{
+ if ( lua["on_timestepping_finished"].valid() )
+ {
+ try {
+ lua["on_timestepping_finished"]();
+ }
+ catch (...)
+ {
+ std::cout << "An error was thrown calling on_timestepping_finished()" << std::endl;
+ }
+ }
+}
+
bool
parameter_loader_base::source_has_changed_state(void) const
{
diff --git a/src/libgmshdg/physical_element.cpp b/src/libgmshdg/physical_element.cpp
index c9190fae938339f837f4051c0981179588681dc2..655c49cffdebad422c532a5e23ec22edf2f5b85d 100644
--- a/src/libgmshdg/physical_element.cpp
+++ b/src/libgmshdg/physical_element.cpp
@@ -1,5 +1,6 @@
#include <iostream>
#include <cassert>
+#include <fstream>
#include "libgmshdg/physical_element.h"
#include "libgmshdg/gmsh_io.h"
diff --git a/src/maxwell/maxwell_postpro.cpp b/src/maxwell/maxwell_postpro.cpp
index 2d162f6c88c5893fb7b0b9ecba87a9779bc9a5a7..7dc612f40e3a8717f46daf5592892b84a1925aec 100644
--- a/src/maxwell/maxwell_postpro.cpp
+++ b/src/maxwell/maxwell_postpro.cpp
@@ -1,4 +1,7 @@
+#include <fstream>
+#include <fcntl.h>
+
#ifdef USE_MPI
#include "common/mpi_helpers.h"
#endif /* USE_MPI */
@@ -37,12 +40,10 @@ compute_error(const model& mod, const solver_state& state, const parameter_loade
const auto& re = e.cell_refelem(pe);
const auto pqps = pe.integration_points();
- const auto dets = pe.determinants();
const auto ofs = e.cell_model_dof_offset(iT);
const auto cbs = re.num_basis_functions();
- assert(pqps.size() == dets.size());
for (size_t iQp = 0; iQp < pqps.size(); iQp++)
{
const auto& pqp = pqps[iQp];
@@ -78,12 +79,10 @@ compute_energy(const model& mod, const solver_state& state, const parameter_load
const auto& re = e.cell_refelem(pe);
const auto pqps = pe.integration_points();
- const auto dets = pe.determinants();
const auto ofs = e.cell_model_dof_offset(iT);
const auto cbs = re.num_basis_functions();
- assert(pqps.size() == dets.size());
const auto bar = pe.barycenter();
auto eps = mpl.epsilon( e.material_tag(), bar );
auto mu = mpl.mu( e.material_tag(), bar );
@@ -113,5 +112,92 @@ validate(const model&, maxwell::solver_state_gpu&,
{}
#endif /* ENABLE_GPU_SOLVER */
+
+void
+compare_at_gauss_points(const model& mod, const solver_state& state, const parameter_loader& mpl)
+{
+ struct qpFval { double Fx, Fy, Fz; };
+
+ size_t point_id = 0;
+ std::ofstream pmap_ofs("pointmap.txt");
+
+ double errorsq_dg_vs_gmshfem = 0.0;
+ double errorsq_dg_vs_ana = 0.0;
+ double errorsq_gmshfem_vs_ana = 0.0;
+ double normsq_dg = 0.0;
+ double normsq_gmshfem = 0.0;
+ double normsq_ana = 0.0;
+ for (auto& e : mod)
+ {
+ std::stringstream ss;
+ ss << "e_" << e.gmsh_dim() << "_" << e.gmsh_tag() << ".dat";
+
+ std::cout << "Loading " << ss.str() << std::endl;
+ int fd = open(ss.str().c_str(), O_RDONLY);
+ if (fd < 0)
+ {
+ std::cout << "Can't open '" << ss.str() << "'" << std::endl;
+ return;
+ }
+ size_t num_vals;
+ read(fd, &num_vals, sizeof(size_t));
+ std::cout << "Reading " << num_vals << " values" << std::endl;
+ std::vector<qpFval> qpFvals(num_vals);
+ read(fd, qpFvals.data(), num_vals*sizeof(qpFval));
+ close(fd);
+
+ std::cout << "Num of cells: " << e.num_cells() << std::endl;
+ for (size_t iT = 0; iT < e.num_cells(); iT++)
+ {
+ const auto& pe = e.cell(iT);
+ const auto& re = e.cell_refelem(pe);
+ const auto ofs = e.cell_model_dof_offset(iT);
+ const auto cbs = re.num_basis_functions();
+ const auto pqps = pe.integration_points();
+
+ const auto gmsh_pos = e.cell_local_index_by_gmsh(iT);
+ const size_t num_qps = pqps.size();
+ const size_t qps_pos_base = num_qps*gmsh_pos;
+ for (size_t iQp = 0; iQp < num_qps; iQp++)
+ {
+ const auto& pqp = pqps[iQp];
+ const vecxd phi = re.basis_functions(iQp);
+ field_values Fh = eval_field(state.emf_curr, ofs, cbs, phi);
+ auto& qpFval = qpFvals.at(qps_pos_base+iQp);
+ auto Exd = Fh.Ex - qpFval.Fx;
+ auto Eyd = Fh.Ey - qpFval.Fy;
+ auto Ezd = Fh.Ez - qpFval.Fz;
+ errorsq_dg_vs_gmshfem += pqp.weight() * (Exd*Exd + Eyd*Eyd + Ezd*Ezd);
+ normsq_dg += pqp.weight() * (Fh.Ex*Fh.Ex + Fh.Ey*Fh.Ey + Fh.Ez*Fh.Ez);
+ normsq_gmshfem += pqp.weight() * (qpFval.Fx*qpFval.Fx + qpFval.Fy*qpFval.Fy + qpFval.Fz*qpFval.Fz);
+
+ auto p = pqp.point();
+ auto Ey_ana = std::sin(M_PI*p.x())*std::sin(M_PI*p.z())*std::cos(2.*M_PI*211.98528005809e6*state.curr_time);
+ normsq_ana += pqp.weight() * Ey_ana * Ey_ana;
+ errorsq_dg_vs_ana += pqp.weight() * ( (Fh.Ey - Ey_ana)*(Fh.Ey - Ey_ana) );
+ errorsq_gmshfem_vs_ana += pqp.weight() * ( (qpFval.Fy - Ey_ana)*(qpFval.Fy - Ey_ana) );
+ pmap_ofs << pe.element_tag() << " " << qps_pos_base+iQp << " " << p.x() << " " << p.y() << " " << p.z() << " " << point_id++ << std::endl;
+ }
+ }
+ }
+
+ double norm_dg_error = 100.*(std::sqrt(normsq_dg) - std::sqrt(normsq_ana))/std::sqrt(normsq_ana);
+ double norm_gmshfem_error = 100.*(std::sqrt(normsq_gmshfem) - std::sqrt(normsq_ana))/std::sqrt(normsq_ana);
+ std::cout << sgr::Bredfg << "Errors and norms:" << sgr::reset << std::endl;
+ std::cout << " Analytic norm: " << std::sqrt(normsq_ana) << std::endl;
+ std::cout << " Computed norm (DG): " << std::sqrt(normsq_dg) << ", off by " << norm_dg_error << " %" << std::endl;
+ std::cout << sgr::BYellowfg;
+ std::cout << " Computed norm (FEM): " << std::sqrt(normsq_gmshfem) << ", off by " << norm_gmshfem_error << " %" << std::endl;
+ std::cout << sgr::reset;
+
+ auto dgvsfem_percent = 100.*std::sqrt(errorsq_dg_vs_gmshfem)/std::sqrt(normsq_gmshfem);
+ std::cout << " DG vs. FEM: " << std::sqrt(errorsq_dg_vs_gmshfem) << ", " << dgvsfem_percent << " %" << std::endl;
+ std::cout << " DG vs. analytic: " << std::sqrt(errorsq_dg_vs_ana) << std::endl;
+ std::cout << sgr::BYellowfg;
+ std::cout << " FEM vs. analytic: " << std::sqrt(errorsq_gmshfem_vs_ana) << std::endl;
+ std::cout << sgr::reset;
+}
+
+
} // namespace maxwell
diff --git a/src/maxwell/maxwell_solver.cpp b/src/maxwell/maxwell_solver.cpp
index 179afbcb2d57edc34abeb1336b0eb06e23b92d42..514d2b73d2694fae782141a464cab9777d66c371 100644
--- a/src/maxwell/maxwell_solver.cpp
+++ b/src/maxwell/maxwell_solver.cpp
@@ -1,3 +1,5 @@
+#include <regex>
+
#ifdef DONT_USE_STD_FILESYSTEM
#include <sys/stat.h>
#include <sys/types.h>
@@ -5,6 +7,9 @@
#include <filesystem>
#endif
+#include <sys/time.h>
+#include <sys/resource.h>
+
#ifdef USE_MPI
#pragma clang diagnostic push
#pragma clang diagnostic ignored "-Weverything"
@@ -191,13 +196,25 @@ void solver_mainloop(const model& mod, State& state, maxwell::parameter_loader&
if (proc_rank == 0) {
#endif /* USE_MPI */
std::cout << " BEGINNING SIMULATION" << std::endl;
- std::cout << "I will do " << num_timesteps << " of " << state.delta_t;
+ std::cout << "I will do " << num_timesteps << " timesteps of " << state.delta_t;
std::cout << "s each." << std::endl;
std::cout << "Time integrator: " << mpl.sim_timeIntegratorName() << std::endl;
#ifdef USE_MPI
}
#endif /* USE_MPI */
+#ifdef USE_MPI
+ size_t my_num_dofs = 6*mod.num_dofs();
+ size_t num_dofs;
+ MPI_Reduce(&my_num_dofs, &num_dofs, 1, MPI_UNSIGNED_LONG_LONG, MPI_SUM, 0, MPI_COMM_WORLD);
+ if (proc_rank == 0)
+ std::cout << "Model has " << num_dofs << " DoFs" << std::endl;
+#else /* USE_MPI */
+ size_t num_dofs = 6*mod.num_dofs();
+ std::cout << "Model has " << num_dofs << " DoFs" << std::endl;
+#endif /* USE_MPI */
+
+
mpl.call_initialization_callback();
#ifdef ENABLE_EXP_GEOMETRIC_DIODE
#ifdef USE_MPI
@@ -208,15 +225,22 @@ void solver_mainloop(const model& mod, State& state, maxwell::parameter_loader&
sol::state& lua = mpl.lua_state();
lua["relmat"] = relmat;
#endif /* ENABLE_EXP_GEOMETRIC_DIODE */
- timecounter tc;
- tc.tic();
+
+ double total_iteration_time = 0.;
+
prepare_sources(mod, state, mpl);
for(size_t i = 0; i < num_timesteps; i++)
{
+ timecounter tc;
+ tc.tic();
+
mpl.call_timestep_callback(i);
timestep(mod, state, mpl, ti);
do_sources(mod, state, mpl);
+ double time = tc.toc();
+ total_iteration_time += time;
+
if ( (silo_output_rate != 0) and (i%silo_output_rate == 0) )
export_fields_to_silo(mod, state, mpl, "");
@@ -224,29 +248,44 @@ void solver_mainloop(const model& mod, State& state, maxwell::parameter_loader&
if ( (cycle_print_rate != 0) and (i%cycle_print_rate == 0) )
{
- double time = tc.toc();
- double dofs_s_proc = (mod.num_dofs()*6*cycle_print_rate)/(time);
+ double dofs_s_proc = (mod.num_dofs()*6)/time;
#ifdef USE_MPI
double dofs_s_total;
MPI_Reduce(&dofs_s_proc, &dofs_s_total, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
if (proc_rank == 0)
{
- std::cout << hidecursor << clrline << "Cycle " << i << ": t = ";
+ std::cout << cr << clrline << "Cycle " << i << ": t = ";
std::cout << state.curr_time << " s" << ", DOFs/s: ";
- std::cout << dofs_s_total << cr << std::flush;
+ std::cout << dofs_s_total << std::flush;
}
#else /* USE_MPI */
- std::cout << hidecursor << clrline << "Cycle " << i << ": t = ";
+ std::cout << cr << clrline << "Cycle " << i << ": t = ";
std::cout << state.curr_time << " s" << ", DOFs/s: ";
- std::cout << dofs_s_proc << cr << std::flush;
+ std::cout << dofs_s_proc << ", time: " << time << " s" << std::flush;
#endif /* USE_MPI */
- tc.tic();
}
}
- #ifdef USE_MPI
+
+#ifdef USE_MPI
+ if (proc_rank == 0)
+#endif
+ std::cout << std::endl;
+
+ mpl.call_finalization_callback();
+
+#ifdef USE_MPI
if (proc_rank == 0)
- #endif /* USE_MPI */
- std::cout << showcursor << std::endl;
+#endif
+ {
+ std::cout << Bon << greenfg;
+ std::cout << "Total iteration time = " << total_iteration_time << " s" << std::endl;
+ std::cout << "Avg iteration time = " << total_iteration_time/num_timesteps << " s" << std::endl;
+#ifdef USE_MPI
+ std::cout << "Avg dofs/s = " << num_timesteps*(mod.num_dofs_world()*6)/total_iteration_time << reset << std::endl;
+#else
+ std::cout << "Avg dofs/s = " << num_timesteps*(mod.num_dofs()*6)/total_iteration_time << reset << std::endl;
+#endif
+ }
}
/****************************************************************************/
@@ -303,6 +342,11 @@ register_lua_usertypes_bystate(maxwell::parameter_loader& mpl,
};
lua["compute_error"] = compute_error_lambda;
+ auto compare_at_gauss_points_lambda = [&]() {
+ compare_at_gauss_points(mod, state, mpl);
+ };
+ lua["compare_at_gauss_points"] = compare_at_gauss_points_lambda;
+
#ifdef ENABLE_EXP_GEOMETRIC_DIODE
lua["integrate_E"] = integrate_electric_field;
#endif /*ENABLE_EXP_GEOMETRIC_DIODE */
@@ -323,6 +367,9 @@ register_lua_usertypes_bystate(maxwell::parameter_loader& mpl,
*/
int main(int argc, char *argv[])
{
+ struct rusage ru_start, ru_end;
+
+ getrusage(RUSAGE_SELF, &ru_start);
#ifdef USE_MPI
MPI_Init(&argc, &argv);
@@ -358,38 +405,48 @@ int main(int argc, char *argv[])
gmsh::initialize();
gmsh::option::setNumber("General.Terminal", 0);
gmsh::open( mpl.sim_gmshmodel() );
+ bool generate = std::regex_match(mpl.sim_gmshmodel(), std::regex(".*\\.geo$") );
+
+ if (generate)
+ {
+ std::cout << "Generating mesh..." << std::flush;
+ gmm::generate( DIMENSION(3) );
+ std::cout << "done" << std::endl;
+ }
+
+ auto [scaled, sf] = mpl.mesh_scalefactor();
+ if (scaled)
+ {
+ std::vector<double> tr = {sf, 0, 0, 0,
+ 0, sf, 0, 0,
+ 0, 0, sf, 0 };
+ gmm::affineTransform(tr);
+ }
+
#ifdef USE_MPI
}
#endif /* USE_MPI */
model mod( mpl.sim_geomorder(), mpl.sim_approxorder() );
- auto [scaled, sf] = mpl.mesh_scalefactor();
#ifdef USE_MPI
if (comm_rank == 0)
- {
- if (scaled)
- mod.generate_mesh(sf);
- else
- mod.generate_mesh();
-
+ {
if (comm_size > 1)
mod.partition_mesh(comm_size);
-
+ std::cout << "Distributing mesh partitions..." << std::flush;
mod.populate_from_gmsh();
+ std::cout << "done" << std::endl;
}
else
{
mod.populate_from_gmsh();
}
#else /* USE_MPI */
- if (scaled)
- mod.build(sf);
- else
- mod.build();
+ mod.populate_from_gmsh();
#endif /* USE_MPI */
-
+
if ( not mpl.validate_model_params(mod) )
@@ -438,5 +495,9 @@ int main(int argc, char *argv[])
MPI_Finalize();
#endif
+ getrusage(RUSAGE_SELF, &ru_end);
+
+ std::cout << "Max RSS: " << (ru_end.ru_maxrss - ru_start.ru_maxrss) << std::endl;
+
return 0;
}
diff --git a/tests/test_basics.cpp b/tests/test_basics.cpp
index ff3818958e5ae379ad76732ea13bff9b741ec0fe..5fcd85317c8bd07edfdd0dfc922cf73baea685bf 100644
--- a/tests/test_basics.cpp
+++ b/tests/test_basics.cpp
@@ -39,6 +39,7 @@ int
test_basics(int geometric_order, int approximation_order)
{
make_geometry(0, 0.1);
+ gmm::generate( DIMENSION(3) );
model mod(geometric_order, approximation_order);
mod.build();
@@ -163,6 +164,8 @@ int
test_normals(int geometric_order, int approximation_order)
{
make_geometry(0, 0.1);
+ gmm::generate( DIMENSION(3) );
+
model mod(geometric_order, approximation_order);
mod.build();
@@ -247,6 +250,8 @@ int
test_normals_2()
{
make_geometry(0, 0.1);
+ gmm::generate( DIMENSION(3) );
+
model mod(1, 1);
mod.build();
diff --git a/tests/test_differentiation.cpp b/tests/test_differentiation.cpp
index 00c7675645c839ec1d58cec1284b4a38df200d32..9c7d56b1aec8070ea1b768898c804ade99097a6f 100644
--- a/tests/test_differentiation.cpp
+++ b/tests/test_differentiation.cpp
@@ -69,6 +69,8 @@ int test_differentiation_convergence(int geometric_order, int approximation_orde
{
double h = sizes[i];
make_geometry(0,h);
+ gmm::generate( DIMENSION(3) );
+
model mod(geometric_order, approximation_order);
mod.build();
diff --git a/tests/test_differentiation_gpu.cpp b/tests/test_differentiation_gpu.cpp
index 607d631385eebdb6074b85c9130714113c134a94..5d60e3d2595ad0cc3b51bb969d7524697470642b 100644
--- a/tests/test_differentiation_gpu.cpp
+++ b/tests/test_differentiation_gpu.cpp
@@ -70,6 +70,8 @@ int test_differentiation_convergence(int geometric_order, int approximation_orde
{
double h = sizes[i];
make_geometry(0,h);
+ gmm::generate( DIMENSION(3) );
+
model mod(geometric_order, approximation_order);
mod.build();
diff --git a/tests/test_lifting.cpp b/tests/test_lifting.cpp
index 609e48552174ecd9344f99cd8c970ced46caa5a4..28f79089a670c4ab69cb6140bd3fb352ba290bb1 100644
--- a/tests/test_lifting.cpp
+++ b/tests/test_lifting.cpp
@@ -81,6 +81,8 @@ int test_lifting(int geometric_order, int approximation_order)
{
double h = sizes[sz];
make_geometry(0,h);
+ gmm::generate( DIMENSION(3) );
+
model mod(geometric_order, approximation_order);
mod.build();
diff --git a/tests/test_lifting_gpu.cpp b/tests/test_lifting_gpu.cpp
index cb82338d6e627079d6572308b737d59c55d9b97a..a88892151eaea4b448da36e9ba54b5460b5321f0 100644
--- a/tests/test_lifting_gpu.cpp
+++ b/tests/test_lifting_gpu.cpp
@@ -82,8 +82,9 @@ int test_lifting(int geometric_order, int approximation_order)
{
double h = sizes[sz];
make_geometry(0,h);
+ gmm::generate( DIMENSION(3) );
model mod(geometric_order, approximation_order);
- mod.build();
+ mod.populate_from_gmsh();
#ifdef WRITE_TEST_OUTPUTS
std::stringstream ss;
diff --git a/tests/test_mass.cpp b/tests/test_mass.cpp
index 72fcbb29b8a19345c12ca6e2412defc657b14653..7d7f90404c889c5c15c3535246a136537f896a83 100644
--- a/tests/test_mass.cpp
+++ b/tests/test_mass.cpp
@@ -39,6 +39,8 @@ int test_mass_convergence(int geometric_order, int approximation_order)
{
double h = sizes[i];
make_geometry(0,h);
+ gmm::generate( DIMENSION(3) );
+
model mod(geometric_order, approximation_order);
mod.build();