diff --git a/doc/lua_api.md b/doc/lua_api.md
index 67c03cca57ed28f74775bc29e649c3cf5fbcab86..409abc65f8e8419bf9f065b493387cc153e20979 100644
--- a/doc/lua_api.md
+++ b/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)
diff --git a/include/param_loader.h b/include/param_loader.h
index 6e2afd77f3c40dfd7850168c9bed186296fd3e28..358c209763871d373283797d8282b727606b2f8d 100644
--- a/include/param_loader.h
+++ b/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; }
};
diff --git a/src/maxwell_solver.cpp b/src/maxwell_solver.cpp
index 60887a56de3783fc1e78895bca943132613e6934..72ac3c5b636055881a8f57b437afcc12a013c5d0 100644
--- a/src/maxwell_solver.cpp
+++ b/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);
}
diff --git a/src/param_loader.cpp b/src/param_loader.cpp
index 5399e34058333bc2b6ba4c313623cabf06f31631..eb9b97127e14d262dd3c981a3ea70db52230f038 100644
--- a/src/param_loader.cpp
+++ b/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
{