Code to handle geometric diode.

doc/lua_api.md

@@ -125,10 +125,18 @@ Currently, only a surface current condition is available:
 None.
 
 ### Callbacks
+- `before_start()`: if defined, the solver calls this function just before beginning the timestepping.
 - `on_timestep(ts)`: if defined, the solver calls this function at every timestep. The solver passes the current timestep number in the parameter `ts`. The function is called **before** the solver state is advanced.
 - `electric_initial_condition(x, y, z)`: if defined, the solver calls this function at the beginning of the simulation in order to initialize the electric field. The parameters `x`, `y` and `z` represent the point at which the field needs to be evaluated. The function has to return a triple specifying the three field components, i.e. `return Ex, Ey, Ez`.
 - `magnetic_initial_condition(x, y, z)`: if defined, the solver calls this function at the beginning of the simulation in order to initialize the magnetic field. The parameters `x`, `y` and `z` represent the point at which the field needs to be evaluated. The function has to return a triple specifying the three field components, i.e. `return Hx, Hy, Hz`.
 
+### Data types
+- `point`: a three dimensional point
+
+### Variables
+- `internal.model`: opaque reference to the internal representation of the GMSH model. Valid only after the solver is fully initialized, i.e. inside the callback functions.
+- `internal.state`: opaque reference to the internal representation of the solver state. Valid only after the solver is fully initialized, i.e. inside the callback functions.
+
 ### Postprocessing
 The export of the fields E, H and J to Silo is controlled using the appropriate sections in the `postpro` table.
 - `postpro["E"].silo_mode` controls the export of the electric field E. Possible values are `"none"` (no export), `"nodal"` (export as a Silo nodal variable) and `"zonal"` (export as a Silo zonal variable) 

include/param_loader.h

@@ -64,6 +64,9 @@ public:
     void            source_was_cleared(void);
 
     void            call_timestep_callback(size_t);
+    void            call_initialization_callback();
+
+    sol::state&     lua_state(void) { return lua; }
 };
 
 

src/maxwell_solver.cpp

@@ -68,21 +68,21 @@ void initialize_solver(const model& mod, State& state, const maxwell::parameter_
     init_matparams(mod, state, mpl);
 }
 
-class line_integrator
+class integration_line
 {
     std::vector<size_t>     element_tags;
     std::vector<point_3d>   sampling_points_phys;
     std::vector<point_3d>   sampling_points_ref;
-    point_3d                increment;
+    point_3d                sampling_increment;
 
 public:
-    line_integrator(const point_3d& start,
+    integration_line(const point_3d& start,
         const point_3d& end, size_t samples)
     {
-        increment = (end - start)/samples;
+        sampling_increment = (end - start)/samples;
         for (size_t i = 0; i < samples; i++)
         {
-            point_3d sp = start + increment*(i+0.5);
+            point_3d sp = start + sampling_increment*(i+0.5);
             sampling_points_phys.push_back(sp);
 
             size_t              etag;
@@ -97,39 +97,95 @@ public:
         }
     }
 
-    double integrate(const model& mod, const vecxd& dofs)
+    size_t sampling_points() const
     {
-        double integral_val = 0.0;
-        for (size_t i = 0; i < sampling_points_phys.size(); i++)
-        {
-            auto [entnum, ofs] = mod.lookup_tag(element_tags.at(i));
-            auto& e = mod.entity_at(entnum);
-            auto& pe = e.cell(ofs);
-            auto& re = e.cell_refelem(pe);
-            auto dof_ofs = e.cell_global_dof_offset(ofs);
-            auto dof_num = re.num_basis_functions();
-            vecxd phi = re.basis_functions( sampling_points_ref[i] );
-            vecxd sol = dofs.segment(dof_ofs, dof_num);
-            integral_val += sol.dot(phi) * increment.x();
-        }
+        assert(element_tags.size() == sampling_points_phys.size());
+        assert(element_tags.size() == sampling_points_ref.size());
+        return element_tags.size();
+    }
+
+    std::tuple<size_t, point_3d, point_3d>
+    operator[](size_t pos)
+    {
+        assert(element_tags.size() == sampling_points_phys.size());
+        assert(element_tags.size() == sampling_points_ref.size());
+        assert(pos < element_tags.size());
+
+        return std::make_tuple(
+            element_tags[pos],
+            sampling_points_phys[pos],
+            sampling_points_ref[pos]
+        );
+    }
 
-        return integral_val;
+    point_3d increment() const
+    {
+        return sampling_increment;
     }
+
 };
 
 double
-integrate(const model& mod, const maxwell::solver_state& state, line_integrator& lint)
+integrate_electric_field(const model& mod, const maxwell::solver_state& state,
+    integration_line& iline)
 {
-    return lint.integrate(mod, state.emf_curr.Ex);
+    double integral_val = 0.0;
+    auto si = iline.increment();
+
+    for (size_t i = 0; i < iline.sampling_points(); i++)
+    {
+        auto [tag, physp, refp] = iline[i];
+        auto [entnum, ofs] = mod.lookup_tag(tag);
+        auto& e = mod.entity_at(entnum);
+        auto& pe = e.cell(ofs);
+        auto& re = e.cell_refelem(pe);
+        auto dof_ofs = e.cell_global_dof_offset(ofs);
+        auto dof_num = re.num_basis_functions();
+        vecxd phi = re.basis_functions( refp );
+        vecxd locEx = state.emf_curr.Ex.segment(dof_ofs, dof_num);
+        vecxd locEy = state.emf_curr.Ey.segment(dof_ofs, dof_num);
+        vecxd locEz = state.emf_curr.Ez.segment(dof_ofs, dof_num);
+        integral_val += locEx.dot(phi) * si.x() +
+                        locEy.dot(phi) * si.y() +
+                        locEz.dot(phi) * si.z();
+    }
+
+    return integral_val;
 }
 
+#ifdef ENABLE_GPU_SOLVER
 double
-integrate(const model& mod, const maxwell::solver_state_gpu& state, line_integrator& lint)
+integrate(const model& mod, const maxwell::solver_state_gpu& state, integration_line& lint)
 {
-    maxwell::field f;
-    f.resize( mod.num_dofs() );
-    state.emf_curr.copyout(f);
-    return lint.integrate(mod, f.Ex);
+    return 0.0;
+}
+#endif
+
+void
+reinit_materials(const model& mod, maxwell::solver_state& state, maxwell::parameter_loader& mpl)
+{
+    for (auto& e : mod)
+    {
+        auto tag = e.gmsh_tag();
+        for (size_t iT = 0; iT < e.num_cells(); iT++)
+        {
+            auto& pe = e.cell(iT);
+            auto& re = e.cell_refelem(pe);
+            auto bar = pe.barycenter();
+            auto epsilon = mpl.epsilon(tag, bar);
+            //auto mu = mpl.mu(tag, bar);
+            auto sigma = mpl.sigma(tag, bar);
+            auto ofs = e.cell_global_dof_offset(iT);
+
+            for (size_t iD = 0; iD < re.num_basis_functions(); iD++)
+            {
+                //state.matparams.inv_epsilon(ofs+iD) = 1./epsilon;
+                //state.matparams.inv_mu(ofs+iD) = 1./mu;
+                state.matparams.sigma(ofs+iD) = sigma;
+                state.matparams.sigma_over_epsilon(ofs+iD) = sigma/epsilon;
+            }
+        }
+    }
 }
 
 template<typename State>
@@ -152,11 +208,12 @@ void test_it(const model& mod, State& state, maxwell::parameter_loader& mpl)
     std::cout << "s each." << std::endl;
     std::cout << "Time integrator: " << mpl.sim_timeIntegratorName() << std::endl;
 
-    point_3d start(-46e-9, 0.0, 38e-9);
-    point_3d end(46e-9, 0.0, 38e-9);
-    line_integrator lint(start, end, 50);
+    mpl.call_initialization_callback();
+
+    auto relmat = [&](){ reinit_materials(mod, state, mpl); };
 
-    std::ofstream ofs("voltage.dat");
+    sol::state& lua = mpl.lua_state();
+    lua["relmat"] = relmat;
 
     timecounter tc;
     tc.tic();
@@ -179,11 +236,43 @@ void test_it(const model& mod, State& state, maxwell::parameter_loader& mpl)
             std::cout << "Cycle " << i << ": t = " << state.curr_time << " s";
             std::cout << ", DOFs/s: " << (mod.num_dofs()*6*cycle_print_rate)/(time) << std::endl;
             tc.tic();
-            ofs << state.curr_time << " " << integrate(mod, state, lint) << std::endl;
         }
     }
 }
 
+void
+register_lua_usertypes(maxwell::parameter_loader& mpl)
+{
+    sol::state& lua = mpl.lua_state();
+    sol::usertype<point_3d> point_3d_ut = lua.new_usertype<point_3d>("point",
+        sol::constructors<point_3d(), point_3d(double, double, double)>()
+        );
+
+    sol::usertype<integration_line> line_integrator_ut = 
+        lua.new_usertype<integration_line>("integration_line",
+            sol::constructors<integration_line(const point_3d&, const point_3d&, size_t)>()
+        );
+}
+
+void
+register_lua_usertypes_bystate(maxwell::parameter_loader& mpl,
+    model& mod, maxwell::solver_state& state)
+{
+    sol::state& lua = mpl.lua_state();
+
+    lua["internal"] = lua.create_table();
+    lua["internal"]["model"] = &mod;
+    lua["internal"]["state"] = &state;
+    lua["integrate_E"] = integrate_electric_field;
+}
+
+#ifdef ENABLE_GPU_SOLVER
+void
+register_lua_usertypes_bystate(maxwell::parameter_loader& mpl,
+    model& mod, maxwell::solver_state_gpu& state)
+{}
+#endif
+
 int main(int argc, const char *argv[])
 {
     gmsh::initialize();
@@ -198,6 +287,8 @@ int main(int argc, const char *argv[])
     }
 
     maxwell::parameter_loader mpl;
+    register_lua_usertypes(mpl);
+
     if ( not mpl.load_file( argv[1] ) )
     {
         std::cout << "Configuration problem, exiting" << std::endl;
@@ -230,6 +321,7 @@ int main(int argc, const char *argv[])
     {
 #ifdef ENABLE_GPU_SOLVER
         maxwell::solver_state_gpu state_g;
+        register_lua_usertypes_bystate(mpl, mod, state_g);
         test_it(mod, state_g, mpl);
 #else
         std::cout << "GPU solver not compiled. Exiting." << std::endl;
@@ -238,6 +330,7 @@ int main(int argc, const char *argv[])
     else
     {
         maxwell::solver_state state_c;
+        register_lua_usertypes_bystate(mpl, mod, state_c);
         test_it(mod, state_c, mpl);
     }
 

src/param_loader.cpp

@@ -211,6 +211,21 @@ parameter_loader_base::call_timestep_callback(size_t timestep)
     }
 }
 
+void
+parameter_loader_base::call_initialization_callback()
+{
+    if ( lua["before_start"].valid() )
+    {
+        try {
+            lua["before_start"]();
+        }
+        catch (...)
+        {
+            std::cout << "An error was thrown calling before_start()" << std::endl;
+        }
+    }
+}
+
 bool
 parameter_loader_base::source_has_changed_state(void) const
 {