diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index 2491767f4dfaa9e988a4960d625f32077d32311d..477e3cc80cf439ae369cb1bd5074e53116eda91b 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -93,8 +93,8 @@ class InterpolateItemValue<OutputType(InputType)> : public PugsFunctionAdapter<O
}
};
-template <size_t Dimension>
-struct GlaceScheme
+template <size_t Dimension, AcousticSolverType solver_type>
+struct AcousticSolverScheme
{
using ConnectivityType = Connectivity<Dimension>;
using MeshType = Mesh<ConnectivityType>;
@@ -103,11 +103,11 @@ struct GlaceScheme
std::shared_ptr<const MeshType> m_mesh;
- GlaceScheme(std::shared_ptr<const IMesh> i_mesh,
- const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
- const FunctionSymbolId& rho_id,
- const FunctionSymbolId& u_id,
- const FunctionSymbolId& p_id)
+ AcousticSolverScheme(std::shared_ptr<const IMesh> i_mesh,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
+ const FunctionSymbolId& rho_id,
+ const FunctionSymbolId& u_id,
+ const FunctionSymbolId& p_id)
: m_mesh{std::dynamic_pointer_cast<const MeshType>(i_mesh)}
{
std::cout << "number of bc descr = " << bc_descriptor_list.size() << '\n';
@@ -160,34 +160,6 @@ struct GlaceScheme
TinyVector<Dimension>)>::template interpolate<ItemType::node>(velocity_id, m_mesh->xr(), node_list);
bc_list.push_back(VelocityBoundaryCondition{node_list, value_list});
-
- // if constexpr (Dimension == 1) {
- // const auto& node_list = ref_face_list.list();
-
- // Array<const TinyVector<Dimension>> value_list = InterpolateItemValue<TinyVector<Dimension>(
- // TinyVector<Dimension>)>::template interpolate<FaceType>(velocity_id, m_mesh->xr(), node_list);
-
- // bc_list.push_back(VelocityBoundaryCondition{node_list, value_list});
- // } else {
- // const auto& face_list = ref_face_list.list();
- // const auto& face_to_node_matrix = m_mesh->connectivity().faceToNodeMatrix();
- // std::set<NodeId> node_set;
- // for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
- // FaceId face_id = face_list[i_face];
- // const auto& face_nodes = face_to_node_matrix[face_id];
- // for (size_t i_node = 0; i_node < face_nodes.size(); ++i_node) {
- // node_set.insert(face_nodes[i_node]);
- // }
- // }
-
- // Array<NodeId> node_list = convert_to_array(node_set);
-
- // Array<const TinyVector<Dimension>> value_list = InterpolateItemValue<TinyVector<Dimension>(
- // TinyVector<Dimension>)>::template interpolate<ItemType::node>(velocity_id, m_mesh->xr(),
- // node_list);
-
- // bc_list.push_back(VelocityBoundaryCondition{node_list, value_list});
- // }
}
}
break;
@@ -249,7 +221,7 @@ struct GlaceScheme
}
unknowns.gammaj().fill(1.4);
- AcousticSolver acoustic_solver(m_mesh, bc_list);
+ AcousticSolver acoustic_solver(m_mesh, bc_list, solver_type);
const double tmax = 0.2;
double t = 0;
@@ -290,12 +262,13 @@ struct GlaceScheme
while ((t < tmax) and (iteration < itermax)) {
MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*m_mesh);
- // vtk_writer.write(m_mesh,
- // {NamedItemValue{"density", rhoj}, NamedItemValue{"velocity", unknowns.uj()},
- // NamedItemValue{"coords", m_mesh->xr()}, NamedItemValue{"xj", mesh_data.xj()},
- // NamedItemValue{"cell_owner", m_mesh->connectivity().cellOwner()},
- // NamedItemValue{"node_owner", m_mesh->connectivity().nodeOwner()}},
- // t);
+ vtk_writer.write(m_mesh,
+ {NamedItemValue{"density", rhoj}, NamedItemValue{"velocity", unknowns.uj()},
+ NamedItemValue{"coords", m_mesh->xr()}, NamedItemValue{"xj", mesh_data.xj()},
+ NamedItemValue{"cell_owner", m_mesh->connectivity().cellOwner()},
+ NamedItemValue{"node_owner", m_mesh->connectivity().nodeOwner()}},
+ t);
+
double dt = 0.4 * acoustic_solver.acoustic_dt(mesh_data.Vj(), cj);
if (t + dt > tmax) {
dt = tmax - t;
@@ -392,35 +365,63 @@ SchemeModule::SchemeModule()
));
- this->_addBuiltinFunction("glace",
- std::make_shared<BuiltinFunctionEmbedder<
- void(std::shared_ptr<const IMesh>,
- const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
- const FunctionSymbolId&, const FunctionSymbolId&, const FunctionSymbolId&)>>(
-
- [](std::shared_ptr<const IMesh> p_mesh,
- const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
- bc_descriptor_list,
- const FunctionSymbolId& rho_id, const FunctionSymbolId& u_id,
- const FunctionSymbolId& p_id) -> void {
- switch (p_mesh->dimension()) {
- case 1: {
- GlaceScheme<1>{p_mesh, bc_descriptor_list, rho_id, u_id, p_id};
- break;
- }
- case 2: {
- GlaceScheme<2>{p_mesh, bc_descriptor_list, rho_id, u_id, p_id};
- break;
- }
- case 3: {
- GlaceScheme<3>{p_mesh, bc_descriptor_list, rho_id, u_id, p_id};
- break;
- }
- default: {
- throw UnexpectedError("invalid mesh dimension");
- }
- }
- }
+ this->_addBuiltinFunction(
+ "glace",
+ std::make_shared<BuiltinFunctionEmbedder<
+ void(std::shared_ptr<const IMesh>, const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
+ const FunctionSymbolId&, const FunctionSymbolId&, const FunctionSymbolId&)>>(
+
+ [](std::shared_ptr<const IMesh> p_mesh,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
+ const FunctionSymbolId& rho_id, const FunctionSymbolId& u_id, const FunctionSymbolId& p_id) -> void {
+ switch (p_mesh->dimension()) {
+ case 1: {
+ AcousticSolverScheme<1, AcousticSolverType::Glace>{p_mesh, bc_descriptor_list, rho_id, u_id, p_id};
+ break;
+ }
+ case 2: {
+ AcousticSolverScheme<2, AcousticSolverType::Glace>{p_mesh, bc_descriptor_list, rho_id, u_id, p_id};
+ break;
+ }
+ case 3: {
+ AcousticSolverScheme<3, AcousticSolverType::Glace>{p_mesh, bc_descriptor_list, rho_id, u_id, p_id};
+ break;
+ }
+ default: {
+ throw UnexpectedError("invalid mesh dimension");
+ }
+ }
+ }
- ));
+ ));
+
+ this->_addBuiltinFunction(
+ "eucclhyd",
+ std::make_shared<BuiltinFunctionEmbedder<
+ void(std::shared_ptr<const IMesh>, const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
+ const FunctionSymbolId&, const FunctionSymbolId&, const FunctionSymbolId&)>>(
+
+ [](std::shared_ptr<const IMesh> p_mesh,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
+ const FunctionSymbolId& rho_id, const FunctionSymbolId& u_id, const FunctionSymbolId& p_id) -> void {
+ switch (p_mesh->dimension()) {
+ case 1: {
+ AcousticSolverScheme<1, AcousticSolverType::Eucclhyd>{p_mesh, bc_descriptor_list, rho_id, u_id, p_id};
+ break;
+ }
+ case 2: {
+ AcousticSolverScheme<2, AcousticSolverType::Eucclhyd>{p_mesh, bc_descriptor_list, rho_id, u_id, p_id};
+ break;
+ }
+ case 3: {
+ AcousticSolverScheme<3, AcousticSolverType::Eucclhyd>{p_mesh, bc_descriptor_list, rho_id, u_id, p_id};
+ break;
+ }
+ default: {
+ throw UnexpectedError("invalid mesh dimension");
+ }
+ }
+ }
+
+ ));
}
diff --git a/src/mesh/MeshData.hpp b/src/mesh/MeshData.hpp
index fa5f2bf63c5e27ca79f0684dee98b2f9c32d249d..aeddc582f09332543df40fbcb11e368f1d7e929e 100644
--- a/src/mesh/MeshData.hpp
+++ b/src/mesh/MeshData.hpp
@@ -32,6 +32,8 @@ class MeshData : public IMeshData
private:
const MeshType& m_mesh;
+ NodeValuePerFace<const Rd> m_Nlr;
+ NodeValuePerFace<const Rd> m_nlr;
NodeValuePerCell<const Rd> m_Cjr;
NodeValuePerCell<const double> m_ljr;
NodeValuePerCell<const Rd> m_njr;
@@ -39,7 +41,6 @@ class MeshData : public IMeshData
CellValue<const Rd> m_xj;
CellValue<const double> m_Vj;
FaceValue<const double> m_ll;
- NodeValuePerFace<const Rd> m_Nlr;
PUGS_INLINE
void
@@ -48,7 +49,7 @@ class MeshData : public IMeshData
if constexpr (Dimension == 1) {
static_assert(Dimension != 1, "ll does not make sense in 1d");
} else {
- const auto& nlr = this->Nlr();
+ const auto& Nlr = this->Nlr();
FaceValue<double> ll{m_mesh.connectivity()};
@@ -59,7 +60,7 @@ class MeshData : public IMeshData
double lenght = 0;
for (size_t i_node = 0; i_node < face_nodes.size(); ++i_node) {
- lenght += l2Norm(nlr(face_id, i_node));
+ lenght += l2Norm(Nlr(face_id, i_node));
}
ll[face_id] = lenght;
@@ -69,6 +70,27 @@ class MeshData : public IMeshData
}
}
+ PUGS_INLINE
+ void
+ _compute_nlr()
+ {
+ if constexpr (Dimension == 1) {
+ static_assert(Dimension != 1, "nlr does not make sense in 1d");
+ } else {
+ const auto& Nlr = this->Nlr();
+
+ NodeValuePerFace<Rd> nlr{m_mesh.connectivity()};
+
+ parallel_for(
+ Nlr.numberOfValues(), PUGS_LAMBDA(size_t lr) {
+ double length = l2Norm(Nlr[lr]);
+ nlr[lr] = 1. / length * Nlr[lr];
+ });
+
+ m_nlr = nlr;
+ }
+ }
+
PUGS_INLINE void
_computeFaceIsobarycenter()
{ // Computes vertices isobarycenter
@@ -368,6 +390,16 @@ class MeshData : public IMeshData
return m_Nlr;
}
+ PUGS_INLINE
+ NodeValuePerFace<const Rd>
+ nlr()
+ {
+ if (not m_nlr.isBuilt()) {
+ this->_compute_nlr();
+ }
+ return m_nlr;
+ }
+
PUGS_INLINE
NodeValuePerCell<const Rd>
Cjr()
diff --git a/src/scheme/AcousticSolver.hpp b/src/scheme/AcousticSolver.hpp
index 88007ac4709f2163bbd0a1737bebbf5633ea074a..1e246f2dbf054715d67e80819078ed84b9441785 100644
--- a/src/scheme/AcousticSolver.hpp
+++ b/src/scheme/AcousticSolver.hpp
@@ -20,6 +20,12 @@
#include <iostream>
+enum class AcousticSolverType
+{
+ Eucclhyd,
+ Glace
+};
+
template <typename MeshType>
class AcousticSolver
{
@@ -47,6 +53,8 @@ class AcousticSolver
const BoundaryConditionList m_boundary_condition_list;
+ const AcousticSolverType m_solver_type;
+
void
_applyPressureBC()
{
@@ -173,11 +181,13 @@ class AcousticSolver
PUGS_INLINE
void
- computeAjr(const CellValue<const double>& rhocj,
- const NodeValuePerCell<const Rd>& Cjr,
- const NodeValuePerCell<const double>& /* ljr */,
- const NodeValuePerCell<const Rd>& njr)
+ _computeGlaceAjr(const CellValue<const double>& rhocj)
{
+ MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*m_mesh);
+
+ const NodeValuePerCell<const Rd> Cjr = mesh_data.Cjr();
+ const NodeValuePerCell<const Rd> njr = mesh_data.njr();
+
parallel_for(
m_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
const size_t& nb_nodes = m_Ajr.numberOfSubValues(j);
@@ -188,6 +198,71 @@ class AcousticSolver
});
}
+ PUGS_INLINE
+ void
+ _computeEucclhydAjr(const CellValue<const double>& rhocj)
+ {
+ MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*m_mesh);
+
+ const NodeValuePerFace<const Rd> Nlr = mesh_data.Nlr();
+ const NodeValuePerFace<const Rd> nlr = mesh_data.nlr();
+
+ const auto& face_to_node_matrix = m_connectivity.faceToNodeMatrix();
+ const auto& cell_to_node_matrix = m_connectivity.cellToNodeMatrix();
+ const auto& cell_to_face_matrix = m_connectivity.cellToFaceMatrix();
+
+ parallel_for(
+ m_Ajr.numberOfValues(), PUGS_LAMBDA(size_t jr) { m_Ajr[jr] = zero; });
+
+ parallel_for(
+ m_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
+ const auto& cell_nodes = cell_to_node_matrix[j];
+
+ const auto& cell_faces = cell_to_face_matrix[j];
+
+ const double& rho_c = rhocj[j];
+
+ for (size_t L = 0; L < cell_faces.size(); ++L) {
+ const FaceId& l = cell_faces[L];
+ const auto& face_nodes = face_to_node_matrix[l];
+
+ auto local_node_number_in_cell = [&](NodeId node_number) {
+ for (size_t i_node = 0; i_node < cell_nodes.size(); ++i_node) {
+ if (node_number == cell_nodes[i_node]) {
+ return i_node;
+ }
+ }
+ return std::numeric_limits<size_t>::max();
+ };
+
+ for (size_t rl = 0; rl < face_nodes.size(); ++rl) {
+ const size_t R = local_node_number_in_cell(face_nodes[rl]);
+ m_Ajr(j, R) += tensorProduct(rho_c * Nlr(l, rl), nlr(l, rl));
+ }
+ }
+ });
+ }
+
+ PUGS_INLINE
+ void
+ computeAjr(const CellValue<const double>& rhocj)
+ {
+ switch (m_solver_type) {
+ case AcousticSolverType::Glace: {
+ this->_computeGlaceAjr(rhocj);
+ break;
+ }
+ case AcousticSolverType::Eucclhyd: {
+ if constexpr (Dimension > 1) {
+ this->_computeEucclhydAjr(rhocj);
+ } else {
+ this->_computeGlaceAjr(rhocj);
+ }
+ break;
+ }
+ }
+ }
+
PUGS_INLINE
const NodeValue<const Rdd>
computeAr(const NodeValuePerCell<const Rdd>& Ajr)
@@ -214,11 +289,13 @@ class AcousticSolver
PUGS_INLINE
const NodeValue<const Rd>
- computeBr(const NodeValuePerCell<Rdd>& Ajr,
- const NodeValuePerCell<const Rd>& Cjr,
- const CellValue<const Rd>& uj,
- const CellValue<const double>& pj)
+ computeBr(const CellValue<const Rd>& uj, const CellValue<const double>& pj)
{
+ MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*m_mesh);
+
+ const NodeValuePerCell<const Rd>& Cjr = mesh_data.Cjr();
+ const NodeValuePerCell<const Rdd>& Ajr = m_Ajr;
+
const auto& node_to_cell_matrix = m_connectivity.nodeToCellMatrix();
const auto& node_local_numbers_in_their_cells = m_connectivity.nodeLocalNumbersInTheirCells();
@@ -246,24 +323,28 @@ class AcousticSolver
this->_applyVelocityBC();
}
- NodeValue<Rd>
- computeUr(const NodeValue<const Rdd>& Ar, const NodeValue<const Rd>& br)
+ void
+ computeUr()
{
+ const NodeValue<const Rdd> Ar = m_Ar;
+ const NodeValue<const Rd> br = m_br;
+
inverse(Ar, m_inv_Ar);
const NodeValue<const Rdd> invAr = m_inv_Ar;
parallel_for(
m_mesh->numberOfNodes(), PUGS_LAMBDA(NodeId r) { m_ur[r] = invAr[r] * br[r]; });
-
- return m_ur;
}
void
- computeFjr(const NodeValuePerCell<Rdd>& Ajr,
- const NodeValue<const Rd>& ur,
- const NodeValuePerCell<const Rd>& Cjr,
- const CellValue<const Rd>& uj,
- const CellValue<const double>& pj)
+ computeFjr(const CellValue<const Rd>& uj, const CellValue<const double>& pj)
{
+ MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*m_mesh);
+
+ const NodeValuePerCell<const Rd> Cjr = mesh_data.Cjr();
+ const NodeValuePerCell<const Rdd> Ajr = m_Ajr;
+
+ const NodeValue<const Rd> ur = m_ur;
+
const auto& cell_to_node_matrix = m_mesh->connectivity().cellToNodeMatrix();
parallel_for(
@@ -288,26 +369,22 @@ class AcousticSolver
computeExplicitFluxes(const CellValue<const double>& rhoj,
const CellValue<const Rd>& uj,
const CellValue<const double>& pj,
- const CellValue<const double>& cj,
- const NodeValuePerCell<const Rd>& Cjr,
- const NodeValuePerCell<const double>& ljr,
- const NodeValuePerCell<const Rd>& njr)
+ const CellValue<const double>& cj)
{
const CellValue<const double> rhocj = computeRhoCj(rhoj, cj);
- computeAjr(rhocj, Cjr, ljr, njr);
+ computeAjr(rhocj);
NodeValuePerCell<const Rdd> Ajr = m_Ajr;
this->computeAr(Ajr);
- this->computeBr(m_Ajr, Cjr, uj, pj);
+ this->computeBr(uj, pj);
this->applyBoundaryConditions();
synchronize(m_Ar);
synchronize(m_br);
- NodeValue<Rd>& ur = m_ur;
- ur = computeUr(m_Ar, m_br);
- computeFjr(m_Ajr, ur, Cjr, uj, pj);
+ computeUr();
+ computeFjr(uj, pj);
}
NodeValue<Rd> m_br;
@@ -320,10 +397,13 @@ class AcousticSolver
CellValue<double> m_Vj_over_cj;
public:
- AcousticSolver(std::shared_ptr<const MeshType> p_mesh, const BoundaryConditionList& bc_list)
+ AcousticSolver(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),
+ m_solver_type(solver_type),
m_br(m_connectivity),
m_Ajr(m_connectivity),
m_Ar(m_connectivity),
@@ -372,12 +452,9 @@ class AcousticSolver
MeshData<Dimension>& mesh_data = MeshDataManager::instance().getMeshData(*m_mesh);
- const CellValue<const double> Vj = mesh_data.Vj();
- const NodeValuePerCell<const Rd> Cjr = mesh_data.Cjr();
- const NodeValuePerCell<const double> ljr = mesh_data.ljr();
- const NodeValuePerCell<const Rd> njr = mesh_data.njr();
+ const CellValue<const double> Vj = mesh_data.Vj();
- computeExplicitFluxes(rhoj, uj, pj, cj, Cjr, ljr, njr);
+ computeExplicitFluxes(rhoj, uj, pj, cj);
const NodeValuePerCell<Rd>& Fjr = m_Fjr;
const NodeValue<const Rd> ur = m_ur;