Better organized examples.

src/maxwell/maxwell_postpro.cpp

+
+#ifdef USE_MPI
+#include <mpi.h>
+#endif /* USE_MPI */
+
+#include "maxwell/maxwell_postpro.h"
+
+namespace maxwell {
+
+field_values
+eval_field(const field& f, size_t ofs, size_t cbs, const vecxd& phi)
+{
+    assert(ofs+cbs <= f.num_dofs);
+    field_values ret;
+    ret.Ex = f.Ex.segment(ofs, cbs).dot(phi);
+    ret.Ey = f.Ey.segment(ofs, cbs).dot(phi);
+    ret.Ez = f.Ez.segment(ofs, cbs).dot(phi);
+    ret.Hx = f.Hx.segment(ofs, cbs).dot(phi);
+    ret.Hy = f.Hy.segment(ofs, cbs).dot(phi);
+    ret.Hz = f.Hz.segment(ofs, cbs).dot(phi);
+    return ret;
+}
+
+field_values
+compute_error(const model& mod, const solver_state& state, const parameter_loader& mpl)
+{
+    field_values err;
+
+    if (not mpl.has_analytical_solution())
+        throw std::invalid_argument("analytical solution not defined in the Lua configuration");
+
+    for (auto& e : mod)
+    {
+        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 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];
+                const vecxd phi = re.basis_functions(iQp);
+
+                field_values Fh = eval_field(state.emf_curr, ofs, cbs, phi);
+
+                auto [E, H] = mpl.eval_analytical_solution(e.material_tag(),
+                    pqp.point(), state.curr_time);
+                field_values F(E,H);
+                field_values Fdiff = F - Fh;
+                err += Fdiff*Fdiff*pqp.weight();
+            }
+        }
+    }
+
+#ifdef USE_MPI
+    MPI_Allreduce(MPI_IN_PLACE, &err, 6, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+#endif /* USE_MPI */
+    return sqrt(err);
+}
+
+field_values
+compute_energy(const model& mod, const solver_state& state, const parameter_loader& mpl)
+{
+    field_values energy;
+
+    for (auto& e : mod)
+    {
+        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 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 );
+
+            field_values mat(eps, eps, eps, mu, mu, mu);
+
+            for (size_t iQp = 0; iQp < pqps.size(); iQp++)
+            {
+                const auto& pqp = pqps[iQp];
+                const vecxd phi = re.basis_functions(iQp);
+                field_values Fh = eval_field(state.emf_curr, ofs, cbs, phi);
+                energy += 0.5*mat*Fh*Fh*pqp.weight();
+            }
+        }
+    }
+
+#ifdef USE_MPI
+    MPI_Allreduce(MPI_IN_PLACE, &energy, 6, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+#endif /* USE_MPI */
+    return energy;
+}
+
+#ifdef ENABLE_GPU_SOLVER
+static void
+validate(const model&, maxwell::solver_state_gpu&,
+    maxwell::parameter_loader&, std::ofstream&)
+{}
+#endif /* ENABLE_GPU_SOLVER */
+
+} // namespace maxwell
\ No newline at end of file

src/maxwell/maxwell_solver.cpp

@@ -19,41 +19,12 @@
 #include "libgmshdg/param_loader.h"
 #include "maxwell/maxwell_interface.h"
 #include "maxwell/maxwell_common.h"
+#include "maxwell/maxwell_postpro.h"
 #include "timecounter.h"
 
 #include "sgr.hpp"
 using namespace sgr;
 
-#if 0
-static void
-make_geometry(int order, double mesh_h)
-{
-    gm::add("difftest");
-    /*
-    std::vector<std::pair<int,int>> objects;
-
-    objects.push_back(
-        std::make_pair(3, gmo::addBox(0.0, 0.0, 0.0, 1.0, 1.0, 0.05) )
-        );
-    objects.push_back(
-        std::make_pair(3, gmo::addBox(0.0, 0.0, 0.05, 1.0, 1.0, 0.05) )
-        );
-
-    std::vector<std::pair<int, int>> tools;
-    gmsh::vectorpair odt;
-    std::vector<gmsh::vectorpair> odtm;
-    gmo::fragment(objects, tools, odt, odtm);
-    */
-    gmo::addBox(0.0, 0.0, 0.0, 1.0, 1.0, 0.1);
-    
-    gmo::synchronize();
-
-    gvp_t vp;
-    gm::getEntities(vp);
-    gmm::setSize(vp, mesh_h);
-}
-#endif
-
 template<typename State>
 void initialize_solver(const model& mod, State& state, const maxwell::parameter_loader& mpl)
 {
@@ -196,265 +167,10 @@ reinit_materials(const model& mod, maxwell::solver_state& state, maxwell::parame
 }
 #endif /* ENABLE_EXP_GEOMETRIC_DIODE */
 
-struct maxwell_6ple {
-    double Ex;
-    double Ey;
-    double Ez;
-    double Hx;
-    double Hy;
-    double Hz;
-
-    maxwell_6ple()
-        : Ex(0.0), Ey(0.0), Ez(0.0),
-          Hx(0.0), Hy(0.0), Hz(0.0)
-    {}
-
-    maxwell_6ple(const vec3d& E, const vec3d& H)
-        : Ex( E(0) ), Ey( E(1) ), Ez( E(2) ),
-          Hx( H(0) ), Hy( H(1) ), Hz( H(2) )
-    {}
-
-    maxwell_6ple(double pEx, double pEy, double pEz, double pHx, double pHy, double pHz)
-        : Ex(pEx), Ey(pEy), Ez(pEz), Hx(pHx), Hy(pHy), Hz(pHz)
-    {}
-
-    maxwell_6ple& operator+=(const maxwell_6ple& other)
-    {
-        Ex += other.Ex;
-        Ey += other.Ey;
-        Ez += other.Ez;
-        Hx += other.Hx;
-        Hy += other.Hy;
-        Hz += other.Hz;
-        return *this;
-    }
-
-    maxwell_6ple operator+(const maxwell_6ple& other) const
-    {
-        maxwell_6ple ret = *this;
-        ret += other;
-        return ret;
-    }
-
-    maxwell_6ple& operator-=(const maxwell_6ple& other)
-    {
-        Ex -= other.Ex;
-        Ey -= other.Ey;
-        Ez -= other.Ez;
-        Hx -= other.Hx;
-        Hy -= other.Hy;
-        Hz -= other.Hz;
-        return *this;
-    }
-
-    maxwell_6ple operator-(const maxwell_6ple& other) const
-    {
-        maxwell_6ple ret = *this;
-        ret -= other;
-        return ret;
-    }
-
-    maxwell_6ple& operator*=(const maxwell_6ple& other)
-    {
-        Ex *= other.Ex;
-        Ey *= other.Ey;
-        Ez *= other.Ez;
-        Hx *= other.Hx;
-        Hy *= other.Hy;
-        Hz *= other.Hz;
-        return *this;
-    }
-
-    maxwell_6ple operator*(const maxwell_6ple& other) const
-    {
-        maxwell_6ple ret = *this;
-        ret *= other;
-        return ret;
-    }
-
-    maxwell_6ple& operator*=(double other)
-    {
-        Ex *= other;
-        Ey *= other;
-        Ez *= other;
-        Hx *= other;
-        Hy *= other;
-        Hz *= other;
-        return *this;
-    }
-
-    maxwell_6ple operator*(double other) const
-    {
-        maxwell_6ple ret = *this;
-        ret *= other;
-        return ret;
-    }
-
-    maxwell_6ple& operator/=(const maxwell_6ple& other)
-    {
-        Ex /= other.Ex;
-        Ey /= other.Ey;
-        Ez /= other.Ez;
-        Hx /= other.Hx;
-        Hy /= other.Hy;
-        Hz /= other.Hz;
-        return *this;
-    }
-
-    maxwell_6ple operator/(const maxwell_6ple& other) const
-    {
-        maxwell_6ple ret = *this;
-        ret /= other;
-        return ret;
-    }
-};
-
-static
-maxwell_6ple operator*(double d, const maxwell_6ple& m)
-{
-    return m*d;
-}
-
-static std::ostream&
-operator<<(std::ostream& os, const maxwell_6ple& m)
-{
-    os << m.Ex << " " << m.Ey << " " << m.Ez << " ";
-    os << m.Hx << " " << m.Hy << " " << m.Hz;
-    return os;
-}
-
-static maxwell_6ple
-sqrt(const maxwell_6ple& m)
-{
-    maxwell_6ple ret;
-    ret.Ex = std::sqrt(m.Ex);
-    ret.Ey = std::sqrt(m.Ey);
-    ret.Ez = std::sqrt(m.Ez);
-    ret.Hx = std::sqrt(m.Hx);
-    ret.Hy = std::sqrt(m.Hy);
-    ret.Hz = std::sqrt(m.Hz);
-    return ret;
-}
-
-static double
-hsum(const maxwell_6ple& m)
-{
-    return m.Ex + m.Ey + m.Ez + m.Hx + m.Hy + m.Hz;
-}
-
-static maxwell_6ple
-eval_field(const maxwell::field& f, size_t ofs, size_t cbs,
-    const vecxd& phi)
-{
-    assert(ofs+cbs <= f.num_dofs);
-    maxwell_6ple ret;
-    ret.Ex = f.Ex.segment(ofs, cbs).dot(phi);
-    ret.Ey = f.Ey.segment(ofs, cbs).dot(phi);
-    ret.Ez = f.Ez.segment(ofs, cbs).dot(phi);
-    ret.Hx = f.Hx.segment(ofs, cbs).dot(phi);
-    ret.Hy = f.Hy.segment(ofs, cbs).dot(phi);
-    ret.Hz = f.Hz.segment(ofs, cbs).dot(phi);
-    return ret;
-}
-
-static maxwell_6ple
-compute_error(const model& mod, maxwell::solver_state& state,
-    maxwell::parameter_loader& mpl)
-{
-    maxwell_6ple err;
-
-    if (not mpl.has_analytical_solution())
-        throw std::invalid_argument("analytical solution not defined in the Lua configuration");
-
-    for (auto& e : mod)
-    {
-        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 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 );
-
-            for (size_t iQp = 0; iQp < pqps.size(); iQp++)
-            {
-                const auto& pqp = pqps[iQp];
-                const vecxd phi = re.basis_functions(iQp);
-
-                maxwell_6ple Fh = eval_field(state.emf_curr, ofs, cbs, phi);
-
-                auto [E, H] = mpl.eval_analytical_solution(e.material_tag(),
-                    pqp.point(), state.curr_time);
-                maxwell_6ple F(E,H);
-                maxwell_6ple Fdiff = F - Fh;
-                err += Fdiff*Fdiff*pqp.weight();
-            }
-        }
-    }
-
-#ifdef USE_MPI
-    MPI_Allreduce(MPI_IN_PLACE, &err, 6, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
-#endif /* USE_MPI */
-    return sqrt(err);
-}
-
-static maxwell_6ple
-compute_energy(const model& mod, const maxwell::solver_state& state,
-    const maxwell::parameter_loader& mpl)
-{
-    maxwell_6ple energy;
-
-    for (auto& e : mod)
-    {
-        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 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 );
-
-            maxwell_6ple mat(eps, eps, eps, mu, mu, mu);
-
-            for (size_t iQp = 0; iQp < pqps.size(); iQp++)
-            {
-                const auto& pqp = pqps[iQp];
-                const vecxd phi = re.basis_functions(iQp);
-                maxwell_6ple Fh = eval_field(state.emf_curr, ofs, cbs, phi);
-                energy += 0.5*mat*Fh*Fh*pqp.weight();
-            }
-        }
-    }
-
-#ifdef USE_MPI
-    MPI_Allreduce(MPI_IN_PLACE, &energy, 6, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
-#endif /* USE_MPI */
-    return energy;
-}
-
-#ifdef ENABLE_GPU_SOLVER
-static void
-validate(const model&, maxwell::solver_state_gpu&,
-    maxwell::parameter_loader&, std::ofstream&)
-{}
-#endif /* ENABLE_GPU_SOLVER */
 
+/****************************************************************************/
+/* Solver driver
+ */
 template<typename State>
 void solver_mainloop(const model& mod, State& state, maxwell::parameter_loader& mpl)
 {
@@ -533,19 +249,25 @@ void solver_mainloop(const model& mod, State& state, maxwell::parameter_loader&
         std::cout << showcursor << std::endl;
 }
 
+/****************************************************************************/
+/* Register Lua usertypes valid for both CPU and GPU
+ */
 static void
 register_lua_usertypes(maxwell::parameter_loader& mpl)
 {
+    using namespace maxwell;
+
     sol::state& lua = mpl.lua_state();
-    sol::usertype<maxwell_6ple> maxwell_6ple_ut = lua.new_usertype<maxwell_6ple>("maxwell_6ple",
-        sol::constructors<maxwell_6ple()>()
+    sol::usertype<field_values> field_values_ut =
+        lua.new_usertype<field_values>("field_values",
+            sol::constructors<field_values()>()
             );
-    maxwell_6ple_ut["Ex"] = &maxwell_6ple::Ex;
-    maxwell_6ple_ut["Ey"] = &maxwell_6ple::Ey;
-    maxwell_6ple_ut["Ez"] = &maxwell_6ple::Ez;
-    maxwell_6ple_ut["Hx"] = &maxwell_6ple::Hx;
-    maxwell_6ple_ut["Hy"] = &maxwell_6ple::Hy;
-    maxwell_6ple_ut["Hz"] = &maxwell_6ple::Hz;
+    field_values_ut["Ex"] = &field_values::Ex;
+    field_values_ut["Ey"] = &field_values::Ey;
+    field_values_ut["Ez"] = &field_values::Ez;
+    field_values_ut["Hx"] = &field_values::Hx;
+    field_values_ut["Hy"] = &field_values::Hy;
+    field_values_ut["Hz"] = &field_values::Hz;
 
 #ifdef ENABLE_EXP_GEOMETRIC_DIODE
     sol::usertype<point_3d> point_3d_ut = lua.new_usertype<point_3d>("point",
@@ -559,6 +281,9 @@ register_lua_usertypes(maxwell::parameter_loader& mpl)
 #endif /* ENABLE_EXP_GEOMETRIC_DIODE */
 }
 
+/****************************************************************************/
+/* Register Lua usertypes valid only for CPU
+ */
 static void
 register_lua_usertypes_bystate(maxwell::parameter_loader& mpl,
     model& mod, maxwell::solver_state& state)
@@ -583,6 +308,9 @@ register_lua_usertypes_bystate(maxwell::parameter_loader& mpl,
 #endif /*ENABLE_EXP_GEOMETRIC_DIODE */
 }
 
+/****************************************************************************/
+/* Register Lua usertypes valid only for GPU
+ */
 #ifdef ENABLE_GPU_SOLVER
 static void
 register_lua_usertypes_bystate(maxwell::parameter_loader& mpl,
@@ -590,6 +318,9 @@ register_lua_usertypes_bystate(maxwell::parameter_loader& mpl,
 {}
 #endif
 
+/****************************************************************************/
+/* Main
+ */
 int main(int argc, char *argv[])
 {
 #ifdef USE_MPI