ExactSolverToro.cpp

Go to the documentation of this file.

///////////////////////////////////////////////////////////////////////////////
//
// File: ExactSolverToro.cpp
//
// For more information, please see: http://www.nektar.info
//
// The MIT License
//
// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
// Department of Aeronautics, Imperial College London (UK), and Scientific
// Computing and Imaging Institute, University of Utah (USA).
//
// Permission is hereby granted, free of charge, to any person obtaining a
// copy of this software and associated documentation files (the "Software"),
// to deal in the Software without restriction, including without limitation
// the rights to use, copy, modify, merge, publish, distribute, sublicense,
// and/or sell copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included
// in all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
// DEALINGS IN THE SOFTWARE.
//
// Description: Exact Riemann solver (Toro 2009).
//
///////////////////////////////////////////////////////////////////////////////
 
#include <CompressibleFlowSolver/RiemannSolvers/ExactSolverToro.h>
 
#define TOL 1e-6
#define NRITER 20
 
namespace Nektar
{
std::string ExactSolverToro::solverName =
    SolverUtils::GetRiemannSolverFactory().RegisterCreatorFunction(
        "ExactToro", ExactSolverToro::create, "Exact Riemann solver");
 
ExactSolverToro::ExactSolverToro(
    const LibUtilities::SessionReaderSharedPtr &pSession)
    : CompressibleSolver(pSession)
{
}
 
/**
 * @brief Use either PVRS, two-rarefaction or two-shock Riemann solvers to
 * calculate an initial pressure for the Newton-Raphson scheme.
 *
 * @param g      Array of calculated gamma values.
 * @param rhoL   Density left state.
 * @param rhoR   Density right state.
 * @param uL     x-velocity component left state.
 * @param uR     x-velocity component right state.
 * @param pL     Pressure component left state.
 * @param pR     Pressure component right state.
 * @param cL     Sound speed component left state.
 * @param cR     Sound speed component right state.
 * @return       Computed initial guess for the Newton-Raphson scheme.
 */
inline NekDouble guessp(NekDouble g[], NekDouble rhoL, NekDouble uL,
                        NekDouble pL, NekDouble cL, NekDouble rhoR,
                        NekDouble uR, NekDouble pR, NekDouble cR)
{
    const NekDouble quser = 2.0;
    NekDouble cup, ppv, pmin, pmax, qmax;
 
    cup  = 0.25 * (rhoL + rhoR) * (cL + cR);
    ppv  = 0.5 * (pL + pR) + 0.5 * (uL - uR) * cup;
    ppv  = std::max(0.0, ppv);
    pmin = std::min(pL, pR);
    pmax = std::max(pL, pR);
    qmax = pmax / pmin;
 
    if (qmax <= quser && pmin <= ppv && ppv <= pmax)
    {
        // Select PVRS Riemann solver.
        return ppv;
    }
    else if (ppv < pmin)
    {
        // Select two-rarefaction Riemann solver.
        NekDouble pq = pow(pL / pR, g[1]);
        NekDouble um =
            (pq * uL / cL + uR / cR + g[4] * (pq - 1.0)) / (pq / cL + 1.0 / cR);
        NekDouble ptL = 1.0 + g[7] * (uL - um) / cL;
        NekDouble ptR = 1.0 + g[7] * (um - uR) / cR;
        return 0.5 * (pL * pow(ptL, g[3]) + pR * pow(ptR, g[3]));
    }
    else
    {
        // Select two-shock Riemann solver with PVRS as estimate.
        NekDouble geL = sqrt((g[5] / rhoL) / (g[6] * pL + ppv));
        NekDouble geR = sqrt((g[5] / rhoR) / (g[6] * pR + ppv));
        return (geL * pL + geR * pR - (uR - uL)) / (geL + geR);
    }
}
 
/**
 * @brief Evaluate pressure functions fL and fR in Newton iteration of
 * Riemann solver (see equation 4.85 and 4.86 of Toro 2009).
 *
 * @param g   Array of gamma parameters.
 * @param p   Pressure at current iteration.
 * @param dk  Density (left or right state)
 * @param pk  Pressure (left or right state)
 * @param ck  Sound speed (left or right state)
 * @param f   Computed pressure (shock).
 * @param fd  Computed pressure (rarefaction).
 */
inline void prefun(NekDouble *g, NekDouble p, NekDouble dk, NekDouble pk,
                   NekDouble ck, NekDouble &f, NekDouble &fd)
{
    if (p <= pk)
    {
        // Rarefaction wave
        NekDouble prat = p / pk;
        f              = g[4] * ck * (pow(prat, g[1]) - 1.0);
        fd             = pow(prat, -g[2]) / (dk * ck);
    }
    else
    {
        // Shock wave
        NekDouble ak  = g[5] / dk;
        NekDouble bk  = g[6] * pk;
        NekDouble qrt = sqrt(ak / (bk + p));
        f             = (p - pk) * qrt;
        fd            = (1.0 - 0.5 * (p - pk) / (bk + p)) * qrt;
    }
}
 
/**
 * @brief Exact Riemann solver for the Euler equations.
 *
 * This algorithm is transcribed from:
 *
 *   "Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical
 *   Introduction", E. F. Toro (3rd edition, 2009).
 *
 * The full Fortran 77 routine can be found at the end of chapter 4 (section
 * 4.9). This transcription is essentially the functions STARPU and SAMPLE
 * glued together, and variable names are kept mostly the same. See the
 * preceding chapter which explains the derivation of the solver. The
 * routines PREFUN and GUESSP are kept separate and are reproduced above.
 *
 * @param rhoL      Density left state.
 * @param rhoR      Density right state.
 * @param rhouL     x-momentum component left state.
 * @param rhouR     x-momentum component right state.
 * @param rhovL     y-momentum component left state.
 * @param rhovR     y-momentum component right state.
 * @param rhowL     z-momentum component left state.
 * @param rhowR     z-momentum component right state.
 * @param EL        Energy left state.
 * @param ER        Energy right state.
 * @param rhof      Computed Riemann flux for density.
 * @param rhouf     Computed Riemann flux for x-momentum component
 * @param rhovf     Computed Riemann flux for y-momentum component
 * @param rhowf     Computed Riemann flux for z-momentum component
 * @param Ef        Computed Riemann flux for energy.
 */
void ExactSolverToro::v_PointSolve(
    NekDouble rhoL, NekDouble rhouL, NekDouble rhovL, NekDouble rhowL,
    NekDouble EL, NekDouble rhoR, NekDouble rhouR, NekDouble rhovR,
    NekDouble rhowR, NekDouble ER, NekDouble &rhof, NekDouble &rhouf,
    NekDouble &rhovf, NekDouble &rhowf, NekDouble &Ef)
{
    static NekDouble gamma = m_params["gamma"]();
 
    // Left and right variables.
    NekDouble uL = rhouL / rhoL;
    NekDouble vL = rhovL / rhoL;
    NekDouble wL = rhowL / rhoL;
    NekDouble uR = rhouR / rhoR;
    NekDouble vR = rhovR / rhoR;
    NekDouble wR = rhowR / rhoR;
 
    // Left and right pressure.
    NekDouble pL =
        (gamma - 1.0) * (EL - 0.5 * (rhouL * uL + rhovL * vL + rhowL * wL));
    NekDouble pR =
        (gamma - 1.0) * (ER - 0.5 * (rhouR * uR + rhovR * vR + rhowR * wR));
 
    // Compute gammas.
    NekDouble g[] = {gamma,
                     (gamma - 1.0) / (2.0 * gamma),
                     (gamma + 1.0) / (2.0 * gamma),
                     2.0 * gamma / (gamma - 1.0),
                     2.0 / (gamma - 1.0),
                     2.0 / (gamma + 1.0),
                     (gamma - 1.0) / (gamma + 1.0),
                     0.5 * (gamma - 1.0),
                     gamma - 1.0};
 
    // Compute sound speeds.
    NekDouble cL = sqrt(gamma * pL / rhoL);
    NekDouble cR = sqrt(gamma * pR / rhoR);
 
    ASSERTL0(g[4] * (cL + cR) > (uR - uL),
             "Vacuum is generated by given data.");
 
    // Guess initial pressure.
    NekDouble pOld  = guessp(g, rhoL, uL, pL, cL, rhoR, uR, pR, cR);
    NekDouble uDiff = uR - uL;
    NekDouble p, fL, fR, fLd, fRd, change;
    int k;
 
    // Newton-Raphson iteration for pressure in star region.
    for (k = 0; k < NRITER; ++k)
    {
        prefun(g, pOld, rhoL, pL, cL, fL, fLd);
        prefun(g, pOld, rhoR, pR, cR, fR, fRd);
        p      = pOld - (fL + fR + uDiff) / (fLd + fRd);
        change = 2 * fabs((p - pOld) / (p + pOld));
 
        if (change <= TOL)
        {
            break;
        }
 
        if (p < 0.0)
        {
            p = TOL;
        }
 
        pOld = p;
    }
 
    ASSERTL0(k < NRITER, "Divergence in Newton-Raphson scheme");
 
    // Compute velocity in star region.
    NekDouble u = 0.5 * (uL + uR + fR - fL);
 
    // -- SAMPLE ROUTINE --
    // The variable S aligns with the Fortran parameter S of the SAMPLE
    // routine, but is hard-coded as 0 (and should be optimised out by the
    // compiler). Since we are using a Godunov scheme we pick this as 0 (see
    // chapter 6 of Toro 2009).
    const NekDouble S = 0.0;
 
    // Computed primitive variables.
    NekDouble outRho, outU, outV, outW, outP;
 
    if (S <= u)
    {
        if (p <= pL)
        {
            // Left rarefaction
            NekDouble shL = uL - cL;
            if (S <= shL)
            {
                // Sampled point is left data state.
                outRho = rhoL;
                outU   = uL;
                outV   = vL;
                outW   = wL;
                outP   = pL;
            }
            else
            {
                NekDouble cmL = cL * pow(p / pL, g[1]);
                NekDouble stL = u - cmL;
 
                if (S > stL)
                {
                    // Sampled point is star left state
                    outRho = rhoL * pow(p / pL, 1.0 / gamma);
                    outU   = u;
                    outV   = vL;
                    outW   = wL;
                    outP   = p;
                }
                else
                {
                    // Sampled point is inside left fan
                    NekDouble c = g[5] * (cL + g[7] * (uL - S));
                    outRho      = rhoL * pow(c / cL, g[4]);
                    outU        = g[5] * (cL + g[7] * uL + S);
                    outV        = vL;
                    outW        = wL;
                    outP        = pL * pow(c / cL, g[3]);
                }
            }
        }
        else
        {
            // Left shock
            NekDouble pmL = p / pL;
            NekDouble SL  = uL - cL * sqrt(g[2] * pmL + g[1]);
            if (S <= SL)
            {
                // Sampled point is left data state
                outRho = rhoL;
                outU   = uL;
                outV   = vL;
                outW   = wL;
                outP   = pL;
            }
            else
            {
                // Sampled point is star left state
                outRho = rhoL * (pmL + g[6]) / (pmL * g[6] + 1.0);
                outU   = u;
                outV   = vL;
                outW   = wL;
                outP   = p;
            }
        }
    }
    else
    {
        if (p > pR)
        {
            // Right shock
            NekDouble pmR = p / pR;
            NekDouble SR  = uR + cR * sqrt(g[2] * pmR + g[1]);
            if (S >= SR)
            {
                // Sampled point is right data state
                outRho = rhoR;
                outU   = uR;
                outV   = vR;
                outW   = wR;
                outP   = pR;
            }
            else
            {
                // Sampled point is star right state
                outRho = rhoR * (pmR + g[6]) / (pmR * g[6] + 1.0);
                outU   = u;
                outV   = vR;
                outW   = wR;
                outP   = p;
            }
        }
        else
        {
            // Right rarefaction
            NekDouble shR = uR + cR;
 
            if (S >= shR)
            {
                // Sampled point is right data state
                outRho = rhoR;
                outU   = uR;
                outV   = vR;
                outW   = wR;
                outP   = pR;
            }
            else
            {
                NekDouble cmR = cR * pow(p / pR, g[1]);
                NekDouble stR = u + cmR;
 
                if (S <= stR)
                {
                    // Sampled point is star right state
                    outRho = rhoR * pow(p / pR, 1.0 / gamma);
                    outU   = u;
                    outV   = vR;
                    outW   = wR;
                    outP   = p;
                }
                else
                {
                    // Sampled point is inside left fan
                    NekDouble c = g[5] * (cR - g[7] * (uR - S));
                    outRho      = rhoR * pow(c / cR, g[4]);
                    outU        = g[5] * (-cR + g[7] * uR + S);
                    outV        = vR;
                    outW        = wR;
                    outP        = pR * pow(c / cR, g[3]);
                }
            }
        }
    }
 
    // Transform computed primitive variables to fluxes.
    rhof  = outRho * outU;
    rhouf = outP + outRho * outU * outU;
    rhovf = outRho * outU * outV;
    rhowf = outRho * outU * outW;
    Ef    = outU *
         (outP / (gamma - 1.0) +
          0.5 * outRho * (outU * outU + outV * outV + outW * outW) + outP);
}
} // namespace Nektar