From 9d36f3d7a505600173747999e0bd2ffd11d7a1f7 Mon Sep 17 00:00:00 2001
From: acrovato <a.crovato@uliege.be>
Date: Fri, 20 Dec 2024 15:43:58 +0100
Subject: [PATCH] Add provisions for second order upwind formulation
---
dart/src/wFluid.cpp | 1 +
dart/src/wFluid.h | 2 +
dart/src/wNewton.cpp | 24 +++++-
dart/src/wPotentialResidual.cpp | 145 ++++++++++++++++++++++++++++++--
dart/src/wPotentialResidual.h | 8 +-
dart/validation/onera.py | 2 +
6 files changed, 172 insertions(+), 10 deletions(-)
diff --git a/dart/src/wFluid.cpp b/dart/src/wFluid.cpp
index f150924..07b018e 100644
--- a/dart/src/wFluid.cpp
+++ b/dart/src/wFluid.cpp
@@ -98,6 +98,7 @@ void Fluid::createMap()
}
// Pair an element and its adjacent elements
adjMap[i] = std::pair<Element *, std::vector<Element *>>(e, vAdjE);
+ adjMap2[e] = i;
i++;
}
}
diff --git a/dart/src/wFluid.h b/dart/src/wFluid.h
index 600c35d..f8c594e 100644
--- a/dart/src/wFluid.h
+++ b/dart/src/wFluid.h
@@ -29,6 +29,7 @@ namespace dart
/**
* @brief Group of elements representing the fluid
* @authors Adrien Crovato
+ * @todo replace adjMap by adjMap2
*/
class DART_API Fluid : public tbox::Group
{
@@ -44,6 +45,7 @@ public:
F0PsPhiInf *phiInf; ///< freestream potential
std::vector<std::pair<tbox::Element *, std::vector<tbox::Element *>>> adjMap; ///< map of adjacent elements
+ std::map<tbox::Element *, size_t> adjMap2; ///< map of adjacent elements
#endif
Fluid(std::shared_ptr<tbox::MshData> _msh, int no, double mInf, int dim, double alpha, double beta);
diff --git a/dart/src/wNewton.cpp b/dart/src/wNewton.cpp
index f4b932f..30672cf 100644
--- a/dart/src/wNewton.cpp
+++ b/dart/src/wNewton.cpp
@@ -293,9 +293,11 @@ void Newton::buildJac(Eigen::SparseMatrix<double, Eigen::RowMajor> &J)
Element *e = p.first;
// Upwind element
Element *eU = p.second[0];
+ // Upwind-upwind element
+ Element *eUU = fluid->adjMap[fluid->adjMap2.at(eU)].second[0];
// Build elementary matrices
- Eigen::MatrixXd Ae1, Ae2;
- std::tie(Ae1, Ae2) = PotentialResidual::buildGradientFlow(*e, *eU, phi, *fluid, muCv, mCOv);
+ Eigen::MatrixXd Ae1, Ae2, Ae3;
+ std::tie(Ae1, Ae2, Ae3) = PotentialResidual::buildGradientFlow(*e, *eU, *eUU, phi, *fluid, muCv, mCOv);
// Assembly (subsonic)
std::deque<Eigen::Triplet<double>> &locT = cmbT.local();
@@ -323,6 +325,20 @@ void Newton::buildJac(Eigen::SparseMatrix<double, Eigen::RowMajor> &J)
locT.push_back(Eigen::Triplet<double>(rows[nodi->row], nodj->row, Ae2(ii, jj)));
}
}
+ // Second order supersonic
+ double ratio = (e->computeGradient(phi, 0)(0) - eU->computeGradient(phi, 0)(0)) / (eU->computeGradient(phi, 0)(0) - eUU->computeGradient(phi, 0)(0));
+ if (ratio >= 0.)
+ {
+ for (size_t ii = 0; ii < e->nodes.size(); ++ii)
+ {
+ Node *nodi = e->nodes[ii];
+ for (size_t jj = 0; jj < eUU->nodes.size(); ++jj)
+ {
+ Node *nodj = eUU->nodes[jj];
+ locT.push_back(Eigen::Triplet<double>(rows[nodi->row], nodj->row, Ae3(ii, jj)));
+ }
+ }
+ }
}
});
tms["0-Jbase"].stop();
@@ -430,8 +446,10 @@ void Newton::buildRes(Eigen::Map<Eigen::VectorXd> &R)
Element *e = p.first;
// Upwind element
Element *eU = p.second[0];
+ // Upwind-upwind element
+ Element *eUU = fluid->adjMap[fluid->adjMap2.at(eU)].second[0];
// Build elementary vector
- Eigen::VectorXd be = PotentialResidual::build(*e, *eU, phi, *fluid, muCv, mCOv);
+ Eigen::VectorXd be = PotentialResidual::build(*e, *eU, *eUU, phi, *fluid, muCv, mCOv);
// Assembly (subsonic)
tbb::spin_mutex::scoped_lock lock(mutex);
for (size_t ii = 0; ii < e->nodes.size(); ++ii)
diff --git a/dart/src/wPotentialResidual.cpp b/dart/src/wPotentialResidual.cpp
index 0a684f1..6ee3654 100644
--- a/dart/src/wPotentialResidual.cpp
+++ b/dart/src/wPotentialResidual.cpp
@@ -52,7 +52,7 @@ Eigen::MatrixXd PotentialResidual::buildFixed(Element const &e, std::vector<doub
*
* b = \int ( (1-mu)*rho + mu*rhoU ) * grad_phi * grad_psi dV
*/
-Eigen::VectorXd PotentialResidual::build(Element const &e, Element const &eU, std::vector<double> const &phi, Fluid const &fluid, double muC, double mCO)
+Eigen::VectorXd PotentialResidual::build(Element const &e, Element const &eU, Element const &eUU, std::vector<double> const &phi, Fluid const &fluid, double muC, double mCO)
{
// Get pre-computed values
Cache &cache = e.getVCache();
@@ -76,6 +76,63 @@ Eigen::VectorXd PotentialResidual::build(Element const &e, Element const &eU, st
// contribution of upwind element
for (size_t k = 0; k < gauss.getN(); ++k)
b += mu * rhoU * (e.getJinv(k) * cache.getDsf(k)).transpose() * e.computeGradient(phi, k) * gauss.getW(k) * e.getDetJ(k);
+
+ // second order upwind
+ double rhoUU = fluid.rho->eval(eUU, phi, 0);
+ double ratio = (e.computeGradient(phi, 0)(0) - eU.computeGradient(phi, 0)(0)) / (eU.computeGradient(phi, 0)(0) - eUU.computeGradient(phi, 0)(0));
+ if (ratio >= 0.)
+ {
+ double alpha = 1.0;
+ if (ratio >= 1.00)
+ alpha = 1.00;
+ else if (ratio >= 0.95)
+ alpha = 0.95;
+ else if (ratio >= 0.90)
+ alpha = 0.90;
+ else if (ratio >= 0.85)
+ alpha = 0.85;
+ else if (ratio >= 0.80)
+ alpha = 0.80;
+ else if (ratio >= 0.75)
+ alpha = 0.75;
+ else if (ratio >= 0.70)
+ alpha = 0.70;
+ else if (ratio >= 0.65)
+ alpha = 0.65;
+ else if (ratio >= 0.60)
+ alpha = 0.60;
+ else if (ratio >= 0.55)
+ alpha = 0.55;
+ else if (ratio >= 0.50)
+ alpha = 0.50;
+ else if (ratio >= 0.45)
+ alpha = 0.45;
+ else if (ratio >= 0.40)
+ alpha = 0.40;
+ else if (ratio >= 0.35)
+ alpha = 0.35;
+ else if (ratio >= 0.30)
+ alpha = 0.30;
+ else if (ratio >= 0.25)
+ alpha = 0.25;
+ else if (ratio >= 0.20)
+ alpha = 0.20;
+ else if (ratio >= 0.15)
+ alpha = 0.15;
+ else if (ratio >= 0.10)
+ alpha = 0.10;
+ else if (ratio >= 0.05)
+ alpha = 0.05;
+ else
+ alpha = 0.025;
+ //alpha = (alpha < 0.5) ? alpha*2 : 1.0;
+
+ for (size_t k = 0; k < gauss.getN(); ++k)
+ {
+ b += alpha * mu * rhoU * (e.getJinv(k) * cache.getDsf(k)).transpose() * e.computeGradient(phi, k) * gauss.getW(k) * e.getDetJ(k);
+ b -= alpha * mu * rhoUU * (e.getJinv(k) * cache.getDsf(k)).transpose() * e.computeGradient(phi, k) * gauss.getW(k) * e.getDetJ(k);
+ }
+ }
}
return b;
}
@@ -89,7 +146,7 @@ Eigen::VectorXd PotentialResidual::build(Element const &e, Element const &eU, st
* + \int ( (1-mu)*rho + mu*rhoU ) * dgrad_phi * grad_psi dV
* + \int (rhoU - rho)*dmu * grad_phi * grad_psi dV
*/
-std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> PotentialResidual::buildGradientFlow(tbox::Element const &e, tbox::Element const &eU, std::vector<double> const &phi, Fluid const &fluid, double muC, double mCO)
+std::tuple<Eigen::MatrixXd, Eigen::MatrixXd, Eigen::MatrixXd> PotentialResidual::buildGradientFlow(Element const &e, Element const &eU, Element const &eUU, std::vector<double> const &phi, Fluid const &fluid, double muC, double mCO)
{
// Get pre-computed values
Cache &cache = e.getVCache();
@@ -111,7 +168,7 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> PotentialResidual::buildGradientFlo
}
// Supersonic contribution
- Eigen::MatrixXd A2;
+ Eigen::MatrixXd A2, A3;
double machC = fluid.mach->eval(e, phi, 0);
double machU = fluid.mach->eval(eU, phi, 0);
double mach = (machC < machU) ? machU : machC;
@@ -128,7 +185,7 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> PotentialResidual::buildGradientFlo
// upwinding
A1 *= 1 - mu;
- A2 = Eigen::MatrixXd::Zero(e.nodes.size(), e.nodes.size());
+ A2 = Eigen::MatrixXd::Zero(e.nodes.size(), eU.nodes.size());
for (size_t k = 0; k < gauss.getN(); ++k)
{
// Gauss point and determinant
@@ -147,8 +204,86 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> PotentialResidual::buildGradientFlo
else
A2 += dmu * (rhoU - fluid.rho->eval(e, phi, k)) * fluid.mach->evalGrad(eU, phi, k) * dSf.transpose() * dPhi * dPhiU.transpose() * dSfU * wdj;
}
+
+ // second order upwind
+ double rhoUU = fluid.rho->eval(eUU, phi, 0);
+ double ratio = (e.computeGradient(phi, 0)(0) - eU.computeGradient(phi, 0)(0)) / (eU.computeGradient(phi, 0)(0) - eUU.computeGradient(phi, 0)(0));
+ if (ratio >= 0.)
+ {
+ double alpha = 1.0;
+ if (ratio >= 1.0)
+ alpha = 1.00;
+ else if (ratio >= 0.95)
+ alpha = 0.95;
+ else if (ratio >= 0.90)
+ alpha = 0.90;
+ else if (ratio >= 0.85)
+ alpha = 0.85;
+ else if (ratio >= 0.80)
+ alpha = 0.80;
+ else if (ratio >= 0.75)
+ alpha = 0.75;
+ else if (ratio >= 0.70)
+ alpha = 0.70;
+ else if (ratio >= 0.65)
+ alpha = 0.65;
+ else if (ratio >= 0.60)
+ alpha = 0.60;
+ else if (ratio >= 0.55)
+ alpha = 0.55;
+ else if (ratio >= 0.50)
+ alpha = 0.50;
+ else if (ratio >= 0.45)
+ alpha = 0.45;
+ else if (ratio >= 0.40)
+ alpha = 0.40;
+ else if (ratio >= 0.35)
+ alpha = 0.35;
+ else if (ratio >= 0.30)
+ alpha = 0.30;
+ else if (ratio >= 0.25)
+ alpha = 0.25;
+ else if (ratio >= 0.20)
+ alpha = 0.20;
+ else if (ratio >= 0.15)
+ alpha = 0.15;
+ else if (ratio >= 0.10)
+ alpha = 0.10;
+ else if (ratio >= 0.05)
+ alpha = 0.05;
+ else
+ alpha = 0.025;
+ //alpha = (alpha < 0.5) ? alpha*2 : 1.0;
+
+ double dRhoUU = fluid.rho->evalGrad(eUU, phi, 0);
+ Eigen::VectorXd dPhiUU = eUU.computeGradient(phi, 0);
+ Eigen::MatrixXd const &dSfUU = eUU.getJinv(0) * eUU.getVCache().getDsf(0);
+
+ A3 = Eigen::MatrixXd::Zero(e.nodes.size(), eUU.nodes.size());
+ for (size_t k = 0; k < gauss.getN(); ++k)
+ {
+ // Gauss point and determinant
+ double wdj = gauss.getW(k) * e.getDetJ(k);
+ // Shape functions and solution gradient
+ Eigen::MatrixXd const &dSf = e.getJinv(k) * cache.getDsf(k);
+ Eigen::VectorXd dPhi = e.computeGradient(phi, k);
+
+ // (mu * rhoU - mu * rhUU) * grad_phi * grad_psi
+ A1 += alpha * mu * rhoU * dSf.transpose() * dSf * wdj;
+ A1 -= alpha * mu * rhoUU * dSf.transpose() * dSf * wdj;
+ // mu * drhoU * grad_phi * grad_psi
+ A2 += alpha * mu * dRhoU * dSf.transpose() * dPhi * dPhiU.transpose() * dSfU * wdj;
+ // mu * drhoUU * grad_phi * grad_psi
+ A3 -= alpha * mu * dRhoUU * dSf.transpose() * dPhi * dPhiUU.transpose() * dSfUU * wdj;
+ // dmu(U) * (rhoU-rhoUU) * grad_phi * grad_psi
+ if (machC >= machU)
+ A1 += alpha * dmu * (rhoU - rhoUU) * fluid.mach->evalGrad(e, phi, k) * dSf.transpose() * dPhi * dPhi.transpose() * dSf * wdj;
+ else
+ A2 += alpha * dmu * (rhoU - rhoUU) * fluid.mach->evalGrad(eU, phi, k) * dSf.transpose() * dPhi * dPhiU.transpose() * dSfU * wdj;
+ }
+ }
}
- return std::make_tuple(A1, A2);
+ return std::make_tuple(A1, A2, A3);
}
/**
diff --git a/dart/src/wPotentialResidual.h b/dart/src/wPotentialResidual.h
index 5007fe1..e5b07cf 100644
--- a/dart/src/wPotentialResidual.h
+++ b/dart/src/wPotentialResidual.h
@@ -29,6 +29,10 @@ namespace dart
* @brief Formulation of nonlinear potential equation residuals
* @todo mu as F0ElMu?
* @todo split stabilization terms from subsonic terms?
+ * @todo check that dR = J
+ * @todo fix SOU oscillations (improve limiter, without limiter, use divU/c (Ducross), use FOU then SOU?)
+ * @todo use aprox. Newton (do not use Ae3)?
+ * @todo clean second order upwind formulation
*/
class DART_API PotentialResidual
{
@@ -36,8 +40,8 @@ public:
// Picard
static Eigen::MatrixXd buildFixed(tbox::Element const &e, std::vector<double> const &phi, Fluid const &fluid);
// Newton
- static Eigen::VectorXd build(tbox::Element const &e, tbox::Element const &eU, std::vector<double> const &phi, Fluid const &fluid, double muC, double mCO);
- static std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> buildGradientFlow(tbox::Element const &e, tbox::Element const &eU, std::vector<double> const &phi, Fluid const &fluid, double muC, double mCO);
+ static Eigen::VectorXd build(tbox::Element const &e, tbox::Element const &eU, tbox::Element const &eUU, std::vector<double> const &phi, Fluid const &fluid, double muC, double mCO);
+ static std::tuple<Eigen::MatrixXd, Eigen::MatrixXd, Eigen::MatrixXd> buildGradientFlow(tbox::Element const &e, tbox::Element const &eU, tbox::Element const &eUU, std::vector<double> const &phi, Fluid const &fluid, double muC, double mCO);
// Adjoint
static std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> buildGradientMesh(tbox::Element const &e, tbox::Element const &eU, std::vector<double> const &phi, Fluid const &fluid, int nDim, double muC, double mCO);
};
diff --git a/dart/validation/onera.py b/dart/validation/onera.py
index cc1c3d3..8b65b86 100644
--- a/dart/validation/onera.py
+++ b/dart/validation/onera.py
@@ -84,6 +84,8 @@ def main():
floU.write_slices(msh.name,[0.01, 0.239, 0.526, 0.777, 0.957, 1.076, 1.136, 1.184], 5)
data = np.loadtxt(f'{msh.name}_slice_2.dat', delimiter=',', skiprows=2)
floD.plot(data[:,0], data[:,-1], {'xlabel': 'x/c', 'ylabel': 'Cp', 'title': 'Pressure distribution at MAC', 'invert': True})
+ data = np.loadtxt(f'{msh.name}_slice_3.dat', delimiter=',', skiprows=2)
+ floD.plot(data[:,0], data[:,-1], {'xlabel': 'x/c', 'ylabel': 'Cp', 'title': 'Pressure distribution at MAC', 'invert': True})
tms['post'].stop()
# display timers
--
GitLab