APE: Retain original APE4 formulation by Ewert and Schröder

The equations that were implemented before had the speed of sound moved into the
divergence term of the p-equation. This implicates an additional source term
that depends on the gradient of c, which was not present. Implementing this
source term has proven to produce instabilities at the element boundaries due to
the unsteady gradients. Accordingly, this implementation modifies the flux
formulation and retains Ewert and Schröders original formulation of the APE 1/4
equations LHS.

This fixes the behavior of the equations at e.g. density or pressure gradients.

solvers/APESolver/EquationSystems/APE.cpp

@@ -175,20 +175,25 @@ void APE::GetFluxVector(
 
         if (i == 0)
         {
-            // F_{adv,p',j} = \gamma p_0 u'_j + p' \bar{u}_j
+            // F_{adv,p',j} = \rho_0 u'_j + p' \bar{u}_j / c^2
             for (int j = 0; j < m_spacedim; ++j)
             {
                 Vmath::Zero(nq, flux[i][j], 1);
 
-                // construct \gamma p_0 u'_j term
-                Vmath::Smul(nq, m_gamma, m_basefield[0], 1, tmp1, 1);
-                Vmath::Vmul(nq, tmp1, 1, physfield[j+1], 1, tmp1, 1);
+                // construct rho_0 u'_j term
+                Vmath::Vmul(nq, m_basefield[1], 1, physfield[j+1], 1, flux[i][j], 1);
 
-                // construct p' \bar{u}_j term
+                // construct p' \bar{u}_j / c^2 term
+                // c^2
+                Vmath::Vdiv(nq, m_basefield[0], 1, m_basefield[1], 1, tmp1, 1);
+                Vmath::Smul(nq, m_gamma, tmp1, 1, tmp1, 1);
+
+                // p' \bar{u}_j / c^2 term
                 Vmath::Vmul(nq, physfield[0], 1, m_basefield[j+2], 1, tmp2, 1);
+                Vmath::Vdiv(nq, tmp2, 1, tmp1, 1, tmp2, 1);
 
-                // add both terms
-                Vmath::Vadd(nq, tmp1, 1, tmp2, 1, flux[i][j], 1);
+                // \rho_0 u'_j + p' \bar{u}_j / c^2
+                Vmath::Vadd(nq, flux[i][j], 1, tmp2, 1, flux[i][j], 1);
             }
         }
         else
@@ -207,8 +212,7 @@ void APE::GetFluxVector(
                     Vmath::Zero(nq, tmp1, 1);
                     for (int k = 0; k < m_spacedim; ++k)
                     {
-                        Vmath::Vmul(nq, physfield[k+1], 1, m_basefield[k+2], 1, tmp2, 1);
-                        Vmath::Vadd(nq, tmp1, 1, tmp2, 1, tmp1, 1);
+                        Vmath::Vvtvp(nq, physfield[k+1], 1, m_basefield[k+2], 1, tmp1, 1, tmp1, 1);
                     }
 
                     // add terms
@@ -228,14 +232,23 @@ void APE::DoOdeRhs(const Array<OneD, const Array<OneD, NekDouble> >&inarray,
                    const NekDouble time)
 {
     int nVariables = inarray.num_elements();
-    int nq         = GetTotPoints();
+    int nq = GetTotPoints();
+    Array<OneD, NekDouble> tmp1(nq);
 
     // WeakDG does not use advVel, so we only provide a dummy array
     Array<OneD, Array<OneD, NekDouble> > advVel(m_spacedim);
     m_advection->Advect(nVariables, m_fields, advVel, inarray, outarray, time);
 
     for (int i = 0; i < nVariables; ++i)
     {
+        if (i ==  0)
+        {
+            // c^2 = gamma*p0/rho0
+            Vmath::Vdiv(nq, m_basefield[0], 1, m_basefield[1], 1, tmp1, 1);
+            Vmath::Smul(nq, m_gamma, tmp1, 1, tmp1, 1);
+            Vmath::Vmul(nq, tmp1, 1, outarray[i], 1, outarray[i], 1);
+        }
+
         Vmath::Neg(nq, outarray[i], 1);
     }
 

solvers/APESolver/RiemannSolvers/LaxFriedrichsSolver.cpp

@@ -95,18 +95,18 @@ void LaxFriedrichsSolver::v_PointSolve(
     // max absolute eigenvalue of the jacobian of F_n1
     NekDouble a_1_max = std::max(std::abs(u0 - c), std::abs(u0 + c));
 
-    NekDouble pFL = gamma*p0*uL + pL*u0;
+    NekDouble pFL = rho0*uL + u0*pL/(c*c);
     NekDouble uFL = pL/rho0 + u0*uL + v0*vL + w0*wL;
     NekDouble vFL = 0;
     NekDouble wFL = 0;
 
-    NekDouble pFR = gamma*p0*uR + pR*u0;
+    NekDouble pFR = rho0*uR + u0*pR/(c*c);
     NekDouble uFR = pR/rho0 + u0*uR + v0*vR + w0*wR;
     NekDouble vFR = 0;
     NekDouble wFR = 0;
 
     // assemble the face-normal fluxes
-    pF = 0.5*(pFL + pFR - a_1_max*(pR - pL));
+    pF = 0.5*(pFL + pFR - a_1_max*(pR/(c*c) - pL/(c*c)));
     uF = 0.5*(uFL + uFR - a_1_max*(uR - uL));
     vF = 0.5*(vFL + vFR - a_1_max*(vR - vL));
     wF = 0.5*(wFL + wFR - a_1_max*(wR - wL));

solvers/APESolver/RiemannSolvers/UpwindSolver.cpp

@@ -86,21 +86,24 @@ void UpwindSolver::v_PointSolve(
     ASSERTL1(CheckParams("Gamma"), "Gamma not defined.");
     const NekDouble &gamma = m_params["Gamma"]();
 
+    // Speed of sound
+    NekDouble c = sqrt(gamma * p0 / rho0);
+
     Array<OneD, NekDouble> characteristic(4);
     Array<OneD, NekDouble> W(2);
     Array<OneD, NekDouble> lambda(2);
 
     // compute the wave speeds
-    lambda[0] = u0 + sqrt(p0*gamma*rho0)/rho0;
-    lambda[1] = u0 - sqrt(p0*gamma*rho0)/rho0;
+    lambda[0] = u0 + c;
+    lambda[1] = u0 - c;
 
     // calculate the caracteristic variables
     //left characteristics
-    characteristic[0] = pL/2 + uL*sqrt(p0*gamma*rho0)/2;
-    characteristic[1] = pL/2 - uL*sqrt(p0*gamma*rho0)/2;
+    characteristic[0] = pL/2 + uL*c*rho0/2;
+    characteristic[1] = pL/2 - uL*c*rho0/2;
     //right characteristics
-    characteristic[2] = pR/2 + uR*sqrt(p0*gamma*rho0)/2;
-    characteristic[3] = pR/2 - uR*sqrt(p0*gamma*rho0)/2;
+    characteristic[2] = pR/2 + uR*c*rho0/2;
+    characteristic[3] = pR/2 - uR*c*rho0/2;
 
     //take left or right value of characteristic variable
     for (int j = 0; j < 2; j++)
@@ -117,15 +120,10 @@ void UpwindSolver::v_PointSolve(
 
     //calculate conservative variables from characteristics
     NekDouble p = W[0] + W[1];
-    NekDouble u = (W[0] - W[1])/sqrt(p0*gamma*rho0);
-    // do not divide by zero
-    if (p0*gamma*rho0 == 0)
-    {
-        u = 0.0;
-    }
+    NekDouble u = (W[0] - W[1])/(c*rho0);
 
     // assemble the fluxes
-    pF = gamma*p0*u + u0*p;
+    pF = rho0*u + u0*p/(c*c);
     uF = p/rho0 + u0*u + v0*vL + w0*wL;
     vF = 0.0;
     wF = 0.0;