Merge branch 'MC/gefpg_lua_binding' into devel

CMakeLists.txt

@@ -43,48 +43,58 @@ if (CMAKE_CXX_COMPILER_ID STREQUAL "AppleClang")
 endif()
 
 ######################################################################
-## Find required libraries
+## Find required libraries: Eigen3
 find_package(Eigen3 REQUIRED)
 set(LINK_LIBS ${LINK_LIBS} Eigen3::Eigen)
 
-find_library(GMSH_LIB gmsh)
-find_path(GMSH_INC gmsh.h)
-if(GMSH_LIB AND GMSH_INC)
-    set(LINK_LIBS ${LINK_LIBS} ${GMSH_LIB})
-    include_directories(${GMSH_INC})
-else()
-  set(GMSH_DIR "/opt/uliege/gmsh" CACHE PATH "GMSH install path")
-  if (DEFINED ENV{GMSH_DIR})
-    set (GMSH_DIR $ENV{GMSH_DIR})
-  endif()
-  if (GMSH_DIR STREQUAL "" OR NOT EXISTS ${GMSH_DIR})
-      message(FATAL_ERROR "GMSH not found")
-  else()
-      set(LINK_LIBS ${LINK_LIBS} gmsh)
-      set(GMSH_INC "${GMSH_DIR}/include")
-      set(GMSH_LIB "${GMSH_DIR}/lib")
-      include_directories(${GMSH_INC})
-      link_directories(${GMSH_LIB})
-  endif()
-endif()
-if(GMSH_INC)
-  file(STRINGS ${GMSH_INC}/gmsh.h GMSH_API_VERSION_STR REGEX "#define GMSH_API_VERSION .*")
+######################################################################
+## Find required libraries: GMSH
+find_package(GMSH REQUIRED)
+set(LINK_LIBS ${LINK_LIBS} GMSH)
+
+file(STRINGS ${GMSH_INCLUDE_DIR}/gmsh.h GMSH_API_VERSION_STR REGEX "#define GMSH_API_VERSION .*")
 string(REGEX MATCH "[\.0-9]+" GMSH_API_VERSION "${GMSH_API_VERSION_STR}")
 if(GMSH_API_VERSION VERSION_LESS "4.9.0")
     add_definitions(-DUSE_INITIAL_4_8_4_API)
 endif()
-endif()
 
+#find_library(GMSH_LIB gmsh)
+#find_path(GMSH_INC gmsh.h)
+#if(GMSH_LIB AND GMSH_INC)
+#    set(LINK_LIBS ${LINK_LIBS} ${GMSH_LIB})
+#    include_directories(${GMSH_INC})
+#else()
+#  set(GMSH_DIR "/opt/uliege/gmsh" CACHE PATH "GMSH install path")
+#  if (DEFINED ENV{GMSH_DIR})
+#    set (GMSH_DIR $ENV{GMSH_DIR})
+#  endif()
+#  if (GMSH_DIR STREQUAL "" OR NOT EXISTS ${GMSH_DIR})
+#      message(FATAL_ERROR "GMSH not found")
+#  else()
+#      set(LINK_LIBS ${LINK_LIBS} gmsh)
+#      set(GMSH_INC "${GMSH_DIR}/include")
+#      set(GMSH_LIB "${GMSH_DIR}/lib")
+#      include_directories(${GMSH_INC})
+#      link_directories(${GMSH_LIB})
+#  endif()
+#endif()
+#if(GMSH_INC)
+#  file(STRINGS ${GMSH_INC}/gmsh.h GMSH_API_VERSION_STR REGEX "#define GMSH_API_VERSION .*")
+#  string(REGEX MATCH "[\.0-9]+" GMSH_API_VERSION "${GMSH_API_VERSION_STR}")
+#  if(GMSH_API_VERSION VERSION_LESS "4.9.0")
+#    add_definitions(-DUSE_INITIAL_4_8_4_API)
+#  endif()
+#endif()
+
+######################################################################
+## Find required libraries: SILO
 find_package(SILO REQUIRED)
-if (SILO_FOUND)
-    add_definitions(-DHAVE_SILO)
-    include_directories("${SILO_INCLUDE_DIRS}")
-    set(LINK_LIBS ${LINK_LIBS} ${SILO_LIBRARIES})
+set(LINK_LIBS ${LINK_LIBS} SILO)
+
 option(OPT_SILO_USE_HDF5 "Use HDF5 driver for silo files" ON)
 if (OPT_SILO_USE_HDF5)
     add_definitions(-DSILO_USE_HDF5)
 endif()
-endif()
 
 find_package(Lua REQUIRED)
 include_directories("${LUA_INCLUDE_DIR}")

cmake/FindGMSH.cmake

+include(FindPackageHandleStandardArgs)
+
+find_path(GMSH_INCLUDE_DIR
+    NAMES gmsh.h
+    HINTS ENV GMSH_ROOT ${GMSH_ROOT}
+    PATH_SUFFIXES include)
+
+find_library(GMSH_LIBRARIES
+    NAMES gmsh
+    HINTS ENV GMSH_ROOT ${GMSH_ROOT}
+    PATH_SUFFIXES lib lib64)
+
+find_package_handle_standard_args(GMSH DEFAULT_MSG GMSH_LIBRARIES GMSH_INCLUDE_DIR)
+
+if (GMSH_FOUND)
+    add_library(GMSH INTERFACE)
+    target_link_libraries(GMSH INTERFACE "${GMSH_LIBRARIES}")
+    target_include_directories(GMSH INTERFACE "${GMSH_INCLUDE_DIR}")
+endif()

cmake/FindSILO.cmake

 include(FindPackageHandleStandardArgs)
 
-
-find_path(SILO_INCLUDE_DIRS
+find_path(SILO_INCLUDE_DIR
     NAMES silo.h
-          PATHS /usr/include /usr/local/include ${SILO_DIR}/include)
+    HINTS ENV SILO_ROOT ${SILO_ROOT}
+    PATH_SUFFIXES include)
 
 find_library(SILO_LIBRARIES
     NAMES silo siloh5
-             PATHS /usr/lib /usr/local/include ${SILO_DIR}/lib)
+    HINTS ENV SILO_ROOT ${SILO_ROOT}
+    PATH_SUFFIXES lib)
+
+find_package_handle_standard_args(SILO DEFAULT_MSG SILO_LIBRARIES SILO_INCLUDE_DIR)
 
-FIND_PACKAGE_HANDLE_STANDARD_ARGS(SILO DEFAULT_MSG SILO_LIBRARIES SILO_INCLUDE_DIRS)
+if (SILO_FOUND)
+    add_library(SILO INTERFACE)
+    target_link_libraries(SILO INTERFACE ${SILO_LIBRARIES})
+    target_include_directories(SILO INTERFACE "${SILO_INCLUDE_DIR}")
+    target_compile_definitions(SILO INTERFACE -DHAVE_SILO)
+endif()
\ No newline at end of file

doc/changelog.md

@@ -29,3 +29,11 @@ Unavailable features:
 - MPI support for the CPU solver
 - Validation API on Lua configuration
 
+## Version v0.4
+- Don't regenerate mesh if input is already `.msh` (7ba32158)
+- Code to allow comparison with GMSH-FEM (a5068e3a)
+- Licensed as AGPLv3 (2c471367)
+- Eigen and Sol symbols hidden by linker script (02ac8c5d)
+- Support GMSH physical groups (be6606d7)
+- Fix performance issues in the evaluation of boundary sources (a247f57e)
+- Don't alloc memory for RK4 auxiliary vectors if RK4 not used (cc387be4)

doc/lua_api.md

 # GMSH DG solver Lua API
 
-The solver employs the Lua programming language for its **configuration**. Lua was chosen first of all for simplicity: it is extremely lightweight and carries almost no dependencies, and this allow to keep the solver small and compact. Secondly, Lua was chosen to *deliberately* limit the possibilities of what the user can do in the configuration files. If configurations become full-fledged programs, it means that something is missing in the solver core or that the solver is being used in the wrong way.
+The solver employs the Lua programming language for its **configuration**. Lua was chosen first of all for simplicity: it is extremely lightweight and carries almost no dependencies, and this allow to keep the solver small and compact. Secondly, Lua was chosen to *deliberately* limit the possibilities of what the user can do in the configuration files. If configurations become full-fledged programs, it means that something is missing in the solver core or that the solver is being used in the wrong way. Third, Lua is *fast*: evaluating user-defined Lua functions is not much slower than native C.
 
 This file documents the Lua API available in the solver. API has a *general* part and a *problem-specific* part.
 
 The *general* part has to do with configuration not related to a specific problem (i.e. timestep, geometry file), whereas the *problem-specific* part configures all the parameter that make sense only on a given problem (i.e. materials and sources). This separation is reflected also in how the configuration is handled internally.
 
 ## API version
-This document describes the API available on the version v0.3 of the solver.
+This document describes the API available on the version v0.4 of the solver.
 
 ## General interface
 
@@ -36,10 +36,15 @@ The parallel solver runs as separate MPI processes. As such, each process loads
 - `parallel.comm_size` (integer): The MPI communicator size
 
 ### Callable functions
+#### Control of the sources
 - `enable_boundary_sources()`: takes a boolean parameter that specifies if boundary sources should be enabled or not
 - `enable_interface_sources()`: takes a boolean parameter that specifies if the sources applied on internal interfaces should be enabled or not
 - `enable_volume_sources()`: takes a boolean parameter that specifies if volumetric sources should be enabled or not
 
+#### Handling of GMSH physical groups
+- `volume_physical_group_entities(pg)`: takes the tag of a volume physical group (`dim == 3` in GMSH parlance) and returns the tags of all the entities included in the physical group. This function is available only after mesh loading, therefore it is usually called inside some callback.
+- `surface_physical_group_entities(pg)`: takes the tag of a surface physical group (`dim == 2` in GMSH parlance) and returns the tags of all the entities included in the physical group. This function is available only after mesh loading, therefore it is usually called inside some callback.
+
 ### Callbacks
 - `on_exit()`: this callback is called at the exit of the program, just before the internal Lua context gets destroyed.
 
@@ -83,10 +88,10 @@ Current sources in volumes can be specified by providing a function that evaluat
 
 ```
 function J_source(tag, x, y, z, t)
-    local Ex = 0
-    local Ey = 0
-    local Ez = math.sin(2*math.pi*freq*t)
-    return Ex, Ey, Ez
+    local Jx = 0
+    local Jy = 0
+    local Jz = math.sin(2*math.pi*freq*t)
+    return Jx, Jy, Jz
 end
 
 sources[4] = J_source
@@ -141,6 +146,7 @@ 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.
+- `on_mesh_loaded()`: if defined, the solver calls this function just after loading the mesh.
 - `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`.
 
@@ -190,4 +196,4 @@ end
 ```
 computes the energy at each timestep and places it in the variable `e`. This variable provides 6 members named `Ex`, `Ey`, `Ez`, `Hx`, `Hy` and `Hz` so, for example, `e.Ex` gives access to the error of the `x` component of the electric field. 
 ### Considerations on the evaluation of sources and the computation on GPU
-Currently sources for the timestep `t+1` are evaluated asynchronously on the CPU while the GPU is computing the timestep `t`. The sources at `t+1` are then uploaded asynchronously from the CPU to the GPU. This usually works relatively well, especially on big domains with high approximation order, however on small domains at low approximation order you can notice slowdowns. For this reason, if the sources do not need to be evaluated along the whole simulation, it is suggested to turn them off at the appropriate timestep by using the `on_timestep()` callback and the appropriate functions described below. Despite of that, I am not sure that implementing repetitive boundary sources will be necessary.
+Currently sources for the timestep `t+1` are evaluated asynchronously on the CPU while the GPU is computing the timestep `t`. The sources at `t+1` are then uploaded asynchronously from the CPU to the GPU. This usually works relatively well, especially on big domains with high approximation order, however on small domains at low approximation order you can notice slowdowns, especially in the GPU code. For this reason, if the sources do not need to be evaluated along the whole simulation, it is suggested to turn them off at the appropriate timestep by using the `on_timestep()` callback and the appropriate functions described above.

include/libgmshdg/connectivity.h

@@ -197,7 +197,7 @@ public:
     }
 
     size_t
-    num_neighbours(const FaceKey& fk)
+    num_neighbours(const FaceKey& fk) const
     {
         size_t nn = face_cell_adj.at(fk).num_neighbours();
         assert(nn == 1 or nn == 2);

include/libgmshdg/gmsh_io.h

@@ -83,6 +83,18 @@ struct boundary_descriptor
     }
 };
 
+struct physical_group_descriptor {
+    int                 dim;
+    int                 tag;
+    std::vector<int>    entity_tags;
+
+    bool
+    operator<(const physical_group_descriptor& other) const
+    {
+        return (dim < other.dim) or ( (dim == other.dim) and (tag < other.tag) );
+    }
+};
+
 #ifdef USE_MPI
 struct comm_descriptor
 {
@@ -106,6 +118,75 @@ std::ostream& operator<<(std::ostream& os, const comm_descriptor& cd)
 
 #define IMPLIES(a, b) ((not(a)) or (b))
 #define IFF(a,b) (IMPLIES(a,b) and IMPLIES(b,a))
+struct surface_identifier
+{
+    int     surface_entity;
+    int     adjacent_volume_entity;
+
+    bool operator<(const surface_identifier& other) const
+    {
+        return (surface_entity < other.surface_entity) or
+                    ((surface_entity == other.surface_entity) and
+                    (adjacent_volume_entity < other.adjacent_volume_entity));
+    }
+};
+
+inline std::ostream&
+operator<<(std::ostream& os, const surface_identifier& si)
+{
+    os << "(" << si.surface_entity << ", " << si.adjacent_volume_entity << ")";
+    return os;
+}
+
+struct surface_descriptor
+{
+    struct face_number
+    {
+        size_t  local;
+        size_t  global;
+    };
+
+    int     surface_entity;
+    int     adjacent_volume_entity;
+#ifdef USE_MPI
+    int      parent_entity;
+#endif /* USE_MPI */
+
+    std::vector<face_number> face_numbers;
+
+    void add_face(size_t local_num, size_t global_num)
+    {
+        face_number fn;
+        fn.local = local_num;
+        fn.global = global_num;
+        face_numbers.push_back( fn );
+    }
+
+    void sort(void)
+    {
+        auto comp = [](const face_number& a, const face_number& b) {
+            return a.global < b.global;
+        };
+
+        std::sort(face_numbers.begin(), face_numbers.end(), comp);
+    }
+
+    int material_tag(void) const
+    {
+#ifdef USE_MPI
+        if (parent_entity != -1)
+            return parent_entity;
+#endif
+        return surface_entity;
+    }
+
+    bool operator<(const surface_descriptor& other) const
+    {
+        return (surface_entity < other.surface_entity) or
+                    ((surface_entity == other.surface_entity) and
+                    (adjacent_volume_entity < other.adjacent_volume_entity));
+    }
+};
 
 class model
 {
@@ -128,6 +209,8 @@ class model
     int                                         num_partitions;
 #endif /* USE_MPI */
 
+    std::vector<physical_group_descriptor>      physical_groups;
+
     /* Map from boundary tag to all the faces of the entity with that tag */
     std::map<int, std::vector<element_key>>     boundary_map;  
 
@@ -135,6 +218,9 @@ class model
     std::vector<boundary_descriptor>            bnd_descriptors;
     element_connectivity<element_key>           conn;               
 
+    std::map<surface_identifier, surface_descriptor>     boundaries;
+    std::map<surface_identifier, surface_descriptor>     interfaces;
+
 #ifdef USE_MPI
     interprocess_connectivity<element_key>      ipconn;
 
@@ -172,7 +258,7 @@ class model
     void    populate_from_gmsh_rank0(void);
 
 #ifdef USE_MPI
-    bool    is_interprocess_boundary(int);
+    bool    is_interprocess_boundary(int) const;
     void    populate_from_gmsh_rankN(void);
     void    import_gmsh_entities_rankN(void);
     int     my_partition(void) const;
@@ -183,6 +269,9 @@ class model
 #endif /* USE_MPI */
 
     void    make_boundary_to_faces_map(void);
+    void    import_physical_groups(void);
+
+    face_type   face_type(size_t, const element_key&) const;
 
 public:
     model();
@@ -196,6 +285,16 @@ public:
     void partition_mesh(int);
     void populate_from_gmsh(void);
 
+    std::vector<int> lookup_physical_group_entities(int dim, int tag);
+
+    void map_boundaries_new(void);
+
+    const std::map<surface_identifier, surface_descriptor>&
+    get_boundaries() const;
+
+    const std::map<surface_identifier, surface_descriptor>&
+    get_interfaces() const;
+
 #ifdef USE_MPI
     const std::map<int, comm_descriptor>& comm_descriptors(void) const {
         return ipc_boundary_comm_table;

include/maxwell/maxwell_common.h

@@ -102,6 +102,99 @@ eval_boundary_sources(const model& mod, const parameter_loader& mpl,
     }
 }
 
+template<typename FunEH>
+matxd
+project_EH_to_face(const entity& e, size_t iF, const FunEH& funEH)
+{
+    auto& pf = e.face(iF);
+    const auto& rf = e.face_refelem(pf);
+    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);
+    matxd ret = matxd::Zero(num_bf, 6);
+
+    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();
+        auto [fE, fH] = funEH(pqp.point());
+        ret.col(0) += pqp.weight() * fE(0) * phi; 
+        ret.col(1) += pqp.weight() * fE(1) * phi; 
+        ret.col(2) += pqp.weight() * fE(2) * phi;
+        ret.col(3) += pqp.weight() * fH(0) * phi; 
+        ret.col(4) += pqp.weight() * fH(1) * phi; 
+        ret.col(5) += pqp.weight() * fH(2) * phi; 
+    }
+
+    Eigen::LDLT<matxd> mass_ldlt(mass);
+    ret = mass_ldlt.solve(ret);
+    return ret;
+}
+
+/* Field must be either 'field' or 'pinned_field'. Too early to use a concept. */
+template<typename State, typename Field>
+void
+eval_boundary_sources_new(const model& mod, const parameter_loader& mpl,
+    State& state, Field& srcfield)
+{
+    auto& bds = mod.get_boundaries();
+    assert( srcfield.num_dofs == mod.num_fluxes() );
+
+    for (auto& [si, sd] : bds)
+    {
+        const auto& e = mod.entity_at(sd.adjacent_volume_entity);
+        auto etag = e.material_tag();
+
+        auto mtag = sd.material_tag();
+
+        switch ( mpl.boundary_type( mtag ) )
+        {
+            case boundary_condition::UNSPECIFIED:
+            case boundary_condition::PEC:
+            case boundary_condition::PMC:
+            case boundary_condition::E_FIELD:
+            case boundary_condition::H_FIELD:
+            case boundary_condition::IMPEDANCE:
+            case boundary_condition::PLANE_WAVE_H:
+            case boundary_condition::SURFACE_CURRENT:
+                    break;
+
+            case boundary_condition::PLANE_WAVE_E: {
+                for (auto& [local, global] : sd.face_numbers)
+                {
+                    auto& pf = e.face(local);
+                    auto bar = pf.barycenter();
+                    auto eps = mpl.epsilon(etag, bar);
+                    auto mu  = mpl.mu(etag, bar);
+                    auto Z = std::sqrt(mu/eps);
+                    vec3d n = state.eds[e.number()].normals.row(local);
+                    
+                    auto fEH = [&](const point_3d& pt) -> std::pair<vec3d, vec3d> {
+                        vec3d E = mpl.eval_boundary_source(mtag, pt, state.curr_time);
+                        vec3d H = n.cross(E)/Z;
+                        return std::make_pair(E, H);
+                    };
+
+                    matxd prj = project_EH_to_face(e, local, fEH);
+                    auto num_bf = prj.rows();
+                    srcfield.Ex.segment(global*num_bf, num_bf) = prj.col(0);
+                    srcfield.Ey.segment(global*num_bf, num_bf) = prj.col(1);
+                    srcfield.Ez.segment(global*num_bf, num_bf) = prj.col(2);
+                    srcfield.Hx.segment(global*num_bf, num_bf) = prj.col(3);
+                    srcfield.Hy.segment(global*num_bf, num_bf) = prj.col(4);
+                    srcfield.Hz.segment(global*num_bf, num_bf) = prj.col(5);
+                }
+            } break; /* case boundary_condition::PLANE_WAVE_E */
+        }
+    }
+}
+
 /* Field must be either 'field' or 'pinned_field'. Too early to use a concept. */
 template<typename State, typename Field>
 void

include/maxwell/maxwell_interface.h

@@ -337,7 +337,7 @@ struct solver_state_gpu
 #endif /* ENABLE_GPU_SOLVER */
 
 #ifndef HIDE_THIS_FROM_NVCC
-void init_from_model(const model&, solver_state&);
+void init_from_model(const model&, solver_state&, const parameter_loader&);
 void init_matparams(const model&, solver_state&, const parameter_loader&);
 //void apply_operator(solver_state&, const field&, field&);
 void export_fields_to_silo(const model&, solver_state&, const parameter_loader&);
@@ -351,7 +351,7 @@ void gather_field(const model&, maxwell::field&, int, MPI_Comm);
 #endif /* USE_MPI */
 
 #ifdef ENABLE_GPU_SOLVER
-void init_from_model(const model&, solver_state_gpu&);
+void init_from_model(const model&, solver_state_gpu&, const parameter_loader&);
 void init_matparams(const model&, solver_state_gpu&, const parameter_loader&);
 //void apply_operator(solver_state_gpu&, const field_gpu&, field_gpu&);
 void export_fields_to_silo(const model&, solver_state_gpu&, const parameter_loader&);

share/validation/params_twomat.lua

+-- Validation test: simulate a plane wave propagating in a parallel plate
+-- waveguide with a discontinuity in materials.
+
+sim.name = "twomat"   -- simulation name
+sim.dt = 1e-12                                  -- timestep size
+sim.timesteps = 10000                             -- num of iterations
+sim.gmsh_model = "twomat.geo"           -- gmsh model filename
+sim.use_gpu = 0                                 -- 0: cpu, 1: gpu
+sim.approx_order = 1                            -- approximation order
+sim.time_integrator = "leapfrog"
+postpro.silo_output_rate = 100
+postpro.cycle_print_rate = 100                  -- console print rate
+
+postpro["E"].mode = "nodal"
+
+materials[1] = {}
+materials[1].epsilon = 1
+materials[1].mu = 1
+materials[1].sigma = 0
+
+materials[2] = {}
+materials[2].epsilon = 4
+materials[2].mu = 1
+materials[2].sigma = 0
+
+-- ** Boundary conditions **
+local PMC_bnds = { 3, 4, 8, 9 }
+for i,v in ipairs(PMC_bnds) do
+    bndconds[v] = {}
+    bndconds[v].kind = "pmc"
+end
+
+bndconds[7] = {}
+bndconds[7].kind = "impedance"
+
+local freq = 3e8
+function source(tag, x, y, z, t)
+    local Ex = 0
+    local Ey = 0
+    local Ez = math.sin(2*math.pi*t*freq)
+    return Ex, Ey, Ez
+end
+
+bndconds[1] = {}
+bndconds[1].kind = "plane_wave_E"
+bndconds[1].source = source
+

share/validation/twomat.geo

+SetFactory("OpenCASCADE");
+Box(1) = {0, 0, 0, 1, 1, 0.1};
+Box(2) = {1, 0, 0, 1, 1, 0.1};
+Coherence;
+
+MeshSize{:} = 0.075;

src/libgmshdg/gmsh_io.cpp

@@ -386,7 +386,7 @@ model::import_gmsh_entities_rankN(void)
 }
 
 bool
-model::is_interprocess_boundary(int tag)
+model::is_interprocess_boundary(int tag) const
 {
     if (num_partitions > 1)
         return comm_map.find(tag) != comm_map.end();
@@ -396,6 +396,55 @@ model::is_interprocess_boundary(int tag)
 
 #endif /* USE_MPI */
 
+void
+model::import_physical_groups(void)
+{
+#ifdef USE_MPI
+    int rank;
+    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+    if (rank == 0) {
+#endif /* USE_MPI */
+        gmsh::vectorpair pgs;
+        gm::getPhysicalGroups(pgs);
+        for (const auto& [dim, tag] : pgs)
+        {
+            std::vector<int> etags;
+            gm::getEntitiesForPhysicalGroup(dim, tag, etags);
+            struct physical_group_descriptor pgd;
+            pgd.dim = dim;
+            pgd.tag = tag;
+            pgd.entity_tags = std::move(etags);
+            physical_groups.push_back( std::move(pgd) );
+        }
+        std::sort(physical_groups.begin(), physical_groups.end());
+#ifdef USE_MPI
+        size_t pgsize = physical_groups.size();
+        MPI_Bcast(&pgsize, 1, MPI_UNSIGNED_LONG_LONG, 0, MPI_COMM_WORLD);
+        for (auto& pg : physical_groups)
+        {
+            MPI_Bcast(&pg.dim, 1, MPI_INT, 0, MPI_COMM_WORLD);
+            MPI_Bcast(&pg.dim, 1, MPI_INT, 0, MPI_COMM_WORLD);
+            size_t esize = pg.entity_tags.size();
+            MPI_Bcast(&esize, 1, MPI_UNSIGNED_LONG_LONG, 0, MPI_COMM_WORLD);
+            MPI_Bcast(pg.entity_tags.data(), esize, MPI_INT, 0, MPI_COMM_WORLD);
+        }
+    } else /* rank != 0 */ {
+        size_t pgsize;
+        MPI_Bcast(&pgsize, 1, MPI_UNSIGNED_LONG_LONG, 0, MPI_COMM_WORLD);
+        physical_groups.resize(pgsize);
+        for (auto& pg : physical_groups)
+        {
+            MPI_Bcast(&pg.dim, 1, MPI_INT, 0, MPI_COMM_WORLD);
+            MPI_Bcast(&pg.dim, 1, MPI_INT, 0, MPI_COMM_WORLD);
+            size_t esize;
+            MPI_Bcast(&esize, 1, MPI_UNSIGNED_LONG_LONG, 0, MPI_COMM_WORLD);
+            pg.entity_tags.resize(esize);
+            MPI_Bcast(pg.entity_tags.data(), esize, MPI_INT, 0, MPI_COMM_WORLD);
+        }
+    }
+#endif /* USE_MPI */
+}
+
 std::vector<model::bfk_t>
 model::get_facekey_tag_pairs(void)
 {
@@ -415,8 +464,24 @@ model::get_facekey_tag_pairs(void)
     return bfk;
 }
 
+face_type
+model::face_type(size_t entity, const element_key& fk) const
+{
+    if (conn.num_neighbours(fk) == 1)
+    {
+#ifdef USE_MPI
+        if ( is_interprocess_boundary(entity) )
+            return face_type::INTERPROCESS_BOUNDARY;
+        else
+#endif /* USE_MPI */
+            return face_type::BOUNDARY;
+    }
+    
+    return face_type::INTERFACE;
+}
+
 void
-model::map_boundaries(void)
+model::map_boundaries_new(void)
 {
     std::vector<bfk_t> bfk = get_facekey_tag_pairs();
 
@@ -441,15 +506,120 @@ model::map_boundaries(void)
             std::cout << " [MPI rank " << rank << "]"; 
 #endif
             std::cout << ". " << reset << std::endl;
-            /*
+        }
+    }
+
+    /* for each entity */
+    for (size_t iE = 0, fbase = 0; iE < num_entities(); iE++)
+    {
+        auto& e = entity_at(iE);
+        /* and for each face */
+        for (size_t iF = 0; iF < e.num_faces(); iF++)
+        {
+            /* 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 not ((*itor).first == fk) )
+                continue;
+            
+            surface_identifier si;
+            si.surface_entity = (*itor).second;
+            si.adjacent_volume_entity = iE;
+
+            using sm_t = std::map<surface_identifier, surface_descriptor>;
+
+            auto add_face = [&](sm_t& m, size_t se, size_t ave) -> void {
+                if ( m.find(si) == m.end() )
+                {
+                    m[si].surface_entity = se;
+                    m[si].adjacent_volume_entity = ave;
 #ifdef USE_MPI
+                    if (num_partitions > 1)
+                    {
+                        auto sitor = surface_to_parent_map.find( se );
+                        if ( sitor != surface_to_parent_map.end() )
+                            m[si].parent_entity = surface_to_parent_map.at( se );
+                    }
+#endif /* USE_MPI */
+                }
+                m[si].add_face(iF, fbase+iF);
+            };
+
+            switch ( face_type(si.surface_entity, fk) )
+            {
+                case face_type::BOUNDARY:
+                    add_face(boundaries, si.surface_entity, si.adjacent_volume_entity);
+                    break;
+
+                case face_type::INTERFACE:
+                    add_face(interfaces, si.surface_entity, si.adjacent_volume_entity);
+                    break;
+
+                default:
+                    break;
+            }
+        }
+
+        fbase += e.num_faces();
+    }
+
+    std::cout << "Boundaries:" << std::endl;
+    for (auto& [si, sd] : boundaries)
+    {
+        std::cout << si << " " << sd.face_numbers.size() << std::endl;
+        sd.sort();
+    }
+
+    std::cout << "Interfaces:" << std::endl;
+    for (auto& [si, sd] : interfaces)
+    {
+        std::cout << si << " " << sd.face_numbers.size() << std::endl;
+        sd.sort();
+    }
+}
+
+const std::map<surface_identifier, surface_descriptor>&
+model::get_boundaries() const
+{
+    return boundaries;
+}
+
+const std::map<surface_identifier, surface_descriptor>&
+model::get_interfaces() const
+{
+    return interfaces;
+}
+
+void
+model::map_boundaries(void)
+{
+    std::vector<bfk_t> bfk = get_facekey_tag_pairs();
+
+    for (size_t i = 1; i < bfk.size(); i++)
+    {
+        if ( bfk[i-1].first == bfk[i].first )
+        {
+            std::cout << redfg << Bon << "WARNING: " << nofg << __FILE__;
+            std::cout << "(" << __LINE__ << "): Face identified by ";
+            std::cout << bfk[i].first << " was found on interfaces ";
+            std::cout << bfk[i-1].second;
+#ifdef USE_MPI
+            if ( is_interprocess_boundary(bfk[i-1].second) )
+                std::cout << " (IP)";
+#endif
+            std::cout << " and " << bfk[i].second;
+#ifdef USE_MPI
+            if ( is_interprocess_boundary(bfk[i].second) )
+                std::cout << " (IP)";
             int rank;
             MPI_Comm_rank(MPI_COMM_WORLD, &rank);
-            if (rank == 0)
+            std::cout << " [MPI rank " << rank << "]"; 
 #endif
-                gmsh::write("crash.msh");
-            throw 42;
-            */
+            std::cout << ". " << reset << std::endl;
         }
     }
 
@@ -524,6 +694,8 @@ model::map_boundaries(void)
     std::cout << bfk.size() << " " << bnd_descriptors.size() << std::endl;
     std::cout << normal << " " << interface << " " << boundary << " " << ipc_boundary << std::endl;
 #endif /* ENABLE_DEBUG_PRINT */
+
+    map_boundaries_new();
 }
 
 #ifdef USE_MPI
@@ -793,6 +965,7 @@ model::populate_from_gmsh()
 #else /* USE_MPI */
     populate_from_gmsh_rank0();
 #endif /* USE_MPI */
+    import_physical_groups();
 }
 
 void
@@ -934,3 +1107,17 @@ model::lookup_tag(size_t tag) const
 {
     return etag_to_entity_offset.at(tag);
 }
+
+std::vector<int>
+model::lookup_physical_group_entities(int dim, int tag)
+{
+    physical_group_descriptor pg;
+    pg.dim = dim;
+    pg.tag = tag;
+    auto elem = std::lower_bound(physical_groups.begin() , physical_groups.end(), pg);
+    if (!(elem == physical_groups.end()) && !(pg < *physical_groups.begin()))
+        return elem->entity_tags;
+
+    return std::vector<int>{};
+}
+

src/libgmshdg/param_loader.cpp

@@ -8,12 +8,33 @@
 
 #include <iostream>
 
+#include "gmsh.h"
+
 #ifdef USE_MPI
 #include "common/mpi_helpers.h"
 #endif /* USE_MPI */
 
 #include "libgmshdg/param_loader.h"
 
+static auto
+lua_getEntitiesForPhysicalGroup(int dim, int tag)
+{
+    std::vector<int> ret;
+#ifdef USE_MPI
+    int rank;
+    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+    if (rank != 0)
+    {
+        std::cout << "[WARNING] getEntitiesForPhysicalGroup() ";
+        std::cout << "can be called only on MPI rank 0. On the other ";
+        std::cout << "ranks you will get an empty list." << std::endl;
+        return sol::as_table(ret);
+    }
+#endif
+    gmsh::model::getEntitiesForPhysicalGroup(dim, tag, ret);
+    return sol::as_table(ret);
+}
+
 parameter_loader_base::parameter_loader_base()
     : m_bnd_sources_active(true), m_ifc_sources_active(true),
       m_vol_sources_active(true)
@@ -65,6 +86,9 @@ parameter_loader_base::init(void)
                      &parameter_loader_base::enable_volume_sources,
                      this);
     
+    lua.set_function("getEntitiesForPhysicalGroup",
+        &lua_getEntitiesForPhysicalGroup);
+
 #ifdef USE_MPI
     int comm_rank, comm_size;
     MPI_Comm_rank(MPI_COMM_WORLD, &comm_rank);

src/maxwell/maxwell_cpu.cpp

@@ -61,7 +61,7 @@ init_matparams(const model& mod, solver_state& state,
 }
 
 void
-init_from_model(const model& mod, solver_state& state)
+init_from_model(const model& mod, solver_state& state, const maxwell::parameter_loader& mpl)
 {
     state.emf_curr.resize( mod.num_dofs() );
     state.emf_next.resize( mod.num_dofs() );
@@ -69,10 +69,14 @@ init_from_model(const model& mod, solver_state& state)
     state.dy.resize( mod.num_dofs() );
     state.dz.resize( mod.num_dofs() );
     
+    if ( mpl.sim_timeIntegrator() == time_integrator_type::RK4 )
+    {
         state.k1.resize( mod.num_dofs() );
         state.k2.resize( mod.num_dofs() );
         state.k3.resize( mod.num_dofs() );
         state.k4.resize( mod.num_dofs() );
+    }
+    
     state.tmp.resize( mod.num_dofs() );
 
     state.jumps.resize( mod.num_fluxes() );
@@ -777,7 +781,7 @@ do_sources(const model& mod, solver_state& state, parameter_loader& mpl)
         eval_volume_sources(mod, mpl, state);
 
     if ( mpl.boundary_sources_enabled() )
-        eval_boundary_sources(mod, mpl, state, state.bndsrcs);
+        eval_boundary_sources_new(mod, mpl, state, state.bndsrcs);
     
     if ( mpl.interface_sources_enabled() )
         eval_interface_sources(mod, mpl, state, state.bndsrcs);
@@ -854,7 +858,7 @@ export_fields_to_silo(const model &mod, maxwell::solver_state &state,
         size_t gdnum = mod.num_dofs_world();
         maxwell::field f;
         f.resize(gdnum);
-        gather_field(mod, state.emf_next, f, 0, MPI_COMM_WORLD);
+        gather_field(mod, state.emf_curr, f, 0, MPI_COMM_WORLD);
 
         if (basename == "")
             basename = "timestep";
@@ -897,6 +901,7 @@ void
 export_fields_to_silo(const model &mod, maxwell::solver_state &state,
                       const maxwell::parameter_loader &mpl, std::string basename)
 {
+    /* This is supposed to be called after swap */
     if (basename == "")
         basename = "timestep";
 
@@ -909,17 +914,17 @@ export_fields_to_silo(const model &mod, maxwell::solver_state &state,
     sdb.write_mesh(state.curr_time, state.curr_timestep);   
 
     auto E_exportmode = mpl.postpro_fieldExportMode("E");
-    export_vector_field(mod, sdb, state.emf_next.Ex, state.emf_next.Ey,
-                state.emf_next.Ez, "E", E_exportmode);
+    export_vector_field(mod, sdb, state.emf_curr.Ex, state.emf_curr.Ey,
+                state.emf_curr.Ez, "E", E_exportmode);
     
     auto H_exportmode = mpl.postpro_fieldExportMode("H");
-    export_vector_field(mod, sdb, state.emf_next.Hx, state.emf_next.Hy,
-                state.emf_next.Hz, "H", H_exportmode);
+    export_vector_field(mod, sdb, state.emf_curr.Hx, state.emf_curr.Hy,
+                state.emf_curr.Hz, "H", H_exportmode);
 
     auto sigma = [&](int tag, const point_3d& pt) -> double { return mpl.sigma(tag, pt); };
     auto J_exportmode = mpl.postpro_fieldExportMode("J");
-    export_vector_field(mod, sdb, state.emf_next.Ex, state.emf_next.Ey,
-                state.emf_next.Ez, sigma, "J", J_exportmode);
+    export_vector_field(mod, sdb, state.emf_curr.Ex, state.emf_curr.Ey,
+                state.emf_curr.Ez, sigma, "J", J_exportmode);
 }
 
 #endif /* USE_MPI */

src/maxwell/maxwell_gpu.cpp

@@ -70,7 +70,61 @@ init_matparams(const model& mod, solver_state_gpu& state,
     state.matparams.bc_coeffs.copyin(bcc.data(), bcc.size());
 }
 
-void init_from_model(const model& mod, solver_state_gpu& state)
+static void
+build_source_compression_tables(const model& mod, solver_state_gpu& state, const parameter_loader& mpl)
+{
+    auto& bds = mod.boundary_descriptors();
+    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& rf = e.face_refelem(pf);
+            auto num_fluxes = rf.num_basis_functions();
+
+            auto gfnum = face_num_base + iF;
+            auto bd = bds[gfnum];
+            if (bd.type == face_type::INTERFACE)
+            {
+                switch ( mpl.interface_type(bd.material_tag()) )
+                {
+                    case interface_condition::E_FIELD:
+                    case interface_condition::SURFACE_CURRENT:
+                        for (size_t i = 0; i < num_fluxes; i++)
+                            state.bndsrcs_decomp_table_cpu.push_back(gfnum*num_fluxes+i);
+                        break;
+                    
+                    default:
+                        break;
+                }
+            }
+
+            if (bd.type == face_type::BOUNDARY)
+            {
+                switch ( mpl.boundary_type(bd.material_tag()) )
+                {
+                    case boundary_condition::PLANE_WAVE_E:
+                        for (size_t i = 0; i < num_fluxes; i++)
+                            state.bndsrcs_decomp_table_cpu.push_back(gfnum*num_fluxes+i);
+                        break;
+                    
+                    default:
+                        break;
+                }
+            }
+        }
+        face_num_base += e.num_faces();
+    }
+
+    state.bndsrcs_decomp_table.copyin( state.bndsrcs_decomp_table_cpu.data(),
+        state.bndsrcs_decomp_table_cpu.size() );
+
+    state.bndsrcs_cpu.resize( state.bndsrcs_decomp_table_cpu.size() );
+    state.bndsrcs_buf.resize( state.bndsrcs_decomp_table_cpu.size() );
+}
+
+void init_from_model(const model& mod, solver_state_gpu& state, const maxwell::parameter_loader& mpl)
 {
     state.emf_curr.resize( mod.num_dofs() );
     state.emf_next.resize( mod.num_dofs() );
@@ -81,10 +135,14 @@ void init_from_model(const model& mod, solver_state_gpu& state)
     state.jumps.resize( mod.num_fluxes() );
     state.fluxes.resize( mod.num_fluxes() );
 
+    if ( mpl.sim_timeIntegrator() == time_integrator_type::RK4 )
+    {
         state.k1.resize( mod.num_dofs() );
         state.k2.resize( mod.num_dofs() );
         state.k3.resize( mod.num_dofs() );
         state.k4.resize( mod.num_dofs() );
+    }
+    
     state.tmp.resize( mod.num_dofs() );
 
     state.Jx_src.resize( mod.num_dofs() );
@@ -109,34 +167,7 @@ void init_from_model(const model& mod, solver_state_gpu& state)
         state.edgs.push_back( std::move(edg) );
     }
 
-    auto& bds = mod.boundary_descriptors();
-    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& rf = e.face_refelem(pf);
-            auto num_fluxes = rf.num_basis_functions();
-
-            auto gfnum = face_num_base + iF;
-            auto bd = bds[gfnum];
-            if (bd.type == face_type::BOUNDARY or bd.type == face_type::INTERFACE)
-            {
-                for (size_t i = 0; i < num_fluxes; i++)
-                {
-                    state.bndsrcs_decomp_table_cpu.push_back(gfnum*num_fluxes+i);
-                }
-            }
-        }
-        face_num_base += e.num_faces();
-    }
-
-    state.bndsrcs_decomp_table.copyin( state.bndsrcs_decomp_table_cpu.data(),
-        state.bndsrcs_decomp_table_cpu.size() );
-
-    state.bndsrcs_cpu.resize( state.bndsrcs_decomp_table_cpu.size() );
-    state.bndsrcs_buf.resize( state.bndsrcs_decomp_table_cpu.size() );
+    build_source_compression_tables(mod, state, mpl);
 
     state.curr_time = 0.0;
     state.curr_timestep = 0;
@@ -466,7 +497,7 @@ prepare_sources(const model& mod, maxwell::solver_state_gpu& state,
     maxwell::parameter_loader& mpl)
 {
     if ( mpl.boundary_sources_enabled() )
-        maxwell::eval_boundary_sources(mod, mpl, state, state.bndsrcs_decomp_cpu);
+        maxwell::eval_boundary_sources_new(mod, mpl, state, state.bndsrcs_decomp_cpu);
 
     if ( mpl.interface_sources_enabled() )
         maxwell::eval_interface_sources(mod, mpl, state, state.bndsrcs_decomp_cpu);
@@ -535,7 +566,7 @@ do_sources(const model& mod, maxwell::solver_state_gpu& state,
         eval_volume_sources(mod, mpl, state);
 
     if ( be )
-        eval_boundary_sources(mod, mpl, state, state.bndsrcs_decomp_cpu);
+        eval_boundary_sources_new(mod, mpl, state, state.bndsrcs_decomp_cpu);
 
     if ( ie )
         eval_interface_sources(mod, mpl, state, state.bndsrcs_decomp_cpu);
@@ -550,6 +581,7 @@ do_sources(const model& mod, maxwell::solver_state_gpu& state,
 void
 swap(solver_state_gpu& state, const parameter_loader& mpl)
 {
+    /* Wait for all async copies to finish */
     checkGPU( gpuDeviceSynchronize() );
 
     auto be = mpl.boundary_sources_enabled();
@@ -561,6 +593,9 @@ swap(solver_state_gpu& state, const parameter_loader& mpl)
     std::swap(state.Jy_src_gpu, state.Jy_src_gpu_buf);
     std::swap(state.Jz_src_gpu, state.Jz_src_gpu_buf);
     std::swap(state.emf_curr, state.emf_next);
+
+    /* Wait for decompression to finish */
+    checkGPU( gpuDeviceSynchronize() );
 }
 
 #ifdef USE_MPI
@@ -574,6 +609,7 @@ void
 export_fields_to_silo(const model& mod, maxwell::solver_state_gpu& state,
     const maxwell::parameter_loader& mpl, std::string basename)
 {
+    /* This is supposed to be called after swap */
     if (basename == "")
         basename = "timestep";
 
@@ -587,7 +623,7 @@ export_fields_to_silo(const model& mod, maxwell::solver_state_gpu& state,
 
     maxwell::field f;
     f.resize( mod.num_dofs() );
-    state.emf_next.copyout(f);
+    state.emf_curr.copyout(f);
 
 
     auto E_exportmode = mpl.postpro_fieldExportMode("E");

src/maxwell/maxwell_solver.cpp

@@ -41,7 +41,7 @@ using namespace sgr;
 template<typename State>
 void initialize_solver(const model& mod, State& state, const maxwell::parameter_loader& mpl)
 {
-    maxwell::init_from_model(mod, state);
+    maxwell::init_from_model(mod, state, mpl);
     
     if ( mpl.initial_Efield_defined() )
     {
@@ -245,6 +245,7 @@ void solver_mainloop(const model& mod, State& state, maxwell::parameter_loader&
         mpl.call_timestep_callback(i);
         timestep(mod, state, mpl, ti);
         do_sources(mod, state, mpl);
+        swap(state, mpl);
         
         double time = tc.toc();
         total_iteration_time += time;
@@ -252,8 +253,6 @@ void solver_mainloop(const model& mod, State& state, maxwell::parameter_loader&
         if ( (silo_output_rate != 0) and (i%silo_output_rate == 0) )
             export_fields_to_silo(mod, state, mpl, "");
 
-        swap(state, mpl);
-
         if ( (cycle_print_rate != 0) and (i%cycle_print_rate == 0) )
         {
             double dofs_s_proc = (mod.num_dofs()*6)/time;
@@ -328,6 +327,24 @@ register_lua_usertypes(maxwell::parameter_loader& mpl)
 #endif /* ENABLE_EXP_GEOMETRIC_DIODE */
 }
 
+/****************************************************************************/
+/* Register GMSH API stuff
+ */
+static void
+register_lua_callbacks(maxwell::parameter_loader& mpl, model& mod)
+{
+    sol::state& lua = mpl.lua_state();
+    auto volume_physical_group_entities = [&](int tag) -> auto {
+        return sol::as_table( mod.lookup_physical_group_entities(3, tag) );
+    };
+    lua["volume_physical_group_entities"] = volume_physical_group_entities;
+
+    auto surface_physical_group_entities = [&](int tag) -> auto {
+        return sol::as_table( mod.lookup_physical_group_entities(2, tag) );
+    };
+    lua["surface_physical_group_entities"] = surface_physical_group_entities;
+}
+
 /****************************************************************************/
 /* Register Lua usertypes valid only for CPU
  */
@@ -455,7 +472,12 @@ int main(int argc, char *argv[])
     mod.populate_from_gmsh();
 #endif /* USE_MPI */
     
+    register_lua_callbacks(mpl, mod);
 
+    sol::state& lua = mpl.lua_state();
+    auto on_mesh_loaded = lua["on_mesh_loaded"];
+    if (on_mesh_loaded.valid())
+        on_mesh_loaded();
 
     if ( not mpl.validate_model_params(mod) )
     {