Sponge test case, physical group handling improvements.

doc/changelog.md

@@ -34,7 +34,10 @@ Unavailable features:
 - Code to allow comparison with GMSH-FEM (a5068e3a)
 - Licensed as AGPLv3 (2c471367)
 - Eigen and Sol symbols hidden by linker script (Only Linux for now) (02ac8c5d)
-- Support GMSH physical groups (be6606d7, ad399724, 0dcb0e72)
+- Support GMSH physical groups (be6606d7, ad399724, 0dcb0e72 and others)
 - Fix performance issues in the evaluation of boundary sources (a247f57e)
 - Don't alloc memory for RK4 auxiliary vectors if RK4 not used (cc387be4)
 - New kernels for direct computation of curls (5b1986be and subsequent).
+- Multi-GPU support (93e017d)
+- Internal profiling infrastructure (c09d790 and subsequent)
+- New kernels to speed up lifting (fe86156)
\ No newline at end of file

doc/lua_api.md

@@ -45,8 +45,8 @@ The parallel solver runs as separate MPI processes. As such, each process loads
 - `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, most probably inside `on_mesh_loaded()` ([see below](#callable-functions-1)).
-- `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, most probably inside `on_mesh_loaded()` ([see below](#callable-functions-1)).
+- `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, most probably inside `on_mesh_loaded()` ([see below](#callable-functions-1)). This function is being replaced by `physical_group_to_volumes()`.
+- `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, most probably inside `on_mesh_loaded()` ([see below](#callable-functions-1)). This function is being replaced by `physical_group_to_surfaces()`.
 - `volume_to_physical_group(entity)`: takes the tag of a volumetric entity and returns the associated physical group. This function is available **only after mesh loading**, therefore it is usually called inside some callback, most probably inside `on_mesh_loaded()` ([see below](#callable-functions-1)). In addition, **exactly one physical group** must be associated to `entity`, otherwise an exception is thrown by the runtime.
 - `volume_to_physical_groups(entity)`: takes the tag of a volumetric entity and returns all the associated physical groups. This function is available **only after mesh loading**, therefore it is usually called inside some callback, most probably inside `on_mesh_loaded()` ([see below](#callable-functions-1)).
 - `surface_to_physical_group(entity)`: takes the tag of a surface entity and returns the associated physical group. This function is available **only after mesh loading**, therefore it is usually called inside some callback, most probably inside `on_mesh_loaded()` ([see below](#callable-functions-1)). In addition, **exactly one physical group** must be associated to `entity`, otherwise an exception is thrown by the runtime.

share/examples/sponge/abc.geo

+SetFactory("OpenCASCADE");
+
+Box(1) = { 1, 1, 0, 1, 1, 0.1 };
+
+Physical Surface("abc", 100) = {1,2,3,4};
+
+MeshSize{ PointsOf{ Volume{1}; } } = 0.11;

share/examples/sponge/params_abc.lua

+sim.dt = 1e-12                          -- timestep size
+sim.timesteps = 20000                   -- num of iterations
+sim.gmsh_model = "abc.geo"
+sim.name = "abc"
+sim.use_gpu = 0                         -- 0: cpu, 1: gpu
+sim.approx_order = 1                    -- approximation order
+sim.time_integrator = "rk4"
+postpro.silo_output_rate = 100
+postpro.cycle_print_rate = 10           -- console print rate
+
+postpro["E"].silo_mode = "nodal"
+postpro["H"].silo_mode = "none"
+postpro["J"].silo_mode = "none"
+
+debug = {};
+debug.dump_cell_ranks = true;
+
+materials.epsilon = function(tag, x, y, z)
+    return 1.0;
+end
+
+materials.mu = function(tag, x, y, z)
+    return 1.0;
+end
+
+materials.sigma = function(tag, x, y, z)
+    return 0.0;
+end
+
+function electric_initial_condition(x,y,z)
+    local rx = x - 1.5;
+    local ry = y - 1.5;
+
+    local Ex = 0;
+    local Ey = 0;
+    local Ez = math.exp( -50*(rx^2 + ry^2) )
+
+    return Ex, Ey, Ez
+end
+
+function setup_boundary_conditions()
+    local abc_surfs = physical_group_to_surfaces(100);
+    for i,v in ipairs(abc_surfs) do
+        bndconds[v] = {};
+        bndconds[v].kind = "impedance";
+    end
+end
+
+function on_mesh_loaded()
+    setup_boundary_conditions()
+end

share/examples/sponge/params_sponge.lua

+sim.dt = 1e-12                      -- timestep size
+sim.timesteps = 20000               -- num of iterations
+sim.gmsh_model = "sponge.geo"
+sim.name = "sponge"
+sim.use_gpu = 0                     -- 0: cpu, 1: gpu
+sim.approx_order = 1                -- approximation order
+sim.time_integrator = "rk4"
+postpro.silo_output_rate = 100
+postpro.cycle_print_rate = 10       -- console print rate
+
+postpro["E"].silo_mode = "nodal"
+postpro["H"].silo_mode = "none"
+postpro["J"].silo_mode = "none"
+
+sigma_max = 0.25;
+
+debug = {};
+debug.dump_cell_ranks = true;
+
+function len(x, y)
+    return math.sqrt(x*x + y*y)
+end
+
+function sigma_profile(r, sigma_max, dr)
+    return sigma_max*(r/dr)^2;
+end
+
+local ep = 1 - sim.dt*sigma_max/const.eps0;
+if ep <= 0 and sim.time_integrator == "rk4" then
+    print("Computation will not be stable: " .. ep)
+end
+
+materials.sigma = function(tag, x, y, z)
+    local dr = 1;
+
+    local pgs = volume_to_physical_groups(tag);
+
+    for i,pg in ipairs(pgs) do
+        if pg == 100 then
+            return 0;
+        end
+
+        -- Bottom left
+        if pg == 101 then
+            local r = len(x-1, y-1)
+            return sigma_profile(r, sigma_max, dr);
+        end
+
+        -- Bottom center
+        if pg == 102 then
+            local r = 1-y
+            return sigma_profile(r, sigma_max, dr);
+        end
+
+        -- Bottom right
+        if pg == 103 then 
+            local r = len(x-2, y-1)
+            return sigma_profile(r, sigma_max, dr);
+        end
+
+        -- Center left
+        if pg == 104 then
+            local r = 1-x
+            return sigma_profile(r, sigma_max, dr);
+        end
+
+        -- Center right
+        if pg == 105 then 
+            local r = x-2
+            return sigma_profile(r, sigma_max, dr);
+        end
+
+        -- Top left
+        if pg == 106 then
+            local r = len(x-1, y-2)
+            return sigma_profile(r, sigma_max, dr);
+        end
+
+        -- Top center
+        if pg == 107 then
+            local r = y-2
+            return sigma_profile(r, sigma_max, dr);
+        end
+
+        -- Top right
+        if pg == 108 then
+            local r = len(x-2, y-2)
+            return sigma_profile(r, sigma_max, dr);
+        end
+    end
+end
+
+materials.epsilon = function(tag, x, y, z)
+    return 1.0;
+end
+
+materials.mu = function(tag, x, y, z)
+    return 1.0;
+end
+
+function electric_initial_condition(x,y,z)
+    local rx = x - 1.5;
+    local ry = y - 1.5;
+
+    local Ex = 0;
+    local Ey = 0;
+    local Ez = math.exp( -50*(rx^2 + ry^2) )
+
+    return Ex, Ey, Ez
+end
+
+

share/examples/sponge/sponge.geo

+SetFactory("OpenCASCADE");
+
+h = 0.10;
+eps = h/10;
+
+Box(1) = { 1, 1, 0, 1, 1, h };    // Center box
+Box(2) = { 0, 0, 0, 1, 1, h };    // Sponge bottom left
+Box(3) = { 1, 0, 0, 1, 1, h };    // Sponge bottom center
+Box(4) = { 2, 0, 0, 1, 1, h };    // Sponge bottom right
+Box(5) = { 0, 1, 0, 1, 1, h };    // Sponge center left
+Box(6) = { 2, 1, 0, 1, 1, h };    // Sponge center right
+Box(7) = { 0, 2, 0, 1, 1, h };    // Sponge top left
+Box(8) = { 1, 2, 0, 1, 1, h };    // Sponge top center
+Box(9) = { 2, 2, 0, 1, 1, h };    // Sponge top right
+
+Coherence;
+
+bnd_bottom() = Surface In BoundingBox{-eps, -eps, -eps, 3+eps, 3+eps, eps};
+bnd_top() = Surface In BoundingBox{-eps, -eps, h-eps, 3+eps, 3+eps, h+eps};
+bnd_left() = Surface In BoundingBox{-eps, -eps, -eps, eps, 3+eps, h+eps};
+bnd_right() = Surface In BoundingBox{3-eps, -eps, -eps, 3+eps, 3+eps, h+eps};
+bnd_front() = Surface In BoundingBox{-eps, -eps, -eps, 3+eps, eps, h+eps};
+bnd_back() = Surface In BoundingBox{-eps, 3-eps, -eps, 3+eps, 3+eps, h+eps};
+
+Physical Volume("domain", 100) = {1};
+Physical Volume("sponge_bl", 101) = {2};
+Physical Volume("sponge_bc", 102) = {3};
+Physical Volume("sponge_br", 103) = {4};
+Physical Volume("sponge_cl", 104) = {5};
+Physical Volume("sponge_cr", 105) = {6};
+Physical Volume("sponge_tl", 106) = {7};
+Physical Volume("sponge_tc", 107) = {8};
+Physical Volume("sponge_tr", 108) = {9};
+
+Physical Volume("all", 110) = {1,2,3,4,5,6,7,8,9};
+
+Physical Surface("bottom", 100) = bnd_bottom();
+Physical Surface("top", 101) = bnd_top();
+Physical Surface("left", 102) = bnd_left();
+Physical Surface("right", 103) = bnd_right();
+Physical Surface("front", 104) = bnd_front();
+Physical Surface("back", 105) = bnd_back();
+
+Mesh.Algorithm3D = 10;
+
+MeshSize{ PointsOf{ Volume{:}; } } = 0.11;

src/maxwell/maxwell_solver.cpp

@@ -346,15 +346,29 @@ register_lua_callbacks(maxwell::parameter_loader& mpl, model& mod)
 {
     sol::state& lua = mpl.lua_state();
     auto volume_physical_group_entities = [&](int tag) -> auto {
+        std::cout << "[GMSH/DG Lua API] volume_physical_group_entities() is ";
+        std::cout << "deprecated, use physical_group_to_volumes()" << std::endl;
         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 {
+        std::cout << "[GMSH/DG Lua API] surface_physical_group_entities() is ";
+        std::cout << "deprecated, use physical_group_to_surfaces()" << std::endl;
         return sol::as_table( mod.lookup_physical_group_entities(2, tag) );
     };
     lua["surface_physical_group_entities"] = surface_physical_group_entities;
 
+    auto physical_group_to_volumes = [&](int tag) -> auto {
+        return sol::as_table( mod.lookup_physical_group_entities(3, tag) );
+    };
+    lua["physical_group_to_volumes"] = physical_group_to_volumes;
+
+    auto physical_group_to_surfaces = [&](int tag) -> auto {
+        return sol::as_table( mod.lookup_physical_group_entities(2, tag) );
+    };
+    lua["physical_group_to_surfaces"] = physical_group_to_surfaces;
+
     auto volume_to_physical_group = [&](int tag) -> auto {
         return mod.lookup_physical_group_for_entity(3, tag);
     };