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/validation/params_resonator_gmshfem_comparison.lua b/share/validation/params_resonator_gmshfem_comparison.lua
new file mode 100644
index 0000000000000000000000000000000000000000..9721a1f6ca9c80a88561a3d9d6dc69da33d7d4c9
--- /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..fc0c04645fca7f0c071c4f07bbab9ca0ad0efd3e 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};
-MeshSize{:} = 0.05;
+Physical Volume(1000) = {1};
+Physical Surface(2000) = {1,2,3,4,5,6};
+
+MeshSize{:} = 0.2;
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..21db7ea7bd0c5a6e19b6bce267ac0a986ead9c56 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,91 @@ 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.Fy*qpFval.Fy);
+
+ 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::BGreenfg;
+ std::cout << " Computed norm (FEM): " << std::sqrt(normsq_gmshfem) << ", off by " << norm_gmshfem_error << " %" << std::endl;
+ std::cout << sgr::reset;
+
+ std::cout << " DG vs. FEM: " << std::sqrt(errorsq_dg_vs_gmshfem) << std::endl;
+ std::cout << " DG vs. analytic: " << std::sqrt(errorsq_dg_vs_ana) << std::endl;
+ std::cout << sgr::BGreenfg;
+ 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 6a3e29a2ff0f5fcbf59fb4086a9d08e086da46ec..3fb6cdf03e4e89711266e2ae2e6d5763e64db79c 100644
--- a/src/maxwell/maxwell_solver.cpp
+++ b/src/maxwell/maxwell_solver.cpp
@@ -191,7 +191,7 @@ 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
@@ -220,15 +220,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, "");
@@ -236,29 +243,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 << ", avg time: " << time/cycle_print_rate << " s" << 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 /* USE_MPI */
- std::cout << showcursor << std::endl;
+#endif
+ std::cout << std::endl;
+
+ mpl.call_finalization_callback();
+
+#ifdef USE_MPI
+ if (proc_rank == 0)
+#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
+ }
}
/****************************************************************************/
@@ -315,6 +337,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 */