This tutorial illustrate how the tools provided by HArDCore can be used to implement a scheme. The examples shown here come from the HHO_fullgradientdiff implementation of the Hybrid High-Order scheme for a diffusion equation.
Overview
The implementation of a scheme starts from a generic class (in src/) that describes the discretisation method, but is oblivious to the particular model to discretise. For example, HHOSpace encodes the degrees of freedom and operators required to implement an HHO scheme: unknowns on the faces and cells (described as coefficients on specific polynomial bases built in the class), local potential reconstruction, gradient and stabilisation, etc.
The scheme itself is implemented in separate files (in Schemes/). These files describe in particular a class, e.g. HArDCore3D::FullGradientDiffusion, that is specific to the model under consideration and the discretisation method. That class has variables to store in particular the system to assemble, and methods to assemble this system. The .hpp file of the scheme also contains particular test cases (choices of unknowns, to fix the source terms and boundary conditions).
A quick look at the .cpp file shows that it starts by defining a series of options that can be passed to the executable, to change some parameters of the simulation without having to re-compile the scheme. For example:
boost::program_options::options_description desc("Allowed options");
desc.add_options()
("help,h", "Display this help message")
("mesh,m", boost::program_options::value<std::string>(), "Set the mesh")
("degree,k", boost::program_options::value<size_t>()->default_value(1), "The polynomial degree of the sequence")
("pthread,p", boost::program_options::value<bool>()->default_value(true), "Use thread-based parallelism")
("solution,s", boost::program_options::value<int>()->default_value(0), "Select the solution")
("export-matrix,e", "Export matrix to Matrix Market format")
("plot", boost::program_options::value<std::string>(), "Save plot of the solution to the given filename")
("solver", boost::program_options::value<std::string>()->default_value("PardisoLU"), "Choice of solver, not case dependent. Options are: PardisoLU, UMFPACK, PaStiXLU, PaStiXLLT, EigenLU, EigenBiCGSTAB (reverts to EigenLU if the selected solver is not available)")
("stabilization-parameter,x", boost::program_options::value<double>(), "Set the stabilization parameter");
boost::program_options::variables_map vm;
boost::program_options::store(boost::program_options::parse_command_line(argc, argv, desc), vm);
boost::program_options::notify(vm);
These options are used to fix some parameters, after which the mesh is build, the boundary conditions are selected (when relevant in the scheme) – with numbering of vertices/edges/faces can be modified, e.g. to facilitate their treatment in the implementation – and the classes describing the method and the scheme are instantiated. Some parameters of these classes are possibly further adjusted depending on the options passed to the executable:
MeshBuilder meshbuilder = MeshBuilder(mesh_file);
std::unique_ptr<Mesh> mesh_ptr = meshbuilder.build_the_mesh();
BC.reorder_faces("start");
bool use_threads = (vm.count(
"pthread") ? vm[
"pthread"].as<
bool>() : true);
std::cout <<
"[main] " << (
use_threads ?
"Parallel execution" :
"Sequential execution") << std:: endl;
(...)
FullGradientDiffusion diff(hho_space, BC,
use_threads);
if(vm.count("stabilization-parameter")) {
diff.stabilizationParameter() = vm["stabilization-parameter"].as<double>();
}
The BoundaryConditions class provides definition of boundary conditions.
Definition BoundaryConditions.hpp:129
bool use_threads
Definition HHO_DiffAdvecReac.hpp:47
size_t K
Definition HHO_DiffAdvecReac.hpp:46
From that point on, the main method to assemble the system are called, the system is solved, and the complete solution is re-assembled (taking into account the boundary conditions, and possibly the static condensation applied when the system was assembled). The solver can be chosen using the HArDCore3D::LinearSolver class, which contains various solvers (depending on the libraries found/made by cmake, see INSTALL_NOTES.txt).
diff.assembleLinearSystem(f, kappa, u, UDir);
(...)
std::string name_solver = vm["solver"].as<std::string>();
LinearSolver<Stokes::SystemMatrixType> solver(name_solver);
std::cout << "[main] Solving the system using " << solver.name() << std::endl;
Eigen::VectorXd uh_solsystem = solver.compute_and_solve(diff.systemMatrix(), diff.systemVector());
Eigen::VectorXd uh = Eigen::VectorXd::Zero(diff.hhospace().dimension());
uh.head(diff.numDirDOFs()) = UDir;
uh.segment(diff.numDirDOFs(), diff.sizeSystem()) = uh_solsystem;
uh.tail(diff.numSCDOFs()) = diff.scVector() + diff.scMatrix() * uh.head(diff.numSkeletalDOFs());
Finally, the output (numerical values or figures) is created.
Space of degrees of freedom
The core class describes the degrees of freedom of the method. There are actually two such classes: HArDCore3D::GlobalDOFSpace for methods that have the same number of degrees of freedom on each mesh entity of the same dimension (all vertices have x numbers of DOFs, all edges have y numbers of DOFs, etc.), and HArDCore3D::VariableDOFSpace for methods with varying number of degree of freedom from mesh entity to mesh entity (vertex 1 has x_1 DOFs, vertex 2 has x_2 DOFs, etc.).
Each of these classes assume that the degrees of freedom are organised from the lowest-dimensional mesh entity to the highest-dimensional ones, and following the (local or global) numbering of these mesh entites. So, for example, if an element \(T\) has vertices \(V_1\), \(V_2\)..., \(V_n\), edges \(E_1\), \(E_2\)..., \(E_r\), and faces \(F_1\), \(F_2\)..., \(F_l\), the degrees of freedom are organised in this order:
- DOFs attached to \(V_1\),
- DOFs attached to \(V_2\),
- ...
- DOFs attached to \(V_n\),
- DOFs attached to \(E_1\),
- DOFs attached to \(E_2\),
- ...
- DOFs attached to \(E_r\),
- DOFs attached to \(F_1\),
- DOFs attached to \(F_2\),
- ...
- DOFs attached to \(F_l\),
- DOFs attached to \(T\).
These classes provide in particular methods to:
- easily access the location (either globally in the mesh, or locally in an element, a face, etc.) of the DOFs attached to particular mesh entities,
- restrict a global vector of DOFs to a particular mesh entity,
- extend a matrix acting, say, on the DOFs of an edge into a matrix acting on the DOFs of a mesh entity and its boundary (filling with 0 for DOFs not involved in the design of this operator),
- list the global numbering of DOFs attached to a given mesh face or element.
Discretisation method
A discretisation method is defined by a choice of DOFs and local operators acting on these DOFs. A class describing a discretisation method is therefore inherited from one of the two classes describing spaces of DOFs, to embed the information relevant to the choice of these degrees of freedom. The class then builds and stores (or provide functions to build on the fly) the relevant local operators.
HArDCore3D::HHOSpace is built on HArDCore3D::GlobalDOFSpace. Its instantiantion therefore fixes how many DOFs each vertex, edge, face and element should have. For a chosen polynomial degree K, the HHO method is based on unknowns that are polynomials of degree K on each face and polynomials of degree K in each element. Passing the argument K to the class as in
therefore results in a HArDCore3D::GlobalDOFSpace designed with zero DOFs on each vertex and edge, \(dim\mathbb{P}^K(\mathbb{R}^2)\) DOFs on each face and \(dim\mathbb{P}^K(\mathbb{R}^3)\) DOFs on each element:
HHOSpace::HHOSpace(
const Mesh & mesh,
size_t K,
bool use_threads, std::ostream & output)
: GlobalDOFSpace(mesh,
0,
0,
PolynomialSpaceDimension<Face>::Poly(
K),
PolynomialSpaceDimension<Cell>::Poly(
K)
),
(...)
Additionally, to give meaning to these DOFs on each mesh entity, the class HArDCore3D::HHOSpace builds and stores bases of the polynomial spaces of degree K on each face and each element (the construction is done here using the multi-threading capacities provided by parallel_for):
std::function<void(size_t, size_t)> construct_all_cell_bases
= [this](size_t start, size_t end)->void
{
for (size_t iT = start; iT < end; iT++) {
this->m_cell_bases[iT].reset( new CellBases(this->_construct_cell_bases(iT)) );
}
};
m_output << "[HHOSpace] Constructing element bases" << std::endl;
parallel_for(mesh.n_cells(), construct_all_cell_bases,
use_threads);
Once a vector of DOFs is known, its corresponding polynomial values on a face/element can therefore be fully described through the coefficients in the vector (the restriction of these to a face/element being accessible via the methods of HArDCore3D::GlobalDOFSpace) and the polynomial basis functions embedded in the instantiation of HArDCore3D::HHOSpace.
Finally, when this class is instantiated, the local operators (potential, gradient, stabilisation) are constructed and stored:
std::function<void(size_t, size_t)> construct_all_operators
= [this](size_t start, size_t end)->void
{
for (size_t iT = start; iT < end; iT++) {
m_operators[iT].reset( new LocalOperators(_compute_operators(iT)) );
}
};
m_output << "[HHOSpace] Constructing operators" << std::endl;
parallel_for(mesh.n_cells(), construct_all_operators,
use_threads);
The construction of these operators is method-dependent, but most often require to set up systems of equations based on matrices obtained by integrating products of polynomial basis functions (these matrices are Gram-matrix like). Such integrals (on edges, faces, elements) can be done using two approaches:
- Standard quadrature rules: weights and nodes are constructed using generate_quadrature_rule(), the basis functions are evaluated at these nodes, and the Gram matrix is built via compute_gram_matrix():
QuadratureRule quad_2k_F = generate_quadrature_rule(F, 2 * degree() );
auto basis_Pk_F_quad = evaluate_quad<Function>::compute(*faceBases(F).Polyk, quad_2k_F);
G = compute_gram_matrix(basis_Pk_F_quad, quad_2k_F);
- Polynomial integration rules: computationally much faster way to compute integrals of polynomial bases, and associated Gram matrix, for certain classes of polynomials (see the mainpage for the list). This uses the GramMatrix() functions.
MonomialCellIntegralsType int_mono_2kp2 = IntegrateCellMonomials(T, 2*degree()+2);
Eigen::MatrixXd MGT = GramMatrix(T, *cellBases(iT).Polykd, int_mono_2kp2);
(calling IntegrateCellMonomials is optional, but avoids re-computing the integrals of the monomials on each cell if they are to be re-used in another GramMatrix() in the same scope).
Local and global assembly
Most numerical schemes are assembled by computing local element-wise contributions, and distributing them in the global system. The local contributions are based on the operators built in the class describing the discretisation, and are therefore very easy to compute from this class. For HArDCore3D::FullGradientDiffusion, this is done by HArDCore3D::FullGradientDiffusion::_compute_local_contribution(), e.g.:
MonomialCellIntegralsType int_mono_2k = IntegrateCellMonomials(T, 2*m_hhospace.degree());
Eigen::MatrixXd mass_Pkd_T = GramMatrix(T, *m_hhospace.cellBases(iT).Polykd, int_mono_2k);
double kappaT = kappa.value(T, T.center_mass());
AT += kappaT * m_hhospace.operators(iT).gradient.transpose() * mass_Pkd_T * m_hhospace.operators(iT).gradient
+ kappaT * m_stab_par * m_hhospace.operators(iT).stabilisation;
QuadratureRule quad_2k_T = generate_quadrature_rule(T, 2 * m_hhospace.degree());
lT.tail(m_hhospace.numLocalDofsCell()) =
integrate(f, evaluate_quad<Function>::compute(*m_hhospace.cellBases(iT).Polyk, quad_2k_T), quad_2k_T);
In theory, the distribution of the local contribution is then straightforward using the local-to-global map available in the class describing the method (which inherits this from the class describing the space of DOFs), see e.g. HArDCore3D::HHOSpace::globalDOFIndices(). This global assembly would look something like this:
for (size_t iT=0; iT<mesh.n_cells(); iT++){
const Cell & T = *mesh.cell(iT);
Eigen::MatrixXd AT;
Eigen::VectorXd lT;
std::tie(AT, lT) = _compute_local_contribution(iT);
std::vector<size_t> IT = m_hhospace.globalDOFIndices(T);
for (size_t i = 0; i < IT.size(); i++){
rhs(IT[i]) += lT(i);
for (size_t j = 0; j < IT.size(); j++){
triplets.emplace_back(IT[i], IT[j], AT(i,j));
}
}
}
Eigen::SparseMatrix<double> GlobalMatrix(N, N);
GlobalMatrix.setFromTriplets(triplets.begin(), triplets.end());
However, to make the assembly more efficient, using the multi-threaded procedure parallel_assembly_system() is recommended. This leads to a slightly more complicated local to global assembly step, based on the definition of a functor that will invoke the corresponding _compute_local_contributions() and is called in multi-threads by parallel_assembly_system(). The function assembleLinearSystem() embeds the definition of the functor and the call to parallel_assembly_system(). The situation can also be complicated when using static condensation, to reduce the size of the globally coupled system; in that case, two systems are assembled at the same time par parallel_assembly_system(): the system on the globally coupled unknowns, and the system used to recover the statically condensed unknowns.
In practice however, applying the multi-threading and the static condensation is relatively easy. For example, the skeleton of HArDCore3D::FullGradientDiffusion::assembleLinearSystem() barely has to be modifed when changing the method or the model. The process goes like this, from the innermost to the outermost procedure:
- _compute_local_contribution() computes the local contributions, element by element,
- These contributions are passed to _assemble_local_contribution() which:
- calls _compute_static_condensation() which prepares the input data to instantiate a HArDCore3D::LocalStaticCondensation class (this input data identifies in particular the DOFs that are statically condensed via a permutation matrix that puts these DOFs at the end, and also contains the global DOF indices for statically condensed and non-statically condensed unknowns).
- this HArDCore3D::LocalStaticCondensation then performs the local static condensation from the local contributions,
- finally, these contributions are distributed in the system corresponding to globally coupled unknowns, and the system used to recover the statically condensed unknowns.
Outputs
Classes corresponding to discretisation methods or schemes often provide specific outputs of interest (e.g. a certain norm, etc.). Additionally, some classes provide methods to compute values of the solution at the vertices of the mesh. In this case, a vtk plot can be generated using HArDCore3D::VtuWriter.
std::cout << "[main] Writing solution to file" << std::endl;
std::string filename = vm["plot"].as<std::string>();
VtuWriter plotdata(mesh_ptr.get());
std::vector<double> exact_u_vertex(mesh_ptr->n_vertices(), 0.);
for (size_t iV=0; iV< mesh_ptr->n_vertices(); iV++){
exact_u_vertex[iV] = u(mesh_ptr->vertex(iV)->coords());
}
std::vector<double> u_vertex = hho_space.computeVertexValues(uh);
plotdata.write_to_vtu(filename + ".vtu",
std::vector<std::vector<double>> {u_vertex, exact_u_vertex},
std::vector<std::string>{"approximate", "exact"});
}
constexpr size_t max_nb_cells_for_plot
Definition sddr-navier-stokes.cpp:31
