From ad399724d242e24c466597cabef39c0fc6ec47eb Mon Sep 17 00:00:00 2001
From: Matteo Cicuttin <datafl4sh@toxicnet.eu>
Date: Mon, 14 Mar 2022 16:08:02 +0100
Subject: [PATCH] Full support for physical entities. Mandatory with MPI.
---
doc/lua_api.md | 1 +
include/libgmshdg/gmsh_io.h | 23 ++++++++++--
include/maxwell/maxwell_common.h | 18 +++++-----
share/validation/params_twomat.lua | 56 ++++++++++++++++++------------
share/validation/twomat.geo | 6 ++++
src/libgmshdg/gmsh_io.cpp | 16 +++++----
src/libgmshdg/param_loader.cpp | 28 ++++-----------
src/maxwell/maxwell_common.cpp | 10 +++---
src/maxwell/maxwell_cpu.cpp | 4 +--
src/maxwell/maxwell_gpu.cpp | 8 ++---
src/maxwell/maxwell_interface.cpp | 21 +++++++++--
src/maxwell/maxwell_postpro.cpp | 4 +--
src/maxwell/maxwell_solver.cpp | 2 +-
13 files changed, 119 insertions(+), 78 deletions(-)
diff --git a/doc/lua_api.md b/doc/lua_api.md
index decb61f..2d01e98 100644
--- a/doc/lua_api.md
+++ b/doc/lua_api.md
@@ -32,6 +32,7 @@ This document describes the API available on the version v0.4 of the solver.
### Parallel execution information (available only if MPI support is enabled)
The parallel solver runs as separate MPI processes. As such, each process loads and executes its own instance of the configuration script. This means that you must beware of global variables or other shared state, because modifications are **not** reflected across processes. Even if the configuration script should in general not depend on the rank of the current process, the following variables are exposed:
+
- `parallel.comm_rank` (integer): The MPI rank of the current solver process
- `parallel.comm_size` (integer): The MPI communicator size
diff --git a/include/libgmshdg/gmsh_io.h b/include/libgmshdg/gmsh_io.h
index c33ac53..e9baad2 100644
--- a/include/libgmshdg/gmsh_io.h
+++ b/include/libgmshdg/gmsh_io.h
@@ -73,6 +73,11 @@ struct boundary_descriptor
#endif /* USE_MPI */
{}
+ int gmsh_tag(void) const
+ {
+ return gmsh_entity;
+ }
+
int material_tag(void) const
{
#ifdef USE_MPI
@@ -95,6 +100,15 @@ struct physical_group_descriptor {
}
};
+inline std::ostream&
+operator<<(std::ostream& os, const physical_group_descriptor& pgd)
+{
+ os << "["<< pgd.dim << ", " << pgd.tag << "] ";
+ for (auto& et : pgd.entity_tags)
+ os << et << " ";
+ return os;
+}
+
#ifdef USE_MPI
struct comm_descriptor
{
@@ -107,8 +121,8 @@ struct comm_descriptor
std::vector<element_key> fks;
};
-inline
-std::ostream& operator<<(std::ostream& os, const comm_descriptor& cd)
+inline std::ostream&
+operator<<(std::ostream& os, const comm_descriptor& cd)
{
os << "tag: " << cd.boundary_tag << ", PM: " << cd.partition_mine << ", PO: " << cd.partition_other;
os << ", RM: " << cd.rank_mine << ", RO: " << cd.rank_other << ", map size: " << cd.dof_mapping.size();
@@ -178,6 +192,11 @@ struct surface_descriptor
std::sort(face_numbers.begin(), face_numbers.end(), comp);
}
+ int gmsh_tag(void) const
+ {
+ return surface_entity;
+ }
+
int material_tag(void) const
{
#ifdef USE_MPI
diff --git a/include/maxwell/maxwell_common.h b/include/maxwell/maxwell_common.h
index 1e902e1..1997435 100644
--- a/include/maxwell/maxwell_common.h
+++ b/include/maxwell/maxwell_common.h
@@ -25,7 +25,7 @@ eval_boundary_sources(const model& mod, const parameter_loader& mpl,
size_t face_num_base = 0;
for (auto& e : mod)
{
- auto etag = e.material_tag();
+ auto etag = e.gmsh_tag();
for (size_t iF = 0; iF < e.num_faces(); iF++)
{
auto gfnum = face_num_base + iF;
@@ -45,7 +45,7 @@ eval_boundary_sources(const model& mod, const parameter_loader& mpl,
vec3d n = state.eds[e.number()].normals.row(iF);
- switch ( mpl.boundary_type(bd.material_tag()) )
+ switch ( mpl.boundary_type(bd.gmsh_tag()) )
{
case boundary_condition::UNSPECIFIED:
break;
@@ -67,10 +67,10 @@ eval_boundary_sources(const model& mod, const parameter_loader& mpl,
case boundary_condition::PLANE_WAVE_E: {
auto fE = [&](const point_3d& pt) -> vec3d {
- return mpl.eval_boundary_source(bd.material_tag(), pt, state.curr_time);
+ return mpl.eval_boundary_source(bd.gmsh_tag(), pt, state.curr_time);
};
auto fH = [&](const point_3d& pt) -> vec3d {
- vec3d H = mpl.eval_boundary_source(bd.material_tag(), pt, state.curr_time);
+ vec3d H = mpl.eval_boundary_source(bd.gmsh_tag(), pt, state.curr_time);
return n.cross(H)/Z;
};
matxd prjE = e.project_on_face(iF, fE);
@@ -149,9 +149,9 @@ eval_boundary_sources_new(const model& mod, const parameter_loader& mpl,
for (auto& [si, sd] : bds)
{
const auto& e = mod.entity_at(sd.adjacent_volume_entity);
- auto etag = e.material_tag();
+ auto etag = e.gmsh_tag();
- auto mtag = sd.material_tag();
+ auto mtag = sd.gmsh_tag();
switch ( mpl.boundary_type( mtag ) )
{
case boundary_condition::UNSPECIFIED:
@@ -219,7 +219,7 @@ eval_interface_sources(const model& mod, const parameter_loader& mpl,
vec3d n = state.eds[e.number()].normals.row(iF);
- switch ( mpl.interface_type(bd.material_tag()) )
+ switch ( mpl.interface_type(bd.gmsh_tag()) )
{
case interface_condition::UNSPECIFIED:
break;
@@ -227,7 +227,7 @@ eval_interface_sources(const model& mod, const parameter_loader& mpl,
case interface_condition::E_FIELD: {
/* This is BS actually... */
auto fE = [&](const point_3d& pt) -> vec3d {
- return mpl.eval_interface_source(bd.material_tag(), pt, state.curr_time);
+ return mpl.eval_interface_source(bd.gmsh_tag(), pt, state.curr_time);
};
matxd prjE = e.project_on_face(iF, fE);
@@ -242,7 +242,7 @@ eval_interface_sources(const model& mod, const parameter_loader& mpl,
case interface_condition::SURFACE_CURRENT: {
auto fJ = [&](const point_3d& pt) -> vec3d {
- vec3d J = mpl.eval_interface_source(bd.material_tag(), pt, state.curr_time);
+ vec3d J = mpl.eval_interface_source(bd.gmsh_tag(), pt, state.curr_time);
return n.cross(J);
};
diff --git a/share/validation/params_twomat.lua b/share/validation/params_twomat.lua
index 67b0160..8be05f6 100644
--- a/share/validation/params_twomat.lua
+++ b/share/validation/params_twomat.lua
@@ -13,35 +13,47 @@ 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"
+function setup_materials()
+ local bulk = volume_physical_group_entities(1);
+ for i,v in ipairs(bulk) do
+ materials[v] = {}
+ materials[v].epsilon = 1
+ materials[v].mu = 1
+ materials[v].sigma = 0
+ end
end
-bndconds[7] = {}
-bndconds[7].kind = "impedance"
-
local freq = 3e8
-function source(tag, x, y, z, t)
+function plane_wave(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
+function setup_boundary_conditions()
+ local source = surface_physical_group_entities(10);
+ for i,v in ipairs(source) do
+ bndconds[v] = {};
+ bndconds[v].kind = "plane_wave_E";
+ bndconds[v].source = plane_wave;
+ end
+
+ local pmc = surface_physical_group_entities(11);
+ for i,v in ipairs(pmc) do
+ bndconds[v] = {};
+ bndconds[v].kind = "pmc";
+ end
+
+ local abc = surface_physical_group_entities(13);
+ for i,v in ipairs(abc) do
+ bndconds[v] = {};
+ bndconds[v].kind = "impedance";
+ end
+end
+
+function on_mesh_loaded()
+ setup_materials();
+ setup_boundary_conditions();
+end
diff --git a/share/validation/twomat.geo b/share/validation/twomat.geo
index 6545506..5ca38ad 100644
--- a/share/validation/twomat.geo
+++ b/share/validation/twomat.geo
@@ -3,4 +3,10 @@ Box(1) = {0, 0, 0, 1, 1, 0.1};
Box(2) = {1, 0, 0, 1, 1, 0.1};
Coherence;
+Physical Volume("bulk", 1) = {1,2};
+Physical Surface("source", 10) = {1};
+Physical Surface("pmc", 11) = {3,4,8,9};
+Physical Surface("pec", 12) = {5,6,10,11};
+Physical Surface("abc", 13) = {7};
+
MeshSize{:} = 0.075;
diff --git a/src/libgmshdg/gmsh_io.cpp b/src/libgmshdg/gmsh_io.cpp
index b67a6b4..7fcedca 100644
--- a/src/libgmshdg/gmsh_io.cpp
+++ b/src/libgmshdg/gmsh_io.cpp
@@ -423,7 +423,7 @@ model::import_physical_groups(void)
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);
+ MPI_Bcast(&pg.tag, 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);
@@ -435,7 +435,9 @@ model::import_physical_groups(void)
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);
+ assert(pg.dim != 0);
+ MPI_Bcast(&pg.tag, 1, MPI_INT, 0, MPI_COMM_WORLD);
+ assert(pg.tag != 0);
size_t esize;
MPI_Bcast(&esize, 1, MPI_UNSIGNED_LONG_LONG, 0, MPI_COMM_WORLD);
pg.entity_tags.resize(esize);
@@ -1114,10 +1116,10 @@ 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;
-
+ auto pg_begin = physical_groups.begin();
+ auto pg_end = physical_groups.end();
+ auto itor = std::lower_bound(pg_begin , pg_end, pg);
+ if (!(itor == pg_end) && !(pg < *pg_begin))
+ return itor->entity_tags;
return std::vector<int>{};
}
-
diff --git a/src/libgmshdg/param_loader.cpp b/src/libgmshdg/param_loader.cpp
index 569a959..0adc570 100644
--- a/src/libgmshdg/param_loader.cpp
+++ b/src/libgmshdg/param_loader.cpp
@@ -16,25 +16,6 @@
#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)
@@ -86,8 +67,13 @@ parameter_loader_base::init(void)
¶meter_loader_base::enable_volume_sources,
this);
- lua.set_function("getEntitiesForPhysicalGroup",
- &lua_getEntitiesForPhysicalGroup);
+ auto exit_abc = [](){
+#ifdef USE_MPI
+ MPI_Finalize();
+#endif /* USE_MPI */
+ exit(0);
+ };
+ lua["exit"] = exit_abc;
#ifdef USE_MPI
int comm_rank, comm_size;
diff --git a/src/maxwell/maxwell_common.cpp b/src/maxwell/maxwell_common.cpp
index bda7f6d..eec2b36 100644
--- a/src/maxwell/maxwell_common.cpp
+++ b/src/maxwell/maxwell_common.cpp
@@ -48,7 +48,7 @@ init_boundary_conditions(const model& mod, const parameter_loader& mpl,
else
{
#endif /* USE_MPI */
- switch ( mpl.boundary_type(bd.material_tag()) )
+ switch ( mpl.boundary_type(bd.gmsh_tag()) )
{
case boundary_condition::UNSPECIFIED:
case boundary_condition::PEC:
@@ -103,7 +103,7 @@ init_lax_milgram(const model& mod, const parameter_loader& mpl,
{
auto& pe = e.cell(iT);
auto bar = pe.barycenter();
- auto tag = e.material_tag();
+ auto tag = e.gmsh_tag();
auto my_eps = mpl.epsilon(tag, bar);
auto my_mu = mpl.mu(tag, bar);
@@ -125,7 +125,7 @@ init_lax_milgram(const model& mod, const parameter_loader& mpl,
auto& ne_e = mod.entity_at(neigh.iE);
auto& ne_pe = ne_e.cell(neigh.iT);
auto ne_bar = ne_pe.barycenter();
- auto ne_tag = ne_e.material_tag();
+ auto ne_tag = ne_e.gmsh_tag();
auto ne_eps = mpl.epsilon(ne_tag, ne_bar);
auto ne_mu = mpl.mu(ne_tag, ne_bar);
@@ -228,7 +228,7 @@ export_vector_field_nodal(const model& mod, silo& sdb, const vecxd& Fx,
auto& pe = e.cell(iT);
auto cgofs = e.cell_world_dof_offset(iT);
auto nodetags = pe.node_tags();
- auto mp = mf(e.material_tag(), pe.barycenter());
+ auto mp = mf(e.gmsh_tag(), pe.barycenter());
for (size_t iD = 0; iD < 4; iD++)
{
size_t tag = nodetags[iD];
@@ -282,7 +282,7 @@ export_vector_field_zonal(const model& mod, silo& sdb, const vecxd& Fx,
auto& pe = e.cell(iT);
auto& re = e.cell_refelem(pe);
vecxd phi_bar = re.basis_functions({1./3., 1./3., 1./3.});
- auto mp = mf(e.material_tag(), pe.barycenter());
+ auto mp = mf(e.gmsh_tag(), pe.barycenter());
auto num_bf = re.num_basis_functions();
auto ofs = e.cell_world_dof_offset(iT);
auto gi = e.cell_world_index_by_gmsh(iT);
diff --git a/src/maxwell/maxwell_cpu.cpp b/src/maxwell/maxwell_cpu.cpp
index b54294b..e6d7782 100644
--- a/src/maxwell/maxwell_cpu.cpp
+++ b/src/maxwell/maxwell_cpu.cpp
@@ -33,7 +33,7 @@ init_matparams(const model& mod, solver_state& state,
for (auto& e : mod)
{
- auto etag = e.material_tag();
+ auto etag = e.gmsh_tag();
for (size_t iT = 0; iT < e.num_cells(); iT++)
{
auto& pe = e.cell(iT);
@@ -787,7 +787,7 @@ eval_volume_sources(const model& mod, const parameter_loader& mpl, solver_state&
{
for (auto& e : mod)
{
- auto etag = e.material_tag();
+ auto etag = e.gmsh_tag();
if ( not mpl.volume_has_source(etag) )
continue;
diff --git a/src/maxwell/maxwell_gpu.cpp b/src/maxwell/maxwell_gpu.cpp
index 5f95498..e13aff2 100644
--- a/src/maxwell/maxwell_gpu.cpp
+++ b/src/maxwell/maxwell_gpu.cpp
@@ -32,7 +32,7 @@ init_matparams(const model& mod, solver_state_gpu& state,
for (auto& e : mod)
{
- auto etag = e.material_tag();
+ auto etag = e.gmsh_tag();
for (size_t iT = 0; iT < e.num_cells(); iT++)
{
auto& pe = e.cell(iT);
@@ -87,7 +87,7 @@ build_source_compression_tables(const model& mod, solver_state_gpu& state, const
auto bd = bds[gfnum];
if (bd.type == face_type::INTERFACE)
{
- switch ( mpl.interface_type(bd.material_tag()) )
+ switch ( mpl.interface_type(bd.gmsh_tag()) )
{
case interface_condition::E_FIELD:
case interface_condition::SURFACE_CURRENT:
@@ -102,7 +102,7 @@ build_source_compression_tables(const model& mod, solver_state_gpu& state, const
if (bd.type == face_type::BOUNDARY)
{
- switch ( mpl.boundary_type(bd.material_tag()) )
+ switch ( mpl.boundary_type(bd.gmsh_tag()) )
{
case boundary_condition::PLANE_WAVE_E:
for (size_t i = 0; i < num_fluxes; i++)
@@ -594,7 +594,7 @@ eval_volume_sources(const model& mod, const parameter_loader& mpl, solver_state_
{
for (auto& e : mod)
{
- auto etag = e.material_tag();
+ auto etag = e.gmsh_tag();
if ( not mpl.volume_has_source(etag) )
continue;
diff --git a/src/maxwell/maxwell_interface.cpp b/src/maxwell/maxwell_interface.cpp
index ab6c81e..32899c4 100644
--- a/src/maxwell/maxwell_interface.cpp
+++ b/src/maxwell/maxwell_interface.cpp
@@ -524,6 +524,20 @@ parameter_loader::parameter_loader()
bool
parameter_loader::validate_materials(const std::string& mpn, int tag) const
{
+#ifdef USE_MPI
+ int rank;
+ MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+ if (rank == 0) {
+#endif /* USE_MPI */
+ int pdim, ptag;
+ gm::getParent(3,tag,pdim,ptag);
+ if (ptag == -1)
+ return true; /* This is not a computational entity, only the parent of a partition. */
+#ifdef USE_MPI
+ } else {
+ return true; /* The script is the same for all ranks, do validation only on rank 0. */
+ }
+#endif /* USE_MPI */
auto mfun = lua[SEC_MATERIALS][mpn];
if ( mfun.valid() )
return true;
@@ -636,12 +650,12 @@ parameter_loader::validate_conditions(const model& mod) const
break;
case face_type::BOUNDARY:
- if (not validate_boundary_condition(mod, bd.material_tag()) )
+ if (not validate_boundary_condition(mod, bd.gmsh_tag()) )
success = false;
break;
case face_type::INTERFACE:
- if (not validate_interface_condition(mod, bd.material_tag()) )
+ if (not validate_interface_condition(mod, bd.gmsh_tag()) )
success = false;
break;
@@ -659,6 +673,7 @@ bool
parameter_loader::validate_model_params(const model& mod) const
{
bool success = true;
+
for (auto& e : mod)
{
auto tag = e.material_tag();
@@ -671,7 +686,7 @@ parameter_loader::validate_model_params(const model& mod) const
if ( not validate_conditions(mod) )
success = false;
-
+
return success;
}
diff --git a/src/maxwell/maxwell_postpro.cpp b/src/maxwell/maxwell_postpro.cpp
index 8490a59..c84baf1 100644
--- a/src/maxwell/maxwell_postpro.cpp
+++ b/src/maxwell/maxwell_postpro.cpp
@@ -92,8 +92,8 @@ compute_energy(const model& mod, const solver_state& state, const parameter_load
const auto cbs = re.num_basis_functions();
const auto bar = pe.barycenter();
- auto eps = mpl.epsilon( e.material_tag(), bar );
- auto mu = mpl.mu( e.material_tag(), bar );
+ auto eps = mpl.epsilon( e.gmsh_tag(), bar );
+ auto mu = mpl.mu( e.gmsh_tag(), bar );
field_values mat(eps, eps, eps, mu, mu, mu);
diff --git a/src/maxwell/maxwell_solver.cpp b/src/maxwell/maxwell_solver.cpp
index 600713f..7c17ab6 100644
--- a/src/maxwell/maxwell_solver.cpp
+++ b/src/maxwell/maxwell_solver.cpp
@@ -157,7 +157,7 @@ reinit_materials(const model& mod, maxwell::solver_state& state, maxwell::parame
{
for (auto& e : mod)
{
- auto tag = e.material_tag();
+ auto tag = e.gmsh_tag();
for (size_t iT = 0; iT < e.num_cells(); iT++)
{
auto& pe = e.cell(iT);
--
GitLab