diff --git a/CMakeLists.txt b/CMakeLists.txt
index 1b198235669b73e991f17b2107f91a28b22eb7db..e3b1b3be6fa2fd70eaf82d98b7c0625ef0c2c74f 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -76,7 +76,7 @@ if (OPT_ENABLE_GPU_SOLVER)
if (CUDAToolkit_FOUND)
set(CMAKE_CUDA_COMPILER ${CUDAToolkit_NVCC_EXECUTABLE})
enable_language(CUDA)
- set(CMAKE_CUDA_STANDARD 14)
+ set(CMAKE_CUDA_STANDARD 17)
set(CMAKE_CUDA_STANDARD_REQUIRED TRUE)
#set(CMAKE_CUDA_FLAGS "${CMAKE_CUDA_FLAGS} -lineinfo --resource-usage -Xptxas -dlcm=ca")
set(CMAKE_CUDA_FLAGS "${CMAKE_CUDA_FLAGS} -lineinfo --resource-usage -prec-div=true")
diff --git a/include/entity.h b/include/entity.h
index 11296407ebe781bf037af16ce785dfa4d3ef78f9..70df6223d9fc821190368363cd348e6508f37c76 100644
--- a/include/entity.h
+++ b/include/entity.h
@@ -66,6 +66,8 @@ public:
void project(const scalar_function&, vecxd&) const;
void project(const vector_function&, vecxd&, vecxd&, vecxd&) const;
+ matxd project_on_face(size_t iF, const vector_function&) const;
+
size_t num_cells(void) const;
size_t num_cell_orientations(void) const;
std::vector<size_t> num_cells_per_orientation(void) const;
diff --git a/include/maxwell_interface.h b/include/maxwell_interface.h
index 0d2910907c8d9fa305f379f1a43b936c2ec6a47e..0c6a57a66078443657ab31645185f66a267cf1e1 100644
--- a/include/maxwell_interface.h
+++ b/include/maxwell_interface.h
@@ -54,6 +54,8 @@ 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;
+
boundary_condition boundary_type(int tag) const;
};
#endif /* HIDE_THIS_FROM_NVCC */
@@ -208,6 +210,8 @@ struct solver_state
field k3;
field k4;
field tmp;
+
+ field bndsrcs;
};
#ifdef ENABLE_GPU_SOLVER
@@ -256,6 +260,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&);
#endif /* HIDE_THIS_FROM_NVCC */
template<typename Function>
diff --git a/src/entity.cpp b/src/entity.cpp
index 73f60ee0dd8a05f6b3cce9c29a5194bd8b59a47b..17c42013e74235d65bba756fe8b56ea2c8e84b48 100644
--- a/src/entity.cpp
+++ b/src/entity.cpp
@@ -165,6 +165,45 @@ entity::project(const vector_function& function, vecxd& pfx, vecxd& pfy, vecxd&
}
}
+matxd
+entity::project_on_face(size_t iF, const vector_function& function) const
+{
+ assert(iF < physical_faces.size());
+
+ const auto& pf = physical_faces[iF];
+ assert( pf.orientation() < reference_faces.size() );
+ const auto& rf = reference_faces[pf.orientation()];
+ const auto pqps = pf.integration_points();
+ const auto dets = pf.determinants();
+
+ assert(pqps.size() == dets.size());
+
+ size_t num_bf = rf.num_basis_functions();
+
+ matxd mass = matxd::Zero(num_bf, num_bf);
+ vecxd rhs_x = vecxd::Zero(num_bf);
+ vecxd rhs_y = vecxd::Zero(num_bf);
+ vecxd rhs_z = vecxd::Zero(num_bf);
+ for (size_t iQp = 0; iQp < pqps.size(); iQp++)
+ {
+ const auto& pqp = pqps[iQp];
+ const vecxd phi = rf.basis_functions(iQp);
+ mass += pqp.weight() * phi * phi.transpose();
+ vec3d fval = function(pqp.point());
+ rhs_x += pqp.weight() * fval(0) * phi;
+ rhs_y += pqp.weight() * fval(1) * phi;
+ rhs_z += pqp.weight() * fval(2) * phi;
+ }
+
+ matxd ret(num_bf, 3);
+
+ Eigen::LDLT<matxd> mass_ldlt(mass);
+ ret.col(0) = mass_ldlt.solve(rhs_x);
+ ret.col(1) = mass_ldlt.solve(rhs_y);
+ ret.col(2) = mass_ldlt.solve(rhs_z);
+ return ret;
+}
+
const reference_element&
entity::cell_refelem(size_t iO) const
{
diff --git a/src/maxwell_common.cpp b/src/maxwell_common.cpp
index d1fd40bd05be8d3ea4a46bae3a236b91b4fc68b4..7a4ff6924c83ff6a8bd8c865c83c45fa4770556d 100644
--- a/src/maxwell_common.cpp
+++ b/src/maxwell_common.cpp
@@ -138,8 +138,8 @@ init_lax_milgram(const model& mod, const parameter_loader& mpl,
parameter_loader::parameter_loader()
{
- lua["const"]["eps-1"] = 8.8541878128e-12;
- lua["const"]["mu-1"] = 4e-7*M_PI;
+ lua["const"]["eps0"] = 8.8541878128e-12;
+ lua["const"]["mu0"] = 4e-7*M_PI;
}
bool
@@ -297,4 +297,13 @@ parameter_loader::boundary_type(int tag) const
return boundary_condition::UNSPECIFIED;
}
+vec3d
+parameter_loader::eval_plane_wave_E(int tag, const point_3d& pt, double t) const
+{
+ vec3d ret;
+ sol::tie(ret(0), ret(1), ret(2)) =
+ lua["bndconds"][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 d47639f050be418546a999922e0de444f550799b..3ce0b419e730f915811fc697888f3cba91ccf195 100644
--- a/src/maxwell_cpu.cpp
+++ b/src/maxwell_cpu.cpp
@@ -48,6 +48,74 @@ init_matparams(const model& mod, solver_state& state,
init_boundary_conditions(mod, mpl, state.matparams.bc_coeffs);
}
+void
+eval_boundary_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::INTERFACE)
+ continue;
+
+ switch ( mpl.boundary_type(bd.gmsh_entity) )
+ {
+ case boundary_condition::E_FIELD:
+ break;
+
+ 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);
+ };
+ auto fH = [&](const point_3d& pt) -> vec3d {
+ vec3d E = mpl.eval_plane_wave_E(bd.gmsh_entity, pt, state.curr_time);
+ return n.cross(E)/Z;
+ };
+ 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);
+
+ matxd prjH = e.project_on_face(iF, fH);
+ srcfield.Hx.segment(gfnum*num_bf, num_bf) = prjH.col(0);
+ srcfield.Hy.segment(gfnum*num_bf, num_bf) = prjH.col(1);
+ srcfield.Hz.segment(gfnum*num_bf, num_bf) = prjH.col(2);
+ } break;
+
+ case boundary_condition::PLANE_WAVE_H:
+ break;
+
+ case boundary_condition::SURFACE_CURRENT:
+ break;
+
+ case boundary_condition::H_FIELD:
+ break;
+
+ default:
+ break;
+ }
+ }
+
+ face_num_base += e.num_faces();
+ }
+}
+
void init_from_model(const model& mod, solver_state& state)
{
@@ -56,7 +124,6 @@ void init_from_model(const model& mod, solver_state& state)
state.dx.resize( mod.num_dofs() );
state.dy.resize( mod.num_dofs() );
state.dz.resize( mod.num_dofs() );
-
state.k1.resize( mod.num_dofs() );
state.k2.resize( mod.num_dofs() );
@@ -67,7 +134,7 @@ void init_from_model(const model& mod, solver_state& state)
state.jumps.resize( mod.num_fluxes() );
state.fluxes.resize( mod.num_fluxes() );
- state.matparams.bc_coeffs = 2*vecxd::Ones( mod.num_fluxes() );
+ state.bndsrcs.resize( mod.num_fluxes() );
for (auto& e : mod)
{
@@ -158,12 +225,12 @@ compute_fluxes_planar(solver_state& state)
auto gofs = ed.flux_base + lofs;
/* Sources are imposed on jumps */
- auto jEx = state.jumps.Ex[gofs];// - boundary_sources.Ex[g_flofs];
- auto jEy = state.jumps.Ey[gofs];// - boundary_sources.Ey[g_flofs];
- auto jEz = state.jumps.Ez[gofs];// - boundary_sources.Ez[g_flofs];
- auto jHx = state.jumps.Hx[gofs];// + boundary_sources.Hx[g_flofs];
- auto jHy = state.jumps.Hy[gofs];// + boundary_sources.Hy[g_flofs];
- auto jHz = state.jumps.Hz[gofs];// + boundary_sources.Hz[g_flofs];
+ auto jEx = state.jumps.Ex[gofs] - state.bndsrcs.Ex[gofs];
+ auto jEy = state.jumps.Ey[gofs] - state.bndsrcs.Ey[gofs];
+ auto jEz = state.jumps.Ez[gofs] - state.bndsrcs.Ez[gofs];
+ auto jHx = state.jumps.Hx[gofs] + state.bndsrcs.Hx[gofs];
+ auto jHy = state.jumps.Hy[gofs] + state.bndsrcs.Hy[gofs];
+ auto jHz = state.jumps.Hz[gofs] + state.bndsrcs.Hz[gofs];
auto ndotE = nx*jEx + ny*jEy + nz*jEz;
auto ndotH = nx*jHx + ny*jHy + nz*jHz;
@@ -246,7 +313,7 @@ compute_rk4_weighted_sum(solver_state& state, const field& in, double dt, field&
void
apply_operator(solver_state& state, const field& curr, field& next)
-{
+{
compute_curls(state, curr, next);
compute_fluxes(state, curr, next);
}
diff --git a/src/maxwell_solver.cpp b/src/maxwell_solver.cpp
index b85bd2478832ec965d67d8f80a9ffeb6edac2e6a..68e43718658239dbdba7bfd003cd84810f765b27 100644
--- a/src/maxwell_solver.cpp
+++ b/src/maxwell_solver.cpp
@@ -64,6 +64,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);
timestep(state);
std::stringstream ss;
ss << "maxwell_" << i << ".silo";
@@ -111,8 +112,8 @@ int main(int argc, const char *argv[])
#ifdef ENABLE_GPU_SOLVER
if ( mpl.sim_usegpu() )
{
- maxwell::solver_state_gpu state_g;
- test_it(mod, state_g, mpl);
+ //maxwell::solver_state_gpu state_g;
+ //test_it(mod, state_g, mpl);
}
else
{