diff --git a/contrib/sgr b/contrib/sgr
index af618523a0c76d5ae4cebde9f5f07c702833d924..79794aae8124bbae25f7b2a13e0f675a37bd7b04 160000
--- a/contrib/sgr
+++ b/contrib/sgr
@@ -1 +1 @@
-Subproject commit af618523a0c76d5ae4cebde9f5f07c702833d924
+Subproject commit 79794aae8124bbae25f7b2a13e0f675a37bd7b04
diff --git a/doc/lua_api.md b/doc/lua_api.md
index c1cf993f9e35e96e749f8fc4d9cdae6a278884ab..9d40e5c71b40e7f115700cc0f07c189bb3825942 100644
--- a/doc/lua_api.md
+++ b/doc/lua_api.md
@@ -145,6 +145,7 @@ None.
### Debug and validation
The solver has some facilities to allow validation and comparison with analytical solutions. Those facilities are accessed via the `debug` table. In the current implementation the `debug` table is undefined by default, so it must be initialized explicitly in the configuration with a statement `debug = {}`. Possible members of the `debug` table are:
- `analytical_solution(tag, x, y, z, t)`: define the analytical solution of the problem. If defined, the numerical solution is compared with the analytical solution at each timestep in the whole domain. The function has to return 6 values, namely `Ex`, `Ey`, `Ez`, `Hx`, `Hy` and `Hz`.
+- `dump_cell_ranks` (bool): if true, and if the solver has MPI support compiled in, add to the Silo output a variable showing the cell-to-rank mapping.
### Postprocessing
The export of the fields E, H and J to Silo is controlled using the appropriate sections in the `postpro` table.
diff --git a/include/maxwell_common.h b/include/maxwell_common.h
index 83bbfe992c797f94d57a3317f25b417fb404dcdb..886e74642a326fe3c8cec53b67d0d8c2f1ce843d 100644
--- a/include/maxwell_common.h
+++ b/include/maxwell_common.h
@@ -165,24 +165,9 @@ eval_interface_sources(const model& mod, const parameter_loader& mpl,
}
}
-void export_E_field_nodal(const model&, silo&, const vecxd&, const vecxd&,
- const vecxd&, const parameter_loader&);
-
-void export_H_field_nodal(const model&, silo&, const vecxd&, const vecxd&,
- const vecxd&, const parameter_loader&);
-
-void export_J_field_nodal(const model&, silo&, const vecxd&, const vecxd&,
- const vecxd&, const parameter_loader&);
-
-void export_E_field_zonal(const model&, silo&, const vecxd&, const vecxd&,
- const vecxd&, const parameter_loader&);
-
-void export_H_field_zonal(const model&, silo&, const vecxd&, const vecxd&,
- const vecxd&, const parameter_loader&);
-
-void export_J_field_zonal(const model&, silo&, const vecxd&, const vecxd&,
- const vecxd&, const parameter_loader&);
-
+#ifdef USE_MPI
+std::vector<double> compute_cell_ranks(const model&);
+#endif /* USE_MPI */
using MaterialFunc = std::function<double(int, const point_3d&)>;
diff --git a/include/maxwell_interface.h b/include/maxwell_interface.h
index 3b24d81c8bb36b6d3aafbb578bd50920de5a69da..4b696dcb929ec212faec1e3e6d9a705def23bae4 100644
--- a/include/maxwell_interface.h
+++ b/include/maxwell_interface.h
@@ -76,6 +76,7 @@ public:
bool has_analytical_solution(void) const;
std::pair<vec3d, vec3d> eval_analytical_solution(int, const point_3d&, double t) const;
+ bool debug_dumpCellRanks(void) const;
boundary_condition boundary_type(int tag) const;
interface_condition interface_type(int tag) const;
diff --git a/src/maxwell_common.cpp b/src/maxwell_common.cpp
index 8c96b3249a35d9f171efa695b992b2c96e2da865..feb067e7c311ed687587202a3505a6f182f7b702 100644
--- a/src/maxwell_common.cpp
+++ b/src/maxwell_common.cpp
@@ -175,6 +175,27 @@ init_lax_milgram(const model& mod, const parameter_loader& mpl,
}
}
+#ifdef USE_MPI
+std::vector<double>
+compute_cell_ranks(const model& mod)
+{
+ /* This is only for debug purposes: make a vector mapping from the
+ * cell number to the rank of the process to which it is assigned */
+ std::vector<double> ret( mod.num_cells_world(), 999999999 );
+ for (auto& etag : mod.world_entities_tags())
+ {
+ auto& e = mod.entity_world_lookup(etag);
+ int er = mod.entity_rank(etag);
+
+ for (size_t iT = 0; iT < e.num_cells(); iT++)
+ {
+ auto gi = e.cell_world_index_by_gmsh(iT);
+ ret[gi] = er;
+ }
+ }
+ return ret;
+}
+#endif /* USE_MPI */
static void
export_vector_field_nodal(const model& mod, silo& sdb, const vecxd& Fx,
@@ -316,256 +337,4 @@ export_vector_field(const model& mod, silo& sdb, const vecxd& Fx, const vecxd& F
}
}
-
-
-void export_E_field_nodal(const model& mod, silo& sdb, const vecxd& Ex,
- const vecxd& Ey, const vecxd& Ez, const parameter_loader&)
-{
- std::vector<double> accum_Ex( sdb.num_nodes(), 0.0 );
- std::vector<double> accum_Ey( sdb.num_nodes(), 0.0 );
- 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_world_dof_offset(iT);
- auto nodetags = pe.node_tags();
- for (size_t iD = 0; iD < 4; iD++)
- {
- size_t tag = nodetags[iD];
- size_t svofs = sdb.node_tag_offset(tag);
- accum_Ex[svofs] += Ex[cgofs+iD];
- accum_Ey[svofs] += Ey[cgofs+iD];
- accum_Ez[svofs] += Ez[cgofs+iD];
- count[svofs]++;
- }
- }
- }
-
- for (size_t i = 0; i < sdb.num_nodes(); i++)
- {
- accum_Ex[i] /= count[i];
- accum_Ey[i] /= count[i];
- accum_Ez[i] /= count[i];
- }
-
- sdb.write_nodal_variable("Ex", accum_Ex);
- sdb.write_nodal_variable("Ey", accum_Ey);
- sdb.write_nodal_variable("Ez", accum_Ez);
- sdb.write_vector_expression("E", "{Ex, Ey, Ez}");
-}
-
-void export_H_field_nodal(const model& mod, silo& sdb, const vecxd& Hx,
- const vecxd& Hy, const vecxd& Hz, const parameter_loader&)
-{
- std::vector<double> accum_Hx( sdb.num_nodes(), 0.0 );
- std::vector<double> accum_Hy( sdb.num_nodes(), 0.0 );
- 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_world_dof_offset(iT);
- auto nodetags = pe.node_tags();
- for (size_t iD = 0; iD < 4; iD++)
- {
- size_t tag = nodetags[iD];
- size_t svofs = sdb.node_tag_offset(tag);
- accum_Hx[svofs] += Hx[cgofs+iD];
- accum_Hy[svofs] += Hy[cgofs+iD];
- accum_Hz[svofs] += Hz[cgofs+iD];
- count[svofs]++;
- }
- }
- }
-
- for (size_t i = 0; i < sdb.num_nodes(); i++)
- {
- accum_Hx[i] /= count[i];
- accum_Hy[i] /= count[i];
- accum_Hz[i] /= count[i];
- }
-
- sdb.write_nodal_variable("Hx", accum_Hx);
- sdb.write_nodal_variable("Hy", accum_Hy);
- sdb.write_nodal_variable("Hz", accum_Hz);
- sdb.write_vector_expression("H", "{Hx, Hy, Hz}");
-}
-
-void export_J_field_nodal(const model& mod, silo& sdb, const vecxd& Ex,
- const vecxd& Ey, const vecxd& Ez, const parameter_loader& mpl)
-{
- std::vector<double> accum_Jx( sdb.num_nodes(), 0.0 );
- std::vector<double> accum_Jy( sdb.num_nodes(), 0.0 );
- 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)
- {
-#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(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++)
- {
- size_t tag = nodetags[iD];
- size_t svofs = sdb.node_tag_offset(tag);
- accum_Jx[svofs] += sigma*Ex[cgofs+iD];
- accum_Jy[svofs] += sigma*Ey[cgofs+iD];
- accum_Jz[svofs] += sigma*Ez[cgofs+iD];
- count[svofs]++;
- }
- }
- }
-
- for (size_t i = 0; i < sdb.num_nodes(); i++)
- {
- accum_Jx[i] /= count[i];
- accum_Jy[i] /= count[i];
- accum_Jz[i] /= count[i];
- }
-
- sdb.write_nodal_variable("Jx", accum_Jx);
- sdb.write_nodal_variable("Jy", accum_Jy);
- sdb.write_nodal_variable("Jz", accum_Jz);
- sdb.write_vector_expression("J", "{Jx, Jy, Jz}");
-}
-
-
-
-void export_E_field_zonal(const model& mod, silo& sdb, const vecxd& Ex,
- const vecxd& Ey, const vecxd& Ez, const parameter_loader&)
-{
- 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() );
-
-#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& 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] = 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);
- }
- }
-
- 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}");
-}
-
-void export_H_field_zonal(const model& mod, silo& sdb, const vecxd& Hx,
- const vecxd& Hy, const vecxd& Hz, const parameter_loader&)
-{
- std::vector<double> out_Hx( mod.num_cells_world() );
- std::vector<double> out_Hy( mod.num_cells_world() );
- std::vector<double> out_Hz( mod.num_cells_world() );
-
-#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& 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_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);
- }
- }
-
- sdb.write_zonal_variable("Hx", out_Hx);
- sdb.write_zonal_variable("Hy", out_Hy);
- sdb.write_zonal_variable("Hz", out_Hz);
- sdb.write_vector_expression("H", "{Hx, Hy, Hz}");
-}
-
-void export_J_field_zonal(const model& mod, silo& sdb, const vecxd& Ex,
- const vecxd& Ey, const vecxd& Ez, const parameter_loader& mpl)
-{
- std::vector<double> out_Jx( mod.num_cells_world() );
- std::vector<double> out_Jy( mod.num_cells_world() );
- std::vector<double> out_Jz( mod.num_cells_world() );
-
-#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& re = e.cell_refelem(pe);
- auto bar = pe.barycenter();
- 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_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);
- }
- }
-
- sdb.write_zonal_variable("Jx", out_Jx);
- sdb.write_zonal_variable("Jy", out_Jy);
- sdb.write_zonal_variable("Jz", out_Jz);
- sdb.write_vector_expression("J", "{Jx, Jy, Jz}");
-}
-
-
} // namespace maxwell
diff --git a/src/maxwell_cpu.cpp b/src/maxwell_cpu.cpp
index 8c5c807ceb8adb5b36c5297c6759042677a31ed6..98aa9d1056bfb99ea3a91a4852cbd3dbe92fe363 100644
--- a/src/maxwell_cpu.cpp
+++ b/src/maxwell_cpu.cpp
@@ -9,6 +9,9 @@
namespace maxwell {
+/**************************************************************************
+ * Solver initialization
+ */
void
init_matparams(const model& mod, solver_state& state,
const parameter_loader& mpl)
@@ -49,10 +52,6 @@ init_matparams(const model& mod, solver_state& state,
init_boundary_conditions(mod, mpl, state.matparams.bc_coeffs);
}
-
-
-
-
void
init_from_model(const model& mod, solver_state& state)
{
@@ -97,6 +96,10 @@ init_from_model(const model& mod, solver_state& state)
state.curr_timestep = 0;
}
+
+/**************************************************************************
+ * Computational kernels: curls
+ */
static void
compute_curls(solver_state& state, const field& curr, field& next)
{
@@ -149,6 +152,9 @@ compute_curls_H(solver_state& state, const field& curr, field& next)
}
+/**************************************************************************
+ * Computational kernels: jumps
+ */
static void
compute_field_jumps(solver_state& state, const field& in)
{
@@ -406,6 +412,9 @@ exchange_boundary_data_H(const model& mod, solver_state& state)
}
#endif /* USE_MPI */
+/**************************************************************************
+ * Computational kernels: fluxes
+ */
static void
compute_fluxes_planar(solver_state& state)
{
@@ -538,6 +547,9 @@ compute_fluxes_planar_H(solver_state& state)
}
}
+/**************************************************************************
+ * Computational kernels: boundary terms
+ */
static void
compute_fluxes(const model& mod, solver_state& state, const field& in, field& out)
{
@@ -598,6 +610,9 @@ compute_fluxes_H(const model& mod, solver_state& state, const field& in, field&
}
}
+/**************************************************************************
+ * Computational kernels: time integration
+ */
static void
compute_euler_update(solver_state& state, const field& y, const field& k, double dt, field& out)
{
@@ -640,6 +655,9 @@ compute_rk4_weighted_sum(solver_state& state, const field& in, double dt, field&
}
}
+/**************************************************************************
+ * Computational kernels: H&W DG operator
+ */
static void
apply_operator(const model& mod, solver_state& state, const field& curr, field& next)
{
@@ -647,6 +665,9 @@ apply_operator(const model& mod, solver_state& state, const field& curr, field&
compute_fluxes(mod, state, curr, next);
}
+/**************************************************************************
+ * Computational kernels: Leapfrog DG
+ */
static void
leapfrog(const model& mod, solver_state& state)
{
@@ -675,6 +696,9 @@ leapfrog(const model& mod, solver_state& state)
}
}
+/**************************************************************************
+ * Computational kernels: Run a whole timestep
+ */
void
timestep(const model& mod, solver_state& state, const parameter_loader&, time_integrator_type ti)
{
@@ -709,6 +733,9 @@ timestep(const model& mod, solver_state& state, const parameter_loader&, time_in
state.curr_timestep += 1;
}
+/**************************************************************************
+ * Sources: volumetric
+ */
void
eval_volume_sources(const model& mod, const parameter_loader& mpl, solver_state& state)
{
@@ -842,104 +869,11 @@ export_fields_to_silo(const model &mod, maxwell::solver_state &state,
auto J_exportmode = mpl.postpro_fieldExportMode("J");
export_vector_field(mod, sdb, f.Ex, f.Ey, f.Ez, sigma, "J", J_exportmode);
- std::vector<double> entranks( mod.num_cells_world(), 999999999 );
- for (auto& etag : mod.world_entities_tags())
- {
- auto& e = mod.entity_world_lookup(etag);
- int er = mod.entity_rank(etag);
-
- for (size_t iT = 0; iT < e.num_cells(); iT++)
- {
- auto gi = e.cell_world_index_by_gmsh(iT);
- entranks[gi] = er;
- }
- }
- sdb.write_zonal_variable("entity_ranks", entranks);
- }
- else
- {
- maxwell::field dummy;
- gather_field(mod, state.emf_next, dummy, 0, MPI_COMM_WORLD);
- }
-}
-
-void
-export_fields_to_silo(const model& mod, maxwell::solver_state& state,
- const maxwell::parameter_loader& mpl)
-{
- int rank;
- MPI_Comm_rank(MPI_COMM_WORLD, &rank);
-
- if (rank == 0)
- {
- size_t gdnum = mod.num_dofs_world();
- maxwell::field f;
- f.resize(gdnum);
- gather_field(mod, state.emf_next, f, 0, MPI_COMM_WORLD);
-
- 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") )
+ if ( mpl.debug_dumpCellRanks() )
{
- 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;
-
- case field_export_mode::ZONAL:
- export_E_field_zonal(mod, sdb, f.Ex, f.Ey, f.Ez, mpl);
- break;
- }
-
- switch( mpl.postpro_fieldExportMode("H") )
- {
- 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;
+ auto entranks = compute_cell_ranks(mod);
+ sdb.write_zonal_variable("entity_ranks", entranks);
}
-
- 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;
- }
-
- std::vector<double> entranks( mod.num_cells_world(), 999999999 );
- for (auto& etag : mod.world_entities_tags())
- {
- auto& e = mod.entity_world_lookup(etag);
- int er = mod.entity_rank(etag);
-
- for (size_t iT = 0; iT < e.num_cells(); iT++)
- {
- auto gi = e.cell_world_index_by_gmsh(iT);
- entranks[gi] = er;
- }
- }
- sdb.write_zonal_variable("entity_ranks", entranks);
}
else
{
@@ -980,67 +914,6 @@ export_fields_to_silo(const model &mod, maxwell::solver_state &state,
state.emf_next.Ez, sigma, "J", J_exportmode);
}
-
-void
-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"))
- {
- case field_export_mode::NONE:
- break;
-
- case field_export_mode::NODAL:
- export_E_field_nodal(mod, sdb, state.emf_next.Ex, state.emf_next.Ey,
- state.emf_next.Ez, mpl);
- break;
-
- case field_export_mode::ZONAL:
- export_E_field_zonal(mod, sdb, state.emf_next.Ex, state.emf_next.Ey,
- state.emf_next.Ez, mpl);
- break;
- }
-
- switch (mpl.postpro_fieldExportMode("H"))
- {
- case field_export_mode::NONE:
- break;
-
- case field_export_mode::NODAL:
- export_H_field_nodal(mod, sdb, state.emf_next.Hx, state.emf_next.Hy,
- state.emf_next.Hz, mpl);
- break;
-
- case field_export_mode::ZONAL:
- export_H_field_zonal(mod, sdb, state.emf_next.Hx, state.emf_next.Hy,
- state.emf_next.Hz, mpl);
- break;
- }
-
- switch (mpl.postpro_fieldExportMode("J"))
- {
- case field_export_mode::NONE:
- break;
-
- case field_export_mode::NODAL:
- export_J_field_nodal(mod, sdb, state.emf_next.Ex, state.emf_next.Ey,
- state.emf_next.Ez, mpl);
- break;
-
- case field_export_mode::ZONAL:
- export_J_field_zonal(mod, sdb, state.emf_next.Ex, state.emf_next.Ey,
- state.emf_next.Ez, mpl);
- break;
- }
-}
#endif /* USE_MPI */
} // namespace maxwell
diff --git a/src/maxwell_gpu.cpp b/src/maxwell_gpu.cpp
index 6756b64bbe831f2cb9c5d12a7dd456aa420c8357..b1d124803e5eea41b6ed7b786ea2e7d9169f42fa 100644
--- a/src/maxwell_gpu.cpp
+++ b/src/maxwell_gpu.cpp
@@ -570,11 +570,14 @@ export_fields_to_silo(const model& mod, maxwell::solver_state_gpu& state,
void
export_fields_to_silo(const model& mod, maxwell::solver_state_gpu& state,
- const maxwell::parameter_loader& mpl)
+ const maxwell::parameter_loader& mpl, std::string basename)
{
+ if (basename == "")
+ basename = "timestep";
+
std::stringstream ss;
- ss << mpl.sim_name() << "/timestep_" << state.curr_timestep << ".silo";
-
+ ss << mpl.sim_name() << "/" << basename << "_" << state.curr_timestep << ".silo";
+
silo sdb;
sdb.create_db(ss.str());
sdb.import_mesh_from_gmsh();
@@ -584,47 +587,16 @@ export_fields_to_silo(const model& mod, maxwell::solver_state_gpu& state,
f.resize( mod.num_dofs() );
state.emf_next.copyout(f);
- 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;
- case field_export_mode::ZONAL:
- export_E_field_zonal(mod, sdb, f.Ex, f.Ey, f.Ez, mpl);
- break;
- }
+ auto E_exportmode = mpl.postpro_fieldExportMode("E");
+ export_vector_field(mod, sdb, f.Ex, f.Ey, f.Ez, "E", E_exportmode);
- switch( mpl.postpro_fieldExportMode("H") )
- {
- case field_export_mode::NONE:
- break;
+ auto H_exportmode = mpl.postpro_fieldExportMode("H");
+ export_vector_field(mod, sdb, f.Hx, f.Hy, f.Hz, "H", H_exportmode);
- 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;
- }
-
- 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;
- }
+ 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, f.Ex, f.Ey, f.Ez, sigma, "J", J_exportmode);
}
#endif /* USE_MPI */
diff --git a/src/maxwell_interface.cpp b/src/maxwell_interface.cpp
index 9bbaa66246b97c1199e5038b7040ff06ed184fbb..8b975b216507297c2efa05edc870b0db0fa5ad2b 100644
--- a/src/maxwell_interface.cpp
+++ b/src/maxwell_interface.cpp
@@ -706,4 +706,14 @@ parameter_loader::eval_analytical_solution(int tag, const point_3d& pt, double t
return std::make_pair(E,H);
}
+bool
+parameter_loader::debug_dumpCellRanks(void) const
+{
+ auto dcr = lua["debug"]["dump_cell_ranks"];
+ if (not dcr.valid())
+ return false;
+
+ return (dcr == true);
+}
+
} // namespace maxwell
diff --git a/src/maxwell_solver.cpp b/src/maxwell_solver.cpp
index 1dca331dec75864435b614612f21785a49d8d89e..63b92d5da852566d28934b8eb34b9f03aaeda18a 100644
--- a/src/maxwell_solver.cpp
+++ b/src/maxwell_solver.cpp
@@ -21,6 +21,9 @@
#include "param_loader.h"
#include "timecounter.h"
+#include "sgr.hpp"
+using namespace sgr;
+
#if 0
static void
make_geometry(int order, double mesh_h)
@@ -227,8 +230,14 @@ maxwell_6ple sqrt(const maxwell_6ple& m)
return ret;
}
+double hsum(const maxwell_6ple& m)
+{
+ return m.Ex + m.Ey + m.Ez + m.Hx + m.Hy + m.Hz;
+}
+
static void
-validate(const model& mod, maxwell::solver_state& state, maxwell::parameter_loader& mpl, std::ofstream& ofs)
+validate(const model& mod, maxwell::solver_state& state,
+ maxwell::parameter_loader& mpl, std::ofstream& ofs)
{
maxwell_6ple err;
maxwell_6ple energy;
@@ -255,7 +264,8 @@ validate(const model& mod, maxwell::solver_state& state, maxwell::parameter_load
{
const auto& pqp = pqps[iQp];
const vecxd phi = re.basis_functions(iQp);
- auto [E, H] = mpl.eval_analytical_solution(e.material_tag(), pqp.point(), state.curr_time);
+ auto [E, H] = mpl.eval_analytical_solution(e.material_tag(),
+ pqp.point(), state.curr_time);
maxwell_6ple Fh;
Fh.Ex = state.emf_next.Ex.segment(ofs, cbs).dot(phi);
@@ -283,15 +293,15 @@ validate(const model& mod, maxwell::solver_state& state, maxwell::parameter_load
}
}
- ofs << energy.Ex << " " << energy.Ey << " " << energy.Ez << " " << energy.Hx;
- ofs << " " << energy.Hy << " " << energy.Hz << " ";
- ofs << energy.Ex + energy.Ey + energy.Ez + energy.Hx + energy.Hy + energy.Hz;// << std::endl;
- ofs << " " << sqrt(err) << std::endl;
+ ofs << energy << " " << hsum(energy) << " " << sqrt(err) << std::endl;
}
-//static void
-//validate(const model& mod, maxwell::solver_state_gpu& state, maxwell::parameter_loader& mpl)
-//{}
+#ifdef ENABLE_GPU_SOLVER
+static void
+validate(const model& mod, maxwell::solver_state_gpu& state,
+ maxwell::parameter_loader& mpl, std::ofstream&)
+{}
+#endif /* ENABLE_GPU_SOLVER */
template<typename State>
void solver_mainloop(const model& mod, State& state, maxwell::parameter_loader& mpl)
@@ -327,6 +337,9 @@ void solver_mainloop(const model& mod, State& state, maxwell::parameter_loader&
mpl.call_initialization_callback();
#ifdef ENABLE_EXP_GEOMETRIC_DIODE
+ #ifdef USE_MPI
+ #error "The geometric diode code is not compatible with MPI yet."
+ #endif
auto relmat = [&](){ reinit_materials(mod, state, mpl); };
sol::state& lua = mpl.lua_state();
@@ -358,16 +371,19 @@ void solver_mainloop(const model& mod, State& state, maxwell::parameter_loader&
MPI_Reduce(&dofs_s_proc, &dofs_s_total, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
if (proc_rank == 0)
{
- std::cout << "Cycle " << i << ": t = " << state.curr_time << " s";
- std::cout << ", DOFs/s: " << dofs_s_total << std::endl;
+ std::cout << hidecursor << clrline << "Cycle " << i << ": t = ";
+ std::cout << state.curr_time << " s" << ", DOFs/s: ";
+ std::cout << dofs_s_total << cr << std::flush;
}
#else /* USE_MPI */
- std::cout << "Cycle " << i << ": t = " << state.curr_time << " s";
- std::cout << ", DOFs/s: " << dofs_s_proc << std::endl;
+ std::cout << hidecursor << clrline << "Cycle " << i << ": t = ";
+ std::cout << state.curr_time << " s" << ", DOFs/s: ";
+ std::cout << dofs_s_proc << cr << std::flush;
#endif /* USE_MPI */
tc.tic();
}
}
+ std::cout << showcursor << std::endl;
}
static void
@@ -483,6 +499,9 @@ int main(int argc, char *argv[])
if ( mpl.sim_usegpu() )
{
#ifdef ENABLE_GPU_SOLVER
+ #ifdef USE_MPI
+ #error "The GPU solver is not compatible with MPI yet."
+ #endif
maxwell::solver_state_gpu state_g;
register_lua_usertypes_bystate(mpl, mod, state_g);
solver_mainloop(mod, state_g, mpl);