Select Git revision
AcousticSolver.cpp
Stéphane Del Pino authored
The code works in 1D/2D but it is very expensive. Probably due to the allocations in the Polynomial class
AcousticSolver.cpp 39.29 KiB
#include <scheme/AcousticSolver.hpp>
#include <analysis/Polynomial1D.hpp>
#include <analysis/Polynomial1DRootComputer.hpp>
#include <language/utils/InterpolateItemValue.hpp>
#include <mesh/ItemValueUtils.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>;
const SolverType m_solver_type;
BoundaryConditionList m_boundary_condition_list;
NodeValue<const Rd> m_ur;
NodeValuePerCell<const Rd> m_Fjr;
NodeValuePerCell<const Rdxd> m_Ajr;
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;
}
NodeValuePerCell<const Rdxd>
_computeEucclhydAjr(const MeshType& mesh, const CellValue<const double>& rhoc)
{
MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
const NodeValuePerFace<const Rd> Nlr = mesh_data.Nlr();
const NodeValuePerFace<const Rd> nlr = mesh_data.nlr();
const auto& face_to_node_matrix = mesh.connectivity().faceToNodeMatrix();
const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
const auto& cell_to_face_matrix = mesh.connectivity().cellToFaceMatrix();
NodeValuePerCell<Rdxd> Ajr{mesh.connectivity()};
parallel_for(
Ajr.numberOfValues(), PUGS_LAMBDA(size_t jr) { Ajr[jr] = zero; });
parallel_for(
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 = rhoc[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]);
Ajr(j, R) += tensorProduct(rho_c * Nlr(l, rl), nlr(l, rl));
}
}
});
return Ajr;
}
NodeValuePerCell<const Rdxd>
_computeAjr(const MeshType& mesh, const CellValue<const double>& rhoc)
{
if constexpr (Dimension == 1) {
return _computeGlaceAjr(mesh, rhoc);
} else {
switch (m_solver_type) {
case SolverType::Glace: {
return _computeGlaceAjr(mesh, rhoc);
}
case SolverType::Eucclhyd: {
return _computeEucclhydAjr(mesh, rhoc);
}
default: {
throw UnexpectedError("invalid solver type");
}
}
}
}
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;
}
double
_computeMaxDensityDt(const MeshType& mesh,
const DiscreteFunctionP0<Dimension, double>& rho,
const double dt_max) const
{
double dt = dt_max;
MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
const NodeValuePerCell<const Rd> Cjr = mesh_data.Cjr();
const CellValue<const double> Vj = mesh_data.Vj();
const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
CellValue<double> dVj{mesh.connectivity()};
parallel_for(
mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
const auto& cell_nodes = cell_to_node_matrix[j];
double dV = 0;
for (size_t R = 0; R < cell_nodes.size(); ++R) {
const NodeId r = cell_nodes[R];
dV += dot(Cjr(j, R), m_ur[r]);
}
dVj[j] = dV;
});
CellValue<double> Gj{mesh.connectivity()};
if constexpr (Dimension == 2) {
for (CellId j = 0; j < mesh.numberOfCells(); ++j) {
const auto& cell_nodes = cell_to_node_matrix[j];
double G = 0;
for (size_t R = 0; R < cell_nodes.size(); ++R) {
const NodeId r = cell_nodes[R];
const NodeId rp1 = cell_nodes[(R + 1) % cell_nodes.size()];
G -= dot(m_ur[rp1], Rd{m_ur[r][1], -m_ur[r][0]});
}
Gj[j] = 0.5 * G;
}
}
for (CellId j = 0; j < mesh.numberOfCells(); ++j) {
Polynomial1D p_tau({1. / rho[j], dVj[j] / (rho[j] * Vj[j])});
if constexpr (Dimension == 2) {
p_tau = p_tau + Polynomial1D({0, 0, Gj[j] / (rho[j] * Vj[j])});
}
Polynomial1DRootComputer computer{p_tau, 0, dt};
std::optional dt_tau = computer.getFirstRoot();
if (dt_tau.has_value()) {
Assert(dt_tau.value() <= dt);
dt = dt_tau.value();
}
}
if (dt != dt_max) {
std::cout << "volume variation imposes time step\n";
}
if (dt < dt_max) {
dt *= 0.95 * dt;
}
return parallel::allReduceMin(dt);
}
double
_computeMaxEnergyDt(const MeshType& mesh,
const DiscreteScalarFunction& rho,
const DiscreteVectorFunction& u,
const DiscreteScalarFunction& E,
const DiscreteScalarFunction& p,
const double dt_max) const
{
double dt = dt_max;
MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
const NodeValuePerCell<const Rd> Cjr = mesh_data.Cjr();
const CellValue<const double> Vj = mesh_data.Vj();
const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
CellValue<double> dVj{mesh.connectivity()};
parallel_for(
mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { const auto& cell_nodes = cell_to_node_matrix[j];
double dV = 0;
for (size_t R = 0; R < cell_nodes.size(); ++R) {
const NodeId r = cell_nodes[R];
dV += dot(Cjr(j, R), m_ur[r]);
}
dVj[j] = dV;
});
CellValue<double> dSj{mesh.connectivity()};
parallel_for(
mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
const auto& cell_nodes = cell_to_node_matrix[j];
double dS = 0;
for (size_t R = 0; R < cell_nodes.size(); ++R) {
const NodeId r = cell_nodes[R];
Rd du = u[j] - m_ur[r];
dS += dot(m_Ajr(j, R) * du, du);
}
dSj[j] = dS;
});
CellValue<Rd> dPj{mesh.connectivity()};
parallel_for(
mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
const auto& cell_nodes = cell_to_node_matrix[j];
Rd dP = zero;
for (size_t R = 0; R < cell_nodes.size(); ++R) {
const NodeId r = cell_nodes[R];
Rd du = u[j] - m_ur[r];
dP += m_Ajr(j, R) * du;
}
dPj[j] = dP;
});
for (CellId j = 0; j < mesh.numberOfCells(); ++j) {
const double inv_mj = 1 / (rho[j] * Vj[j]);
Polynomial1D p_e(
{E[j] - 0.5 * dot(u[j], u[j]), inv_mj * (dSj[j] - p[j] * dVj[j]), -inv_mj * inv_mj * dot(dPj[j], dPj[j])});
Polynomial1DRootComputer computer{p_e, 0, dt};
std::optional dt_e = computer.getFirstRoot();
if (dt_e.has_value()) {
Assert(dt_e.value() <= dt);
dt = dt_e.value();
}
}
if (dt < dt_max) {
dt *= 0.95 * dt;
}
if (dt != dt_max) {
std::cout << "internal energy variation imposes time step\n";
}
return parallel::allReduceMin(dt); }
double
_computeMaxEntropyDt(const MeshType& mesh,
const DiscreteScalarFunction& rho,
const DiscreteVectorFunction& u,
const DiscreteScalarFunction& p,
const double dt_max) const
{
double dt = dt_max;
MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
const NodeValuePerCell<const Rd> Cjr = mesh_data.Cjr();
const CellValue<const double> Vj = mesh_data.Vj();
const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
CellValue<double> dVj{mesh.connectivity()};
parallel_for(
mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
const auto& cell_nodes = cell_to_node_matrix[j];
double dV = 0;
for (size_t R = 0; R < cell_nodes.size(); ++R) {
const NodeId r = cell_nodes[R];
dV += dot(Cjr(j, R), m_ur[r]);
}
dVj[j] = dV;
});
CellValue<double> Gj{mesh.connectivity()};
if constexpr (Dimension == 2) {
for (CellId j = 0; j < mesh.numberOfCells(); ++j) {
const auto& cell_nodes = cell_to_node_matrix[j];
double G = 0;
for (size_t R = 0; R < cell_nodes.size(); ++R) {
const NodeId r = cell_nodes[R];
const NodeId rp1 = cell_nodes[(R + 1) % cell_nodes.size()];
G -= dot(m_ur[rp1], Rd{m_ur[r][1], -m_ur[r][0]});
}
Gj[j] = 0.5 * G;
}
}
CellValue<double> dSj{mesh.connectivity()};
parallel_for(
mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
const auto& cell_nodes = cell_to_node_matrix[j];
double dS = 0;
for (size_t R = 0; R < cell_nodes.size(); ++R) {
const NodeId r = cell_nodes[R];
const Rd du = u[j] - m_ur[r];
dS += dot(m_Ajr(j, R) * du, du);
}
dSj[j] = dS;
});
CellValue<Rd> dPj{mesh.connectivity()};
parallel_for(
mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
const auto& cell_nodes = cell_to_node_matrix[j];
Rd dP = zero;
for (size_t R = 0; R < cell_nodes.size(); ++R) {
const NodeId r = cell_nodes[R];
Rd du = u[j] - m_ur[r];
dP += m_Ajr(j, R) * du;
}
dPj[j] = dP;
});
#warning fixed gamma value
const double gamma = 1.4;
for (CellId j = 0; j < mesh.numberOfCells(); ++j) {
const double inv_mj = 1 / (rho[j] * Vj[j]);
const double inv_Vj = 1 / Vj[j];
Polynomial1D DT({0, 1});
std::vector<double> delta_S = {dSj[j], -0.5 * inv_mj * dot(dPj[j], dPj[j])};
if constexpr (Dimension == 2) {
delta_S[1] += p[j] * Gj[j];
}
Polynomial1D p_S0(delta_S);
std::vector<double> delta_V{dVj[j] /*, dG */};
if constexpr (Dimension == 2) {
delta_V.push_back(Gj[j]);
}
Polynomial1D p_deltaV(delta_V);
Polynomial1D p_S1({dSj[j] - p[j] * dVj[j], -0.5 * inv_mj * dot(dPj[j], dPj[j])});
Polynomial1D p_S = //
p_S0 //
+ (inv_Vj * DT) * (((gamma - 1) * p_S1) * p_deltaV + (gamma - 2) * p[j] * p_deltaV * p_deltaV) //
+ (inv_Vj * DT) * (inv_Vj * DT) * (((gamma - 1) * (gamma - 2)) * p_S1) * p_deltaV * p_deltaV;
Polynomial1DRootComputer computer{p_S, 0, dt};
std::optional dt_S = computer.getFirstRoot();
if (dt_S.has_value()) {
Assert(dt_S.value() <= dt);
dt = dt_S.value();
}
}
if (dt != dt_max) {
std::cout << "entropy variation imposes time step\n";
}
return parallel::allReduceMin(dt);
}
AcousticSolver(SolverType solver_type,
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_solver_type{solver_type}, 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->_computeAjr(mesh, rhoc);
m_Ajr = Ajr;
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 += dot(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));
}
std::tuple<double,
std::shared_ptr<const IMesh>,
std::shared_ptr<const DiscreteScalarFunction>,
std::shared_ptr<const DiscreteVectorFunction>,
std::shared_ptr<const DiscreteScalarFunction>>
apply_computing_entropic_dt(const double& dt_max,
const std::shared_ptr<const MeshType>& mesh,
const DiscreteScalarFunction& rho,
const DiscreteVectorFunction& u,
const DiscreteScalarFunction& E,
const DiscreteScalarFunction& p) const
{
const auto& cell_to_node_matrix = mesh->connectivity().cellToNodeMatrix();
double dt = dt_max;
dt = _computeMaxDensityDt(*mesh, rho, dt);
dt = _computeMaxEnergyDt(*mesh, rho, u, E, p, dt);
dt = _computeMaxEntropyDt(*mesh, rho, u, p, dt);
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 += dot(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 {dt, 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<double,
std::shared_ptr<const IMesh>,
std::shared_ptr<const IDiscreteFunction>,
std::shared_ptr<const IDiscreteFunction>,
std::shared_ptr<const IDiscreteFunction>>
apply_computing_entropic_dt(const double& dt_max,
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>& p) 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_computing_entropic_dt(dt_max, 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),
*std::dynamic_pointer_cast<const DiscreteScalarFunction>(p));
}
AcousticSolver(SolverType solver_type,
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{solver_type,
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(
SolverType solver_type,
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>>(solver_type, i_mesh, rho, c, u, p, bc_descriptor_list);
break;
}
case 2: {
m_acoustic_solver = std::make_unique<AcousticSolver<2>>(solver_type, i_mesh, rho, c, u, p, bc_descriptor_list);
break;
}
case 3: {
m_acoustic_solver = std::make_unique<AcousticSolver<3>>(solver_type, i_mesh, rho, c, u, p, bc_descriptor_list);
break;
}
default: {
throw UnexpectedError("invalid mesh dimension");
}
}
}