Parametrizing the solver with Lua: plane wave BCs on CPU.

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")

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;

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>

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
 {

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

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)
 {
@@ -57,7 +125,6 @@ void init_from_model(const model& mod, solver_state& state)
     state.dy.resize( mod.num_dofs() );
     state.dz.resize( mod.num_dofs() );
     
-    
     state.k1.resize( mod.num_dofs() );
     state.k2.resize( mod.num_dofs() );
     state.k3.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;

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
     {