From a1ca0ecf39267c32633655d175cb5aa3db91ae2a Mon Sep 17 00:00:00 2001
From: Matteo Cicuttin <datafl4sh@toxicnet.eu>
Date: Wed, 5 May 2021 17:47:05 +0200
Subject: [PATCH] Parametrizing the solver with Lua: initial conditions.
---
CMakeLists.txt | 1 +
include/connectivity.h | 8 ++
include/gmsh_io.h | 24 +++++
include/param_loader.h | 49 ++++++++++
src/gmsh_io.cpp | 106 ++++++++++++++++++---
src/maxwell_solver.cpp | 63 ++++++++-----
src/param_loader.cpp | 192 +++++++++++++++++++++++++++++++++++++++
src/physical_element.cpp | 3 +
8 files changed, 412 insertions(+), 34 deletions(-)
create mode 100644 include/param_loader.h
create mode 100644 src/param_loader.cpp
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 0c07aab..99ca31b 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -256,6 +256,7 @@ set(LIBGMSHDG_SOURCES src/gmsh_io.cpp
src/maxwell_cpu.cpp
src/maxwell_gpu.cpp
src/maxwell_kernels_cuda.cu
+ src/param_loader.cpp
)
add_library(gmshdg SHARED ${LIBGMSHDG_SOURCES})
diff --git a/include/connectivity.h b/include/connectivity.h
index 352d8eb..16be7ca 100644
--- a/include/connectivity.h
+++ b/include/connectivity.h
@@ -164,6 +164,14 @@ public:
return std::make_pair(neighbour_descriptor(), false); /* Silence the compiler */
}
+ size_t
+ num_neighbours(const FaceKey& fk)
+ {
+ size_t nn = face_cell_adj.at(fk).num_neighbours();
+ assert(nn == 1 or nn == 2);
+ return nn;
+ }
+
std::set<ConnData>
neighbours(size_t iE, size_t iT) const
{
diff --git a/include/gmsh_io.h b/include/gmsh_io.h
index c569808..24434c0 100644
--- a/include/gmsh_io.h
+++ b/include/gmsh_io.h
@@ -26,6 +26,23 @@ std::string basis_grad_name(int);
using face_key = std::array<size_t, 3>;
+enum class face_type
+{
+ NONE,
+ INTERFACE,
+ BOUNDARY
+};
+
+struct boundary_descriptor
+{
+ face_type type;
+ int gmsh_entity;
+
+ boundary_descriptor()
+ : type(face_type::NONE), gmsh_entity(0)
+ {}
+};
+
class model
{
int geometric_order;
@@ -33,6 +50,7 @@ class model
std::vector<entity> entities;
std::map<size_t, std::vector<element_key>> boundary_map;
+ std::vector<boundary_descriptor> bnd_descriptors;
element_connectivity<element_key> conn;
void map_boundaries(void);
@@ -54,6 +72,11 @@ public:
return conn;
}
+ const std::vector<boundary_descriptor>& boundary_descriptors() const
+ {
+ return bnd_descriptors;
+ }
+
std::vector<element_key> get_bnd(size_t which)
{
return boundary_map.at(which);
@@ -62,6 +85,7 @@ public:
size_t num_dofs() const;
size_t num_fluxes() const;
size_t num_cells() const;
+ size_t num_faces() const;
size_t num_entities() const;
std::vector<entity>::const_iterator begin() const;
diff --git a/include/param_loader.h b/include/param_loader.h
new file mode 100644
index 0000000..e0f15d9
--- /dev/null
+++ b/include/param_loader.h
@@ -0,0 +1,49 @@
+#pragma once
+
+#include "sol/sol.hpp"
+#include "gmsh_io.h"
+
+class parameter_loader
+{
+protected:
+ sol::state lua;
+
+private:
+ void init(void);
+ bool validate_simulation_params(void) const;
+
+public:
+ parameter_loader();
+
+ bool load_file(const std::string&);
+
+ int approx_order(void) const;
+ int geom_order(void) const;
+ std::string simulation_name(void) const;
+ double delta_t(void) const;
+ size_t timesteps(void) const;
+ std::string gmsh_model_filename(void) const;
+ bool use_gpu(void) const;
+ size_t silo_output_rate(void) const;
+};
+
+class maxwell_parameter_loader : public parameter_loader
+{
+ void init(void);
+
+ bool validate_materials(const std::string&, int);
+
+public:
+ maxwell_parameter_loader();
+
+ bool validate_model_params(const model&);
+ double epsilon(int, const point_3d&) const;
+ double mu(int, const point_3d&) const;
+ double sigma(int, const point_3d&) const;
+
+ bool initial_Efield_defined(void) const;
+ vec3d Efield(const point_3d&) const;
+
+ bool initial_Hfield_defined(void) const;
+ vec3d Hfield(const point_3d&) const;
+};
\ No newline at end of file
diff --git a/src/gmsh_io.cpp b/src/gmsh_io.cpp
index 24ec863..8ec1a0b 100644
--- a/src/gmsh_io.cpp
+++ b/src/gmsh_io.cpp
@@ -75,9 +75,31 @@ model::remesh()
{
gmm::generate( DIMENSION(3) );
gmm::setOrder( geometric_order );
+
+ gvp_t gmsh_entities;
+
+ /* Create a table mapping the tag of a boundary to the list of
+ * its face keys. This must happen before import_gmsh_entities()
+ * because otherwise you get spurious 2D entities (and this in
+ * turn because the constructor of entity calls addDiscreteEntity()
+ * to get all the element faces) .*/
+ std::map<size_t, std::vector<face_key>> boundaries_elemTags;
+ gm::getEntities(gmsh_entities, DIMENSION(2));
+ for (auto [dim, tag] : gmsh_entities)
+ {
+ std::vector<int> elemTypes;
+ gmm::getElementTypes(elemTypes, dim, tag);
+ assert(elemTypes.size() == 1);
+ for (auto& elemType : elemTypes)
+ {
+ element_key_factory ekf(DIMENSION(2), tag, elemType);
+ boundary_map[tag] = ekf.get_keys();
+ }
+ }
+
entities.clear();
import_gmsh_entities();
- map_boundaries();
+ map_boundaries(); /* This must happen after import_gmsh_entities(). */
}
void
@@ -153,6 +175,16 @@ model::num_cells(void) const
return ret;
}
+size_t
+model::num_faces(void) const
+{
+ size_t ret = 0;
+ for (const auto& e : entities)
+ ret += e.num_faces();
+
+ return ret;
+}
+
size_t
model::num_entities(void) const
{
@@ -216,20 +248,68 @@ model::import_gmsh_entities(void)
void
model::map_boundaries(void)
{
- gvp_t entities;
- /* Create a table mapping the tag of a boundary to the list of
- * its face keys */
- std::map<size_t, std::vector<face_key>> boundaries_elemTags;
- gm::getEntities(entities, DIMENSION(2));
- for (auto [dim, tag] : entities)
+ /* Make a vector mapping element_key to entity tag */
+ using bfk_t = std::pair<element_key, int>;
+ std::vector<bfk_t> bfk;
+ for (auto& [tag, keys] : boundary_map)
+ for (auto& k : keys)
+ bfk.push_back( std::make_pair(k, tag) );
+
+ /* Sort it for lookups */
+ auto comp = [](const bfk_t& a, const bfk_t& b) -> bool {
+ return a.first < b.first;
+ };
+ std::sort(bfk.begin(), bfk.end(), comp);
+
+ bnd_descriptors.resize( num_faces() );
+ size_t fbase = 0;
+ /* For each entity */
+ for (auto& e : entities)
{
- std::vector<int> elemTypes;
- gmm::getElementTypes(elemTypes, dim, tag);
- assert(elemTypes.size() == 1);
- for (auto& elemType : elemTypes)
+ /* and for each face */
+ for (size_t iF = 0; iF < e.num_faces(); iF++)
{
- element_key_factory ekf(DIMENSION(2), tag, elemType);
- boundary_map[tag] = ekf.get_keys();
+ /* get the element key */
+ element_key fk(e.face(iF));
+ auto lbcomp = [](const bfk_t& a, const element_key& b) -> bool {
+ return a.first < b;
+ };
+ /* and lookup it in the vector we created before. */
+ auto itor = std::lower_bound(bfk.begin(), bfk.end(), fk, lbcomp);
+ if ( (itor == bfk.end()) or ((*itor).first != fk) )
+ continue;
+
+ /* If found */
+ boundary_descriptor bd;
+ /* determine if it is an interface or a boundary */
+ if (conn.num_neighbours(fk) == 1)
+ bd.type = face_type::BOUNDARY;
+ else
+ bd.type = face_type::INTERFACE;
+ /* and store the gmsh tag in the descriptor. */
+ bd.gmsh_entity = (*itor).second;
+
+ /* Finally, put the descriptor in the global array of faces. */
+ bnd_descriptors.at(fbase + iF) = bd;
+ }
+ fbase += e.num_faces();
+ }
+
+ /*
+ size_t normal = 0;
+ size_t interface = 0;
+ size_t boundary = 0;
+ for (auto& bd : boundary_descriptors)
+ {
+ switch (bd.type)
+ {
+ case face_type::NONE: normal++; break;
+ case face_type::INTERFACE: interface++; break;
+ case face_type::BOUNDARY: boundary++; break;
}
}
+
+ std::cout << bfk.size() << " " << boundary_descriptors.size() << std::endl;
+ std::cout << normal << " " << interface << " " << boundary << std::endl;
+ */
}
\ No newline at end of file
diff --git a/src/maxwell_solver.cpp b/src/maxwell_solver.cpp
index e8c6208..a820347 100644
--- a/src/maxwell_solver.cpp
+++ b/src/maxwell_solver.cpp
@@ -1,6 +1,7 @@
#include "gmsh.h"
#include "silo_output.hpp"
#include "maxwell_interface.h"
+#include "param_loader.h"
static void
make_geometry(int order, double mesh_h)
@@ -29,34 +30,37 @@ make_geometry(int order, double mesh_h)
}
template<typename State>
-void test_it(const model& mod, State& state)
+void test_it(const model& mod, State& state, const maxwell_parameter_loader& mpl)
{
- auto E = [](const point_3d& pt) -> vec3d {
- vec3d ret;
- ret(0) = 0;
- ret(1) = 0;
- ret(2) = std::sin(M_PI*pt.x()) * std::sin(M_PI*pt.y());
- return ret;
- };
-
maxwell::init_from_model(mod, state);
- maxwell::init_E_field(mod, state, E);
+
+ if ( mpl.initial_Efield_defined() )
+ {
+ auto E = [&](const point_3d& pt) -> vec3d { return mpl.Efield(pt); };
+ maxwell::init_E_field(mod, state, E);
+ }
+
+ if ( mpl.initial_Hfield_defined() )
+ {
+ auto H = [&](const point_3d& pt) -> vec3d { return mpl.Hfield(pt); };
+ maxwell::init_H_field(mod, state, H);
+ }
- state.delta_t = 0.001;
+ state.delta_t = mpl.delta_t();
+ auto num_timesteps = mpl.timesteps();
+ auto silo_output_rate = mpl.silo_output_rate();
- for(size_t i = 0; i < 10000; i++)
+ for(size_t i = 0; i < num_timesteps; i++)
{
timestep(state);
std::stringstream ss;
ss << "maxwell_" << i << ".silo";
- if (i%10 == 0)
+ if (i%silo_output_rate == 0)
maxwell::export_to_silo(mod, state, ss.str());
std::swap(state.emf_curr, state.emf_next);
}
-
-
}
int main(int argc, const char *argv[])
@@ -64,13 +68,30 @@ int main(int argc, const char *argv[])
gmsh::initialize();
make_geometry(1, 0.1);
- model mod(1,3);
-
- maxwell::solver_state state_c;
- test_it(mod, state_c);
+ maxwell_parameter_loader mpl;
+ if ( not mpl.load_file("params.lua") )
+ {
+ std::cout << "Configuration problem, exiting" << std::endl;
+ return 1;
+ }
+
+ model mod( mpl.geom_order(), mpl.approx_order() );
+ if ( not mpl.validate_model_params(mod) )
+ {
+ std::cout << "Configuration problem, exiting" << std::endl;
+ return 1;
+ }
- //maxwell::solver_state_gpu state_g;
- //test_it(mod, state_g);
+ if ( mpl.use_gpu() )
+ {
+ maxwell::solver_state_gpu state_g;
+ test_it(mod, state_g, mpl);
+ }
+ else
+ {
+ maxwell::solver_state state_c;
+ test_it(mod, state_c, mpl);
+ }
//gmsh::finalize();
diff --git a/src/param_loader.cpp b/src/param_loader.cpp
new file mode 100644
index 0000000..aa89fea
--- /dev/null
+++ b/src/param_loader.cpp
@@ -0,0 +1,192 @@
+#include <iostream>
+#include "param_loader.h"
+
+parameter_loader::parameter_loader()
+{
+ lua.open_libraries(sol::lib::base, sol::lib::math);
+ init();
+}
+
+void
+parameter_loader::init(void)
+{
+ lua["const"] = lua.create_table();
+ lua["const"]["eps0"] = 8.8541878128e-12;
+ lua["const"]["mu0"] = 4e-7*M_PI;
+
+ lua["sim"] = lua.create_table();
+ lua["sim"]["approx_order"] = 1;
+ lua["sim"]["geom_order"] = 1;
+ lua["sim"]["use_gpu"] = 0;
+ //lua["sim"]["name"];
+ //lua["sim"]["dt"];
+ //lua["sim"]["timesteps"];
+ //lua["sim"]["gmsh_model"];
+
+ lua["materials"] = lua.create_table();
+ lua["bndconds"] = lua.create_table();
+ lua["ifaceconds"] = lua.create_table();
+
+ lua["postpro"] = lua.create_table();
+ lua["postpro"]["silo_output_rate"] = 0;
+}
+
+bool
+parameter_loader::validate_simulation_params(void) const
+{
+ bool success = true;
+
+ auto sim_name = lua["sim"]["name"];
+ if (not sim_name.valid())
+ {
+ std::cout << "[CONFIG] 'sim.name' not set" << std::endl;
+ success = false;
+ }
+
+ auto sim_dt = lua["sim"]["dt"];
+ if (not sim_dt.valid())
+ {
+ std::cout << "[CONFIG] 'sim.dt' not set" << std::endl;
+ success = false;
+ }
+
+ auto sim_timesteps = lua["sim"]["timesteps"];
+ if (not sim_timesteps.valid())
+ {
+ std::cout << "[CONFIG] 'sim.timesteps' not set" << std::endl;
+ success = false;
+ }
+
+ auto sim_gmshmodel = lua["sim"]["gmsh_model"];
+ if (not sim_gmshmodel.valid())
+ {
+ std::cout << "[CONFIG] 'sim.gmsh_model' not set" << std::endl;
+ success = false;
+ }
+
+ return success;
+}
+
+bool
+parameter_loader::load_file(const std::string& fn)
+{
+ auto script = lua.script_file(fn);
+ if (not script.valid())
+ return false;
+
+ return validate_simulation_params();
+}
+
+int
+parameter_loader::approx_order(void) const
+{
+ return lua["sim"]["approx_order"];
+}
+
+int
+parameter_loader::geom_order(void) const
+{
+ return lua["sim"]["geom_order"];
+}
+
+double
+parameter_loader::delta_t(void) const
+{
+ return lua["sim"]["dt"];
+}
+
+size_t
+parameter_loader::timesteps(void) const
+{
+ return lua["sim"]["timesteps"];
+}
+
+std::string
+parameter_loader::simulation_name(void) const
+{
+ return lua["sim"]["name"];
+}
+
+bool
+parameter_loader::use_gpu(void) const
+{
+ return lua["sim"]["use_gpu"];
+}
+
+size_t
+parameter_loader::silo_output_rate(void) const
+{
+ return lua["postpro"]["silo_output_rate"];
+}
+
+maxwell_parameter_loader::maxwell_parameter_loader()
+{}
+
+bool
+maxwell_parameter_loader::validate_materials(const std::string& mpn, int tag)
+{
+ auto mfun = lua["materials"][mpn];
+ if ( mfun.valid() )
+ return true;
+
+ auto mparams = lua["materials"][tag];
+ if ( mparams.valid() )
+ {
+ auto matparams_mpn = lua["materials"][tag][mpn];
+ if ( matparams_mpn.valid() )
+ return true;
+ }
+
+ std::cout << "[CONFIG] 'materials[" << tag << "]." << mpn << "' not defined and ";
+ std::cout << "'materials." << mpn << "(tag,x,y,z)' not present." << std::endl;
+ return false;
+}
+
+bool
+maxwell_parameter_loader::validate_model_params(const model& mod)
+{
+ bool success = true;
+ for (auto& e : mod)
+ {
+ auto tag = e.gmsh_tag();
+ bool eps_valid = validate_materials("epsilon", tag);
+ bool mu_valid = validate_materials("mu", tag);
+ bool sigma_valid = validate_materials("sigma", tag);
+ if ( (not eps_valid) or (not mu_valid) or (not sigma_valid) )
+ success = false;
+ }
+
+ return success;
+}
+
+bool
+maxwell_parameter_loader::initial_Efield_defined(void) const
+{
+ auto eic = lua["electric_initial_condition"];
+ return eic.valid();
+}
+
+vec3d
+maxwell_parameter_loader::Efield(const point_3d& pt) const
+{
+ vec3d ret;
+ sol::tie(ret(0), ret(1), ret(2)) =
+ lua["electric_initial_condition"](pt.x(), pt.y(), pt.z());
+ return ret;
+}
+
+bool
+maxwell_parameter_loader::initial_Hfield_defined(void) const
+{
+ auto mic = lua["magnetic_initial_condition"];
+ return mic.valid();
+}
+
+vec3d
+maxwell_parameter_loader::Hfield(const point_3d& pt) const
+{
+ vec3d ret;
+ sol::tie(ret(0), ret(1), ret(2)) =
+ lua["magnetic_initial_condition"](pt.x(), pt.y(), pt.z());
+ return ret;
+}
\ No newline at end of file
diff --git a/src/physical_element.cpp b/src/physical_element.cpp
index 5859b58..b871a2e 100644
--- a/src/physical_element.cpp
+++ b/src/physical_element.cpp
@@ -312,6 +312,9 @@ element_key_factory::element_key_factory(int dim, int tag, int etype)
ekeys.push_back( ek );
}
std::sort(ekeys.begin(), ekeys.end());
+
+ /* Not sure why there are duplicates */
+ ekeys.erase( std::unique(ekeys.begin(), ekeys.end()), ekeys.end() );
}
std::vector<element_key>
--
GitLab