From ecec18f3807f0df415ec2fc706ead8a1e5ee76c3 Mon Sep 17 00:00:00 2001
From: Matteo Cicuttin <datafl4sh@toxicnet.eu>
Date: Thu, 6 May 2021 18:18:47 +0200
Subject: [PATCH] Parametrizing the solver with Lua: interface conditions on
CPU.
---
include/maxwell_interface.h | 16 +++++++--
src/maxwell_common.cpp | 42 ++++++++++++++++++++++-
src/maxwell_cpu.cpp | 68 +++++++++++++++++++++++++++++++++++--
src/maxwell_solver.cpp | 1 +
4 files changed, 121 insertions(+), 6 deletions(-)
diff --git a/include/maxwell_interface.h b/include/maxwell_interface.h
index 0c6a57a..27d62b4 100644
--- a/include/maxwell_interface.h
+++ b/include/maxwell_interface.h
@@ -34,6 +34,14 @@ enum class boundary_condition
SURFACE_CURRENT
};
+enum class interface_condition
+{
+ UNSPECIFIED,
+ E_FIELD,
+ H_FIELD,
+ SURFACE_CURRENT
+};
+
class parameter_loader : public parameter_loader_base
{
void init(void);
@@ -54,9 +62,11 @@ public:
bool initial_Hfield_defined(void) const;
vec3d initial_Hfield(const point_3d&) const;
- vec3d eval_plane_wave_E(int, const point_3d&, double) const;
+ vec3d eval_boundary_source(int, const point_3d&, double) const;
+ vec3d eval_interface_source(int, const point_3d&, double) const;
- boundary_condition boundary_type(int tag) const;
+ boundary_condition boundary_type(int tag) const;
+ interface_condition interface_type(int tag) const;
};
#endif /* HIDE_THIS_FROM_NVCC */
@@ -243,7 +253,6 @@ struct solver_state_gpu
#endif /* ENABLE_GPU_SOLVER */
#ifndef HIDE_THIS_FROM_NVCC
-
void init_from_model(const model&, solver_state&);
void init_matparams(const model&, solver_state&, const parameter_loader&);
void apply_operator(solver_state&, const field&, field&);
@@ -261,6 +270,7 @@ void timestep(solver_state_gpu&);
void init_boundary_conditions(const model&, const parameter_loader&, vecxd&);
void init_lax_milgram(const model&, const parameter_loader&, vecxd&, vecxd&, vecxd&, vecxd&);
void eval_boundary_sources(const model&, const parameter_loader&, solver_state&, field&);
+void eval_interface_sources(const model&, const parameter_loader&, solver_state&, field&);
#endif /* HIDE_THIS_FROM_NVCC */
template<typename Function>
diff --git a/src/maxwell_common.cpp b/src/maxwell_common.cpp
index 7a4ff69..9f5dff7 100644
--- a/src/maxwell_common.cpp
+++ b/src/maxwell_common.cpp
@@ -297,8 +297,39 @@ parameter_loader::boundary_type(int tag) const
return boundary_condition::UNSPECIFIED;
}
+interface_condition
+parameter_loader::interface_type(int tag) const
+{
+ auto bnd_data = lua["ifaceconds"][tag];
+ if (not bnd_data.valid())
+ return interface_condition::UNSPECIFIED;
+
+ auto kind = bnd_data["kind"];
+ if (not kind.valid())
+ {
+ std::cout << "[CONFIG] warning: 'kind' not specified on interface ";
+ std::cout << tag << std::endl;
+ return interface_condition::UNSPECIFIED;
+ }
+
+ std::string kind_str = kind;
+
+ if (kind_str == "E_field")
+ return interface_condition::E_FIELD;
+
+ if (kind_str == "H_field")
+ return interface_condition::H_FIELD;
+
+ if (kind_str == "surface_current")
+ return interface_condition::SURFACE_CURRENT;
+
+ std::cout << "[CONFIG] warning: 'kind' invalid on boundary ";
+ std::cout << tag << std::endl;
+ return interface_condition::UNSPECIFIED;
+}
+
vec3d
-parameter_loader::eval_plane_wave_E(int tag, const point_3d& pt, double t) const
+parameter_loader::eval_boundary_source(int tag, const point_3d& pt, double t) const
{
vec3d ret;
sol::tie(ret(0), ret(1), ret(2)) =
@@ -306,4 +337,13 @@ parameter_loader::eval_plane_wave_E(int tag, const point_3d& pt, double t) const
return ret;
}
+vec3d
+parameter_loader::eval_interface_source(int tag, const point_3d& pt, double t) const
+{
+ vec3d ret;
+ sol::tie(ret(0), ret(1), ret(2)) =
+ lua["ifaceconds"][tag]["source"](tag, pt.x(), pt.y(), pt.z(), t);
+ return ret;
+}
+
} // namespace maxwell
\ No newline at end of file
diff --git a/src/maxwell_cpu.cpp b/src/maxwell_cpu.cpp
index 3ce0b41..c8ffddb 100644
--- a/src/maxwell_cpu.cpp
+++ b/src/maxwell_cpu.cpp
@@ -80,10 +80,10 @@ eval_boundary_sources(const model& mod, const parameter_loader& mpl,
case boundary_condition::PLANE_WAVE_E: {
auto fE = [&](const point_3d& pt) -> vec3d {
- return mpl.eval_plane_wave_E(bd.gmsh_entity, pt, state.curr_time);
+ return mpl.eval_boundary_source(bd.gmsh_entity, pt, state.curr_time);
};
auto fH = [&](const point_3d& pt) -> vec3d {
- vec3d E = mpl.eval_plane_wave_E(bd.gmsh_entity, pt, state.curr_time);
+ vec3d E = mpl.eval_boundary_source(bd.gmsh_entity, pt, state.curr_time);
return n.cross(E)/Z;
};
matxd prjE = e.project_on_face(iF, fE);
@@ -116,6 +116,70 @@ eval_boundary_sources(const model& mod, const parameter_loader& mpl,
}
}
+void
+eval_interface_sources(const model& mod, const parameter_loader& mpl,
+ solver_state& state, field& srcfield)
+{
+ auto& bds = mod.boundary_descriptors();
+ assert( srcfield.num_dofs == mod.num_fluxes() );
+
+ size_t face_num_base = 0;
+ for (auto& e : mod)
+ {
+ for (size_t iF = 0; iF < e.num_faces(); iF++)
+ {
+ auto& pf = e.face(iF);
+ auto bar = pf.barycenter();
+ auto eps = mpl.epsilon(e.gmsh_tag(), bar);
+ auto mu = mpl.mu(e.gmsh_tag(), bar);
+ auto Z = std::sqrt(mu/eps);
+
+ vec3d n = state.eds[e.number()].normals.row(iF);
+
+ auto gfnum = face_num_base + iF;
+ auto bd = bds[gfnum];
+ if (bd.type == face_type::NONE or bd.type == face_type::BOUNDARY)
+ continue;
+
+ switch ( mpl.interface_type(bd.gmsh_entity) )
+ {
+ case interface_condition::E_FIELD: {
+ /* This is BS actually... */
+ auto fE = [&](const point_3d& pt) -> vec3d {
+ return mpl.eval_interface_source(bd.gmsh_entity, pt, state.curr_time);
+ };
+
+ matxd prjE = e.project_on_face(iF, fE);
+ auto num_bf = prjE.rows();
+ srcfield.Ex.segment(gfnum*num_bf, num_bf) = prjE.col(0);
+ srcfield.Ey.segment(gfnum*num_bf, num_bf) = prjE.col(1);
+ srcfield.Ez.segment(gfnum*num_bf, num_bf) = prjE.col(2);
+ } break;
+
+ case interface_condition::SURFACE_CURRENT: {
+ auto fJ = [&](const point_3d& pt) -> vec3d {
+ vec3d J = mpl.eval_interface_source(bd.gmsh_entity, pt, state.curr_time);
+ return n.cross(J);
+ };
+
+ matxd prjJ = e.project_on_face(iF, fJ);
+ auto num_bf = prjJ.rows();
+ srcfield.Hx.segment(gfnum*num_bf, num_bf) = prjJ.col(0);
+ srcfield.Hy.segment(gfnum*num_bf, num_bf) = prjJ.col(1);
+ srcfield.Hz.segment(gfnum*num_bf, num_bf) = prjJ.col(2);
+ } break;
+
+ case interface_condition::H_FIELD:
+ break;
+
+ default:
+ break;
+ }
+ }
+
+ face_num_base += e.num_faces();
+ }
+}
void init_from_model(const model& mod, solver_state& state)
{
diff --git a/src/maxwell_solver.cpp b/src/maxwell_solver.cpp
index 68e4371..5e3fdfe 100644
--- a/src/maxwell_solver.cpp
+++ b/src/maxwell_solver.cpp
@@ -65,6 +65,7 @@ void test_it(const model& mod, State& state, const maxwell::parameter_loader& mp
for(size_t i = 0; i < num_timesteps; i++)
{
maxwell::eval_boundary_sources(mod, mpl, state, state.bndsrcs);
+ maxwell::eval_interface_sources(mod, mpl, state, state.bndsrcs);
timestep(state);
std::stringstream ss;
ss << "maxwell_" << i << ".silo";
--
GitLab