Nektar++
CompressibleSolver.cpp
Go to the documentation of this file.
1///////////////////////////////////////////////////////////////////////////////
2//
3// File: CompressibleSolver.cpp
4//
5// For more information, please see: http://www.nektar.info
6//
7// The MIT License
8//
9// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10// Department of Aeronautics, Imperial College London (UK), and Scientific
11// Computing and Imaging Institute, University of Utah (USA).
12//
13// Permission is hereby granted, free of charge, to any person obtaining a
14// copy of this software and associated documentation files (the "Software"),
15// to deal in the Software without restriction, including without limitation
16// the rights to use, copy, modify, merge, publish, distribute, sublicense,
17// and/or sell copies of the Software, and to permit persons to whom the
18// Software is furnished to do so, subject to the following conditions:
19//
20// The above copyright notice and this permission notice shall be included
21// in all copies or substantial portions of the Software.
22//
23// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29// DEALINGS IN THE SOFTWARE.
30//
31// Description: Compressible Riemann solver.
32//
33///////////////////////////////////////////////////////////////////////////////
34
35#include "CompressibleSolver.h"
36#include <boost/algorithm/string/predicate.hpp>
37
38namespace Nektar
39{
40
43 : RiemannSolver(pSession), m_pointSolve(true)
44{
45 m_requiresRotation = true;
46
47 // Create equation of state object
48 std::string eosType;
49 pSession->LoadSolverInfo("EquationOfState", eosType, "IdealGas");
50 m_eos = GetEquationOfStateFactory().CreateInstance(eosType, pSession);
51 // Check if using ideal gas
52 m_idealGas = boost::iequals(eosType, "IdealGas");
53}
54
56 : RiemannSolver(), m_pointSolve(true), m_idealGas(true)
57{
58 m_requiresRotation = true;
59}
60
61/**
62 *
63 */
65 const int nDim, const Array<OneD, const Array<OneD, NekDouble>> &Fwd,
66 const Array<OneD, const Array<OneD, NekDouble>> &Bwd,
68{
69 if (m_pointSolve)
70 {
71 size_t expDim = nDim;
72 NekDouble rhouf{}, rhovf{};
73
74 if (expDim == 1)
75 {
76 for (size_t i = 0; i < Fwd[0].size(); ++i)
77 {
78 v_PointSolve(Fwd[0][i], Fwd[1][i], 0.0, 0.0, Fwd[2][i],
79 Bwd[0][i], Bwd[1][i], 0.0, 0.0, Bwd[2][i],
80 flux[0][i], flux[1][i], rhouf, rhovf, flux[2][i]);
81 }
82 }
83 else if (expDim == 2)
84 {
85 for (size_t i = 0; i < Fwd[0].size(); ++i)
86 {
87 v_PointSolve(Fwd[0][i], Fwd[1][i], Fwd[2][i], 0.0, Fwd[3][i],
88 Bwd[0][i], Bwd[1][i], Bwd[2][i], 0.0, Bwd[3][i],
89 flux[0][i], flux[1][i], flux[2][i], rhovf,
90 flux[3][i]);
91 }
92 }
93 else if (expDim == 3)
94 {
95 for (size_t i = 0; i < Fwd[0].size(); ++i)
96 {
97 v_PointSolve(Fwd[0][i], Fwd[1][i], Fwd[2][i], Fwd[3][i],
98 Fwd[4][i], Bwd[0][i], Bwd[1][i], Bwd[2][i],
99 Bwd[3][i], Bwd[4][i], flux[0][i], flux[1][i],
100 flux[2][i], flux[3][i], flux[4][i]);
101 }
102 }
103 }
104 else
105 {
106 v_ArraySolve(Fwd, Bwd, flux);
107 }
108}
109
110/**
111 *
112 */
114 NekDouble rhoL, NekDouble pL, NekDouble eL, [[maybe_unused]] NekDouble HL,
115 [[maybe_unused]] NekDouble srL, NekDouble rhoR, NekDouble pR, NekDouble eR,
116 [[maybe_unused]] NekDouble HR, [[maybe_unused]] NekDouble srR,
117 NekDouble HRoe, NekDouble URoe2, [[maybe_unused]] NekDouble srLR)
118{
119 static NekDouble gamma = m_params["gamma"]();
120 NekDouble cRoe;
121 if (m_idealGas)
122 {
123 cRoe = sqrt((gamma - 1.0) * (HRoe - 0.5 * URoe2));
124 }
125 else
126 {
127 // Calculate static enthalpy of left and right states
128 NekDouble hL = eL + pL / rhoL;
129 NekDouble hR = eR + pR / rhoR;
130
131 // Get partial derivatives of P(rho,e)
132 NekDouble dpdeL = m_eos->GetDPDe_rho(rhoL, eL);
133 NekDouble dpdeR = m_eos->GetDPDe_rho(rhoR, eR);
134 NekDouble dpdrhoL = m_eos->GetDPDrho_e(rhoL, eL);
135 NekDouble dpdrhoR = m_eos->GetDPDrho_e(rhoR, eR);
136
137 // Evaluate chi and kappa parameters
138 NekDouble chiL = dpdrhoL - eL / rhoL * dpdeL;
139 NekDouble kappaL = dpdeL / rhoL;
140 NekDouble chiR = dpdrhoR - eR / rhoR * dpdeR;
141 NekDouble kappaR = dpdeR / rhoR;
142
143 //
144 // Calculate interface speed of sound using procedure from
145 // Vinokur, M.; Montagné, J.-L. "Generalized Flux-Vector
146 // Splitting and Roe Average for an Equilibrium Real Gas",
147 // JCP (1990).
148 //
149
150 // Calculate averages
151 NekDouble avgChi = 0.5 * (chiL + chiR);
152 NekDouble avgKappa = 0.5 * (kappaL + kappaR);
153 NekDouble avgKappaH = 0.5 * (kappaL * hL + kappaR * hR);
154
155 // Calculate jumps
156 NekDouble deltaP = pR - pL;
157 NekDouble deltaRho = rhoR - rhoL;
158 NekDouble deltaRhoe = rhoR * eR - rhoL * eL;
159
160 // Evaluate dP: equation (64) from Vinokur-Montagné
161 NekDouble dP = deltaP - avgChi * deltaRho - avgKappa * deltaRhoe;
162 // s (eq 66)
163 NekDouble s = avgChi + avgKappaH;
164 // D (eq 65)
165 NekDouble D = (s * deltaRho) * (s * deltaRho) + deltaP * deltaP;
166 // chiRoe and kappaRoe (eq 66)
167 NekDouble chiRoe, kappaRoe;
168 NekDouble fac = D - deltaP * deltaRho;
170 {
171 chiRoe = (D * avgChi + s * s * deltaRho * dP) / fac;
172 kappaRoe = D * avgKappa / fac;
173 }
174 else
175 {
176 chiRoe = avgChi;
177 kappaRoe = avgKappa;
178 }
179 // Speed of sound (eq 53)
180 cRoe = sqrt(chiRoe + kappaRoe * (HRoe - 0.5 * URoe2));
181 }
182 return cRoe;
183}
184
185} // namespace Nektar
virtual void v_PointSolve(ND rhoL, ND rhouL, ND rhovL, ND rhowL, ND EL, ND rhoR, ND rhouR, ND rhovR, ND rhowR, ND ER, ND &rhof, ND &rhouf, ND &rhovf, ND &rhowf, ND &Ef)
ND GetRoeSoundSpeed(ND rhoL, ND pL, ND eL, ND HL, ND srL, ND rhoR, ND pR, ND eR, ND HR, ND srR, ND HRoe, ND URoe2, ND srLR)
CompressibleSolver()
Programmatic ctor.
EquationOfStateSharedPtr m_eos
virtual void v_ArraySolve(const Array< OneD, const Array< OneD, ND > > &Fwd, const Array< OneD, const Array< OneD, ND > > &Bwd, Array< OneD, Array< OneD, ND > > &flux)
void v_Solve(const int nDim, const Array< OneD, const Array< OneD, ND > > &Fwd, const Array< OneD, const Array< OneD, ND > > &Bwd, Array< OneD, Array< OneD, ND > > &flux) override
tBaseSharedPtr CreateInstance(tKey idKey, tParam... args)
Create an instance of the class referred to by idKey.
The RiemannSolver class provides an abstract interface under which solvers for various Riemann proble...
Definition: RiemannSolver.h:58
bool m_requiresRotation
Indicates whether the Riemann solver requires a rotation to be applied to the velocity fields.
std::map< std::string, RSParamFuncType > m_params
Map of parameter function types.
std::shared_ptr< SessionReader > SessionReaderSharedPtr
static const NekDouble kNekZeroTol
EquationOfStateFactory & GetEquationOfStateFactory()
Declaration of the equation of state factory singleton.
double NekDouble
scalarT< T > abs(scalarT< T > in)
Definition: scalar.hpp:289
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:285