Allow to export data for comparison with other GMSH-based solvers.

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; }
 };

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
 

share/validation/params_resonator_gmshfem_comparison.lua

+-- 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

share/validation/resonator.geo

 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;

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
 {

src/libgmshdg/physical_element.cpp

 #include <iostream>
 #include <cassert>
+#include <fstream>
 
 #include "libgmshdg/physical_element.h"
 #include "libgmshdg/gmsh_io.h"

src/maxwell/maxwell_postpro.cpp

 
+#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
 

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 */