From bf18a1885076646111b52d9279866f8f9383d994 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?St=C3=A9phane=20Del=20Pino?= <stephane.delpino44@gmail.com>
Date: Mon, 8 Mar 2021 07:58:44 +0100
Subject: [PATCH] Add a single time step version of the Glace scheme
This allows to write the time loop directly in the pgs file.
---
.../algorithms/AcousticSolverAlgorithm.cpp | 10 +-
src/language/modules/SchemeModule.cpp | 126 ++++
src/mesh/MeshData.hpp | 39 +
src/scheme/AcousticSolver.cpp | 713 ++++++++++++++++++
src/scheme/AcousticSolver.hpp | 61 +-
src/scheme/CMakeLists.txt | 5 +-
src/scheme/DiscreteFunctionP0.hpp | 13 +-
src/scheme/DiscreteFunctionUtils.cpp | 115 +++
src/scheme/DiscreteFunctionUtils.hpp | 43 ++
src/scheme/PerfectGas.cpp | 248 ++++++
src/scheme/PerfectGas.hpp | 32 +
11 files changed, 1390 insertions(+), 15 deletions(-)
create mode 100644 src/scheme/AcousticSolver.cpp
create mode 100644 src/scheme/DiscreteFunctionUtils.cpp
create mode 100644 src/scheme/DiscreteFunctionUtils.hpp
create mode 100644 src/scheme/PerfectGas.cpp
create mode 100644 src/scheme/PerfectGas.hpp
diff --git a/src/language/algorithms/AcousticSolverAlgorithm.cpp b/src/language/algorithms/AcousticSolverAlgorithm.cpp
index edce9e025..0af71d5a5 100644
--- a/src/language/algorithms/AcousticSolverAlgorithm.cpp
+++ b/src/language/algorithms/AcousticSolverAlgorithm.cpp
@@ -24,7 +24,7 @@ AcousticSolverAlgorithm<Dimension>::AcousticSolverAlgorithm(
std::shared_ptr<const MeshType> mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
- typename AcousticSolver<MeshType>::BoundaryConditionList bc_list;
+ typename AcousticSolverOld<MeshType>::BoundaryConditionList bc_list;
{
constexpr ItemType FaceType = [] {
if constexpr (Dimension > 1) {
@@ -39,7 +39,7 @@ AcousticSolverAlgorithm<Dimension>::AcousticSolverAlgorithm(
switch (bc_descriptor->type()) {
case IBoundaryConditionDescriptor::Type::symmetry: {
- using SymmetryBoundaryCondition = typename AcousticSolver<MeshType>::SymmetryBoundaryCondition;
+ using SymmetryBoundaryCondition = typename AcousticSolverOld<MeshType>::SymmetryBoundaryCondition;
const SymmetryBoundaryConditionDescriptor& sym_bc_descriptor =
dynamic_cast<const SymmetryBoundaryConditionDescriptor&>(*bc_descriptor);
@@ -60,7 +60,7 @@ AcousticSolverAlgorithm<Dimension>::AcousticSolverAlgorithm(
const DirichletBoundaryConditionDescriptor& dirichlet_bc_descriptor =
dynamic_cast<const DirichletBoundaryConditionDescriptor&>(*bc_descriptor);
if (dirichlet_bc_descriptor.name() == "velocity") {
- using VelocityBoundaryCondition = typename AcousticSolver<MeshType>::VelocityBoundaryCondition;
+ using VelocityBoundaryCondition = typename AcousticSolverOld<MeshType>::VelocityBoundaryCondition;
for (size_t i_ref_face_list = 0;
i_ref_face_list < mesh->connectivity().template numberOfRefItemList<FaceType>(); ++i_ref_face_list) {
@@ -80,7 +80,7 @@ AcousticSolverAlgorithm<Dimension>::AcousticSolverAlgorithm(
}
}
} else if (dirichlet_bc_descriptor.name() == "pressure") {
- using PressureBoundaryCondition = typename AcousticSolver<MeshType>::PressureBoundaryCondition;
+ using PressureBoundaryCondition = typename AcousticSolverOld<MeshType>::PressureBoundaryCondition;
for (size_t i_ref_face_list = 0;
i_ref_face_list < mesh->connectivity().template numberOfRefItemList<FaceType>(); ++i_ref_face_list) {
@@ -141,7 +141,7 @@ AcousticSolverAlgorithm<Dimension>::AcousticSolverAlgorithm(
}
unknowns.gammaj().fill(1.4);
- AcousticSolver acoustic_solver(mesh, bc_list, solver_type);
+ AcousticSolverOld acoustic_solver(mesh, bc_list, solver_type);
const double tmax = 0.2;
double t = 0;
diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index c23d3794e..212464f79 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -4,9 +4,11 @@
#include <language/utils/BuiltinFunctionEmbedder.hpp>
#include <language/utils/TypeDescriptor.hpp>
#include <mesh/Mesh.hpp>
+#include <scheme/AcousticSolver.hpp>
#include <scheme/DirichletBoundaryConditionDescriptor.hpp>
#include <scheme/DiscreteFunctionDescriptorP0.hpp>
#include <scheme/DiscreteFunctionInterpoler.hpp>
+#include <scheme/DiscreteFunctionUtils.hpp>
#include <scheme/FourierBoundaryConditionDescriptor.hpp>
#include <scheme/FreeBoundaryConditionDescriptor.hpp>
#include <scheme/IBoundaryConditionDescriptor.hpp>
@@ -18,6 +20,8 @@
#include <scheme/NumberedBoundaryDescriptor.hpp>
#include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
+#include <scheme/PerfectGas.hpp>
+
#include <memory>
SchemeModule::SchemeModule()
@@ -201,5 +205,127 @@ SchemeModule::SchemeModule()
}
}
+ ));
+
+ this
+ ->_addBuiltinFunction("perfect_gas_epsilon_from_rho_p_gamma",
+ std::make_shared<BuiltinFunctionEmbedder<
+ std::shared_ptr<const IDiscreteFunction>(const std::shared_ptr<const IDiscreteFunction>&,
+ const std::shared_ptr<const IDiscreteFunction>&,
+ const std::shared_ptr<const IDiscreteFunction>&)>>(
+
+ [](const std::shared_ptr<const IDiscreteFunction>& rho,
+ const std::shared_ptr<const IDiscreteFunction>& p,
+ const std::shared_ptr<const IDiscreteFunction>& gamma)
+ -> std::shared_ptr<const IDiscreteFunction> {
+ return perfect_gas::epsilonFromRhoPGamma(rho, p, gamma);
+ }
+
+ ));
+
+ this
+ ->_addBuiltinFunction("perfect_gas_p_from_rho_epsilon_gamma",
+ std::make_shared<BuiltinFunctionEmbedder<
+ std::shared_ptr<const IDiscreteFunction>(const std::shared_ptr<const IDiscreteFunction>&,
+ const std::shared_ptr<const IDiscreteFunction>&,
+ const std::shared_ptr<const IDiscreteFunction>&)>>(
+
+ [](const std::shared_ptr<const IDiscreteFunction>& rho,
+ const std::shared_ptr<const IDiscreteFunction>& epsilon,
+ const std::shared_ptr<const IDiscreteFunction>& gamma)
+ -> std::shared_ptr<const IDiscreteFunction> {
+ return perfect_gas::pFromRhoEpsilonGamma(rho, epsilon, gamma);
+ }
+
+ ));
+
+ this
+ ->_addBuiltinFunction("perfect_gas_sound_speed",
+ std::make_shared<BuiltinFunctionEmbedder<
+ std::shared_ptr<const IDiscreteFunction>(const std::shared_ptr<const IDiscreteFunction>&,
+ const std::shared_ptr<const IDiscreteFunction>&,
+ const std::shared_ptr<const IDiscreteFunction>&)>>(
+
+ [](const std::shared_ptr<const IDiscreteFunction>& rho,
+ const std::shared_ptr<const IDiscreteFunction>& p,
+ const std::shared_ptr<const IDiscreteFunction>& gamma)
+ -> std::shared_ptr<const IDiscreteFunction> {
+ return perfect_gas::soundSpeed(rho, p, gamma);
+ }
+
+ ));
+
+ this
+ ->_addBuiltinFunction("E_from_epsilon_u",
+ std::make_shared<BuiltinFunctionEmbedder<
+ std::shared_ptr<const IDiscreteFunction>(const std::shared_ptr<const IDiscreteFunction>&,
+ const std::shared_ptr<const IDiscreteFunction>&)>>(
+
+ [](const std::shared_ptr<const IDiscreteFunction>& epsilon,
+ const std::shared_ptr<const IDiscreteFunction>& u)
+ -> std::shared_ptr<const IDiscreteFunction> {
+ return perfect_gas::totalEnergyFromEpsilonU(epsilon, u);
+ }
+
+ ));
+
+ this
+ ->_addBuiltinFunction("epsilon_from_E_u",
+ std::make_shared<BuiltinFunctionEmbedder<
+ std::shared_ptr<const IDiscreteFunction>(const std::shared_ptr<const IDiscreteFunction>&,
+ const std::shared_ptr<const IDiscreteFunction>&)>>(
+
+ [](const std::shared_ptr<const IDiscreteFunction>& E,
+ const std::shared_ptr<const IDiscreteFunction>& u)
+ -> std::shared_ptr<const IDiscreteFunction> {
+ return perfect_gas::epsilonFromTotalEnergyU(E, u);
+ }
+
+ ));
+
+ this->_addBuiltinFunction(
+ "glace_solver",
+ std::make_shared<BuiltinFunctionEmbedder<std::tuple<
+ std::shared_ptr<const IMesh>,
+ std::shared_ptr<const IDiscreteFunction>,
+ std::shared_ptr<const IDiscreteFunction>,
+ std::shared_ptr<const IDiscreteFunction>>(const std::shared_ptr<const IDiscreteFunction>&,
+ const std::shared_ptr<const IDiscreteFunction>&,
+ const std::shared_ptr<const IDiscreteFunction>&,
+ const std::shared_ptr<const IDiscreteFunction>&,
+ const std::shared_ptr<const IDiscreteFunction>&,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
+ const double&)>>(
+
+ [](const std::shared_ptr<const IDiscreteFunction>& rho, const std::shared_ptr<const IDiscreteFunction>& u,
+ const std::shared_ptr<const IDiscreteFunction>& E, const std::shared_ptr<const IDiscreteFunction>& c,
+ const std::shared_ptr<const IDiscreteFunction>& p,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list, const double& dt)
+ -> std::tuple<std::shared_ptr<const IMesh>, std::shared_ptr<const IDiscreteFunction>,
+ std::shared_ptr<const IDiscreteFunction>,
+ std::shared_ptr<const IDiscreteFunction>> {
+ return AcousticSolverHandler{rho, c, u, p, bc_descriptor_list}.solver().apply(dt, rho, u, E);
+ }
+
+ ));
+
+ this
+ ->_addBuiltinFunction("lagrangian",
+ std::make_shared<BuiltinFunctionEmbedder<
+ std::shared_ptr<const IDiscreteFunction>(const std::shared_ptr<const IMesh>&,
+ const std::shared_ptr<const IDiscreteFunction>&)>>(
+
+ [](const std::shared_ptr<const IMesh>& mesh,
+ const std::shared_ptr<const IDiscreteFunction>& v)
+ -> std::shared_ptr<const IDiscreteFunction> { return shallowCopy(mesh, v); }
+
+ ));
+
+ this->_addBuiltinFunction("acoustic_dt",
+ std::make_shared<
+ BuiltinFunctionEmbedder<double(const std::shared_ptr<const IDiscreteFunction>&)>>(
+
+ [](const std::shared_ptr<const IDiscreteFunction>& c) -> double { return acoustic_dt(c); }
+
));
}
diff --git a/src/mesh/MeshData.hpp b/src/mesh/MeshData.hpp
index c82e9ce06..99cc8aeeb 100644
--- a/src/mesh/MeshData.hpp
+++ b/src/mesh/MeshData.hpp
@@ -41,6 +41,7 @@ class MeshData : public IMeshData
CellValue<const Rd> m_cell_centroid;
CellValue<const Rd> m_cell_iso_barycenter;
CellValue<const double> m_Vj;
+ CellValue<const double> m_sum_r_ljr;
FaceValue<const double> m_ll;
PUGS_INLINE
@@ -435,6 +436,34 @@ class MeshData : public IMeshData
}
}
+ PUGS_INLINE
+ void
+ _computeSumOverRljr()
+ {
+ CellValue<double> sum_r_ljr(m_mesh.connectivity());
+
+ if constexpr (Dimension == 1) {
+ parallel_for(
+ m_mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { sum_r_ljr[j] = 2; });
+ } else {
+ auto ljr = this->ljr();
+
+ const auto& cell_to_node_matrix = m_mesh.connectivity().cellToNodeMatrix();
+
+ parallel_for(
+ m_mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
+ const auto& cell_nodes = cell_to_node_matrix[j];
+ double sum = 0;
+ for (size_t r = 0; r < cell_nodes.size(); ++r) {
+ sum += ljr(j, r);
+ }
+ sum_r_ljr[j] = sum;
+ });
+ }
+
+ m_sum_r_ljr = sum_r_ljr;
+ }
+
void
_checkCellVolume()
{
@@ -563,6 +592,16 @@ class MeshData : public IMeshData
return m_Vj;
}
+ PUGS_INLINE
+ CellValue<const double>
+ sumOverRLjr()
+ {
+ if (not m_sum_r_ljr.isBuilt()) {
+ this->_computeSumOverRljr();
+ }
+ return m_sum_r_ljr;
+ }
+
private:
// MeshData **must** be constructed through MeshDataManager
friend class MeshDataManager;
diff --git a/src/scheme/AcousticSolver.cpp b/src/scheme/AcousticSolver.cpp
new file mode 100644
index 000000000..05a8adaf5
--- /dev/null
+++ b/src/scheme/AcousticSolver.cpp
@@ -0,0 +1,713 @@
+#include <scheme/AcousticSolver.hpp>
+
+#include <mesh/MeshNodeBoundary.hpp>
+#include <scheme/DirichletBoundaryConditionDescriptor.hpp>
+#include <scheme/DiscreteFunctionP0.hpp>
+#include <scheme/DiscreteFunctionUtils.hpp>
+#include <scheme/IBoundaryConditionDescriptor.hpp>
+#include <scheme/IDiscreteFunction.hpp>
+#include <scheme/IDiscreteFunctionDescriptor.hpp>
+#include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
+
+#include <variant>
+#include <vector>
+
+template <size_t Dimension>
+double
+acoustic_dt(const DiscreteFunctionP0<Dimension, double>& c)
+{
+ const std::shared_ptr<const Mesh<Connectivity<Dimension>>> mesh =
+ std::dynamic_pointer_cast<const Mesh<Connectivity<Dimension>>>(c.mesh());
+
+ const auto Vj = MeshDataManager::instance().getMeshData(*mesh).Vj();
+ const auto Sj = MeshDataManager::instance().getMeshData(*mesh).sumOverRLjr();
+
+ CellValue<double> local_dt{mesh->connectivity()};
+ parallel_for(
+ mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { local_dt[j] = 2 * Vj[j] / (Sj[j] * c[j]); });
+
+ return min(local_dt);
+}
+
+double
+acoustic_dt(const std::shared_ptr<const IDiscreteFunction>& c)
+{
+ if ((c->descriptor().type() != DiscreteFunctionType::P0) or (c->dataType() != ASTNodeDataType::double_t)) {
+ throw NormalError("invalid discrete function type");
+ }
+
+ std::shared_ptr mesh = c->mesh();
+
+ switch (mesh->dimension()) {
+ case 1: {
+ return acoustic_dt(dynamic_cast<const DiscreteFunctionP0<1, double>&>(*c));
+ }
+ case 2: {
+ return acoustic_dt(dynamic_cast<const DiscreteFunctionP0<2, double>&>(*c));
+ }
+ case 3: {
+ return acoustic_dt(dynamic_cast<const DiscreteFunctionP0<3, double>&>(*c));
+ }
+ default: {
+ throw UnexpectedError("invalid mesh dimension");
+ }
+ }
+}
+
+template <size_t Dimension>
+class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler::IAcousticSolver
+{
+ private:
+ using Rdxd = TinyMatrix<Dimension>;
+ using Rd = TinyVector<Dimension>;
+
+ using MeshType = Mesh<Connectivity<Dimension>>;
+ using MeshDataType = MeshData<Dimension>;
+
+ using DiscreteScalarFunction = DiscreteFunctionP0<Dimension, double>;
+ using DiscreteVectorFunction = DiscreteFunctionP0<Dimension, Rd>;
+
+ class PressureBoundaryCondition;
+ class SymmetryBoundaryCondition;
+ class VelocityBoundaryCondition;
+
+ using BoundaryCondition =
+ std::variant<PressureBoundaryCondition, SymmetryBoundaryCondition, VelocityBoundaryCondition>;
+
+ using BoundaryConditionList = std::vector<BoundaryCondition>;
+
+ BoundaryConditionList m_boundary_condition_list;
+
+ NodeValue<const Rd> m_ur;
+ NodeValuePerCell<const Rd> m_Fjr;
+
+ CellValue<const double>
+ _getRhoC(const DiscreteScalarFunction& rho, const DiscreteScalarFunction& c)
+ {
+ Assert(rho.mesh() == c.mesh());
+
+ std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(rho.mesh());
+
+ CellValue<double> rhoc{mesh->connectivity()};
+ parallel_for(
+ mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { rhoc[cell_id] = rho[cell_id] * c[cell_id]; });
+
+ return rhoc;
+ }
+
+ NodeValuePerCell<const Rdxd>
+ _computeGlaceAjr(const MeshType& mesh, const CellValue<const double>& rhoc)
+ {
+ MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+
+ const NodeValuePerCell<const TinyVector<Dimension>> Cjr = mesh_data.Cjr();
+ const NodeValuePerCell<const TinyVector<Dimension>> njr = mesh_data.njr();
+
+ NodeValuePerCell<Rdxd> Ajr{mesh.connectivity()};
+
+ parallel_for(
+ mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
+ const size_t& nb_nodes = Ajr.numberOfSubValues(j);
+ const double& rhoc_j = rhoc[j];
+ for (size_t r = 0; r < nb_nodes; ++r) {
+ Ajr(j, r) = tensorProduct(rhoc_j * Cjr(j, r), njr(j, r));
+ }
+ });
+
+ return Ajr;
+ }
+
+ NodeValue<Rdxd>
+ _computeAr(const MeshType& mesh, const NodeValuePerCell<const Rdxd>& Ajr)
+ {
+ const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix();
+ const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells();
+
+ NodeValue<Rdxd> Ar{mesh.connectivity()};
+
+ parallel_for(
+ mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) {
+ Rdxd sum = zero;
+ const auto& node_to_cell = node_to_cell_matrix[r];
+ const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemValues(r);
+
+ for (size_t j = 0; j < node_to_cell.size(); ++j) {
+ const CellId J = node_to_cell[j];
+ const unsigned int R = node_local_number_in_its_cells[j];
+ sum += Ajr(J, R);
+ }
+ Ar[r] = sum;
+ });
+
+ return Ar;
+ }
+
+ NodeValue<Rd>
+ _computeBr(const Mesh<Connectivity<Dimension>>& mesh,
+ const NodeValuePerCell<const Rdxd>& Ajr,
+ const DiscreteVectorFunction& u,
+ const DiscreteScalarFunction& p)
+ {
+ MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+
+ const NodeValuePerCell<const Rd>& Cjr = mesh_data.Cjr();
+
+ const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix();
+ const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells();
+
+ NodeValue<Rd> b{mesh.connectivity()};
+
+ parallel_for(
+ mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) {
+ const auto& node_to_cell = node_to_cell_matrix[r];
+ const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemValues(r);
+
+ Rd br = zero;
+ for (size_t j = 0; j < node_to_cell.size(); ++j) {
+ const CellId J = node_to_cell[j];
+ const unsigned int R = node_local_number_in_its_cells[j];
+ br += Ajr(J, R) * u[J] + p[J] * Cjr(J, R);
+ }
+
+ b[r] = br;
+ });
+
+ return b;
+ }
+
+ BoundaryConditionList
+ _getBCList(const std::shared_ptr<const MeshType>& mesh,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list)
+ {
+ BoundaryConditionList bc_list;
+
+ constexpr ItemType FaceType = [] {
+ if constexpr (Dimension > 1) {
+ return ItemType::face;
+ } else {
+ return ItemType::node;
+ }
+ }();
+
+ for (const auto& bc_descriptor : bc_descriptor_list) {
+ bool is_valid_boundary_condition = true;
+
+ switch (bc_descriptor->type()) {
+ case IBoundaryConditionDescriptor::Type::symmetry: {
+ const SymmetryBoundaryConditionDescriptor& sym_bc_descriptor =
+ dynamic_cast<const SymmetryBoundaryConditionDescriptor&>(*bc_descriptor);
+ for (size_t i_ref_face_list = 0;
+ i_ref_face_list < mesh->connectivity().template numberOfRefItemList<FaceType>(); ++i_ref_face_list) {
+ const auto& ref_face_list = mesh->connectivity().template refItemList<FaceType>(i_ref_face_list);
+ const RefId& ref = ref_face_list.refId();
+ if (ref == sym_bc_descriptor.boundaryDescriptor()) {
+ SymmetryBoundaryCondition{MeshFlatNodeBoundary<Dimension>(mesh, ref_face_list)};
+ bc_list.push_back(SymmetryBoundaryCondition{MeshFlatNodeBoundary<Dimension>(mesh, ref_face_list)});
+ }
+ }
+ is_valid_boundary_condition = true;
+ break;
+ }
+
+ case IBoundaryConditionDescriptor::Type::dirichlet: {
+ const DirichletBoundaryConditionDescriptor& dirichlet_bc_descriptor =
+ dynamic_cast<const DirichletBoundaryConditionDescriptor&>(*bc_descriptor);
+ if (dirichlet_bc_descriptor.name() == "velocity") {
+ for (size_t i_ref_face_list = 0;
+ i_ref_face_list < mesh->connectivity().template numberOfRefItemList<FaceType>(); ++i_ref_face_list) {
+ const auto& ref_face_list = mesh->connectivity().template refItemList<FaceType>(i_ref_face_list);
+ const RefId& ref = ref_face_list.refId();
+ if (ref == dirichlet_bc_descriptor.boundaryDescriptor()) {
+ MeshNodeBoundary<Dimension> mesh_node_boundary{mesh, ref_face_list};
+
+ const FunctionSymbolId velocity_id = dirichlet_bc_descriptor.rhsSymbolId();
+
+ const auto& node_list = mesh_node_boundary.nodeList();
+
+ Array<const TinyVector<Dimension>> value_list = InterpolateItemValue<TinyVector<Dimension>(
+ TinyVector<Dimension>)>::template interpolate<ItemType::node>(velocity_id, mesh->xr(), node_list);
+
+ bc_list.push_back(VelocityBoundaryCondition{node_list, value_list});
+ }
+ }
+ } else if (dirichlet_bc_descriptor.name() == "pressure") {
+ for (size_t i_ref_face_list = 0;
+ i_ref_face_list < mesh->connectivity().template numberOfRefItemList<FaceType>(); ++i_ref_face_list) {
+ const auto& ref_face_list = mesh->connectivity().template refItemList<FaceType>(i_ref_face_list);
+ const RefId& ref = ref_face_list.refId();
+ if (ref == dirichlet_bc_descriptor.boundaryDescriptor()) {
+ const auto& face_list = ref_face_list.list();
+
+ const FunctionSymbolId pressure_id = dirichlet_bc_descriptor.rhsSymbolId();
+
+ Array<const double> face_values = [&] {
+ if constexpr (Dimension == 1) {
+ return InterpolateItemValue<double(
+ TinyVector<Dimension>)>::template interpolate<FaceType>(pressure_id, mesh->xr(), face_list);
+ } else {
+ MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
+
+ return InterpolateItemValue<double(
+ TinyVector<Dimension>)>::template interpolate<FaceType>(pressure_id, mesh_data.xl(), face_list);
+ }
+ }();
+
+ bc_list.push_back(PressureBoundaryCondition{face_list, face_values});
+ }
+ }
+ } else {
+ is_valid_boundary_condition = false;
+ }
+ break;
+ }
+ default: {
+ is_valid_boundary_condition = false;
+ }
+ }
+ if (not is_valid_boundary_condition) {
+ std::ostringstream error_msg;
+ error_msg << *bc_descriptor << " is an invalid boundary condition for acoustic solver";
+ throw NormalError(error_msg.str());
+ }
+ }
+
+ return bc_list;
+ }
+
+ void _applyPressureBC(const MeshType& mesh, NodeValue<Rd>& br);
+ void _applySymmetryBC(NodeValue<Rdxd>& Ar, NodeValue<Rd>& br);
+ void _applyVelocityBC(NodeValue<Rdxd>& Ar, NodeValue<Rd>& br);
+
+ void
+ _applyBoundaryConditions(const MeshType& mesh, NodeValue<Rdxd>& Ar, NodeValue<Rd>& br)
+ {
+ this->_applyPressureBC(mesh, br);
+ this->_applySymmetryBC(Ar, br);
+ this->_applyVelocityBC(Ar, br);
+ }
+
+ NodeValue<const Rd>
+ _computeUr(const MeshType& mesh, const NodeValue<Rdxd>& Ar, const NodeValue<Rd>& br)
+ {
+ NodeValue<Rd> u{mesh.connectivity()};
+ parallel_for(
+ mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) { u[r] = inverse(Ar[r]) * br[r]; });
+
+ return u;
+ }
+
+ NodeValuePerCell<Rd>
+ _computeFjr(const MeshType& mesh,
+ const NodeValuePerCell<const Rdxd>& Ajr,
+ const NodeValue<const Rd>& ur,
+ const DiscreteVectorFunction& u,
+ const DiscreteScalarFunction& p)
+ {
+ MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+
+ const NodeValuePerCell<const Rd> Cjr = mesh_data.Cjr();
+
+ const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
+
+ NodeValuePerCell<Rd> F{mesh.connectivity()};
+ parallel_for(
+ mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
+ const auto& cell_nodes = cell_to_node_matrix[j];
+
+ for (size_t r = 0; r < cell_nodes.size(); ++r) {
+ F(j, r) = Ajr(j, r) * (u[j] - ur[cell_nodes[r]]) + p[j] * Cjr(j, r);
+ }
+ });
+
+ return F;
+ }
+
+ AcousticSolver(const std::shared_ptr<const MeshType>& p_mesh,
+ const DiscreteScalarFunction& rho,
+ const DiscreteScalarFunction& c,
+ const DiscreteVectorFunction& u,
+ const DiscreteScalarFunction& p,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list)
+ : m_boundary_condition_list{this->_getBCList(p_mesh, bc_descriptor_list)}
+ {
+ const MeshType& mesh = *p_mesh;
+
+ CellValue<const double> rhoc = this->_getRhoC(rho, c);
+ NodeValuePerCell<const Rdxd> Ajr = this->_computeGlaceAjr(mesh, rhoc);
+
+ NodeValue<Rdxd> Ar = this->_computeAr(mesh, Ajr);
+ NodeValue<Rd> br = this->_computeBr(mesh, Ajr, u, p);
+
+ this->_applyBoundaryConditions(mesh, Ar, br);
+
+ synchronize(Ar);
+ synchronize(br);
+
+ m_ur = this->_computeUr(mesh, Ar, br);
+ m_Fjr = this->_computeFjr(mesh, Ajr, m_ur, u, p);
+ }
+
+ public:
+ std::tuple<std::shared_ptr<const IMesh>,
+ std::shared_ptr<const DiscreteFunctionP0<Dimension, double>>,
+ std::shared_ptr<const DiscreteFunctionP0<Dimension, Rd>>,
+ std::shared_ptr<const DiscreteFunctionP0<Dimension, double>>>
+ apply(const double& dt,
+ const std::shared_ptr<const MeshType>& mesh,
+ const DiscreteFunctionP0<Dimension, double>& rho,
+ const DiscreteFunctionP0<Dimension, Rd>& u,
+ const DiscreteFunctionP0<Dimension, double>& E) const
+ {
+ const auto& cell_to_node_matrix = mesh->connectivity().cellToNodeMatrix();
+
+ NodeValue<Rd> new_xr = copy(mesh->xr());
+ parallel_for(
+ mesh->numberOfNodes(), PUGS_LAMBDA(NodeId r) { new_xr[r] += dt * m_ur[r]; });
+
+ std::shared_ptr<const MeshType> new_mesh = std::make_shared<MeshType>(mesh->shared_connectivity(), new_xr);
+
+ CellValue<const double> Vj = MeshDataManager::instance().getMeshData(*mesh).Vj();
+
+ CellValue<double> new_rho = copy(rho.cellValues());
+ CellValue<Rd> new_u = copy(u.cellValues());
+ CellValue<double> new_E = copy(E.cellValues());
+
+ parallel_for(
+ mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
+ const auto& cell_nodes = cell_to_node_matrix[j];
+
+ Rd momentum_fluxes = zero;
+ double energy_fluxes = 0;
+ for (size_t R = 0; R < cell_nodes.size(); ++R) {
+ const NodeId r = cell_nodes[R];
+ momentum_fluxes += m_Fjr(j, R);
+ energy_fluxes += (m_Fjr(j, R), m_ur[r]);
+ }
+ const double dt_over_Mj = dt / (rho[j] * Vj[j]);
+ new_u[j] -= dt_over_Mj * momentum_fluxes;
+ new_E[j] -= dt_over_Mj * energy_fluxes;
+ });
+
+ CellValue<const double> new_Vj = MeshDataManager::instance().getMeshData(*new_mesh).Vj();
+
+ parallel_for(
+ mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { new_rho[j] *= Vj[j] / new_Vj[j]; });
+
+ return {new_mesh, std::make_shared<DiscreteScalarFunction>(new_mesh, new_rho),
+ std::make_shared<DiscreteVectorFunction>(new_mesh, new_u),
+ std::make_shared<DiscreteScalarFunction>(new_mesh, new_E)};
+ }
+
+ std::tuple<std::shared_ptr<const IMesh>,
+ std::shared_ptr<const IDiscreteFunction>,
+ std::shared_ptr<const IDiscreteFunction>,
+ std::shared_ptr<const IDiscreteFunction>>
+ apply(const double& dt,
+ const std::shared_ptr<const IDiscreteFunction>& rho,
+ const std::shared_ptr<const IDiscreteFunction>& u,
+ const std::shared_ptr<const IDiscreteFunction>& E) const
+ {
+ std::shared_ptr i_mesh = getCommonMesh({rho, u, E});
+ if (not i_mesh) {
+ throw NormalError("discrete functions are not defined on the same mesh");
+ }
+
+ if (not checkDiscretizationType({rho, u, E}, DiscreteFunctionType::P0)) {
+ throw NormalError("acoustic solver expects P0 functions");
+ }
+
+ return this->apply(dt, std::dynamic_pointer_cast<const MeshType>(i_mesh),
+ *std::dynamic_pointer_cast<const DiscreteScalarFunction>(rho),
+ *std::dynamic_pointer_cast<const DiscreteVectorFunction>(u),
+ *std::dynamic_pointer_cast<const DiscreteScalarFunction>(E));
+ }
+
+ AcousticSolver(const std::shared_ptr<const IMesh>& mesh,
+ const std::shared_ptr<const IDiscreteFunction>& rho,
+ const std::shared_ptr<const IDiscreteFunction>& c,
+ const std::shared_ptr<const IDiscreteFunction>& u,
+ const std::shared_ptr<const IDiscreteFunction>& p,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list)
+ : AcousticSolver{std::dynamic_pointer_cast<const Mesh<Connectivity<Dimension>>>(mesh),
+ *std::dynamic_pointer_cast<const DiscreteScalarFunction>(rho),
+ *std::dynamic_pointer_cast<const DiscreteScalarFunction>(c),
+ *std::dynamic_pointer_cast<const DiscreteVectorFunction>(u),
+ *std::dynamic_pointer_cast<const DiscreteScalarFunction>(p),
+ bc_descriptor_list}
+ {}
+
+ AcousticSolver() = default;
+ AcousticSolver(AcousticSolver&&) = default;
+ ~AcousticSolver() = default;
+};
+
+template <size_t Dimension>
+void
+AcousticSolverHandler::AcousticSolver<Dimension>::_applyPressureBC(const MeshType& mesh, NodeValue<Rd>& br)
+{
+ for (const auto& boundary_condition : m_boundary_condition_list) {
+ std::visit(
+ [&](auto&& bc) {
+ using T = std::decay_t<decltype(bc)>;
+ if constexpr (std::is_same_v<PressureBoundaryCondition, T>) {
+ MeshData<Dimension>& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+ if constexpr (Dimension == 1) {
+ const NodeValuePerCell<const Rd> Cjr = mesh_data.Cjr();
+
+ const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix();
+ const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells();
+
+ const auto& node_list = bc.faceList();
+ const auto& value_list = bc.valueList();
+ parallel_for(
+ node_list.size(), PUGS_LAMBDA(size_t i_node) {
+ const NodeId node_id = node_list[i_node];
+ const auto& node_cell_list = node_to_cell_matrix[node_id];
+ Assert(node_cell_list.size() == 1);
+
+ CellId node_cell_id = node_cell_list[0];
+ size_t node_local_number_in_cell = node_local_numbers_in_their_cells(node_id, 0);
+
+ br[node_id] -= value_list[i_node] * Cjr(node_cell_id, node_local_number_in_cell);
+ });
+ } else {
+ const NodeValuePerFace<const Rd> Nlr = mesh_data.Nlr();
+
+ const auto& face_to_cell_matrix = mesh.connectivity().faceToCellMatrix();
+ const auto& face_to_node_matrix = mesh.connectivity().faceToNodeMatrix();
+ const auto& face_local_numbers_in_their_cells = mesh.connectivity().faceLocalNumbersInTheirCells();
+ const auto& face_cell_is_reversed = mesh.connectivity().cellFaceIsReversed();
+
+ const auto& face_list = bc.faceList();
+ const auto& value_list = bc.valueList();
+ for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
+ const FaceId face_id = face_list[i_face];
+ const auto& face_cell_list = face_to_cell_matrix[face_id];
+ Assert(face_cell_list.size() == 1);
+
+ CellId face_cell_id = face_cell_list[0];
+ size_t face_local_number_in_cell = face_local_numbers_in_their_cells(face_id, 0);
+
+ const double sign = face_cell_is_reversed(face_cell_id, face_local_number_in_cell) ? -1 : 1;
+
+ const auto& face_nodes = face_to_node_matrix[face_id];
+
+ for (size_t i_node = 0; i_node < face_nodes.size(); ++i_node) {
+ NodeId node_id = face_nodes[i_node];
+ br[node_id] -= sign * value_list[i_face] * Nlr(face_id, i_node);
+ }
+ }
+ }
+ }
+ },
+ boundary_condition);
+ }
+}
+
+template <size_t Dimension>
+void
+AcousticSolverHandler::AcousticSolver<Dimension>::_applySymmetryBC(NodeValue<Rdxd>& Ar, NodeValue<Rd>& br)
+{
+ for (const auto& boundary_condition : m_boundary_condition_list) {
+ std::visit(
+ [&](auto&& bc) {
+ using T = std::decay_t<decltype(bc)>;
+ if constexpr (std::is_same_v<SymmetryBoundaryCondition, T>) {
+ const Rd& n = bc.outgoingNormal();
+
+ const Rdxd I = identity;
+ const Rdxd nxn = tensorProduct(n, n);
+ const Rdxd P = I - nxn;
+
+ const Array<const NodeId>& node_list = bc.nodeList();
+ parallel_for(
+ bc.numberOfNodes(), PUGS_LAMBDA(int r_number) {
+ const NodeId r = node_list[r_number];
+
+ Ar[r] = P * Ar[r] * P + nxn;
+ br[r] = P * br[r];
+ });
+ }
+ },
+ boundary_condition);
+ }
+}
+
+template <size_t Dimension>
+void
+AcousticSolverHandler::AcousticSolver<Dimension>::_applyVelocityBC(NodeValue<Rdxd>& Ar, NodeValue<Rd>& br)
+{
+ for (const auto& boundary_condition : m_boundary_condition_list) {
+ std::visit(
+ [&](auto&& bc) {
+ using T = std::decay_t<decltype(bc)>;
+ if constexpr (std::is_same_v<VelocityBoundaryCondition, T>) {
+ const auto& node_list = bc.nodeList();
+ const auto& value_list = bc.valueList();
+
+ parallel_for(
+ node_list.size(), PUGS_LAMBDA(size_t i_node) {
+ NodeId node_id = node_list[i_node];
+ const auto& value = value_list[i_node];
+
+ Ar[node_id] = identity;
+ br[node_id] = value;
+ });
+ }
+ },
+ boundary_condition);
+ }
+}
+
+template <size_t Dimension>
+class AcousticSolverHandler::AcousticSolver<Dimension>::PressureBoundaryCondition
+{
+ private:
+ const Array<const double> m_value_list;
+ const Array<const FaceId> m_face_list;
+
+ public:
+ const Array<const FaceId>&
+ faceList() const
+ {
+ return m_face_list;
+ }
+
+ const Array<const double>&
+ valueList() const
+ {
+ return m_value_list;
+ }
+
+ PressureBoundaryCondition(const Array<const FaceId>& face_list, const Array<const double>& value_list)
+ : m_value_list{value_list}, m_face_list{face_list}
+ {}
+
+ ~PressureBoundaryCondition() = default;
+};
+
+template <>
+class AcousticSolverHandler::AcousticSolver<1>::PressureBoundaryCondition
+{
+ private:
+ const Array<const double> m_value_list;
+ const Array<const NodeId> m_face_list;
+
+ public:
+ const Array<const NodeId>&
+ faceList() const
+ {
+ return m_face_list;
+ }
+
+ const Array<const double>&
+ valueList() const
+ {
+ return m_value_list;
+ }
+
+ PressureBoundaryCondition(const Array<const NodeId>& face_list, const Array<const double>& value_list)
+ : m_value_list{value_list}, m_face_list{face_list}
+ {}
+
+ ~PressureBoundaryCondition() = default;
+};
+
+template <size_t Dimension>
+class AcousticSolverHandler::AcousticSolver<Dimension>::VelocityBoundaryCondition
+{
+ private:
+ const Array<const TinyVector<Dimension>> m_value_list;
+ const Array<const NodeId> m_node_list;
+
+ public:
+ const Array<const NodeId>&
+ nodeList() const
+ {
+ return m_node_list;
+ }
+
+ const Array<const TinyVector<Dimension>>&
+ valueList() const
+ {
+ return m_value_list;
+ }
+
+ VelocityBoundaryCondition(const Array<const NodeId>& node_list, const Array<const TinyVector<Dimension>>& value_list)
+ : m_value_list{value_list}, m_node_list{node_list}
+ {}
+
+ ~VelocityBoundaryCondition() = default;
+};
+
+template <size_t Dimension>
+class AcousticSolverHandler::AcousticSolver<Dimension>::SymmetryBoundaryCondition
+{
+ public:
+ using Rd = TinyVector<Dimension, double>;
+
+ private:
+ const MeshFlatNodeBoundary<Dimension> m_mesh_flat_node_boundary;
+
+ public:
+ const Rd&
+ outgoingNormal() const
+ {
+ return m_mesh_flat_node_boundary.outgoingNormal();
+ }
+
+ size_t
+ numberOfNodes() const
+ {
+ return m_mesh_flat_node_boundary.nodeList().size();
+ }
+
+ const Array<const NodeId>&
+ nodeList() const
+ {
+ return m_mesh_flat_node_boundary.nodeList();
+ }
+
+ SymmetryBoundaryCondition(const MeshFlatNodeBoundary<Dimension>& mesh_flat_node_boundary)
+ : m_mesh_flat_node_boundary(mesh_flat_node_boundary)
+ {
+ ;
+ }
+
+ ~SymmetryBoundaryCondition() = default;
+};
+
+AcousticSolverHandler::AcousticSolverHandler(
+ const std::shared_ptr<const IDiscreteFunction>& rho,
+ const std::shared_ptr<const IDiscreteFunction>& c,
+ const std::shared_ptr<const IDiscreteFunction>& u,
+ const std::shared_ptr<const IDiscreteFunction>& p,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list)
+{
+ std::shared_ptr i_mesh = getCommonMesh({rho, c, u, p});
+ if (not i_mesh) {
+ throw NormalError("discrete functions are not defined on the same mesh");
+ }
+
+ if (not checkDiscretizationType({rho, c, u, p}, DiscreteFunctionType::P0)) {
+ throw NormalError("acoustic solver expects P0 functions");
+ }
+
+ switch (i_mesh->dimension()) {
+ case 1: {
+ m_acoustic_solver = std::make_unique<AcousticSolver<1>>(i_mesh, rho, c, u, p, bc_descriptor_list);
+ break;
+ }
+ case 2: {
+ m_acoustic_solver = std::make_unique<AcousticSolver<2>>(i_mesh, rho, c, u, p, bc_descriptor_list);
+ break;
+ }
+ case 3: {
+ m_acoustic_solver = std::make_unique<AcousticSolver<3>>(i_mesh, rho, c, u, p, bc_descriptor_list);
+ break;
+ }
+ default: {
+ throw UnexpectedError("invalid mesh dimension");
+ }
+ }
+}
diff --git a/src/scheme/AcousticSolver.hpp b/src/scheme/AcousticSolver.hpp
index 3e3151d25..6618361cc 100644
--- a/src/scheme/AcousticSolver.hpp
+++ b/src/scheme/AcousticSolver.hpp
@@ -21,8 +21,53 @@
#include <iostream>
+class IDiscreteFunction;
+class IBoundaryConditionDescriptor;
+
+double acoustic_dt(const std::shared_ptr<const IDiscreteFunction>& c);
+
+class AcousticSolverHandler
+{
+ private:
+ struct IAcousticSolver
+ {
+ virtual std::tuple<std::shared_ptr<const IMesh>,
+ std::shared_ptr<const IDiscreteFunction>,
+ std::shared_ptr<const IDiscreteFunction>,
+ std::shared_ptr<const IDiscreteFunction>>
+ apply(const double& dt,
+ const std::shared_ptr<const IDiscreteFunction>& rho,
+ const std::shared_ptr<const IDiscreteFunction>& u,
+ const std::shared_ptr<const IDiscreteFunction>& E) const = 0;
+
+ IAcousticSolver() = default;
+ IAcousticSolver(IAcousticSolver&&) = default;
+ IAcousticSolver& operator=(IAcousticSolver&&) = default;
+
+ virtual ~IAcousticSolver() = default;
+ };
+
+ template <size_t Dimension>
+ class AcousticSolver;
+
+ std::unique_ptr<IAcousticSolver> m_acoustic_solver;
+
+ public:
+ const IAcousticSolver&
+ solver() const
+ {
+ return *m_acoustic_solver;
+ }
+
+ AcousticSolverHandler(const std::shared_ptr<const IDiscreteFunction>& rho,
+ const std::shared_ptr<const IDiscreteFunction>& c,
+ const std::shared_ptr<const IDiscreteFunction>& u,
+ const std::shared_ptr<const IDiscreteFunction>& p,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list);
+};
+
template <typename MeshType>
-class AcousticSolver
+class AcousticSolverOld
{
public:
class PressureBoundaryCondition;
@@ -392,9 +437,9 @@ class AcousticSolver
CellValue<double> m_Vj_over_cj;
public:
- AcousticSolver(std::shared_ptr<const MeshType> p_mesh,
- const BoundaryConditionList& bc_list,
- const AcousticSolverType solver_type)
+ AcousticSolverOld(std::shared_ptr<const MeshType> p_mesh,
+ const BoundaryConditionList& bc_list,
+ const AcousticSolverType solver_type)
: m_mesh(p_mesh),
m_connectivity(m_mesh->connectivity()),
m_boundary_condition_list(bc_list),
@@ -490,7 +535,7 @@ class AcousticSolver
};
template <typename MeshType>
-class AcousticSolver<MeshType>::PressureBoundaryCondition
+class AcousticSolverOld<MeshType>::PressureBoundaryCondition
{
private:
const Array<const double> m_value_list;
@@ -517,7 +562,7 @@ class AcousticSolver<MeshType>::PressureBoundaryCondition
};
template <>
-class AcousticSolver<Mesh<Connectivity<1>>>::PressureBoundaryCondition
+class AcousticSolverOld<Mesh<Connectivity<1>>>::PressureBoundaryCondition
{
private:
const Array<const double> m_value_list;
@@ -544,7 +589,7 @@ class AcousticSolver<Mesh<Connectivity<1>>>::PressureBoundaryCondition
};
template <typename MeshType>
-class AcousticSolver<MeshType>::VelocityBoundaryCondition
+class AcousticSolverOld<MeshType>::VelocityBoundaryCondition
{
private:
const Array<const TinyVector<Dimension>> m_value_list;
@@ -571,7 +616,7 @@ class AcousticSolver<MeshType>::VelocityBoundaryCondition
};
template <typename MeshType>
-class AcousticSolver<MeshType>::SymmetryBoundaryCondition
+class AcousticSolverOld<MeshType>::SymmetryBoundaryCondition
{
public:
static constexpr size_t Dimension = MeshType::Dimension;
diff --git a/src/scheme/CMakeLists.txt b/src/scheme/CMakeLists.txt
index 1467393f2..b7522b9bb 100644
--- a/src/scheme/CMakeLists.txt
+++ b/src/scheme/CMakeLists.txt
@@ -2,4 +2,7 @@
add_library(
PugsScheme
- DiscreteFunctionInterpoler.cpp)
+ AcousticSolver.cpp
+ DiscreteFunctionInterpoler.cpp
+ DiscreteFunctionUtils.cpp
+ PerfectGas.cpp)
diff --git a/src/scheme/DiscreteFunctionP0.hpp b/src/scheme/DiscreteFunctionP0.hpp
index bf97079b7..c51e8e09d 100644
--- a/src/scheme/DiscreteFunctionP0.hpp
+++ b/src/scheme/DiscreteFunctionP0.hpp
@@ -28,7 +28,7 @@ class DiscreteFunctionP0 : public IDiscreteFunction
return ast_node_data_type_from<DataType>;
}
- CellValue<DataType>
+ const CellValue<DataType>&
cellValues() const
{
return m_cell_values;
@@ -63,6 +63,17 @@ class DiscreteFunctionP0 : public IDiscreteFunction
mesh_data.xj());
}
+ DiscreteFunctionP0(const std::shared_ptr<const MeshType>& mesh) : m_mesh{mesh}, m_cell_values{mesh->connectivity()}
+ {
+ ;
+ }
+
+ DiscreteFunctionP0(const std::shared_ptr<const MeshType>& mesh, const CellValue<DataType>& cell_value)
+ : m_mesh{mesh}, m_cell_values{cell_value}
+ {
+ ;
+ }
+
DiscreteFunctionP0(const DiscreteFunctionP0&) noexcept = default;
DiscreteFunctionP0(DiscreteFunctionP0&&) noexcept = default;
diff --git a/src/scheme/DiscreteFunctionUtils.cpp b/src/scheme/DiscreteFunctionUtils.cpp
new file mode 100644
index 000000000..56453a746
--- /dev/null
+++ b/src/scheme/DiscreteFunctionUtils.cpp
@@ -0,0 +1,115 @@
+#include <scheme/DiscreteFunctionUtils.hpp>
+
+#include <mesh/Connectivity.hpp>
+#include <mesh/IMesh.hpp>
+#include <mesh/Mesh.hpp>
+#include <scheme/DiscreteFunctionP0.hpp>
+
+template <size_t Dimension, typename DataType>
+std::shared_ptr<const IDiscreteFunction>
+shallowCopy(const std::shared_ptr<const Mesh<Connectivity<Dimension>>>& mesh,
+ const std::shared_ptr<const DiscreteFunctionP0<Dimension, DataType>>& discrete_function)
+{
+ return std::make_shared<DiscreteFunctionP0<Dimension, DataType>>(mesh, discrete_function->cellValues());
+}
+
+template <size_t Dimension>
+std::shared_ptr<const IDiscreteFunction>
+shallowCopy(const std::shared_ptr<const Mesh<Connectivity<Dimension>>>& mesh,
+ const std::shared_ptr<const IDiscreteFunction>& discrete_function)
+{
+ const std::shared_ptr function_mesh =
+ std::dynamic_pointer_cast<const Mesh<Connectivity<Dimension>>>(discrete_function->mesh());
+
+ if (mesh->shared_connectivity() != function_mesh->shared_connectivity()) {
+ throw NormalError("incompatible connectivities");
+ }
+
+ switch (discrete_function->descriptor().type()) {
+ case DiscreteFunctionType::P0: {
+ switch (discrete_function->dataType()) {
+ case ASTNodeDataType::double_t: {
+ return shallowCopy(mesh,
+ std::dynamic_pointer_cast<const DiscreteFunctionP0<Dimension, double>>(discrete_function));
+ }
+ case ASTNodeDataType::vector_t: {
+ switch (discrete_function->dataType().dimension()) {
+ case 1: {
+ return shallowCopy(mesh, std::dynamic_pointer_cast<const DiscreteFunctionP0<Dimension, TinyVector<1>>>(
+ discrete_function));
+ }
+ case 2: {
+ return shallowCopy(mesh, std::dynamic_pointer_cast<const DiscreteFunctionP0<Dimension, TinyVector<2>>>(
+ discrete_function));
+ }
+ case 3: {
+ return shallowCopy(mesh, std::dynamic_pointer_cast<const DiscreteFunctionP0<Dimension, TinyVector<3>>>(
+ discrete_function));
+ }
+ default: {
+ throw UnexpectedError("invalid data vector dimension: " +
+ std::to_string(discrete_function->dataType().dimension()));
+ }
+ }
+ }
+ case ASTNodeDataType::matrix_t: {
+ if (discrete_function->dataType().nbRows() != discrete_function->dataType().nbColumns()) {
+ throw UnexpectedError(
+ "invalid data matrix dimensions: " + std::to_string(discrete_function->dataType().nbRows()) + "x" +
+ std::to_string(discrete_function->dataType().nbColumns()));
+ }
+ switch (discrete_function->dataType().nbRows()) {
+ case 1: {
+ return shallowCopy(mesh, std::dynamic_pointer_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<1>>>(
+ discrete_function));
+ }
+ case 2: {
+ return shallowCopy(mesh, std::dynamic_pointer_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<2>>>(
+ discrete_function));
+ }
+ case 3: {
+ return shallowCopy(mesh, std::dynamic_pointer_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<3>>>(
+ discrete_function));
+ }
+ default: {
+ throw UnexpectedError(
+ "invalid data matrix dimensions: " + std::to_string(discrete_function->dataType().nbRows()) + "x" +
+ std::to_string(discrete_function->dataType().nbColumns()));
+ }
+ }
+ }
+ default: {
+ throw UnexpectedError("invalid kind of P0 function: invalid data type");
+ }
+ }
+ }
+ default: {
+ throw NormalError("invalid discretization type");
+ }
+ }
+}
+
+std::shared_ptr<const IDiscreteFunction>
+shallowCopy(const std::shared_ptr<const IMesh>& mesh, const std::shared_ptr<const IDiscreteFunction>& discrete_function)
+{
+ if (mesh == discrete_function->mesh()) {
+ return discrete_function;
+ } else if (mesh->dimension() != discrete_function->mesh()->dimension()) {
+ throw NormalError("incompatible mesh dimensions");
+ }
+
+ switch (mesh->dimension()) {
+ case 1: {
+ return shallowCopy(std::dynamic_pointer_cast<const Mesh<Connectivity<1>>>(mesh), discrete_function);
+ }
+ case 2: {
+ return shallowCopy(std::dynamic_pointer_cast<const Mesh<Connectivity<2>>>(mesh), discrete_function);
+ }
+ case 3: {
+ return shallowCopy(std::dynamic_pointer_cast<const Mesh<Connectivity<3>>>(mesh), discrete_function);
+ }
+ default: {
+ throw UnexpectedError("invalid mesh dimension");
+ }
+ }
+}
diff --git a/src/scheme/DiscreteFunctionUtils.hpp b/src/scheme/DiscreteFunctionUtils.hpp
new file mode 100644
index 000000000..1eea93a7e
--- /dev/null
+++ b/src/scheme/DiscreteFunctionUtils.hpp
@@ -0,0 +1,43 @@
+#ifndef DISCRETE_FUNCTION_UTILS_HPP
+#define DISCRETE_FUNCTION_UTILS_HPP
+
+#include <scheme/DiscreteFunctionType.hpp>
+#include <scheme/IDiscreteFunction.hpp>
+#include <scheme/IDiscreteFunctionDescriptor.hpp>
+
+#include <vector>
+
+PUGS_INLINE
+bool
+checkDiscretizationType(const std::vector<std::shared_ptr<const IDiscreteFunction>>& discrete_function_list,
+ const DiscreteFunctionType& discrete_function_type)
+{
+ for (auto discrete_function : discrete_function_list) {
+ if (discrete_function->descriptor().type() != discrete_function_type) {
+ return false;
+ }
+ }
+ return true;
+}
+
+PUGS_INLINE
+std::shared_ptr<const IMesh>
+getCommonMesh(const std::vector<std::shared_ptr<const IDiscreteFunction>>& discrete_function_list)
+{
+ std::shared_ptr<const IMesh> i_mesh;
+ for (auto discrete_function : discrete_function_list) {
+ if (not i_mesh.use_count()) {
+ i_mesh = discrete_function->mesh();
+ } else {
+ if (i_mesh != discrete_function->mesh()) {
+ return {};
+ }
+ }
+ }
+ return i_mesh;
+}
+
+std::shared_ptr<const IDiscreteFunction> shallowCopy(const std::shared_ptr<const IMesh>& mesh,
+ const std::shared_ptr<const IDiscreteFunction>& discrete_function);
+
+#endif // DISCRETE_FUNCTION_UTILS_HPP
diff --git a/src/scheme/PerfectGas.cpp b/src/scheme/PerfectGas.cpp
new file mode 100644
index 000000000..b371dea3c
--- /dev/null
+++ b/src/scheme/PerfectGas.cpp
@@ -0,0 +1,248 @@
+#include <scheme/PerfectGas.hpp>
+
+#include <scheme/DiscreteFunctionP0.hpp>
+#include <scheme/DiscreteFunctionUtils.hpp>
+
+namespace perfect_gas
+{
+template <size_t Dimension, typename DataType>
+std::shared_ptr<IDiscreteFunction>
+soundSpeed(const DiscreteFunctionP0<Dimension, DataType>& rho,
+ const DiscreteFunctionP0<Dimension, DataType>& p,
+ const DiscreteFunctionP0<Dimension, DataType>& gamma)
+{
+ using MeshType = Mesh<Connectivity<Dimension>>;
+ std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(rho.mesh());
+
+ DiscreteFunctionP0<Dimension, DataType> sound_speed{mesh};
+
+ parallel_for(
+ mesh->numberOfCells(),
+ PUGS_LAMBDA(CellId cell_id) { sound_speed[cell_id] = std::sqrt(gamma[cell_id] * p[cell_id] / rho[cell_id]); });
+
+ return std::make_shared<DiscreteFunctionP0<Dimension, DataType>>(sound_speed);
+}
+
+std::shared_ptr<IDiscreteFunction>
+soundSpeed(const std::shared_ptr<const IDiscreteFunction>& rho,
+ const std::shared_ptr<const IDiscreteFunction>& p,
+ const std::shared_ptr<const IDiscreteFunction>& gamma)
+{
+ std::shared_ptr<const IMesh> mesh = getCommonMesh({rho, p, gamma});
+ if (not mesh.use_count()) {
+ throw NormalError("rho, p and gamma are not defined on the same mesh");
+ }
+
+ switch (mesh->dimension()) {
+ case 1: {
+ return soundSpeed(dynamic_cast<const DiscreteFunctionP0<1, double>&>(*rho),
+ dynamic_cast<const DiscreteFunctionP0<1, double>&>(*p),
+ dynamic_cast<const DiscreteFunctionP0<1, double>&>(*gamma));
+ }
+ case 2: {
+ return soundSpeed(dynamic_cast<const DiscreteFunctionP0<2, double>&>(*rho),
+ dynamic_cast<const DiscreteFunctionP0<2, double>&>(*p),
+ dynamic_cast<const DiscreteFunctionP0<2, double>&>(*gamma));
+ }
+ case 3: {
+ return soundSpeed(dynamic_cast<const DiscreteFunctionP0<3, double>&>(*rho),
+ dynamic_cast<const DiscreteFunctionP0<3, double>&>(*p),
+ dynamic_cast<const DiscreteFunctionP0<3, double>&>(*gamma));
+ }
+ default: {
+ throw UnexpectedError("invalid mesh dimension");
+ }
+ }
+}
+
+template <size_t Dimension, typename DataType>
+std::shared_ptr<IDiscreteFunction>
+epsilonFromRhoPGamma(const DiscreteFunctionP0<Dimension, DataType>& rho,
+ const DiscreteFunctionP0<Dimension, DataType>& p,
+ const DiscreteFunctionP0<Dimension, DataType>& gamma)
+{
+ using MeshType = Mesh<Connectivity<Dimension>>;
+ std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(rho.mesh());
+
+ DiscreteFunctionP0<Dimension, DataType> epsilon{mesh};
+
+ parallel_for(
+ mesh->numberOfCells(),
+ PUGS_LAMBDA(CellId cell_id) { epsilon[cell_id] = p[cell_id] / ((gamma[cell_id] - 1) * rho[cell_id]); });
+
+ return std::make_shared<DiscreteFunctionP0<Dimension, DataType>>(epsilon);
+}
+
+std::shared_ptr<IDiscreteFunction>
+epsilonFromRhoPGamma(const std::shared_ptr<const IDiscreteFunction>& rho,
+ const std::shared_ptr<const IDiscreteFunction>& p,
+ const std::shared_ptr<const IDiscreteFunction>& gamma)
+{
+ std::shared_ptr<const IMesh> mesh = getCommonMesh({rho, p, gamma});
+ if (not mesh.use_count()) {
+ throw NormalError("rho, p and gamma are not defined on the same mesh");
+ }
+
+ switch (mesh->dimension()) {
+ case 1: {
+ return epsilonFromRhoPGamma(dynamic_cast<const DiscreteFunctionP0<1, double>&>(*rho),
+ dynamic_cast<const DiscreteFunctionP0<1, double>&>(*p),
+ dynamic_cast<const DiscreteFunctionP0<1, double>&>(*gamma));
+ }
+ case 2: {
+ return epsilonFromRhoPGamma(dynamic_cast<const DiscreteFunctionP0<2, double>&>(*rho),
+ dynamic_cast<const DiscreteFunctionP0<2, double>&>(*p),
+ dynamic_cast<const DiscreteFunctionP0<2, double>&>(*gamma));
+ }
+ case 3: {
+ return epsilonFromRhoPGamma(dynamic_cast<const DiscreteFunctionP0<3, double>&>(*rho),
+ dynamic_cast<const DiscreteFunctionP0<3, double>&>(*p),
+ dynamic_cast<const DiscreteFunctionP0<3, double>&>(*gamma));
+ }
+ default: {
+ throw UnexpectedError("invalid mesh dimension");
+ }
+ }
+}
+
+template <size_t Dimension, typename DataType>
+std::shared_ptr<IDiscreteFunction>
+pFromRhoEpsilonGamma(const DiscreteFunctionP0<Dimension, DataType>& rho,
+ const DiscreteFunctionP0<Dimension, DataType>& epsilon,
+ const DiscreteFunctionP0<Dimension, DataType>& gamma)
+{
+ using MeshType = Mesh<Connectivity<Dimension>>;
+ std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(rho.mesh());
+
+ DiscreteFunctionP0<Dimension, DataType> p{mesh};
+
+ parallel_for(
+ mesh->numberOfCells(),
+ PUGS_LAMBDA(CellId cell_id) { p[cell_id] = ((gamma[cell_id] - 1) * rho[cell_id]) * epsilon[cell_id]; });
+
+ return std::make_shared<DiscreteFunctionP0<Dimension, DataType>>(p);
+}
+
+std::shared_ptr<IDiscreteFunction>
+pFromRhoEpsilonGamma(const std::shared_ptr<const IDiscreteFunction>& rho,
+ const std::shared_ptr<const IDiscreteFunction>& epsilon,
+ const std::shared_ptr<const IDiscreteFunction>& gamma)
+{
+ std::shared_ptr<const IMesh> mesh = getCommonMesh({rho, epsilon, gamma});
+ if (not mesh.use_count()) {
+ throw NormalError("rho, epsilon and gamma are not defined on the same mesh");
+ }
+
+ switch (mesh->dimension()) {
+ case 1: {
+ return pFromRhoEpsilonGamma(dynamic_cast<const DiscreteFunctionP0<1, double>&>(*rho),
+ dynamic_cast<const DiscreteFunctionP0<1, double>&>(*epsilon),
+ dynamic_cast<const DiscreteFunctionP0<1, double>&>(*gamma));
+ }
+ case 2: {
+ return pFromRhoEpsilonGamma(dynamic_cast<const DiscreteFunctionP0<2, double>&>(*rho),
+ dynamic_cast<const DiscreteFunctionP0<2, double>&>(*epsilon),
+ dynamic_cast<const DiscreteFunctionP0<2, double>&>(*gamma));
+ }
+ case 3: {
+ return pFromRhoEpsilonGamma(dynamic_cast<const DiscreteFunctionP0<3, double>&>(*rho),
+ dynamic_cast<const DiscreteFunctionP0<3, double>&>(*epsilon),
+ dynamic_cast<const DiscreteFunctionP0<3, double>&>(*gamma));
+ }
+ default: {
+ throw UnexpectedError("invalid mesh dimension");
+ }
+ }
+}
+
+template <size_t Dimension, typename DataType>
+std::shared_ptr<IDiscreteFunction>
+totalEnergyFromEpsilonU(const DiscreteFunctionP0<Dimension, DataType>& epsilon,
+ const DiscreteFunctionP0<Dimension, TinyVector<Dimension, DataType>>& u)
+{
+ using MeshType = Mesh<Connectivity<Dimension>>;
+ std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(epsilon.mesh());
+
+ DiscreteFunctionP0<Dimension, DataType> E{mesh};
+
+ parallel_for(
+ mesh->numberOfCells(),
+ PUGS_LAMBDA(CellId cell_id) { E[cell_id] = epsilon[cell_id] + 0.5 * (u[cell_id], u[cell_id]); });
+
+ return std::make_shared<DiscreteFunctionP0<Dimension, DataType>>(epsilon);
+}
+
+std::shared_ptr<IDiscreteFunction>
+totalEnergyFromEpsilonU(const std::shared_ptr<const IDiscreteFunction>& epsilon,
+ const std::shared_ptr<const IDiscreteFunction>& u)
+{
+ std::shared_ptr<const IMesh> mesh = getCommonMesh({epsilon, u});
+ if (not mesh.use_count()) {
+ throw NormalError("epsilon and u are not defined on the same mesh");
+ }
+
+ switch (mesh->dimension()) {
+ case 1: {
+ return totalEnergyFromEpsilonU(dynamic_cast<const DiscreteFunctionP0<1, double>&>(*epsilon),
+ dynamic_cast<const DiscreteFunctionP0<1, TinyVector<1, double>>&>(*u));
+ }
+ case 2: {
+ return totalEnergyFromEpsilonU(dynamic_cast<const DiscreteFunctionP0<2, double>&>(*epsilon),
+ dynamic_cast<const DiscreteFunctionP0<2, TinyVector<2, double>>&>(*u));
+ }
+ case 3: {
+ return totalEnergyFromEpsilonU(dynamic_cast<const DiscreteFunctionP0<3, double>&>(*epsilon),
+ dynamic_cast<const DiscreteFunctionP0<3, TinyVector<3, double>>&>(*u));
+ }
+ default: {
+ throw UnexpectedError("invalid mesh dimension");
+ }
+ }
+}
+
+template <size_t Dimension, typename DataType>
+std::shared_ptr<IDiscreteFunction>
+epsilonFromTotalEnergyU(const DiscreteFunctionP0<Dimension, DataType>& E,
+ const DiscreteFunctionP0<Dimension, TinyVector<Dimension, DataType>>& u)
+{
+ using MeshType = Mesh<Connectivity<Dimension>>;
+ std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(u.mesh());
+
+ DiscreteFunctionP0<Dimension, DataType> epsilon{mesh};
+
+ parallel_for(
+ mesh->numberOfCells(),
+ PUGS_LAMBDA(CellId cell_id) { epsilon[cell_id] = E[cell_id] - 0.5 * (u[cell_id], u[cell_id]); });
+
+ return std::make_shared<DiscreteFunctionP0<Dimension, DataType>>(epsilon);
+}
+
+std::shared_ptr<IDiscreteFunction>
+epsilonFromTotalEnergyU(const std::shared_ptr<const IDiscreteFunction>& E,
+ const std::shared_ptr<const IDiscreteFunction>& u)
+{
+ std::shared_ptr<const IMesh> mesh = getCommonMesh({E, u});
+ if (not mesh.use_count()) {
+ throw NormalError("E and u are not defined on the same mesh");
+ }
+
+ switch (mesh->dimension()) {
+ case 1: {
+ return epsilonFromTotalEnergyU(dynamic_cast<const DiscreteFunctionP0<1, double>&>(*E),
+ dynamic_cast<const DiscreteFunctionP0<1, TinyVector<1, double>>&>(*u));
+ }
+ case 2: {
+ return epsilonFromTotalEnergyU(dynamic_cast<const DiscreteFunctionP0<2, double>&>(*E),
+ dynamic_cast<const DiscreteFunctionP0<2, TinyVector<2, double>>&>(*u));
+ }
+ case 3: {
+ return epsilonFromTotalEnergyU(dynamic_cast<const DiscreteFunctionP0<3, double>&>(*E),
+ dynamic_cast<const DiscreteFunctionP0<3, TinyVector<3, double>>&>(*u));
+ }
+ default: {
+ throw UnexpectedError("invalid mesh dimension");
+ }
+ }
+}
+
+} // namespace perfect_gas
diff --git a/src/scheme/PerfectGas.hpp b/src/scheme/PerfectGas.hpp
new file mode 100644
index 000000000..fb7b9d212
--- /dev/null
+++ b/src/scheme/PerfectGas.hpp
@@ -0,0 +1,32 @@
+#ifndef PERFECT_GAS_HPP
+#define PERFECT_GAS_HPP
+
+#include <scheme/IDiscreteFunction.hpp>
+
+namespace perfect_gas
+{
+std::shared_ptr<IDiscreteFunction> soundSpeed(const std::shared_ptr<const IDiscreteFunction>& rho,
+ const std::shared_ptr<const IDiscreteFunction>& p,
+ const std::shared_ptr<const IDiscreteFunction>& gamma);
+
+std::shared_ptr<IDiscreteFunction> epsilonFromRhoPGamma(const std::shared_ptr<const IDiscreteFunction>& rho,
+ const std::shared_ptr<const IDiscreteFunction>& p,
+ const std::shared_ptr<const IDiscreteFunction>& gamma);
+
+std::shared_ptr<IDiscreteFunction> pFromRhoEpsilonGamma(const std::shared_ptr<const IDiscreteFunction>& rho,
+ const std::shared_ptr<const IDiscreteFunction>& epsilon,
+ const std::shared_ptr<const IDiscreteFunction>& gamma);
+
+// This is a temporary function
+// Todo: TO REMOVE
+std::shared_ptr<IDiscreteFunction> totalEnergyFromEpsilonU(const std::shared_ptr<const IDiscreteFunction>& epsilon,
+ const std::shared_ptr<const IDiscreteFunction>& u);
+
+// This is a temporary function
+// Todo: TO REMOVE
+std::shared_ptr<IDiscreteFunction> epsilonFromTotalEnergyU(const std::shared_ptr<const IDiscreteFunction>& E,
+ const std::shared_ptr<const IDiscreteFunction>& u);
+
+} // namespace perfect_gas
+
+#endif // PERFECT_GAS_HPP
--
GitLab