Parametrizing the solver with Lua: initial conditions.

CMakeLists.txt

@@ -256,6 +256,7 @@ set(LIBGMSHDG_SOURCES src/gmsh_io.cpp
                       src/maxwell_cpu.cpp
                       src/maxwell_gpu.cpp
                       src/maxwell_kernels_cuda.cu
+                      src/param_loader.cpp
                     )
 
 add_library(gmshdg SHARED ${LIBGMSHDG_SOURCES})

include/connectivity.h

@@ -164,6 +164,14 @@ public:
         return std::make_pair(neighbour_descriptor(), false); /* Silence the compiler */
     }
 
+    size_t
+    num_neighbours(const FaceKey& fk)
+    {
+        size_t nn = face_cell_adj.at(fk).num_neighbours();
+        assert(nn == 1 or nn == 2);
+        return nn;
+    }
+
     std::set<ConnData>
     neighbours(size_t iE, size_t iT) const
     {

include/gmsh_io.h

@@ -26,6 +26,23 @@ std::string basis_grad_name(int);
 
 using face_key = std::array<size_t, 3>;
 
+enum class face_type
+{
+    NONE,
+    INTERFACE,
+    BOUNDARY
+};
+
+struct boundary_descriptor
+{
+    face_type   type;
+    int         gmsh_entity;
+
+    boundary_descriptor()
+        : type(face_type::NONE), gmsh_entity(0)
+    {}
+};
+
 class model
 {
     int                             geometric_order;
@@ -33,6 +50,7 @@ class model
 
     std::vector<entity>                         entities;
     std::map<size_t, std::vector<element_key>>  boundary_map;
+    std::vector<boundary_descriptor>            bnd_descriptors;
     element_connectivity<element_key>           conn;
 
     void    map_boundaries(void);
@@ -54,6 +72,11 @@ public:
         return conn;
     }
 
+    const std::vector<boundary_descriptor>& boundary_descriptors() const
+    {
+        return bnd_descriptors;
+    }
+
     std::vector<element_key> get_bnd(size_t which)
     {
         return boundary_map.at(which);
@@ -62,6 +85,7 @@ public:
     size_t  num_dofs() const;
     size_t  num_fluxes() const;
     size_t  num_cells() const;
+    size_t  num_faces() const;
     size_t  num_entities() const;
 
     std::vector<entity>::const_iterator begin() const;

include/param_loader.h

+#pragma once
+
+#include "sol/sol.hpp"
+#include "gmsh_io.h"
+
+class parameter_loader
+{
+protected:
+    sol::state lua;
+
+private:
+    void init(void);
+    bool validate_simulation_params(void) const;
+
+public:
+    parameter_loader();
+
+    bool            load_file(const std::string&);
+
+    int             approx_order(void) const;
+    int             geom_order(void) const;
+    std::string     simulation_name(void) const;
+    double          delta_t(void) const;
+    size_t          timesteps(void) const;
+    std::string     gmsh_model_filename(void) const;
+    bool            use_gpu(void) const;
+    size_t          silo_output_rate(void) const;
+};
+
+class maxwell_parameter_loader : public parameter_loader
+{
+    void init(void);
+
+    bool    validate_materials(const std::string&, int);
+
+public:
+    maxwell_parameter_loader();
+
+    bool    validate_model_params(const model&);
+    double  epsilon(int, const point_3d&) const;
+    double  mu(int, const point_3d&) const;
+    double  sigma(int, const point_3d&) const;
+
+    bool    initial_Efield_defined(void) const;
+    vec3d   Efield(const point_3d&) const;
+
+    bool    initial_Hfield_defined(void) const;
+    vec3d   Hfield(const point_3d&) const;
+};
\ No newline at end of file

src/gmsh_io.cpp

@@ -75,9 +75,31 @@ model::remesh()
 {
     gmm::generate( DIMENSION(3) );
     gmm::setOrder( geometric_order );
+
+    gvp_t gmsh_entities;
+
+    /* Create a table mapping the tag of a boundary to the list of
+     * its face keys. This must happen before import_gmsh_entities()
+     * because otherwise you get spurious 2D entities (and this in
+     * turn because the constructor of entity calls addDiscreteEntity()
+     * to get all the element faces) .*/
+    std::map<size_t, std::vector<face_key>>   boundaries_elemTags;
+    gm::getEntities(gmsh_entities, DIMENSION(2));
+    for (auto [dim, tag] : gmsh_entities)
+    {
+        std::vector<int> elemTypes;
+        gmm::getElementTypes(elemTypes, dim, tag);
+        assert(elemTypes.size() == 1);
+        for (auto& elemType : elemTypes)
+        {
+            element_key_factory ekf(DIMENSION(2), tag, elemType);
+            boundary_map[tag] = ekf.get_keys();
+        }
+    }
+
     entities.clear();
     import_gmsh_entities();
-    map_boundaries();
+    map_boundaries(); /* This must happen after import_gmsh_entities(). */
 }
 
 void
@@ -153,6 +175,16 @@ model::num_cells(void) const
     return ret;
 }
 
+size_t
+model::num_faces(void) const
+{
+    size_t ret = 0;
+    for (const auto& e : entities)
+        ret += e.num_faces();
+
+    return ret;
+}
+
 size_t
 model::num_entities(void) const
 {
@@ -216,20 +248,68 @@ model::import_gmsh_entities(void)
 void
 model::map_boundaries(void)
 {
-    gvp_t entities;
-    /* Create a table mapping the tag of a boundary to the list of
-     * its face keys */
-    std::map<size_t, std::vector<face_key>>   boundaries_elemTags;
-    gm::getEntities(entities, DIMENSION(2));
-    for (auto [dim, tag] : entities)
+    /* Make a vector mapping element_key to entity tag */
+    using bfk_t = std::pair<element_key, int>;
+    std::vector<bfk_t> bfk;
+    for (auto& [tag, keys] : boundary_map)
+        for (auto& k : keys)
+            bfk.push_back( std::make_pair(k, tag) );
+
+    /* Sort it for lookups */
+    auto comp = [](const bfk_t& a, const bfk_t& b) -> bool {
+        return a.first < b.first;
+    };
+    std::sort(bfk.begin(), bfk.end(), comp);
+
+    bnd_descriptors.resize( num_faces() );
+    size_t fbase = 0;
+    /* For each entity */
+    for (auto& e : entities)
     {
-        std::vector<int> elemTypes;
-        gmm::getElementTypes(elemTypes, dim, tag);
-        assert(elemTypes.size() == 1);
-        for (auto& elemType : elemTypes)
+        /* and for each face */
+        for (size_t iF = 0; iF < e.num_faces(); iF++)
         {
-            element_key_factory ekf(DIMENSION(2), tag, elemType);
-            boundary_map[tag] = ekf.get_keys();
+            /* get the element key */
+            element_key fk(e.face(iF));
+            auto lbcomp = [](const bfk_t& a, const element_key& b) -> bool {
+                return a.first < b;
+            };
+            /* and lookup it in the vector we created before. */
+            auto itor = std::lower_bound(bfk.begin(), bfk.end(), fk, lbcomp);
+            if ( (itor == bfk.end()) or ((*itor).first != fk) )
+                continue;
+            
+            /* If found */
+            boundary_descriptor bd;
+            /* determine if it is an interface or a boundary */
+            if (conn.num_neighbours(fk) == 1)
+                bd.type = face_type::BOUNDARY;
+            else
+                bd.type = face_type::INTERFACE;
+            /* and store the gmsh tag in the descriptor. */
+            bd.gmsh_entity = (*itor).second;
+
+            /* Finally, put the descriptor in the global array of faces. */
+            bnd_descriptors.at(fbase + iF) = bd;
+        }
+        fbase += e.num_faces();
+    }
+
+    /*
+    size_t normal = 0;
+    size_t interface = 0;
+    size_t boundary = 0;
+    for (auto& bd : boundary_descriptors)
+    {
+        switch (bd.type)
+        {
+            case face_type::NONE: normal++; break;
+            case face_type::INTERFACE: interface++; break;
+            case face_type::BOUNDARY: boundary++; break;
         }
     }
+
+    std::cout << bfk.size() << " " << boundary_descriptors.size() << std::endl;
+    std::cout << normal << " " << interface << " " << boundary << std::endl;
+    */
 }
\ No newline at end of file

src/maxwell_solver.cpp

 #include "gmsh.h"
 #include "silo_output.hpp"
 #include "maxwell_interface.h"
+#include "param_loader.h"
 
 static void
 make_geometry(int order, double mesh_h)
@@ -29,34 +30,37 @@ make_geometry(int order, double mesh_h)
 }
 
 template<typename State>
-void test_it(const model& mod, State& state)
+void test_it(const model& mod, State& state, const maxwell_parameter_loader& mpl)
 {
-    auto E = [](const point_3d& pt) -> vec3d {
-        vec3d ret;
-        ret(0) = 0;
-        ret(1) = 0;
-        ret(2) = std::sin(M_PI*pt.x()) * std::sin(M_PI*pt.y());
-        return ret;
-    };
-
     maxwell::init_from_model(mod, state);
-    maxwell::init_E_field(mod, state, E);
+    
+    if ( mpl.initial_Efield_defined() )
+    {
+        auto E = [&](const point_3d& pt) -> vec3d { return mpl.Efield(pt); };
+        maxwell::init_E_field(mod, state, E);
+    }
+
+    if ( mpl.initial_Hfield_defined() )
+    {
+        auto H = [&](const point_3d& pt) -> vec3d { return mpl.Hfield(pt); };
+        maxwell::init_H_field(mod, state, H);
+    }
 
-    state.delta_t = 0.001;
+    state.delta_t = mpl.delta_t();
+    auto num_timesteps = mpl.timesteps();
+    auto silo_output_rate = mpl.silo_output_rate();
 
-    for(size_t i = 0; i < 10000; i++)
+    for(size_t i = 0; i < num_timesteps; i++)
     {
         timestep(state);
         std::stringstream ss;
         ss << "maxwell_" << i << ".silo";
 
-        if (i%10 == 0)
+        if (i%silo_output_rate == 0)
             maxwell::export_to_silo(mod, state, ss.str());
         
         std::swap(state.emf_curr, state.emf_next);
     }
-
-    
 }
 
 int main(int argc, const char *argv[])
@@ -64,13 +68,30 @@ int main(int argc, const char *argv[])
     gmsh::initialize();
     make_geometry(1, 0.1);
 
-    model mod(1,3);
-    
-    maxwell::solver_state state_c;
-    test_it(mod, state_c);
+    maxwell_parameter_loader mpl;
+    if ( not mpl.load_file("params.lua") )
+    {
+        std::cout << "Configuration problem, exiting" << std::endl;
+        return 1;
+    }
+
+    model mod( mpl.geom_order(), mpl.approx_order() );
+    if ( not mpl.validate_model_params(mod) )
+    {
+        std::cout << "Configuration problem, exiting" << std::endl;
+        return 1;
+    }
 
-    //maxwell::solver_state_gpu state_g;
-    //test_it(mod, state_g);
+    if ( mpl.use_gpu() )
+    {
+        maxwell::solver_state_gpu state_g;
+        test_it(mod, state_g, mpl);
+    }
+    else
+    {
+        maxwell::solver_state state_c;
+        test_it(mod, state_c, mpl);
+    }
 
     //gmsh::finalize();
 

src/param_loader.cpp

+#include <iostream>
+#include "param_loader.h"
+
+parameter_loader::parameter_loader()
+{
+    lua.open_libraries(sol::lib::base, sol::lib::math);
+    init();
+}
+
+void
+parameter_loader::init(void)
+{
+    lua["const"] = lua.create_table();
+    lua["const"]["eps0"] = 8.8541878128e-12;
+    lua["const"]["mu0"] = 4e-7*M_PI;
+
+    lua["sim"] = lua.create_table();
+    lua["sim"]["approx_order"] = 1;
+    lua["sim"]["geom_order"] = 1;
+    lua["sim"]["use_gpu"] = 0;
+    //lua["sim"]["name"];
+    //lua["sim"]["dt"];
+    //lua["sim"]["timesteps"];
+    //lua["sim"]["gmsh_model"];
+
+    lua["materials"] = lua.create_table();
+    lua["bndconds"] = lua.create_table();
+    lua["ifaceconds"] = lua.create_table();
+
+    lua["postpro"] = lua.create_table();
+    lua["postpro"]["silo_output_rate"] = 0;
+}
+
+bool
+parameter_loader::validate_simulation_params(void) const
+{
+    bool success = true;
+
+    auto sim_name = lua["sim"]["name"];
+    if (not sim_name.valid())
+    {
+        std::cout << "[CONFIG] 'sim.name' not set" << std::endl;
+        success = false;
+    } 
+
+    auto sim_dt = lua["sim"]["dt"];
+    if (not sim_dt.valid())
+    {
+        std::cout << "[CONFIG] 'sim.dt' not set" << std::endl;
+        success = false;
+    } 
+
+    auto sim_timesteps = lua["sim"]["timesteps"];
+    if (not sim_timesteps.valid())
+    {
+        std::cout << "[CONFIG] 'sim.timesteps' not set" << std::endl;
+        success = false;
+    } 
+
+    auto sim_gmshmodel = lua["sim"]["gmsh_model"];
+    if (not sim_gmshmodel.valid())
+    {
+        std::cout << "[CONFIG] 'sim.gmsh_model' not set" << std::endl;
+        success = false;
+    } 
+
+    return success;
+}
+
+bool
+parameter_loader::load_file(const std::string& fn)
+{
+    auto script = lua.script_file(fn);
+    if (not script.valid())
+        return false;
+
+    return validate_simulation_params();
+}
+
+int
+parameter_loader::approx_order(void) const
+{
+    return lua["sim"]["approx_order"];
+}
+
+int
+parameter_loader::geom_order(void) const
+{
+    return lua["sim"]["geom_order"];
+}
+
+double
+parameter_loader::delta_t(void) const
+{
+    return lua["sim"]["dt"];
+}
+
+size_t
+parameter_loader::timesteps(void) const
+{
+    return lua["sim"]["timesteps"];
+}
+
+std::string
+parameter_loader::simulation_name(void) const
+{
+    return lua["sim"]["name"];
+}
+
+bool
+parameter_loader::use_gpu(void) const
+{
+    return lua["sim"]["use_gpu"];
+}
+
+size_t
+parameter_loader::silo_output_rate(void) const
+{
+    return lua["postpro"]["silo_output_rate"];
+}
+
+maxwell_parameter_loader::maxwell_parameter_loader()
+{}
+
+bool
+maxwell_parameter_loader::validate_materials(const std::string& mpn, int tag)
+{
+    auto mfun = lua["materials"][mpn];
+    if ( mfun.valid() )
+        return true;
+
+    auto mparams = lua["materials"][tag];
+    if ( mparams.valid() )
+    {
+        auto matparams_mpn = lua["materials"][tag][mpn];
+        if ( matparams_mpn.valid() )
+            return true;
+    }
+
+    std::cout << "[CONFIG] 'materials[" << tag << "]." << mpn << "' not defined and ";
+    std::cout << "'materials." << mpn << "(tag,x,y,z)' not present." << std::endl;
+    return false;
+}
+
+bool
+maxwell_parameter_loader::validate_model_params(const model& mod)
+{
+    bool success = true;
+    for (auto& e : mod)
+    {
+        auto tag = e.gmsh_tag();
+        bool eps_valid = validate_materials("epsilon", tag);
+        bool mu_valid = validate_materials("mu", tag);
+        bool sigma_valid = validate_materials("sigma", tag);
+        if ( (not eps_valid) or (not mu_valid) or (not sigma_valid) )
+            success = false;
+    }
+
+    return success;
+}
+
+bool
+maxwell_parameter_loader::initial_Efield_defined(void) const
+{
+    auto eic = lua["electric_initial_condition"];
+    return eic.valid();
+}
+
+vec3d
+maxwell_parameter_loader::Efield(const point_3d& pt) const
+{
+    vec3d ret;
+    sol::tie(ret(0), ret(1), ret(2)) =
+        lua["electric_initial_condition"](pt.x(), pt.y(), pt.z());
+    return ret;
+}
+
+bool
+maxwell_parameter_loader::initial_Hfield_defined(void) const
+{
+    auto mic = lua["magnetic_initial_condition"];
+    return mic.valid();
+}
+
+vec3d
+maxwell_parameter_loader::Hfield(const point_3d& pt) const
+{
+    vec3d ret;
+    sol::tie(ret(0), ret(1), ret(2)) =
+        lua["magnetic_initial_condition"](pt.x(), pt.y(), pt.z());
+    return ret;
+}
\ No newline at end of file

src/physical_element.cpp

@@ -312,6 +312,9 @@ element_key_factory::element_key_factory(int dim, int tag, int etype)
         ekeys.push_back( ek );
     }
     std::sort(ekeys.begin(), ekeys.end());
+
+    /* Not sure why there are duplicates */
+    ekeys.erase( std::unique(ekeys.begin(), ekeys.end()), ekeys.end() );
 }
 
 std::vector<element_key>