diff --git a/src/analysis/Polynomial1D.hpp b/src/analysis/Polynomial1D.hpp
index df62a381b31dd8449df7ffb34c2eaa53175a014a..db033a4e83fb4b297909421ecb10613921501a2a 100644
--- a/src/analysis/Polynomial1D.hpp
+++ b/src/analysis/Polynomial1D.hpp
@@ -4,6 +4,8 @@
#include <utils/Array.hpp>
#include <utils/Exceptions.hpp>
+#include <algorithm>
+
class [[nodiscard]] Polynomial1D
{
private:
@@ -22,17 +24,18 @@ class [[nodiscard]] Polynomial1D
PUGS_INLINE
friend Polynomial1D _simplify(Polynomial1D && p)
{
- size_t real_degree = p._getRealDegree();
-
- if (real_degree != p.degree()) {
- Polynomial1D q(real_degree);
- for (size_t i = 0; i <= real_degree; ++i) {
- q.coefficient(i) = p.coefficient(i);
- }
- return q;
- } else {
- return std::move(p);
- }
+ return std::move(p);
+ // size_t real_degree = p._getRealDegree();
+
+ // if (real_degree != p.degree()) {
+ // Polynomial1D q(real_degree);
+ // for (size_t i = 0; i <= real_degree; ++i) {
+ // q.coefficient(i) = p.coefficient(i);
+ // }
+ // return q;
+ // } else {
+ // return std::move(p);
+ // }
}
public:
@@ -41,27 +44,27 @@ class [[nodiscard]] Polynomial1D
bool has_written = false;
for (size_t i = 0; i < p.m_coefficients.size(); ++i) {
const double& coef = p.m_coefficients[i];
- if (coef != 0) {
- if (coef > 0) {
- if (has_written) {
- os << " + ";
- }
- os << coef;
+ // if (coef != 0) {
+ if (coef > 0) {
+ if (has_written) {
+ os << " + ";
+ }
+ os << coef;
+ } else {
+ if (has_written) {
+ os << " - " << -coef;
} else {
- if (has_written) {
- os << " - " << -coef;
- } else {
- os << coef;
- }
+ os << coef;
}
- if (i > 0) {
- os << "*x";
- if (i > 1) {
- os << '^' << i;
- }
+ }
+ if (i > 0) {
+ os << "*x";
+ if (i > 1) {
+ os << '^' << i;
}
- has_written = true;
}
+ has_written = true;
+ // }
}
if (not has_written) {
os << 0;
@@ -169,9 +172,15 @@ class [[nodiscard]] Polynomial1D
}
}
- for (size_t j = q_degree; j <= remaining.degree(); ++j) {
- remaining.m_coefficients[j] = 0;
+ if (q_degree == 0) {
+ remaining.m_coefficients.resize(1);
+ remaining.m_coefficients[0] = 0;
+ } else {
+ remaining.m_coefficients.resize(std::max(q_degree - 1, 1ul));
}
+ // for (size_t j = q_degree; j <= remaining.degree(); ++j) {
+ // remaining.m_coefficients[j] = 0;
+ // }
return _simplify(std::move(remaining));
}
diff --git a/src/analysis/Polynomial1DRootComputer.cpp b/src/analysis/Polynomial1DRootComputer.cpp
index c86d2ad3117a734979d41b9161cef3d7e8a4177e..a2c54f88c7bc2c7e87a2be7c096fe3f12c897a6e 100644
--- a/src/analysis/Polynomial1DRootComputer.cpp
+++ b/src/analysis/Polynomial1DRootComputer.cpp
@@ -13,7 +13,7 @@ Polynomial1DRootComputer::_computeRootBetween(const PositionSignChanges& left_bo
double left_value = p0(left);
int cpt = 0;
- while (right - left > 1E-12 * interval_size) {
+ while (right - left > 1E-15 * interval_size) {
double middle = 0.5 * (right + left);
double middle_value = p0(middle);
if (middle_value * left_value < 0) {
diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index ef4df609f50d032d360242e2d1020d5bbe2ed518..1863043263f8f584bdd16a138bd39821ab477ac4 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -311,30 +311,17 @@ SchemeModule::SchemeModule()
));
- this->_addBuiltinFunction("legacy_entropy_dt",
+ this->_addBuiltinFunction("acoustic_pg_entropy_dt",
std::function(
[](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<const SubItemValuePerItemVariant>& Ajr,
- const std::shared_ptr<const ItemValueVariant>& ur, const double& dt_max) -> double {
- return acoustic_entropy_dt(rho, u, E, p, Ajr, ur, dt_max);
- }
-
- ));
-
- this->_addBuiltinFunction("acoustic_entropy_dt",
- std::function(
-
- [](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>& epsilon,
const std::shared_ptr<const IDiscreteFunction>& p,
+ const std::shared_ptr<const IDiscreteFunction>& gamma,
const std::shared_ptr<const SubItemValuePerItemVariant>& Ajr,
const std::shared_ptr<const ItemValueVariant>& ur, const double& dt_max) -> double {
- return acoustic_entropy_dt(rho, u, E, p, Ajr, ur, dt_max);
+ return acoustic_pg_entropy_dt(rho, u, epsilon, p, gamma, Ajr, ur, dt_max);
}
));
diff --git a/src/scheme/AcousticSolver.cpp b/src/scheme/AcousticSolver.cpp
index 09698e685833781c802b777396ea6a80f72f6145..f63eef8aef7b3b99c31ba280b5f9cb9efd120cf4 100644
--- a/src/scheme/AcousticSolver.cpp
+++ b/src/scheme/AcousticSolver.cpp
@@ -122,11 +122,11 @@ acoustic_rho_dt(const DiscreteFunctionP0<Dimension, double>& rho,
double dt = dt_max;
for (CellId j = 0; j < mesh.numberOfCells(); ++j) {
- Polynomial1D p_tau({1. / rho[j], dVj[j] / (rho[j] * Vj[j])});
+ Polynomial1D q_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])});
+ q_tau = q_tau + Polynomial1D({0, 0, Gj[j] / (rho[j] * Vj[j])});
}
- Polynomial1DRootComputer computer{p_tau, 0, dt};
+ Polynomial1DRootComputer computer{q_tau, 0, dt};
std::optional dt_tau = computer.getFirstRoot();
if (dt_tau.has_value()) {
Assert(dt_tau.value() <= dt);
@@ -134,15 +134,11 @@ acoustic_rho_dt(const DiscreteFunctionP0<Dimension, double>& rho,
}
}
- // if (dt != dt_max) {
- // std::cout << "volume variation imposes time step\n";
- // }
-
- // if (dt < dt_max) {
- // dt *= 0.95 * dt;
- // }
+ if (dt < dt_max) {
+ dt *= 0.95;
+ }
- return dt;
+ return parallel::allReduceMin(dt);
}
template <size_t Dimension>
@@ -237,138 +233,9 @@ acoustic_epsilon_dt(const DiscreteFunctionP0<Dimension, double>& rho,
}
}
- return dt;
-}
-
-template <size_t Dimension>
-double
-legacy_entropy_dt(const DiscreteFunctionP0<Dimension, double>& rho,
- const DiscreteFunctionP0<Dimension, TinyVector<Dimension>>& u,
- const DiscreteFunctionP0<Dimension, double>& E,
- const DiscreteFunctionP0<Dimension, double>& p,
- const std::shared_ptr<const SubItemValuePerItemVariant>& Ajr_variant,
- const std::shared_ptr<const ItemValueVariant>& ur_variant,
- const double& dt_max)
-{
- double dt = dt_max;
-
- using Rd = TinyVector<Dimension>;
- using Rdxd = TinyMatrix<Dimension>;
- using MeshType = Mesh<Connectivity<Dimension>>;
- using MeshDataType = MeshData<Dimension>;
-
- const MeshType& mesh = dynamic_cast<const MeshType&>(*rho.mesh());
-
- const NodeValue<const Rd> ur = ur_variant->get<NodeValue<const Rd>>();
- const NodeValuePerCell<const Rdxd> Ajr = Ajr_variant->get<NodeValuePerCell<const Rdxd>>();
- 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), 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(ur[rp1], Rd{ur[r][1], -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] - ur[r];
- dS += dot(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] - ur[r];
- dP += 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]};
- 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 << "energy variation imposes time step\n";
+ dt *= 0.95;
}
return parallel::allReduceMin(dt);
@@ -376,13 +243,14 @@ legacy_entropy_dt(const DiscreteFunctionP0<Dimension, double>& rho,
template <size_t Dimension>
double
-acoustic_entropy_dt(const DiscreteFunctionP0<Dimension, double>& rho,
- const DiscreteFunctionP0<Dimension, TinyVector<Dimension>>& u,
- const DiscreteFunctionP0<Dimension, double>& E,
- const DiscreteFunctionP0<Dimension, double>& p,
- const std::shared_ptr<const SubItemValuePerItemVariant>& Ajr_variant,
- const std::shared_ptr<const ItemValueVariant>& ur_variant,
- const double& dt_max)
+acoustic_pg_entropy_dt(const DiscreteFunctionP0<Dimension, double>& rho,
+ const DiscreteFunctionP0<Dimension, TinyVector<Dimension>>& u,
+ const DiscreteFunctionP0<Dimension, double>& epsilon,
+ const DiscreteFunctionP0<Dimension, double>& p,
+ const DiscreteFunctionP0<Dimension, double>& gamma,
+ const std::shared_ptr<const SubItemValuePerItemVariant>& Ajr_variant,
+ const std::shared_ptr<const ItemValueVariant>& ur_variant,
+ const double& dt_max)
{
double dt = dt_max;
@@ -419,7 +287,7 @@ acoustic_entropy_dt(const DiscreteFunctionP0<Dimension, double>& rho,
});
CellValue<double> Gj{mesh.connectivity()};
-
+ Gj.fill(0);
if constexpr (Dimension == 2) {
for (CellId j = 0; j < mesh.numberOfCells(); ++j) {
const auto& cell_nodes = cell_to_node_matrix[j];
@@ -470,32 +338,51 @@ acoustic_entropy_dt(const DiscreteFunctionP0<Dimension, double>& rho,
dPj[j] = dP;
});
-#warning fixed gamma value
- const double gamma = 1.4;
+ DiscreteFunctionP0<Dimension, const double> tauj_gamma_1 = pow(rho, 1 - gamma);
+ DiscreteFunctionP0<Dimension, const double> tauj_gamma_2 = tauj_gamma_1 * rho;
+ std::vector<Polynomial1D> q_Sj(mesh.numberOfCells(), Polynomial1D(std::vector<double>{0}));
+
+ Polynomial1D DT({0, 1});
for (CellId j = 0; j < mesh.numberOfCells(); ++j) {
- const double inv_mj = 1 / (rho[j] * Vj[j]);
- const double inv_Vj = 1 / Vj[j];
+ if (dSj[j] == 0)
+ continue;
+ const double mj = (rho[j] * Vj[j]);
+ const double inv_mj = 1 / mj;
+
+ const double delta_S = dSj[j] * tauj_gamma_1[j];
+ const double mj_pi_Gj = [&]() -> double {
+ if constexpr (Dimension == 1) {
+ return 0;
+ } else if constexpr (Dimension == 2) {
+ return mj * p[j] * Gj[j];
+ } else {
+ throw NotImplementedError("3d");
+ }
+ }();
+
+ const double dPj2 = dot(dPj[j], dPj[j]);
+
+ Polynomial1D dS1({dSj[j] - p[j] * dVj[j], -0.5 * inv_mj * dPj2});
- Polynomial1D DT({0, 1});
- std::vector<double> delta_S = {dSj[j], -0.5 * inv_mj * dot(dPj[j], dPj[j])};
+ Polynomial1D q_tau({1. / rho[j], inv_mj * dVj[j]});
if constexpr (Dimension == 2) {
- delta_S[1] += p[j] * Gj[j];
+ q_tau = q_tau + Polynomial1D({0, 0, inv_mj * Gj[j]});
}
- Polynomial1D p_S0(delta_S);
- std::vector<double> delta_V{dVj[j]};
- if constexpr (Dimension == 2) {
- delta_V.push_back(Gj[j]);
+ Polynomial1D S1 = Polynomial1D(std::vector<double>{(mj_pi_Gj - 0.5 * dPj2) * tauj_gamma_1[j]}) +
+ ((gamma[j] - 1) * (dVj[j] /* +DT*G[j] */) * tauj_gamma_2[j]) * dS1;
+
+ if (gamma[j] < 2) {
+ const double tauj_min = std::min(q_tau(0), q_tau(dt_max));
+ S1 = S1 + Polynomial1D(std::vector<double>{0.5 * (gamma[j] - 1) * (gamma[j] - 2) * (dVj[j] /* +DT*G[j] */) *
+ (dVj[j] /* +DT*G[j] */) * std::pow(tauj_min, gamma[j] - 3)}) *
+ (Polynomial1D(std::vector<double>{epsilon[j]}) + (inv_mj * DT) * dS1);
}
- 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 = Polynomial1D(std::vector<double>{delta_S}) + (inv_mj * DT) * S1;
- 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;
+ q_Sj[j] = p_S;
Polynomial1DRootComputer computer{p_S, 0, dt};
std::optional dt_S = computer.getFirstRoot();
@@ -582,90 +469,48 @@ acoustic_epsilon_dt(const std::shared_ptr<const IDiscreteFunction>& rho,
}
double
-legacy_entropy_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<const IDiscreteFunction>& p,
- const std::shared_ptr<const SubItemValuePerItemVariant>& Ajr_variant,
- const std::shared_ptr<const ItemValueVariant>& ur_variant,
- const double& dt_max)
+acoustic_pg_entropy_dt(const std::shared_ptr<const IDiscreteFunction>& rho,
+ const std::shared_ptr<const IDiscreteFunction>& u,
+ const std::shared_ptr<const IDiscreteFunction>& epsilon,
+ const std::shared_ptr<const IDiscreteFunction>& p,
+ const std::shared_ptr<const IDiscreteFunction>& gamma,
+ const std::shared_ptr<const SubItemValuePerItemVariant>& Ajr_variant,
+ const std::shared_ptr<const ItemValueVariant>& ur_variant,
+ const double& dt_max)
{
- if (not checkDiscretizationType({rho, u, E, p}, DiscreteFunctionType::P0)) {
+ if (not checkDiscretizationType({rho, u, gamma, p}, DiscreteFunctionType::P0)) {
throw NormalError("acoustic solver expects P0 functions");
}
- std::shared_ptr mesh = getCommonMesh({rho, u, E, p});
+ std::shared_ptr mesh = getCommonMesh({rho, u, gamma, p});
switch (mesh->dimension()) {
case 1: {
constexpr size_t Dimension = 1;
- return legacy_entropy_dt(dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*rho),
- dynamic_cast<const DiscreteFunctionP0<Dimension, TinyVector<Dimension>>&>(*u),
- dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*E),
- dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p), Ajr_variant, ur_variant,
- dt_max);
+ return acoustic_pg_entropy_dt(dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*rho),
+ dynamic_cast<const DiscreteFunctionP0<Dimension, TinyVector<Dimension>>&>(*u),
+ dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*epsilon),
+ dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p),
+ dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*gamma), Ajr_variant,
+ ur_variant, dt_max);
}
case 2: {
constexpr size_t Dimension = 2;
- return legacy_entropy_dt(dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*rho),
- dynamic_cast<const DiscreteFunctionP0<Dimension, TinyVector<Dimension>>&>(*u),
- dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*E),
- dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p), Ajr_variant, ur_variant,
- dt_max);
+ return acoustic_pg_entropy_dt(dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*rho),
+ dynamic_cast<const DiscreteFunctionP0<Dimension, TinyVector<Dimension>>&>(*u),
+ dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*epsilon),
+ dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p),
+ dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*gamma), Ajr_variant,
+ ur_variant, dt_max);
}
case 3: {
constexpr size_t Dimension = 3;
- return legacy_entropy_dt(dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*rho),
- dynamic_cast<const DiscreteFunctionP0<Dimension, TinyVector<Dimension>>&>(*u),
- dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*E),
- dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p), Ajr_variant, ur_variant,
- dt_max);
- }
- default: {
- throw UnexpectedError("invalid mesh dimension");
- }
- }
-}
-
-double
-acoustic_entropy_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<const IDiscreteFunction>& p,
- const std::shared_ptr<const SubItemValuePerItemVariant>& Ajr_variant,
- const std::shared_ptr<const ItemValueVariant>& ur_variant,
- const double& dt_max)
-{
- if (not checkDiscretizationType({rho, u, E, p}, DiscreteFunctionType::P0)) {
- throw NormalError("acoustic solver expects P0 functions");
- }
-
- std::shared_ptr mesh = getCommonMesh({rho, u, E, p});
-
- switch (mesh->dimension()) {
- case 1: {
- constexpr size_t Dimension = 1;
- return acoustic_entropy_dt(dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*rho),
- dynamic_cast<const DiscreteFunctionP0<Dimension, TinyVector<Dimension>>&>(*u),
- dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*E),
- dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p), Ajr_variant, ur_variant,
- dt_max);
- }
- case 2: {
- constexpr size_t Dimension = 2;
- return acoustic_entropy_dt(dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*rho),
- dynamic_cast<const DiscreteFunctionP0<Dimension, TinyVector<Dimension>>&>(*u),
- dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*E),
- dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p), Ajr_variant, ur_variant,
- dt_max);
- }
- case 3: {
- constexpr size_t Dimension = 3;
- return acoustic_entropy_dt(dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*rho),
- dynamic_cast<const DiscreteFunctionP0<Dimension, TinyVector<Dimension>>&>(*u),
- dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*E),
- dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p), Ajr_variant, ur_variant,
- dt_max);
+ return acoustic_pg_entropy_dt(dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*rho),
+ dynamic_cast<const DiscreteFunctionP0<Dimension, TinyVector<Dimension>>&>(*u),
+ dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*epsilon),
+ dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p),
+ dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*gamma), Ajr_variant,
+ ur_variant, dt_max);
}
default: {
throw UnexpectedError("invalid mesh dimension");
@@ -1083,9 +928,23 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler
});
CellValue<const double> new_Vj = MeshDataManager::instance().getMeshData(*new_mesh).Vj();
+ NodeValuePerCell<const Rd> Cjr = MeshDataManager::instance().getMeshData(mesh).Cjr();
parallel_for(
- mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { new_rho[j] *= Vj[j] / new_Vj[j]; });
+ mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
+ const auto& cell_nodes = cell_to_node_matrix[j];
+
+ const double mj = rho[j] * Vj[j];
+
+ new_rho[j] = mj / new_Vj[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), ur[r]);
+ }
+
+ new_rho[j] = 1. / (1. / rho[j] + dt / mj * dV);
+ });
return {new_mesh, std::make_shared<DiscreteScalarFunction>(new_mesh, new_rho),
std::make_shared<DiscreteVectorFunction>(new_mesh, new_u),
diff --git a/src/scheme/AcousticSolver.hpp b/src/scheme/AcousticSolver.hpp
index a00ce4e143f1d66fe8d851b9b2c2cc5a8422d39e..1fc8e03557542d5e9861db7c7cf511f188ded586 100644
--- a/src/scheme/AcousticSolver.hpp
+++ b/src/scheme/AcousticSolver.hpp
@@ -27,27 +27,21 @@ double acoustic_epsilon_dt(const std::shared_ptr<const IDiscreteFunction>& rho,
double acoustic_entropy_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<const IDiscreteFunction>& epsilon,
const std::shared_ptr<const IDiscreteFunction>& p,
+ const std::shared_ptr<const IDiscreteFunction>& gamma,
const std::shared_ptr<const SubItemValuePerItemVariant>& Ajr_variant,
const std::shared_ptr<const ItemValueVariant>& ur_variant,
const double& dt_max);
-double legacy_entropy_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<const IDiscreteFunction>& p,
- const std::shared_ptr<const SubItemValuePerItemVariant>& Ajr_variant,
- const std::shared_ptr<const ItemValueVariant>& ur_variant,
- const double& dt_max);
-
-double acoustic_entropy_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<const IDiscreteFunction>& p,
- const std::shared_ptr<const SubItemValuePerItemVariant>& Ajr_variant,
- const std::shared_ptr<const ItemValueVariant>& ur_variant,
- const double& dt_max);
+double acoustic_pg_entropy_dt(const std::shared_ptr<const IDiscreteFunction>& rho,
+ const std::shared_ptr<const IDiscreteFunction>& u,
+ const std::shared_ptr<const IDiscreteFunction>& epsilon,
+ const std::shared_ptr<const IDiscreteFunction>& p,
+ const std::shared_ptr<const IDiscreteFunction>& gamma,
+ const std::shared_ptr<const SubItemValuePerItemVariant>& Ajr_variant,
+ const std::shared_ptr<const ItemValueVariant>& ur_variant,
+ const double& dt_max);
class AcousticSolverHandler
{