Skip to content
Snippets Groups Projects
Commit f86e9cce authored by Matteo Cicuttin's avatar Matteo Cicuttin
Browse files

Got rid of the old Silo epxort code.

parent 9284b9e5
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
sgr @ 79794aae
Compare af618523...79794aae
Subproject commit af618523a0c76d5ae4cebde9f5f07c702833d924 Subproject commit 79794aae8124bbae25f7b2a13e0f675a37bd7b04
doc/lua_api.md
+ 1
− 0
View file @ f86e9cce
...@@ -145,6 +145,7 @@ None. ...@@ -145,6 +145,7 @@ None.
### Debug and validation ### 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: 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`. - `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 ### Postprocessing
The export of the fields E, H and J to Silo is controlled using the appropriate sections in the `postpro` table. The export of the fields E, H and J to Silo is controlled using the appropriate sections in the `postpro` table.
......
include/maxwell_common.h
+ 3
− 18
View file @ f86e9cce
...@@ -165,24 +165,9 @@ eval_interface_sources(const model& mod, const parameter_loader& mpl, ...@@ -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&, #ifdef USE_MPI
const vecxd&, const parameter_loader&); std::vector<double> compute_cell_ranks(const model&);
#endif /* USE_MPI */
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&);
using MaterialFunc = std::function<double(int, const point_3d&)>; using MaterialFunc = std::function<double(int, const point_3d&)>;
......
include/maxwell_interface.h
+ 1
− 0
View file @ f86e9cce
...@@ -76,6 +76,7 @@ public: ...@@ -76,6 +76,7 @@ public:
bool has_analytical_solution(void) const; bool has_analytical_solution(void) const;
std::pair<vec3d, vec3d> eval_analytical_solution(int, const point_3d&, double t) 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; boundary_condition boundary_type(int tag) const;
interface_condition interface_type(int tag) const; interface_condition interface_type(int tag) const;
......
src/maxwell_common.cpp
+ 21
− 252
View file @ f86e9cce
...@@ -175,6 +175,27 @@ init_lax_milgram(const model& mod, const parameter_loader& mpl, ...@@ -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 static void
export_vector_field_nodal(const model& mod, silo& sdb, const vecxd& Fx, 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 ...@@ -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 } // namespace maxwell
src/maxwell_cpu.cpp
+ 34
− 161
View file @ f86e9cce
...@@ -9,6 +9,9 @@ ...@@ -9,6 +9,9 @@
namespace maxwell { namespace maxwell {
/**************************************************************************
* Solver initialization
*/
void void
init_matparams(const model& mod, solver_state& state, init_matparams(const model& mod, solver_state& state,
const parameter_loader& mpl) const parameter_loader& mpl)
...@@ -49,10 +52,6 @@ init_matparams(const model& mod, solver_state& state, ...@@ -49,10 +52,6 @@ init_matparams(const model& mod, solver_state& state,
init_boundary_conditions(mod, mpl, state.matparams.bc_coeffs); init_boundary_conditions(mod, mpl, state.matparams.bc_coeffs);
} }
void void
init_from_model(const model& mod, solver_state& state) init_from_model(const model& mod, solver_state& state)
{ {
...@@ -97,6 +96,10 @@ 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; state.curr_timestep = 0;
} }
/**************************************************************************
* Computational kernels: curls
*/
static void static void
compute_curls(solver_state& state, const field& curr, field& next) 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) ...@@ -149,6 +152,9 @@ compute_curls_H(solver_state& state, const field& curr, field& next)
} }
/**************************************************************************
* Computational kernels: jumps
*/
static void static void
compute_field_jumps(solver_state& state, const field& in) compute_field_jumps(solver_state& state, const field& in)
{ {
...@@ -406,6 +412,9 @@ exchange_boundary_data_H(const model& mod, solver_state& state) ...@@ -406,6 +412,9 @@ exchange_boundary_data_H(const model& mod, solver_state& state)
} }
#endif /* USE_MPI */ #endif /* USE_MPI */
/**************************************************************************
* Computational kernels: fluxes
*/
static void static void
compute_fluxes_planar(solver_state& state) compute_fluxes_planar(solver_state& state)
{ {
...@@ -538,6 +547,9 @@ compute_fluxes_planar_H(solver_state& state) ...@@ -538,6 +547,9 @@ compute_fluxes_planar_H(solver_state& state)
} }
} }
/**************************************************************************
* Computational kernels: boundary terms
*/
static void static void
compute_fluxes(const model& mod, solver_state& state, const field& in, field& out) 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& ...@@ -598,6 +610,9 @@ compute_fluxes_H(const model& mod, solver_state& state, const field& in, field&
} }
} }
/**************************************************************************
* Computational kernels: time integration
*/
static void static void
compute_euler_update(solver_state& state, const field& y, const field& k, double dt, field& out) 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& ...@@ -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 static void
apply_operator(const model& mod, solver_state& state, const field& curr, field& next) 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& ...@@ -647,6 +665,9 @@ apply_operator(const model& mod, solver_state& state, const field& curr, field&
compute_fluxes(mod, state, curr, next); compute_fluxes(mod, state, curr, next);
} }
/**************************************************************************
* Computational kernels: Leapfrog DG
*/
static void static void
leapfrog(const model& mod, solver_state& state) leapfrog(const model& mod, solver_state& state)
{ {
...@@ -675,6 +696,9 @@ 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 void
timestep(const model& mod, solver_state& state, const parameter_loader&, time_integrator_type ti) 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 ...@@ -709,6 +733,9 @@ timestep(const model& mod, solver_state& state, const parameter_loader&, time_in
state.curr_timestep += 1; state.curr_timestep += 1;
} }
/**************************************************************************
* Sources: volumetric
*/
void void
eval_volume_sources(const model& mod, const parameter_loader& mpl, solver_state& state) 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, ...@@ -842,104 +869,11 @@ export_fields_to_silo(const model &mod, maxwell::solver_state &state,
auto J_exportmode = mpl.postpro_fieldExportMode("J"); auto J_exportmode = mpl.postpro_fieldExportMode("J");
export_vector_field(mod, sdb, f.Ex, f.Ey, f.Ez, sigma, "J", J_exportmode); export_vector_field(mod, sdb, f.Ex, f.Ey, f.Ez, sigma, "J", J_exportmode);
std::vector<double> entranks( mod.num_cells_world(), 999999999 ); if ( mpl.debug_dumpCellRanks() )
for (auto& etag : mod.world_entities_tags())
{ {
auto& e = mod.entity_world_lookup(etag); auto entranks = compute_cell_ranks(mod);
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); 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") )
{
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;
}
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 else
{ {
...@@ -980,67 +914,6 @@ export_fields_to_silo(const model &mod, maxwell::solver_state &state, ...@@ -980,67 +914,6 @@ export_fields_to_silo(const model &mod, maxwell::solver_state &state,
state.emf_next.Ez, sigma, "J", J_exportmode); 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 */ #endif /* USE_MPI */
} // namespace maxwell } // namespace maxwell
src/maxwell_gpu.cpp
+ 13
− 41
View file @ f86e9cce
...@@ -570,10 +570,13 @@ export_fields_to_silo(const model& mod, maxwell::solver_state_gpu& state, ...@@ -570,10 +570,13 @@ export_fields_to_silo(const model& mod, maxwell::solver_state_gpu& state,
void void
export_fields_to_silo(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) const maxwell::parameter_loader& mpl, std::string basename)
{ {
if (basename == "")
basename = "timestep";
std::stringstream ss; std::stringstream ss;
ss << mpl.sim_name() << "/timestep_" << state.curr_timestep << ".silo"; ss << mpl.sim_name() << "/" << basename << "_" << state.curr_timestep << ".silo";
silo sdb; silo sdb;
sdb.create_db(ss.str()); sdb.create_db(ss.str());
...@@ -584,47 +587,16 @@ export_fields_to_silo(const model& mod, maxwell::solver_state_gpu& state, ...@@ -584,47 +587,16 @@ export_fields_to_silo(const model& mod, maxwell::solver_state_gpu& state,
f.resize( mod.num_dofs() ); f.resize( mod.num_dofs() );
state.emf_next.copyout(f); 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: auto E_exportmode = mpl.postpro_fieldExportMode("E");
export_E_field_zonal(mod, sdb, f.Ex, f.Ey, f.Ez, mpl); export_vector_field(mod, sdb, f.Ex, f.Ey, f.Ez, "E", E_exportmode);
break;
}
switch( mpl.postpro_fieldExportMode("H") )
{
case field_export_mode::NONE:
break;
case field_export_mode::NODAL: auto H_exportmode = mpl.postpro_fieldExportMode("H");
export_H_field_nodal(mod, sdb, f.Hx, f.Hy, f.Hz, mpl); export_vector_field(mod, sdb, f.Hx, f.Hy, f.Hz, "H", H_exportmode);
break;
case field_export_mode::ZONAL: auto sigma = [&](int tag, const point_3d& pt) -> double { return mpl.sigma(tag, pt); };
export_H_field_zonal(mod, sdb, f.Hx, f.Hy, f.Hz, mpl); auto J_exportmode = mpl.postpro_fieldExportMode("J");
break; export_vector_field(mod, sdb, f.Ex, f.Ey, f.Ez, sigma, "J", J_exportmode);
}
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;
}
} }
#endif /* USE_MPI */ #endif /* USE_MPI */
......
src/maxwell_interface.cpp
+ 10
− 0
View file @ f86e9cce
...@@ -706,4 +706,14 @@ parameter_loader::eval_analytical_solution(int tag, const point_3d& pt, double t ...@@ -706,4 +706,14 @@ parameter_loader::eval_analytical_solution(int tag, const point_3d& pt, double t
return std::make_pair(E,H); 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 } // namespace maxwell
src/maxwell_solver.cpp
+ 32
− 13
View file @ f86e9cce
...@@ -21,6 +21,9 @@ ...@@ -21,6 +21,9 @@
#include "param_loader.h" #include "param_loader.h"
#include "timecounter.h" #include "timecounter.h"
#include "sgr.hpp"
using namespace sgr;
#if 0 #if 0
static void static void
make_geometry(int order, double mesh_h) make_geometry(int order, double mesh_h)
...@@ -227,8 +230,14 @@ maxwell_6ple sqrt(const maxwell_6ple& m) ...@@ -227,8 +230,14 @@ maxwell_6ple sqrt(const maxwell_6ple& m)
return ret; return ret;
} }
double hsum(const maxwell_6ple& m)
{
return m.Ex + m.Ey + m.Ez + m.Hx + m.Hy + m.Hz;
}
static void 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 err;
maxwell_6ple energy; maxwell_6ple energy;
...@@ -255,7 +264,8 @@ validate(const model& mod, maxwell::solver_state& state, maxwell::parameter_load ...@@ -255,7 +264,8 @@ validate(const model& mod, maxwell::solver_state& state, maxwell::parameter_load
{ {
const auto& pqp = pqps[iQp]; const auto& pqp = pqps[iQp];
const vecxd phi = re.basis_functions(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; maxwell_6ple Fh;
Fh.Ex = state.emf_next.Ex.segment(ofs, cbs).dot(phi); 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 ...@@ -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 << " " << hsum(energy) << " " << sqrt(err) << std::endl;
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;
} }
//static void #ifdef ENABLE_GPU_SOLVER
//validate(const model& mod, maxwell::solver_state_gpu& state, maxwell::parameter_loader& mpl) static void
//{} validate(const model& mod, maxwell::solver_state_gpu& state,
maxwell::parameter_loader& mpl, std::ofstream&)
{}
#endif /* ENABLE_GPU_SOLVER */
template<typename State> template<typename State>
void solver_mainloop(const model& mod, State& state, maxwell::parameter_loader& mpl) 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& ...@@ -327,6 +337,9 @@ void solver_mainloop(const model& mod, State& state, maxwell::parameter_loader&
mpl.call_initialization_callback(); mpl.call_initialization_callback();
#ifdef ENABLE_EXP_GEOMETRIC_DIODE #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); }; auto relmat = [&](){ reinit_materials(mod, state, mpl); };
sol::state& lua = mpl.lua_state(); sol::state& lua = mpl.lua_state();
...@@ -358,16 +371,19 @@ void solver_mainloop(const model& mod, State& state, maxwell::parameter_loader& ...@@ -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); MPI_Reduce(&dofs_s_proc, &dofs_s_total, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
if (proc_rank == 0) if (proc_rank == 0)
{ {
std::cout << "Cycle " << i << ": t = " << state.curr_time << " s"; std::cout << hidecursor << clrline << "Cycle " << i << ": t = ";
std::cout << ", DOFs/s: " << dofs_s_total << std::endl; std::cout << state.curr_time << " s" << ", DOFs/s: ";
std::cout << dofs_s_total << cr << std::flush;
} }
#else /* USE_MPI */ #else /* USE_MPI */
std::cout << "Cycle " << i << ": t = " << state.curr_time << " s"; std::cout << hidecursor << clrline << "Cycle " << i << ": t = ";
std::cout << ", DOFs/s: " << dofs_s_proc << std::endl; std::cout << state.curr_time << " s" << ", DOFs/s: ";
std::cout << dofs_s_proc << cr << std::flush;
#endif /* USE_MPI */ #endif /* USE_MPI */
tc.tic(); tc.tic();
} }
} }
std::cout << showcursor << std::endl;
} }
static void static void
...@@ -483,6 +499,9 @@ int main(int argc, char *argv[]) ...@@ -483,6 +499,9 @@ int main(int argc, char *argv[])
if ( mpl.sim_usegpu() ) if ( mpl.sim_usegpu() )
{ {
#ifdef ENABLE_GPU_SOLVER #ifdef ENABLE_GPU_SOLVER
#ifdef USE_MPI
#error "The GPU solver is not compatible with MPI yet."
#endif
maxwell::solver_state_gpu state_g; maxwell::solver_state_gpu state_g;
register_lua_usertypes_bystate(mpl, mod, state_g); register_lua_usertypes_bystate(mpl, mod, state_g);
solver_mainloop(mod, state_g, mpl); solver_mainloop(mod, state_g, mpl);
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment