Comments (3)
pjb236
commented on August 19, 2024
@fnrizzi I have 1d Roe flux working: it was able to solve the Sod Shock tube correctly. I've added a test for the flux function, but we don't have a test where we solve a 1D Euler problem with Roe Flux. Should I add a new test or modify one of the existing tests, such as the 1D sod cases, to use the Roe flux instead? Perhaps I could add separate sub-directories in 1d_sod for Roe flux tests like the weno3 and weno5 subdirectories?
I'll start adding 2D and 3D Roe fluxes tomorrow.
from pressio-demoapps.
pjb236
commented on August 19, 2024
@fnrizzi Good news: I've implemented 2D Roe flux and it works for the 2D normal shock and 2D double Mach reflection with WENO3. Bad news: I can't get it to work for WENO5. The branch is up to date if you want to try, it is set up to make a WENO5 test for the double Mach reflection that uses the Roe Flux. Looks like some of the boundaries are unstable after one timestep.
from pressio-demoapps.
fnrizzi
commented on August 19, 2024
template<class T, typename sc_t>
void eeRoeFluxThreeDof(T & F,
const T & qL,
const T & qR,
const sc_t gamma)
{
constexpr auto one = static_cast<sc_t>(1);
constexpr auto two = static_cast<sc_t>(2);
constexpr auto half = one/two;
constexpr sc_t es_ef{0.1}; // Entropy fix parameter
constexpr sc_t es{1.e-30};
sc_t FL[3];
sc_t FR[3];
sc_t K1[3];
sc_t K2[3];
sc_t K3[3];
// Compute Fluxes
// left
const auto rL = qL(0);
const auto uL = qL(1)/(rL + es);
const auto pL = (gamma-1)*( qL(2) - half*rL*(uL*uL) );
const auto HL = ( qL(2) + pL ) / rL;
FL[0] = rL*uL;
FL[1] = rL*uL*uL + pL;
FL[2] = rL*uL*HL;
// right
const auto rR = qR(0);
const auto uR = qR(1)/(rR + es);
const auto pR = (gamma-1)*( qR(2) - half*rR*(uR*uR) );
const auto HR = ( qR(2) + pR ) / rR;
FR[0] = rR*uR;
FR[1] = rR*uR*uR + pR;
FR[2] = rR*uR*HR;
// Compute Roe averaged variables
const auto rLsqrt = std::sqrt(rL);
const auto rRsqrt = std::sqrt(rR);
const auto rSqrtSum = rLsqrt + rRsqrt;
const auto uTilde = (rLsqrt * uL + rRsqrt * uR) / rSqrtSum;
const auto HTilde = (rLsqrt * HL + rRsqrt * HR) / rSqrtSum;
const auto aTilde = std::sqrt((gamma-1)*(HTilde - half * uTilde * uTilde));
// Compute Eigenvalues
auto lam1 = uTilde - aTilde;
auto lam2 = uTilde;
auto lam3 = uTilde + aTilde;
// Entropy fix (San and Kara: Evaluation of Riemann flux solvers for WENO reconstruction
// schemes: Kelvin–Helmholtz instability), Computers & Fluids, 117 (2015), 24-41. Eqn. 35
const auto tol = 2 * es_ef * aTilde;
if (std::abs(lam1) < tol) lam1 = lam1 * lam1 / (two * tol) + es_ef * aTilde;
if (std::abs(lam2) < tol) lam2 = lam2 * lam2 / (two * tol) + es_ef * aTilde;
if (std::abs(lam3) < tol) lam3 = lam3 * lam3 / (two * tol) + es_ef * aTilde;
// Compute Eigenvectors
K1[0] = one;
K1[1] = uTilde - aTilde;
K1[2] = HTilde - uTilde * aTilde;
K2[0] = one;
K2[1] = uTilde;
K2[2] = half * uTilde * uTilde;
K3[0] = one;
K3[1] = uTilde + aTilde;
K3[2] = HTilde + uTilde * aTilde;
// Compute Eigenvector weights
const auto dq1 = qR[0] - qL[0];
const auto dq2 = qR[1] - qL[1];
const auto dq3 = qR[2] - qL[2];
const auto alpha2 = (gamma-1) / (aTilde * aTilde) * (dq1 * (HTilde - uTilde * uTilde) + uTilde * dq2 - dq3);
const auto alpha1 = (half / aTilde) * (dq1 * (uTilde + aTilde) - dq2 - aTilde * alpha2);
const auto alpha3 = dq1 - alpha1 - alpha2;
// Compute Flux
F(0) = half * (FL[0] + FR[0] - alpha1 * std::abs(lam1) * K1[0] - alpha2 * std::abs(lam2) * K2[0] - alpha3 * std::abs(lam3) * K3[0] );
F(1) = half * (FL[1] + FR[1] - alpha1 * std::abs(lam1) * K1[1] - alpha2 * std::abs(lam2) * K2[1] - alpha3 * std::abs(lam3) * K3[1] );
F(2) = half * (FL[2] + FR[2] - alpha1 * std::abs(lam1) * K1[2] - alpha2 * std::abs(lam2) * K2[2] - alpha3 * std::abs(lam3) * K3[2] );
}
template<class T, class normal_t, typename sc_t>
void eeRoeFluxFourDof(T & F,
const T & qL,
const T & qR,
const normal_t & n,
const sc_t gamma)
{
constexpr auto one = static_cast<sc_t>(1);
constexpr auto two = static_cast<sc_t>(2);
constexpr auto half = one/two;
constexpr sc_t es_ef{0.1}; // Entropy fix parameter
constexpr sc_t es{1.e-30};
sc_t FL[4];
sc_t FR[4];
sc_t K1[4];
sc_t K2[4];
sc_t K3[4];
sc_t K4[4];
// Compute Fluxes
// left
const auto rL = qL(0);
const auto uL = qL(1)/(rL + es);
const auto vL = qL(2)/(rL + es);
const auto unL = uL*n[0] + vL*n[1];
const auto pL = (gamma-1)*( qL(3) - half*rL*(uL*uL+vL*vL) );
const auto HL = ( qL(3) + pL ) / rL;
FL[0] = rL*unL;
FL[1] = rL*unL*uL + pL*n[0];
FL[2] = rL*unL*vL + pL*n[1];
FL[3] = rL*unL*HL;
// right
const auto rR = qR(0);
const auto uR = qR(1)/(rR + es);
const auto vR = qR(2)/(rR + es);
const auto unR = uR*n[0] + vR*n[1];
const auto pR = (gamma-1)*( qR(3) - half*rR*(uR*uR+vR*vR) );
const auto HR = ( qR(3) + pR ) / rR;
FR[0] = rR*unR;
FR[1] = rR*unR*uR + pR*n[0];
FR[2] = rR*unR*vR + pR*n[1];
FR[3] = rR*unR*HR;
// Compute Roe averaged variables
const auto rLsqrt = std::sqrt(rL);
const auto rRsqrt = std::sqrt(rR);
const auto rSqrtSum = rLsqrt + rRsqrt;
const auto uTilde = (rLsqrt * uL + rRsqrt * uR) / rSqrtSum;
const auto vTilde = (rLsqrt * vL + rRsqrt * vR) / rSqrtSum;
const auto unTilde = (rLsqrt * unL + rRsqrt * unR) / rSqrtSum;
const auto HTilde = (rLsqrt * HL + rRsqrt * HR) / rSqrtSum;
const auto VTilde2 = uTilde * uTilde + vTilde * vTilde;
const auto aTilde = std::sqrt((gamma-1)*(HTilde - half * VTilde2));
// Compute Eigenvalues
auto lam1 = unTilde - aTilde;
auto lam2 = unTilde;
auto lam3 = unTilde;
auto lam4 = unTilde + aTilde;
// Entropy fix (San and Kara: Evaluation of Riemann flux solvers for WENO reconstruction
// schemes: Kelvin–Helmholtz instability), Computers & Fluids, 117 (2015), 24-41. Eqn. 35
const auto tol = two * es_ef * aTilde;
if (std::abs(lam1) < tol) lam1 = lam1 * lam1 / (two * tol) + es_ef * aTilde;
if (std::abs(lam2) < tol) lam2 = lam2 * lam2 / (two * tol) + es_ef * aTilde;
if (std::abs(lam3) < tol) lam3 = lam3 * lam3 / (two * tol) + es_ef * aTilde;
if (std::abs(lam4) < tol) lam4 = lam4 * lam4 / (two * tol) + es_ef * aTilde;
// Compute Eigenvectors
K1[0] = one;
K1[1] = uTilde - aTilde*n[0];
K1[2] = vTilde - aTilde*n[1];
K1[3] = HTilde - unTilde * aTilde;
K2[0] = one;
K2[1] = uTilde;
K2[2] = vTilde;
K2[3] = half * VTilde2;
K3[0] = 0.0;
K3[1] = n[1];
K3[2] = n[0];
K3[3] = n[0]*vTilde + n[1]*uTilde;
K4[0] = one;
K4[1] = uTilde + aTilde*n[0];
K4[2] = vTilde + aTilde*n[1];
K4[3] = HTilde + unTilde * aTilde;
// Compute Eigenvector weights
const auto dq1 = qR[0] - qL[0];
const auto dq2 = qR[1] - qL[1];
const auto dq3 = qR[2] - qL[2];
const auto dq4 = qR[3] - qL[3];
const auto alpha3 = (dq3*n[0] - dq2*n[1] + dq1*n[1]*uTilde - dq1*n[0]*vTilde)/(n[0]*n[0] - n[1]*n[1]);
const auto alpha2 = (gamma-1) / (aTilde * aTilde) * (dq1 * (HTilde - VTilde2) + uTilde * dq2 -dq4 +dq3*vTilde);
const auto b1 = dq1-alpha2;
const auto b2 = dq2-alpha2*uTilde-alpha3*n[1];
const auto b3 = dq3 - alpha2*vTilde-alpha3*n[0];
const auto alpha1 = (half / aTilde) * ((b1*(uTilde+aTilde*n[0])-b2)*n[0] + (b1*(vTilde+aTilde*n[1])-b3)*n[1]);
const auto alpha4 = dq1 - alpha1 - alpha2;
// Compute Flux
F(0) = half * (FL[0] + FR[0] - alpha1 * std::abs(lam1) * K1[0] - alpha2 * std::abs(lam2) * K2[0] - alpha3 * std::abs(lam3) * K3[0] - alpha4 * std::abs(lam4) * K4[0]);
F(1) = half * (FL[1] + FR[1] - alpha1 * std::abs(lam1) * K1[1] - alpha2 * std::abs(lam2) * K2[1] - alpha3 * std::abs(lam3) * K3[1] - alpha4 * std::abs(lam4) * K4[1]);
F(2) = half * (FL[2] + FR[2] - alpha1 * std::abs(lam1) * K1[2] - alpha2 * std::abs(lam2) * K2[2] - alpha3 * std::abs(lam3) * K3[2] - alpha4 * std::abs(lam4) * K4[2]);
F(3) = half * (FL[3] + FR[3] - alpha1 * std::abs(lam1) * K1[3] - alpha2 * std::abs(lam2) * K2[3] - alpha3 * std::abs(lam3) * K3[3] - alpha4 * std::abs(lam4) * K4[3]);
}since code has changed quite a bit, putting these here and we also need at some point to add Jacobians
from pressio-demoapps.
Related Issues (20)
- add radial advection problem HOT 1
- doc: verify commands for building the mesh for each problem
- 1D Burgers problem
- Documentation 2D Shallow water HOT 1
- burgers2d: erro in doc HOT 1
- 1D Burgers' equation HOT 6
- possible typo in documentation HOT 1
- ci too slow in debug mode: need to make tests smaller
- test auto add
- schwarz: semantics of problem classes
- schwarz: make ccu mesh default constructible, assignable and movable
- schwarz: `bool isFullyPeriodic()` should not be a method inside the ccu mesh class
- schwarz: expose nested alias for the connectivity graph
- schwarz: support custom BCs in a practical, less intrusive way
- cleanup headers includes
- make `numDofPerCell` NOT a template in the problem classes
- add method to api adapter to query numDofPerCells
- add method to set BCFunctor internal pointer HOT 2
- add dummy icIdentifier to Swe2d create_problem_eigen for consistent Schwarz subdomain initialization
- add first problem for c++ to doc
Recommend Projects
-
React
A declarative, efficient, and flexible JavaScript library for building user interfaces.
-
Vue.js
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
-
Typescript
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
-
TensorFlow
An Open Source Machine Learning Framework for Everyone
-
Django
The Web framework for perfectionists with deadlines.
-
Laravel
A PHP framework for web artisans
-
D3
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
-
Recommend Topics
-
javascript
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
-
web
Some thing interesting about web. New door for the world.
-
server
A server is a program made to process requests and deliver data to clients.
-
Machine learning
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
-
Visualization
Some thing interesting about visualization, use data art
-
Game
Some thing interesting about game, make everyone happy.
Recommend Org
-
Facebook
We are working to build community through open source technology. NB: members must have two-factor auth.
-
Microsoft
Open source projects and samples from Microsoft.
-
Google
Google ❤️ Open Source for everyone.
-
Alibaba
Alibaba Open Source for everyone
-
D3
Data-Driven Documents codes.
-
Tencent
China tencent open source team.
from pressio-demoapps.