From 4765d09fecf2d809f452285a1db8ea9fcbd672c2 Mon Sep 17 00:00:00 2001
From: Matteo Cicuttin <datafl4sh@toxicnet.eu>
Date: Mon, 30 Aug 2021 15:43:49 +0200
Subject: [PATCH] Reworked Silo export for MPI.
---
CMakeLists.txt | 15 +++++--
include/maxwell_interface.h | 6 +--
src/maxwell_common.cpp | 69 +++++++++++++++++++++++---------
src/maxwell_cpu.cpp | 80 +++++++++++++++++++++++--------------
src/maxwell_gpu.cpp | 10 +++--
src/maxwell_solver.cpp | 35 ++++++++--------
6 files changed, 137 insertions(+), 78 deletions(-)
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 81dbb89..105a8c9 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -275,6 +275,13 @@ option(OPT_DISABLE_TEXTURE_CACHE "On ROCm 4.1 texture stuff is broken, disable i
if (OPT_DISABLE_TEXTURE_CACHE)
add_definitions(-DDISABLE_TEXTURE_CACHE)
endif()
+
+######################################################################
+## Enable code to handle geometric diode simulation for UniUD project
+option(OPT_ENABLE_EXP_GEOMETRIC_DIODE "Enable Geometric Diode stuff [EXPERIMENTAL][NO GPU]")
+if (OPT_ENABLE_EXP_GEOMETRIC_DIODE)
+ add_definitions(-DENABLE_EXP_GEOMETRIC_DIODE)
+endif()
######################################################################
## Configure targets
include_directories(contrib/sgr)
@@ -312,11 +319,11 @@ endif()
########## mpi testing stuff
-add_executable(testpart src/testpart.cpp)
-target_link_libraries(testpart gmshdg)
+#add_executable(testpart src/testpart.cpp)
+#target_link_libraries(testpart gmshdg)
-add_executable(testmod src/testmod.cpp)
-target_link_libraries(testmod gmshdg)
+#add_executable(testmod src/testmod.cpp)
+#target_link_libraries(testmod gmshdg)
########## maxwell solver
set(MAXWELL_SOURCES src/maxwell_interface.cpp
diff --git a/include/maxwell_interface.h b/include/maxwell_interface.h
index fc686a3..2c344cc 100644
--- a/include/maxwell_interface.h
+++ b/include/maxwell_interface.h
@@ -325,8 +325,7 @@ struct solver_state_gpu
void init_from_model(const model&, solver_state&);
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&, const solver_state&, const parameter_loader&);
-void export_fields_to_silo_2(const model&, solver_state&, const parameter_loader&);
+void export_fields_to_silo(const model&, solver_state&, const parameter_loader&);
void timestep(solver_state&, const parameter_loader&, time_integrator_type);
void prepare_sources(const model&, maxwell::solver_state&, maxwell::parameter_loader&);
void do_sources(const model&, maxwell::solver_state&, maxwell::parameter_loader&);
@@ -339,8 +338,7 @@ void gather_field(const model&, maxwell::field&, int, MPI_Comm);
void init_from_model(const model&, solver_state_gpu&);
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&, const solver_state_gpu&, const parameter_loader&);
-void export_fields_to_silo_2(const model&, solver_state_gpu&, const parameter_loader&);
+void export_fields_to_silo(const model&, solver_state_gpu&, const parameter_loader&);
void timestep(solver_state_gpu&, const parameter_loader&, time_integrator_type);
void prepare_sources(const model&, maxwell::solver_state_gpu&, maxwell::parameter_loader&);
void do_sources(const model&, maxwell::solver_state_gpu&, maxwell::parameter_loader&);
diff --git a/src/maxwell_common.cpp b/src/maxwell_common.cpp
index 3a8bf17..ba1e8db 100644
--- a/src/maxwell_common.cpp
+++ b/src/maxwell_common.cpp
@@ -151,12 +151,18 @@ void export_E_field_nodal(const model& mod, silo& sdb, const vecxd& Ex,
std::vector<double> accum_Ez( sdb.num_nodes(), 0.0 );
std::vector<size_t> count( sdb.num_nodes(), 0 );
+#ifdef USE_MPI
+ for (auto& etag : mod.world_entities_tags())
+ {
+ auto& e = mod.entity_world_lookup(etag);
+#else /* USE_MPI */
for (auto& e : mod)
{
+#endif /* USE_MPI */
for (size_t iT = 0; iT < e.num_cells(); iT++)
{
auto& pe = e.cell(iT);
- auto cgofs = e.cell_model_dof_offset(iT);
+ auto cgofs = e.cell_world_dof_offset(iT);
auto nodetags = pe.node_tags();
for (size_t iD = 0; iD < 4; iD++)
{
@@ -191,12 +197,18 @@ void export_H_field_nodal(const model& mod, silo& sdb, const vecxd& Hx,
std::vector<double> accum_Hz( sdb.num_nodes(), 0.0 );
std::vector<size_t> count( sdb.num_nodes(), 0 );
+#ifdef USE_MPI
+ for (auto& etag : mod.world_entities_tags())
+ {
+ auto& e = mod.entity_world_lookup(etag);
+#else /* USE_MPI */
for (auto& e : mod)
{
+#endif /* USE_MPI */
for (size_t iT = 0; iT < e.num_cells(); iT++)
{
auto& pe = e.cell(iT);
- auto cgofs = e.cell_model_dof_offset(iT);
+ auto cgofs = e.cell_world_dof_offset(iT);
auto nodetags = pe.node_tags();
for (size_t iD = 0; iD < 4; iD++)
{
@@ -231,15 +243,20 @@ void export_J_field_nodal(const model& mod, silo& sdb, const vecxd& Ex,
std::vector<double> accum_Jz( sdb.num_nodes(), 0.0 );
std::vector<size_t> count( sdb.num_nodes(), 0 );
+#ifdef USE_MPI
+ for (auto& etag : mod.world_entities_tags())
+ {
+ auto& e = mod.entity_world_lookup(etag);
+#else /* USE_MPI */
for (auto& e : mod)
{
- auto etag = e.material_tag();
+#endif /* USE_MPI */
for (size_t iT = 0; iT < e.num_cells(); iT++)
{
auto& pe = e.cell(iT);
auto bar = pe.barycenter();
- auto sigma = mpl.sigma(etag, bar);
- auto cgofs = e.cell_model_dof_offset(iT);
+ auto sigma = mpl.sigma(e.material_tag(), bar);
+ auto cgofs = e.cell_world_dof_offset(iT);
auto nodetags = pe.node_tags();
for (size_t iD = 0; iD < 4; iD++)
{
@@ -275,17 +292,22 @@ void export_E_field_zonal(const model& mod, silo& sdb, const vecxd& Ex,
std::vector<double> out_Ey( mod.num_cells() );
std::vector<double> out_Ez( mod.num_cells() );
- for (size_t iE = 0; iE < mod.num_entities(); iE++)
+#ifdef USE_MPI
+ for (auto& etag : mod.world_entities_tags())
+ {
+ auto& e = mod.entity_world_lookup(etag);
+#else /* USE_MPI */
+ for (auto& e : mod)
{
- auto& e = mod.entity_at(iE);
+#endif /* USE_MPI */
for (size_t iT = 0; iT < e.num_cells(); iT++)
{
auto& pe = e.cell(iT);
auto& re = e.cell_refelem(pe);
vecxd phi_bar = re.basis_functions({1./3., 1./3., 1./3.});
auto num_bf = re.num_basis_functions();
- auto ofs = e.cell_model_dof_offset(iT);
- auto gi = e.cell_model_index_by_gmsh(iT);
+ auto ofs = e.cell_world_dof_offset(iT);
+ auto gi = e.cell_world_index_by_gmsh(iT);
out_Ex[gi] = Ex.segment(ofs, num_bf).dot(phi_bar);
out_Ey[gi] = Ey.segment(ofs, num_bf).dot(phi_bar);
out_Ez[gi] = Ez.segment(ofs, num_bf).dot(phi_bar);
@@ -305,17 +327,22 @@ void export_H_field_zonal(const model& mod, silo& sdb, const vecxd& Hx,
std::vector<double> out_Hy( mod.num_cells() );
std::vector<double> out_Hz( mod.num_cells() );
- for (size_t iE = 0; iE < mod.num_entities(); iE++)
+#ifdef USE_MPI
+ for (auto& etag : mod.world_entities_tags())
{
- auto& e = mod.entity_at(iE);
+ auto& e = mod.entity_world_lookup(etag);
+#else /* USE_MPI */
+ for (auto& e : mod)
+ {
+#endif /* USE_MPI */
for (size_t iT = 0; iT < e.num_cells(); iT++)
{
auto& pe = e.cell(iT);
auto& re = e.cell_refelem(pe);
vecxd phi_bar = re.basis_functions({1./3., 1./3., 1./3.});
auto num_bf = re.num_basis_functions();
- auto ofs = e.cell_model_dof_offset(iT);
- auto gi = e.cell_model_index_by_gmsh(iT);
+ auto ofs = e.cell_world_dof_offset(iT);
+ auto gi = e.cell_world_index_by_gmsh(iT);
out_Hx[gi] = Hx.segment(ofs, num_bf).dot(phi_bar);
out_Hy[gi] = Hy.segment(ofs, num_bf).dot(phi_bar);
out_Hz[gi] = Hz.segment(ofs, num_bf).dot(phi_bar);
@@ -335,20 +362,24 @@ void export_J_field_zonal(const model& mod, silo& sdb, const vecxd& Ex,
std::vector<double> out_Jy( mod.num_cells() );
std::vector<double> out_Jz( mod.num_cells() );
- for (size_t iE = 0; iE < mod.num_entities(); iE++)
+#ifdef USE_MPI
+ for (auto& etag : mod.world_entities_tags())
+ {
+ auto& e = mod.entity_world_lookup(etag);
+#else /* USE_MPI */
+ for (auto& e : mod)
{
- auto& e = mod.entity_at(iE);
- auto etag = e.material_tag();
+#endif /* USE_MPI */
for (size_t iT = 0; iT < e.num_cells(); iT++)
{
auto& pe = e.cell(iT);
auto& re = e.cell_refelem(pe);
auto bar = pe.barycenter();
- auto sigma = mpl.sigma(etag, bar);
+ auto sigma = mpl.sigma(e.material_tag(), bar);
vecxd phi_bar = re.basis_functions({1./3., 1./3., 1./3.});
auto num_bf = re.num_basis_functions();
- auto ofs = e.cell_model_dof_offset(iT);
- auto gi = e.cell_model_index_by_gmsh(iT);
+ auto ofs = e.cell_world_dof_offset(iT);
+ auto gi = e.cell_world_index_by_gmsh(iT);
out_Jx[gi] = sigma*Ex.segment(ofs, num_bf).dot(phi_bar);
out_Jy[gi] = sigma*Ey.segment(ofs, num_bf).dot(phi_bar);
out_Jz[gi] = sigma*Ez.segment(ofs, num_bf).dot(phi_bar);
diff --git a/src/maxwell_cpu.cpp b/src/maxwell_cpu.cpp
index 0491151..2af0267 100644
--- a/src/maxwell_cpu.cpp
+++ b/src/maxwell_cpu.cpp
@@ -609,7 +609,7 @@ gather_field(const model& mod, maxwell::field& in, maxwell::field& out,
auto& e = mod.entity_world_lookup(tag);
auto base = e.dof_base_world();
auto len = e.num_dofs();
- std::cout << "GR: " << tag << " " << len << std::endl;
+ //std::cout << "GR: " << tag << " " << len << std::endl;
receive_field(out, mod.entity_rank(tag), tag, base, len, comm);
}
@@ -630,14 +630,14 @@ gather_field(const model& mod, maxwell::field& in, maxwell::field& out,
{
for (auto& e : mod)
{
- std::cout << "GS: " << e.gmsh_tag() << " " << e.num_dofs() << std::endl;
+ //std::cout << "GS: " << e.gmsh_tag() << " " << e.num_dofs() << std::endl;
send_field(in, root, e.gmsh_tag(), e.dof_base(), e.num_dofs(), comm);
}
}
}
void
-export_fields_to_silo_2(const model& mod, maxwell::solver_state& state,
+export_fields_to_silo(const model& mod, maxwell::solver_state& state,
const maxwell::parameter_loader& mpl)
{
int rank;
@@ -658,31 +658,47 @@ export_fields_to_silo_2(const model& mod, maxwell::solver_state& state,
sdb.import_mesh_from_gmsh();
sdb.write_mesh(state.curr_time, state.curr_timestep);
- std::vector<double> out_Ex( mod.num_cells_world() );
- std::vector<double> out_Ey( mod.num_cells_world() );
- std::vector<double> out_Ez( mod.num_cells_world() );
+ switch( mpl.postpro_fieldExportMode("E") )
+ {
+ case field_export_mode::NONE:
+ break;
+
+ case field_export_mode::NODAL:
+ export_E_field_nodal(mod, sdb, f.Ex, f.Ey, f.Ez, mpl);
+ break;
- for (auto& tag : mod.world_entities_tags())
+ case field_export_mode::ZONAL:
+ export_E_field_zonal(mod, sdb, f.Ex, f.Ey, f.Ez, mpl);
+ break;
+ }
+
+ switch( mpl.postpro_fieldExportMode("H") )
{
- auto& e = mod.entity_world_lookup(tag);
- for (size_t iT = 0; iT < e.num_cells(); iT++)
- {
- auto& pe = e.cell(iT);
- auto& re = e.cell_refelem(pe);
- vecxd phi_bar = re.basis_functions({1./3., 1./3., 1./3.});
- auto num_bf = re.num_basis_functions();
- auto ofs = e.cell_world_dof_offset(iT);
- auto gi = e.cell_world_index_by_gmsh(iT);
- out_Ex[gi] = f.Ex.segment(ofs, num_bf).dot(phi_bar);
- out_Ey[gi] = f.Ey.segment(ofs, num_bf).dot(phi_bar);
- out_Ez[gi] = f.Ez.segment(ofs, num_bf).dot(phi_bar);
- }
+ case field_export_mode::NONE:
+ break;
+
+ case field_export_mode::NODAL:
+ export_H_field_nodal(mod, sdb, f.Hx, f.Hy, f.Hz, mpl);
+ break;
+
+ case field_export_mode::ZONAL:
+ export_H_field_zonal(mod, sdb, f.Hx, f.Hy, f.Hz, mpl);
+ break;
}
- sdb.write_zonal_variable("Ex", out_Ex);
- sdb.write_zonal_variable("Ey", out_Ey);
- sdb.write_zonal_variable("Ez", out_Ez);
- sdb.write_vector_expression("E", "{Ex, Ey, Ez}");
+ switch( mpl.postpro_fieldExportMode("J") )
+ {
+ case field_export_mode::NONE:
+ break;
+
+ case field_export_mode::NODAL:
+ export_J_field_nodal(mod, sdb, f.Ex, f.Ey, f.Ez, mpl);
+ break;
+
+ case field_export_mode::ZONAL:
+ export_J_field_zonal(mod, sdb, f.Ex, f.Ey, f.Ez, mpl);
+ break;
+ }
}
else
{
@@ -690,21 +706,22 @@ export_fields_to_silo_2(const model& mod, maxwell::solver_state& state,
gather_field(mod, state.emf_next, dummy, 0, MPI_COMM_WORLD);
}
}
-#endif /* USE_MPI */
+
+#else /* USE_MPI */
void
-export_fields_to_silo(const model& mod, const maxwell::solver_state& state,
- const maxwell::parameter_loader& mpl)
+export_fields_to_silo(const model &mod, maxwell::solver_state &state,
+ const maxwell::parameter_loader &mpl)
{
std::stringstream ss;
ss << mpl.sim_name() << "/timestep_" << state.curr_timestep << ".silo";
-
+
silo sdb;
sdb.create_db(ss.str());
sdb.import_mesh_from_gmsh();
sdb.write_mesh(state.curr_time, state.curr_timestep);
- switch( mpl.postpro_fieldExportMode("E") )
+ switch (mpl.postpro_fieldExportMode("E"))
{
case field_export_mode::NONE:
break;
@@ -720,7 +737,7 @@ export_fields_to_silo(const model& mod, const maxwell::solver_state& state,
break;
}
- switch( mpl.postpro_fieldExportMode("H") )
+ switch (mpl.postpro_fieldExportMode("H"))
{
case field_export_mode::NONE:
break;
@@ -736,7 +753,7 @@ export_fields_to_silo(const model& mod, const maxwell::solver_state& state,
break;
}
- switch( mpl.postpro_fieldExportMode("J") )
+ switch (mpl.postpro_fieldExportMode("J"))
{
case field_export_mode::NONE:
break;
@@ -752,5 +769,6 @@ export_fields_to_silo(const model& mod, const maxwell::solver_state& state,
break;
}
}
+#endif /* USE_MPI */
} // namespace maxwell
diff --git a/src/maxwell_gpu.cpp b/src/maxwell_gpu.cpp
index 00ff575..2314859 100644
--- a/src/maxwell_gpu.cpp
+++ b/src/maxwell_gpu.cpp
@@ -33,7 +33,7 @@ init_matparams(const model& mod, solver_state_gpu& state,
auto v_epsilon = mpl.epsilon(etag, bar);
auto v_mu = mpl.mu(etag, bar);
auto v_sigma = mpl.sigma(etag, bar);
- auto ofs = e.cell_global_dof_offset(iT);
+ auto ofs = e.cell_model_dof_offset(iT);
for (size_t iD = 0; iD < re.num_basis_functions(); iD++)
{
@@ -555,18 +555,21 @@ swap(solver_state_gpu& state, const parameter_loader& mpl)
std::swap(state.emf_curr, state.emf_next);
}
+#ifdef USE_MPI
void
gather_field(const model& mod, const maxwell::field_gpu& in, maxwell::field_gpu& f,
int root, MPI_Comm comm)
{}
void
-export_fields_to_silo_2(const model& mod, maxwell::solver_state_gpu& state,
+export_fields_to_silo(const model& mod, maxwell::solver_state_gpu& state,
const maxwell::parameter_loader& mpl)
{}
+#else /* USE_MPI */
+
void
-export_fields_to_silo(const model& mod, const maxwell::solver_state_gpu& state,
+export_fields_to_silo(const model& mod, maxwell::solver_state_gpu& state,
const maxwell::parameter_loader& mpl)
{
std::stringstream ss;
@@ -623,5 +626,6 @@ export_fields_to_silo(const model& mod, const maxwell::solver_state_gpu& state,
break;
}
}
+#endif /* USE_MPI */
} // namespace maxwell
diff --git a/src/maxwell_solver.cpp b/src/maxwell_solver.cpp
index a4d7480..ebe6249 100644
--- a/src/maxwell_solver.cpp
+++ b/src/maxwell_solver.cpp
@@ -68,6 +68,7 @@ void initialize_solver(const model& mod, State& state, const maxwell::parameter_
init_matparams(mod, state, mpl);
}
+#ifdef ENABLE_EXP_GEOMETRIC_DIODE
class integration_line
{
std::vector<size_t> element_tags;
@@ -187,6 +188,7 @@ reinit_materials(const model& mod, maxwell::solver_state& state, maxwell::parame
}
}
}
+#endif /* ENABLE_EXP_GEOMETRIC_DIODE */
template<typename State>
void test_it(const model& mod, State& state, maxwell::parameter_loader& mpl)
@@ -205,12 +207,12 @@ void test_it(const model& mod, State& state, maxwell::parameter_loader& mpl)
std::cout << "Time integrator: " << mpl.sim_timeIntegratorName() << std::endl;
mpl.call_initialization_callback();
-
+#ifdef ENABLE_EXP_GEOMETRIC_DIODE
auto relmat = [&](){ reinit_materials(mod, state, mpl); };
sol::state& lua = mpl.lua_state();
lua["relmat"] = relmat;
-
+#endif /* ENABLE_EXP_GEOMETRIC_DIODE */
timecounter tc;
tc.tic();
prepare_sources(mod, state, mpl);
@@ -220,13 +222,8 @@ void test_it(const model& mod, State& state, maxwell::parameter_loader& mpl)
timestep(state, mpl, ti);
do_sources(mod, state, mpl);
-#ifdef USE_MPI
- if ( (silo_output_rate != 0) and (i%silo_output_rate == 0) )
- export_fields_to_silo_2(mod, state, mpl);
-#else
- if ( (silo_output_rate != 0) and (i%silo_output_rate == 0) )
- export_fields_to_silo(mod, state, mpl);
-#endif /* USE_MPI */
+ if ( (silo_output_rate != 0) and (i%silo_output_rate == 0) )
+ export_fields_to_silo(mod, state, mpl);
swap(state, mpl);
@@ -243,6 +240,7 @@ void test_it(const model& mod, State& state, maxwell::parameter_loader& mpl)
void
register_lua_usertypes(maxwell::parameter_loader& mpl)
{
+#ifdef ENABLE_EXP_GEOMETRIC_DIODE
sol::state& lua = mpl.lua_state();
sol::usertype<point_3d> point_3d_ut = lua.new_usertype<point_3d>("point",
sol::constructors<point_3d(), point_3d(double, double, double)>()
@@ -252,18 +250,21 @@ register_lua_usertypes(maxwell::parameter_loader& mpl)
lua.new_usertype<integration_line>("integration_line",
sol::constructors<integration_line(const point_3d&, const point_3d&, size_t)>()
);
+#endif /* ENABLE_EXP_GEOMETRIC_DIODE */
}
void
register_lua_usertypes_bystate(maxwell::parameter_loader& mpl,
model& mod, maxwell::solver_state& state)
{
+#ifdef ENABLE_EXP_GEOMETRIC_DIODE
sol::state& lua = mpl.lua_state();
lua["internal"] = lua.create_table();
lua["internal"]["model"] = &mod;
lua["internal"]["state"] = &state;
lua["integrate_E"] = integrate_electric_field;
+#endif /*ENABLE_EXP_GEOMETRIC_DIODE */
}
#ifdef ENABLE_GPU_SOLVER
@@ -273,15 +274,15 @@ register_lua_usertypes_bystate(maxwell::parameter_loader& mpl,
{}
#endif
-int main(int argc, const char *argv[])
+int main(int argc, char *argv[])
{
#ifdef USE_MPI
MPI_Init(&argc, &argv);
- int rank, size;
- MPI_Comm_rank(MPI_COMM_WORLD, &rank);
- MPI_Comm_size(MPI_COMM_WORLD, &size);
- std::cout << "Running with MPI: " << rank+1 << "/" << size << std::endl;
+ int comm_rank, comm_size;
+ MPI_Comm_rank(MPI_COMM_WORLD, &comm_rank);
+ MPI_Comm_size(MPI_COMM_WORLD, &comm_size);
+ std::cout << "Running with MPI: " << comm_rank+1 << "/" << comm_size << std::endl;
#endif
if (argc != 2)
@@ -312,11 +313,11 @@ int main(int argc, const char *argv[])
model mod( mpl.sim_geomorder(), mpl.sim_approxorder() );
#ifdef USE_MPI
- if (rank == 0)
+ if (comm_rank == 0)
{
mod.generate_mesh();
- if (size > 1)
- mod.partition_mesh(size);
+ if (comm_size > 1)
+ mod.partition_mesh(comm_size);
mod.populate_from_gmsh();
}
--
GitLab