#include <ExactSolverToro.h>

Inheritance diagram for Nektar::ExactSolverToro:

Collaboration diagram for Nektar::ExactSolverToro:

Static Public Member Functions
static RiemannSolverSharedPtr	create ()

Static Public Attributes
static std::string	solverName

Protected Member Functions
	ExactSolverToro ()
virtual void	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)
	Exact Riemann solver for the Euler equations.
Protected Member Functions inherited from Nektar::CompressibleSolver
	CompressibleSolver ()
virtual void	v_Solve (const int nDim, const Array< OneD, const Array< OneD, NekDouble > > &Fwd, const Array< OneD, const Array< OneD, NekDouble > > &Bwd, Array< OneD, Array< OneD, NekDouble > > &flux)
virtual void	v_ArraySolve (const Array< OneD, const Array< OneD, NekDouble > > &Fwd, const Array< OneD, const Array< OneD, NekDouble > > &Bwd, Array< OneD, Array< OneD, NekDouble > > &flux)
virtual void	v_PointSolveVisc (NekDouble rhoL, NekDouble rhouL, NekDouble rhovL, NekDouble rhowL, NekDouble EL, NekDouble EpsL, NekDouble rhoR, NekDouble rhouR, NekDouble rhovR, NekDouble rhowR, NekDouble ER, NekDouble EpsR, NekDouble &rhof, NekDouble &rhouf, NekDouble &rhovf, NekDouble &rhowf, NekDouble &Ef, NekDouble &Epsf)
Protected Member Functions inherited from Nektar::SolverUtils::RiemannSolver
SOLVER_UTILS_EXPORT	RiemannSolver ()
void	GenerateRotationMatrices (const Array< OneD, const Array< OneD, NekDouble > > &normals)
	Generate rotation matrices for 3D expansions.
void	FromToRotation (Array< OneD, const NekDouble > &from, Array< OneD, const NekDouble > &to, NekDouble *mat)
	A function for creating a rotation matrix that rotates a vector from into another vector to.
SOLVER_UTILS_EXPORT void	rotateToNormal (const Array< OneD, const Array< OneD, NekDouble > > &inarray, const Array< OneD, const Array< OneD, NekDouble > > &normals, const Array< OneD, const Array< OneD, NekDouble > > &vecLocs, Array< OneD, Array< OneD, NekDouble > > &outarray)
	Rotate a vector field to trace normal.
SOLVER_UTILS_EXPORT void	rotateFromNormal (const Array< OneD, const Array< OneD, NekDouble > > &inarray, const Array< OneD, const Array< OneD, NekDouble > > &normals, const Array< OneD, const Array< OneD, NekDouble > > &vecLocs, Array< OneD, Array< OneD, NekDouble > > &outarray)
	Rotate a vector field from trace normal.
SOLVER_UTILS_EXPORT bool	CheckScalars (std::string name)
	Determine whether a scalar has been defined in m_scalars.
SOLVER_UTILS_EXPORT bool	CheckVectors (std::string name)
	Determine whether a vector has been defined in m_vectors.
SOLVER_UTILS_EXPORT bool	CheckParams (std::string name)
	Determine whether a parameter has been defined in m_params.
SOLVER_UTILS_EXPORT bool	CheckAuxScal (std::string name)
	Determine whether a scalar has been defined in m_auxScal.
SOLVER_UTILS_EXPORT bool	CheckAuxVec (std::string name)
	Determine whether a vector has been defined in m_auxVec.

Additional Inherited Members
Public Member Functions inherited from Nektar::SolverUtils::RiemannSolver
SOLVER_UTILS_EXPORT void	Solve (const int nDim, const Array< OneD, const Array< OneD, NekDouble > > &Fwd, const Array< OneD, const Array< OneD, NekDouble > > &Bwd, Array< OneD, Array< OneD, NekDouble > > &flux)
	Perform the Riemann solve given the forwards and backwards spaces.
template<typename FuncPointerT , typename ObjectPointerT >
void	SetScalar (std::string name, FuncPointerT func, ObjectPointerT obj)
void	SetScalar (std::string name, RSScalarFuncType fp)
template<typename FuncPointerT , typename ObjectPointerT >
void	SetVector (std::string name, FuncPointerT func, ObjectPointerT obj)
void	SetVector (std::string name, RSVecFuncType fp)
template<typename FuncPointerT , typename ObjectPointerT >
void	SetParam (std::string name, FuncPointerT func, ObjectPointerT obj)
void	SetParam (std::string name, RSParamFuncType fp)
template<typename FuncPointerT , typename ObjectPointerT >
void	SetAuxScal (std::string name, FuncPointerT func, ObjectPointerT obj)
template<typename FuncPointerT , typename ObjectPointerT >
void	SetAuxVec (std::string name, FuncPointerT func, ObjectPointerT obj)
std::map< std::string, RSScalarFuncType > &	GetScalars ()
std::map< std::string, RSVecFuncType > &	GetVectors ()
std::map< std::string, RSParamFuncType > &	GetParams ()
Public Attributes inherited from Nektar::SolverUtils::RiemannSolver
int	m_spacedim
Protected Attributes inherited from Nektar::CompressibleSolver
bool	m_pointSolve

Detailed Description

Definition at line 43 of file ExactSolverToro.h.

Constructor & Destructor Documentation

Nektar::ExactSolverToro::ExactSolverToro ( )

protected

Definition at line 49 of file ExactSolverToro.cpp.

Referenced by create().

                                     : CompressibleSolver()
    {
    }

Member Function Documentation

static RiemannSolverSharedPtr Nektar::ExactSolverToro::create ( )

inlinestatic

Definition at line 46 of file ExactSolverToro.h.

References ExactSolverToro().

        {
            return RiemannSolverSharedPtr(
                new ExactSolverToro());
        }

void Nektar::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
	)

protectedvirtual

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.

Parameters

rhoL	Density left state.
rhoR	Density right state.
rhouL	x-momentum component left state.
rhouR	x-momentum component right state.
rhovL	y-momentum component left state.
rhovR	y-momentum component right state.
rhowL	z-momentum component left state.
rhowR	z-momentum component right state.
EL	Energy left state.
ER	Energy right state.
rhof	Computed Riemann flux for density.
rhouf	Computed Riemann flux for x-momentum component
rhovf	Computed Riemann flux for y-momentum component
rhowf	Computed Riemann flux for z-momentum component
Ef	Computed Riemann flux for energy.

Reimplemented from Nektar::CompressibleSolver.

Definition at line 170 of file ExactSolverToro.cpp.

References ASSERTL0, Nektar::guessp(), Nektar::SolverUtils::RiemannSolver::m_params, NRITER, Nektar::prefun(), and TOL.

    {
        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);
    }

Member Data Documentation

std::string Nektar::ExactSolverToro::solverName

static

Initial value:

        SolverUtils::GetRiemannSolverFactory().RegisterCreatorFunction(
            "ExactToro",
            ExactSolverToro::create,
            "Exact Riemann solver")

Definition at line 52 of file ExactSolverToro.h.

Static Public Member Functions

Static Public Attributes

Protected Member Functions

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation