From f10638eb3c6e5efdb5b4f95871266d7b1aa64faa Mon Sep 17 00:00:00 2001
From: acrovato <a.crovato@uliege.be>
Date: Tue, 11 Jan 2022 16:16:44 +0100
Subject: [PATCH] Change F0El. Functions can be clamped or not. Clamping is
performed on density base, and applied to density, Mach and pressure.
---
dart/src/dart.h | 10 +-
dart/src/wF0El.cpp | 258 ++++++++++++++++++++++++++++++--------------
dart/src/wF0El.h | 149 ++++++++++++++++++-------
dart/src/wFluid.cpp | 12 +--
4 files changed, 300 insertions(+), 129 deletions(-)
diff --git a/dart/src/dart.h b/dart/src/dart.h
index 2afc459..ad5628c 100644
--- a/dart/src/dart.h
+++ b/dart/src/dart.h
@@ -69,12 +69,16 @@ class Adjoint;
class F0Ps;
class F0El;
class F0ElC;
+class F0ElRhoBase;
class F0ElRho;
-class F0ElRhoL;
+class F0ElRhoClamp;
+class F0ElRhoLin;
class F0ElMach;
-class F0ElMachL;
+class F0ElMachClamp;
+class F0ElMachLin;
class F0ElCp;
-class F0ElCpL;
+class F0ElCpClamp;
+class F0ElCpLin;
class F1El;
class F1ElVi;
class F1Ct;
diff --git a/dart/src/wF0El.cpp b/dart/src/wF0El.cpp
index 8a2deef..e17c9d1 100644
--- a/dart/src/wF0El.cpp
+++ b/dart/src/wF0El.cpp
@@ -28,153 +28,249 @@ void F0ElC::update(double _v)
/**
* @brief Evaluate the constant
*/
-double F0ElC::eval(Element const &e, std::vector<double> const &u, size_t k) const
+double F0ElC::eval(Element const &e, std::vector<double> const &phi, size_t k) const
{
return v;
}
/**
* @brief Evaluate the gradient of the constant
*/
-double F0ElC::evalGrad(Element const &e, std::vector<double> const &u, size_t k) const
+double F0ElC::evalGrad(Element const &e, std::vector<double> const &phi, size_t k) const
{
return 0;
}
/**
- * @brief Evaluate the nonlinear density (constant over element)
+ * @brief Precompute values for Padé approximation
+ *
+ * The input _mC2 is the square of the critical Mach number (used to compute
+ * the critical velocity).
*/
-double F0ElRho::eval(Element const &e, std::vector<double> const &u, size_t k) const
+F0ElRhoBase::F0ElRhoBase(double _mInf, double _mC2) : F0El(), gamma(1.4), mInf(_mInf)
{
- double gradU = e.computeGradient(u, k).norm(); // norm of velocity
-
- if (gradU < gradUC) // true density
- {
- // particularized: pow(1 + 0.5 * (gamma - 1) * mInf * mInf * (1 - gradU2), 1 / (gamma - 1));
- double a = 1 + 0.5 * (gamma - 1) * mInf * mInf * (1 - gradU * gradU);
- return sqrt(a * a * a * a * a);
- }
+ uC = sqrt(_mC2 * (1 / (mInf * mInf) + (gamma - 1) / 2) / (1 + (gamma - 1) / 2 * _mC2)); // critical velocity
+ rC = 1 + 0.5 * (gamma - 1) * mInf * mInf * (1 - uC * uC); // r(uC)
+ beta = ((gamma - 1) * mInf * mInf * uC) / rC; // dr/du(uC) / -r(uC)
+}
+/**
+ * @brief Evaluate the base of the isentropic density
+ *
+ * r = [1 + (gamma-1)/2*Mi^2*(1-u^2)], u <= uC
+ * r = r(uC) / (1-(dr/du(uC)/f(uC))*(u/uC-1)), u > uC
+ */
+double F0ElRhoBase::eval(Element const &e, std::vector<double> const &phi, size_t k) const
+{
+ double u = e.computeGradient(phi, k).norm(); // norm of velocity
+ if (u <= uC)
+ return 1 + 0.5 * (gamma - 1) * mInf * mInf * (1 - u * u);
else
- {
- double a = 1 + 0.5 * (gamma - 1) * mInf * mInf * (1 - gradUC * gradUC);
- double beta = mInf * mInf * gradUC / a;
- return sqrt(a * a * a * a * a) / (1 + beta * (gradU / gradUC - 1)); // limited density
- }
+ return rC / (1 + beta * (u / uC - 1));
}
/**
- * @brief Evaluate the nonlinear density derivative (constant over element)
+ * @brief Evaluate the gradient of the base of the isentropic density
+ *
+ * dr = -(gamma-1)*Mi^2 * u, u <= uC
+ * dr = -r(uC) / (1-(dr/du(uC)/f(uC))*(u/uC-1))^2 * (dr/du(uC)/f(uC)) / uC, u > uC
+ * Computed value has to be multiplied by u to recover actual gradient
*/
-double F0ElRho::evalGrad(Element const &e, std::vector<double> const &u, size_t k) const
+double F0ElRhoBase::evalGrad(Element const &e, std::vector<double> const &phi, size_t k) const
{
- double gradU = e.computeGradient(u, k).norm(); // norm of velocity
+ double u = e.computeGradient(phi, k).norm(); // norm of velocity
+ if (u <= uC)
+ return -(gamma - 1) * mInf * mInf;
+ else
+ return -rC / ((1 + beta * (u / uC - 1)) * (1 + beta * (u / uC - 1))) * beta / uC / u;
+}
- if (gradU < gradUC) // true density
- {
- //particularized: -mInf * mInf * pow(rho, 2 - gamma);
- double a = 1 + 0.5 * (gamma - 1) * mInf * mInf * (1 - gradU * gradU);
- return -mInf * mInf * sqrt(a * a * a); // has to be multiplied by grad u to recover true derivative
- }
- else // limited density
- {
- double a = 1 + 0.5 * (gamma - 1) * mInf * mInf * (1 - gradUC * gradUC);
- double beta = mInf * mInf * gradUC / a;
- return -sqrt(a * a * a * a * a) / ((1 + beta * (gradU / gradUC - 1)) * (1 + beta * (gradU / gradUC - 1))) * beta / gradUC / gradU; // has to be multiplied by grad u to recover true derivative
- }
+/**
+ * @brief Evaluate the isentropic density
+ *
+ * rho = (1 + (gamma-1)/2*Mi^2*(1-u^2))^(1/(gamma-1))
+ */
+double F0ElRho::eval(Element const &e, std::vector<double> const &phi, size_t k) const
+{
+ double u = e.computeGradient(phi, k).norm(); // norm of velocity
+ double r = 1 + 0.5 * (gamma - 1) * mInf * mInf * (1 - u * u);
+ return sqrt(r * r * r * r * r); // particularized to gamma=1.4
+}
+/**
+ * @brief Evaluate the isentropic density gradient
+ *
+ * drho = -Mi^2 * rho^(2-gamma) * u
+ * Computed value has to be multiplied by u to recover actual gradient
+ */
+double F0ElRho::evalGrad(Element const &e, std::vector<double> const &phi, size_t k) const
+{
+ double u = e.computeGradient(phi, k).norm(); // norm of velocity
+ double r = 1 + 0.5 * (gamma - 1) * mInf * mInf * (1 - u * u);
+ return -mInf * mInf * sqrt(r * r * r); // particularized to gamma=1.4
}
/**
- * @brief Evaluate the linear density (constant over element)
+ * @brief Evaluate the (clamped) isentropic density
+ *
+ * rho = r^(1/(gamma-1))
*/
-double F0ElRhoL::eval(Element const &e, std::vector<double> const &u, size_t k) const
+double F0ElRhoClamp::eval(Element const &e, std::vector<double> const &phi, size_t k) const
+{
+ double r = F0ElRhoBase::eval(e, phi, k);
+ return sqrt(r * r * r * r * r); // particularized to gamma=1.4
+}
+/**
+ * @brief Evaluate the (clamped) isentropic density gradient
+ *
+ * drho = 1/(gamma-1) * r^(2-gamma) * dr
+ * Computed value has to be multiplied by u to recover actual gradient
+ */
+double F0ElRhoClamp::evalGrad(Element const &e, std::vector<double> const &phi, size_t k) const
+{
+ double r = F0ElRhoBase::eval(e, phi, k);
+ double dr = F0ElRhoBase::evalGrad(e, phi, k);
+ return 1 / (gamma - 1) * sqrt(r * r * r) * dr; // particularized to gamma=1.4
+}
+
+/**
+ * @brief Evaluate the linear density
+ */
+double F0ElRhoLin::eval(Element const &e, std::vector<double> const &phi, size_t k) const
{
return 1.;
}
/**
- * @brief Evaluate the linear density derivative (constant over element)
+ * @brief Evaluate the linear density gradient
*/
-double F0ElRhoL::evalGrad(Element const &e, std::vector<double> const &u, size_t k) const
+double F0ElRhoLin::evalGrad(Element const &e, std::vector<double> const &phi, size_t k) const
{
return 0.;
}
/**
- * @brief Evaluate the nonlinear Mach number (constant over element)
+ * @brief Evaluate the isentropic Mach number
+ *
+ * M = u / sqrt(1/mi^2 + (gamma-1)/2*(1-u^2))
+ */
+double F0ElMach::eval(Element const &e, std::vector<double> const &phi, size_t k) const
+{
+ double u2 = e.computeGradient(phi, k).squaredNorm(); // velocity squared
+ double a2 = 1 / (mInf * mInf) + 0.5 * (gamma - 1) * (1 - u2); // speed of sound squared
+ return sqrt(u2) / sqrt(a2);
+}
+/**
+ * @brief Evaluate the isentropic Mach number gradient
+ *
+ * dM = (1/(u*a) + (gamma-1)/2*u/a^3) * u
+ * Computed value has to be multiplied by u to recover actual gradient
*/
-double F0ElMach::eval(Element const &e, std::vector<double> const &u, size_t k) const
+double F0ElMach::evalGrad(Element const &e, std::vector<double> const &phi, size_t k) const
{
- Eigen::VectorXd gradU = e.computeGradient(u, k); // velocity
- double gradU2 = gradU.squaredNorm(); // velocity squared
- double a2 = 1 / (mInf * mInf) + 0.5 * (gamma - 1) * (1 - gradU2); // speed of sound squared
+ double u2 = e.computeGradient(phi, k).squaredNorm(); // velocity squared
+ double a2 = 1 / (mInf * mInf) + 0.5 * (gamma - 1) * (1 - u2); // speed of sound squared
+ return 1 / sqrt(u2 * a2) + 0.5 * (gamma - 1) * sqrt(u2) / sqrt(a2 * a2 * a2);
+}
- return sqrt(gradU2) / sqrt(a2);
+/**
+ * @brief Evaluate the (clamped) isentropic Mach number
+ *
+ * M = u / sqrt(r/Mi^2)
+ */
+double F0ElMachClamp::eval(Element const &e, std::vector<double> const &phi, size_t k) const
+{
+ double u = e.computeGradient(phi, k).norm(); // norm of velocity
+ double r = F0ElRhoBase::eval(e, phi, k);
+ return mInf * u / sqrt(r);
}
/**
- * @brief Evaluate the nonlinear Mach number derivative (constant over element)
+ * @brief Evaluate the (clamped) isentropic Mach number gradient
+ *
+ * dM = 1 / sqrt(r/Mi^2) + u / sqrt(r^3/Mi^2) * dr
+ * Computed value has to be multiplied by u to recover actual gradient
*/
-double F0ElMach::evalGrad(Element const &e, std::vector<double> const &u, size_t k) const
+double F0ElMachClamp::evalGrad(Element const &e, std::vector<double> const &phi, size_t k) const
{
- Eigen::VectorXd gradU = e.computeGradient(u, k); // velocity
- double gradU2 = gradU.squaredNorm(); // velocity squared
- double a2 = 1 / (mInf * mInf) + 0.5 * (gamma - 1) * (1 - gradU2);
-
- return 1 / sqrt(gradU2 * a2) + 0.5 * (gamma - 1) * sqrt(gradU2) / sqrt(a2 * a2 * a2); // has to be multiplied by grad u to recover true derivative
+ double u = e.computeGradient(phi, k).norm(); // norm of velocity
+ double r = F0ElRhoBase::eval(e, phi, k);
+ double dr = F0ElRhoBase::evalGrad(e, phi, k);
+ return mInf * (1 / sqrt(r) / u - 0.5 * u / sqrt(r * r * r) * dr);
}
/**
- * @brief Evaluate the linear Mach number (constant over element)
+ * @brief Evaluate the linear Mach number
*/
-double F0ElMachL::eval(Element const &e, std::vector<double> const &u, size_t k) const
+double F0ElMachLin::eval(Element const &e, std::vector<double> const &phi, size_t k) const
{
return 0.;
}
/**
- * @brief Evaluate the linear Mach number derivative (constant over element)
+ * @brief Evaluate the linear Mach number gradient
*/
-double F0ElMachL::evalGrad(Element const &e, std::vector<double> const &u, size_t k) const
+double F0ElMachLin::evalGrad(Element const &e, std::vector<double> const &phi, size_t k) const
{
return 0.;
}
/**
- * @brief Evaluate the nonlinear pressure coefficient (constant over element)
+ * @brief Evaluate the isentropic pressure coefficient
+ *
+ * Cp = 2/(gamma*Mi^2) * ((1 + (gamma-1)/2*Mi^2*(1-u^2))^(gamma/(gamma-1)) - 1)
*/
-double F0ElCp::eval(Element const &e, std::vector<double> const &u, size_t k) const
+double F0ElCp::eval(Element const &e, std::vector<double> const &phi, size_t k) const
{
- Eigen::VectorXd gradU = e.computeGradient(u, k); // velocity
- double gradU2 = gradU.squaredNorm(); // velocity squared
- //particularized: 2 / (gamma * mInf * mInf) * (pow(1 + 0.5 * (gamma - 1) * mInf * mInf * (1 - gradU2), gamma / (gamma - 1)) - 1);
- double a = 1 + 0.5 * (gamma - 1) * mInf * mInf * (1 - gradU2);
-
- if (a <= 0)
- return -2 / (gamma * mInf * mInf); // limit the pressure coefficient when vacuum conditions are reached
- else
- return 2 / (gamma * mInf * mInf) * (sqrt(a * a * a * a * a * a * a) - 1);
+ double u2 = e.computeGradient(phi, k).squaredNorm(); // velocity squared
+ double r = 1 + 0.5 * (gamma - 1) * mInf * mInf * (1 - u2);
+ return 2 / (gamma * mInf * mInf) * (sqrt(r * r * r * r * r * r * r) - 1); // particularized to gamma=1.4
}
/**
- * @brief Evaluate the nonlinear pressure coefficient derivative (constant over element)
+ * @brief Evaluate the isentropic pressure coefficient gradient
+ *
+ * dCp = -2*(1 + (gamma-1)/2*Mi^2*(1-u^2))^(1/(gamma-1)) * u
+ * Computed value has to multiplied by u to recover actual gradient
*/
-double F0ElCp::evalGrad(Element const &e, std::vector<double> const &u, size_t k) const
+double F0ElCp::evalGrad(Element const &e, std::vector<double> const &phi, size_t k) const
{
- Eigen::VectorXd gradU = e.computeGradient(u, k); // velocity
- double gradU2 = gradU.squaredNorm(); // velocity squared
- //particularized: -2 * pow(1 + 0.5 * (gamma - 1) * mInf * mInf * (1 - gradU2), 1 / (gamma - 1));
- double a = 1 + 0.5 * (gamma - 1) * mInf * mInf * (1 - gradU2);
+ double u2 = e.computeGradient(phi, k).squaredNorm(); // velocity squared
+ double r = 1 + 0.5 * (gamma - 1) * mInf * mInf * (1 - u2);
+ return -2 * sqrt(r * r * r * r * r); // particularized to gamma=1.4
+}
- if (a <= 0)
- return 0; // limit the gradient when vacuum conditions are reached
- else
- return -2 * sqrt(a * a * a * a * a); // has to be multiplied by grad u to recover true derivative
+/**
+ * @brief Evaluate the (clamped) isentropic pressure coefficient
+ *
+ * Cp = 2/(gamma*Mi^2) * (r^(gamma/(gamma-1)) - 1)
+ */
+double F0ElCpClamp::eval(Element const &e, std::vector<double> const &phi, size_t k) const
+{
+ double r = F0ElRhoBase::eval(e, phi, k);
+ return 2 / (gamma * mInf * mInf) * (sqrt(r * r * r * r * r * r * r) - 1); // particularized to gamma=1.4
+}
+/**
+ * @brief Evaluate the (clamped) isentropic pressure coefficient gradient
+ *
+ * dCp = 2/((gamma-1)*Mi^2) * r^(1/(gamma-1)) * dr
+ * Computed value has to multiplied by u to recover actual gradient
+ */
+double F0ElCpClamp::evalGrad(Element const &e, std::vector<double> const &phi, size_t k) const
+{
+ double r = F0ElRhoBase::eval(e, phi, k);
+ double dr = F0ElRhoBase::evalGrad(e, phi, k);
+ return 2 / ((gamma - 1) * mInf * mInf) * sqrt(r * r * r * r * r) * dr; // particularized to gamma=1.4
}
/**
- * @brief Evaluate the linear pressure coefficient (constant over element)
+ * @brief Evaluate the linear pressure coefficient
+ *
+ * Cp = 1-u^2
*/
-double F0ElCpL::eval(Element const &e, std::vector<double> const &u, size_t k) const
+double F0ElCpLin::eval(Element const &e, std::vector<double> const &phi, size_t k) const
{
- Eigen::VectorXd gradU = e.computeGradient(u, k); // velocity
- double gradU2 = gradU.squaredNorm(); // velocity squared
- return 1 - gradU2;
+ return 1 - e.computeGradient(phi, k).squaredNorm();
}
-double F0ElCpL::evalGrad(Element const &e, std::vector<double> const &u, size_t k) const
+/**
+ * @brief Evaluate the linear pressure coefficient gradient
+ *
+ * dCp = -2*u
+ * Computed value has to multiplied by u to recover actual gradient
+ */
+double F0ElCpLin::evalGrad(Element const &e, std::vector<double> const &phi, size_t k) const
{
- return -2.; // has to be multiplied by grad u to recover true derivative
+ return -2.;
}
\ No newline at end of file
diff --git a/dart/src/wF0El.h b/dart/src/wF0El.h
index fc4850c..4f4a774 100644
--- a/dart/src/wF0El.h
+++ b/dart/src/wF0El.h
@@ -34,9 +34,9 @@ public:
F0El() : Observer() {}
virtual ~F0El() {}
- virtual double eval(tbox::Element const &e, std::vector<double> const &u,
+ virtual double eval(tbox::Element const &e, std::vector<double> const &phi,
size_t k) const = 0;
- virtual double evalGrad(tbox::Element const &e, std::vector<double> const &u,
+ virtual double evalGrad(tbox::Element const &e, std::vector<double> const &phi,
size_t k) const = 0;
};
@@ -51,55 +51,93 @@ public:
F0ElC(double _v) : F0El(), v(_v) {}
virtual void update(double _v) override;
- virtual double eval(tbox::Element const &e, std::vector<double> const &u,
+ virtual double eval(tbox::Element const &e, std::vector<double> const &phi,
size_t k) const override;
- virtual double evalGrad(tbox::Element const &e, std::vector<double> const &u,
+ virtual double evalGrad(tbox::Element const &e, std::vector<double> const &phi,
size_t k) const override;
};
/**
- * @brief Nonlinear density
+ * @brief Base of isentropic density
*
- * The density is limited (Padé approximation) for large value of velocity and
- * particularized for diatomic gas. The input _mC2 is the square of the critical
- * Mach number used to compute the critical velocity.
+ * This function computes the base of the isentropic density, and is then used
+ * to compute the actual density, Mach number and pressure coefficient. This
+ * function is clamped for large velocities using a Padé approximation and
+ * particularized to diatomic gas.
+ */
+class DART_API F0ElRhoBase : public F0El
+{
+protected:
+ double gamma; ///< heat capacity ratio (diatomic gas only)
+ double mInf; ///< freestream Mach number
+private:
+ double uC; ///< critical velocity
+ double rC; ///< critical density base
+ double beta; ///< ratio for Padé
+
+public:
+ F0ElRhoBase(double _mInf, double _mC2 = 5);
+
+ virtual double eval(tbox::Element const &e, std::vector<double> const &phi,
+ size_t k) const override;
+ virtual double evalGrad(tbox::Element const &e, std::vector<double> const &phi,
+ size_t k) const override;
+};
+
+/**
+ * @brief Isentropic density
+ *
+ * The density is particularized to diatomic gas.
*/
class DART_API F0ElRho : public F0El
{
- double gamma; ///< heat capacity ratio (diatomic gas only)
- double mInf; ///< freestream Mach number
- double gradUC; ///< critical velocity
+ double gamma; ///< heat capacity ratio (diatomic gas only)
+ double mInf; ///< freestream Mach number
public:
- F0ElRho(double _mInf, double _mC2 = 5) : F0El(), gamma(1.4), mInf(_mInf)
- {
- gradUC = sqrt(_mC2 * (1 / (mInf * mInf) + (gamma - 1) / 2) / (1 + (gamma - 1) / 2 * _mC2));
- }
+ F0ElRho(double _mInf) : F0El(), gamma(1.4), mInf(_mInf) {}
+
+ virtual double eval(tbox::Element const &e, std::vector<double> const &phi,
+ size_t k) const override;
+ virtual double evalGrad(tbox::Element const &e, std::vector<double> const &phi,
+ size_t k) const override;
+};
+
+/**
+ * @brief Isentropic clamped density
+ *
+ * The density is clamped (Padé approximation) for large values of velocity and
+ * particularized to diatomic gas.
+ */
+class DART_API F0ElRhoClamp : public F0ElRhoBase
+{
+public:
+ F0ElRhoClamp(double _mInf, double _mC2 = 5) : F0ElRhoBase(_mInf, _mC2) {}
- virtual double eval(tbox::Element const &e, std::vector<double> const &u,
+ virtual double eval(tbox::Element const &e, std::vector<double> const &phi,
size_t k) const override;
- virtual double evalGrad(tbox::Element const &e, std::vector<double> const &u,
+ virtual double evalGrad(tbox::Element const &e, std::vector<double> const &phi,
size_t k) const override;
};
/**
* @brief Linear density (constant)
*/
-class DART_API F0ElRhoL : public F0El
+class DART_API F0ElRhoLin : public F0El
{
public:
- F0ElRhoL() : F0El() {}
+ F0ElRhoLin() : F0El() {}
- virtual double eval(tbox::Element const &e, std::vector<double> const &u,
+ virtual double eval(tbox::Element const &e, std::vector<double> const &phi,
size_t k) const override;
- virtual double evalGrad(tbox::Element const &e, std::vector<double> const &u,
+ virtual double evalGrad(tbox::Element const &e, std::vector<double> const &phi,
size_t k) const override;
};
/**
- * @brief Nonlinear Mach number
+ * @brief Isentropic Mach number
*
- * Particularized for diatomic gas
+ * The Mach number is particularized to diatomic gas.
*/
class DART_API F0ElMach : public F0El
{
@@ -109,30 +147,47 @@ class DART_API F0ElMach : public F0El
public:
F0ElMach(double _mInf) : F0El(), gamma(1.4), mInf(_mInf) {}
- virtual double eval(tbox::Element const &e, std::vector<double> const &u,
+ virtual double eval(tbox::Element const &e, std::vector<double> const &phi,
+ size_t k) const override;
+ virtual double evalGrad(tbox::Element const &e, std::vector<double> const &phi,
+ size_t k) const override;
+};
+
+/**
+ * @brief Isentropic clamped Mach number
+ *
+ * The Mach number is clamped (Padé approximation) for large values of velocity
+ * and particularized to diatomic gas.
+ */
+class DART_API F0ElMachClamp : public F0ElRhoBase
+{
+public:
+ F0ElMachClamp(double _mInf, double _mC2 = 5) : F0ElRhoBase(_mInf, _mC2) {}
+
+ virtual double eval(tbox::Element const &e, std::vector<double> const &phi,
size_t k) const override;
- virtual double evalGrad(tbox::Element const &e, std::vector<double> const &u,
+ virtual double evalGrad(tbox::Element const &e, std::vector<double> const &phi,
size_t k) const override;
};
/**
* @brief Linear Mach number (constant)
*/
-class DART_API F0ElMachL : public F0El
+class DART_API F0ElMachLin : public F0El
{
public:
- F0ElMachL() : F0El() {}
+ F0ElMachLin() : F0El() {}
- virtual double eval(tbox::Element const &e, std::vector<double> const &u,
+ virtual double eval(tbox::Element const &e, std::vector<double> const &phi,
size_t k) const override;
- virtual double evalGrad(tbox::Element const &e, std::vector<double> const &u,
+ virtual double evalGrad(tbox::Element const &e, std::vector<double> const &phi,
size_t k) const override;
};
/**
- * @brief Nonlinear pressure coefficient
+ * @brief Isentropic pressure coefficient
*
- * Particularized for diatomic gas
+ * The pressure coefficient is particularized to diatomic gas.
*/
class DART_API F0ElCp : public F0El
{
@@ -142,27 +197,43 @@ class DART_API F0ElCp : public F0El
public:
F0ElCp(double _mInf) : F0El(), gamma(1.4), mInf(_mInf) {}
- virtual double eval(tbox::Element const &e, std::vector<double> const &u,
+ virtual double eval(tbox::Element const &e, std::vector<double> const &phi,
size_t k) const override;
- virtual double evalGrad(tbox::Element const &e, std::vector<double> const &u,
+ virtual double evalGrad(tbox::Element const &e, std::vector<double> const &phi,
size_t k) const override;
};
/**
- * @brief Linear pressure coefficient
+ * @brief Isentropic clamped pressure coefficient
+ *
+ * The pressure coefficient is clamped (Padé approximation) for large values of
+ * velocity and particularized to diatomic gas.
*/
-class DART_API F0ElCpL : public F0El
+class DART_API F0ElCpClamp : public F0ElRhoBase
{
+public:
+ F0ElCpClamp(double _mInf, double _mC2 = 5) : F0ElRhoBase(_mInf, _mC2) {}
+ virtual double eval(tbox::Element const &e, std::vector<double> const &phi,
+ size_t k) const override;
+ virtual double evalGrad(tbox::Element const &e, std::vector<double> const &phi,
+ size_t k) const override;
+};
+
+/**
+ * @brief Linear pressure coefficient
+ */
+class DART_API F0ElCpLin : public F0El
+{
public:
- F0ElCpL() : F0El() {}
+ F0ElCpLin() : F0El() {}
- virtual double eval(tbox::Element const &e, std::vector<double> const &u,
+ virtual double eval(tbox::Element const &e, std::vector<double> const &phi,
size_t k) const override;
- virtual double evalGrad(tbox::Element const &e, std::vector<double> const &u,
+ virtual double evalGrad(tbox::Element const &e, std::vector<double> const &phi,
size_t k) const override;
};
} // namespace dart
-#endif //WF0EL_H
+#endif // WF0EL_H
diff --git a/dart/src/wFluid.cpp b/dart/src/wFluid.cpp
index 1713078..a209f4b 100644
--- a/dart/src/wFluid.cpp
+++ b/dart/src/wFluid.cpp
@@ -54,15 +54,15 @@ void Fluid::initFun(double mInf, int dim, double alpha, double beta)
{
if (mInf == 0)
{
- rho = new F0ElRhoL();
- mach = new F0ElMachL();
- cP = new F0ElCpL();
+ rho = new F0ElRhoLin();
+ mach = new F0ElMachLin();
+ cP = new F0ElCpLin();
}
else
{
- rho = new F0ElRho(mInf);
- mach = new F0ElMach(mInf);
- cP = new F0ElCp(mInf);
+ rho = new F0ElRhoClamp(mInf);
+ mach = new F0ElMachClamp(mInf);
+ cP = new F0ElCpClamp(mInf);
}
phiInf = new F0PsPhiInf(dim, alpha, beta);
}
--
GitLab