Nektar++
IsentropicVortex.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File IsentropicVortex.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: Euler equations for isentropic vortex
32 //
33 ///////////////////////////////////////////////////////////////////////////////
34 
35 #include <boost/core/ignore_unused.hpp>
36 
38 
39 using namespace std;
40 
41 namespace Nektar
42 {
43 string IsentropicVortex::className =
45  "IsentropicVortex", IsentropicVortex::create,
46  "Euler equations for the isentropic vortex test case.");
47 
48 IsentropicVortex::IsentropicVortex(
51  : UnsteadySystem(pSession, pGraph), EulerCFE(pSession, pGraph)
52 {
53 }
54 
55 /**
56  * @brief Destructor for EulerCFE class.
57  */
59 {
60 }
61 
62 /**
63  * @brief Print out a summary with some relevant information.
64  */
66 {
68  SolverUtils::AddSummaryItem(s, "Problem Type", "Isentropic Vortex");
69 }
70 
71 /**
72  * @brief Set the initial conditions.
73  */
75  bool dumpInitialConditions,
76  const int domain)
77 {
78  boost::ignore_unused(domain);
79 
80  int nTotQuadPoints = GetTotPoints();
81  Array<OneD, NekDouble> x(nTotQuadPoints);
82  Array<OneD, NekDouble> y(nTotQuadPoints);
83  Array<OneD, NekDouble> z(nTotQuadPoints);
85 
86  m_fields[0]->GetCoords(x, y, z);
87 
88  for (int i = 0; i < m_spacedim + 2; ++i)
89  {
90  u[i] = Array<OneD, NekDouble>(nTotQuadPoints);
91  }
92 
93  EvaluateIsentropicVortex(x, y, z, u, initialtime);
94 
95  // Forward transform to fill the coefficient space
96  for (int i = 0; i < m_fields.size(); ++i)
97  {
98  Vmath::Vcopy(nTotQuadPoints, u[i], 1, m_fields[i]->UpdatePhys(), 1);
99  m_fields[i]->SetPhysState(true);
100  m_fields[i]->FwdTrans(m_fields[i]->GetPhys(),
101  m_fields[i]->UpdateCoeffs());
102  }
103 
104  if (dumpInitialConditions && m_checksteps)
105  {
106  // Dump initial conditions to file
108  }
109 }
110 
111 /**
112  * @brief Get the exact solutions for isentropic vortex and Ringleb
113  * flow problems.
114  */
116  Array<OneD, NekDouble> &outfield,
117  const NekDouble time)
118 {
119  int nTotQuadPoints = GetTotPoints();
120  Array<OneD, NekDouble> x(nTotQuadPoints);
121  Array<OneD, NekDouble> y(nTotQuadPoints);
122  Array<OneD, NekDouble> z(nTotQuadPoints);
124 
125  m_fields[0]->GetCoords(x, y, z);
126 
127  for (int i = 0; i < m_spacedim + 2; ++i)
128  {
129  u[i] = Array<OneD, NekDouble>(nTotQuadPoints);
130  }
131 
132  EvaluateIsentropicVortex(x, y, z, u, time);
133 
134  Vmath::Vcopy(nTotQuadPoints, u[field], 1, outfield, 1);
135 }
136 
140  NekDouble time, const int o)
141 {
142  boost::ignore_unused(z);
143 
144  int nq = x.size();
145 
146  // Flow parameters
147  const NekDouble x0 = 5.0;
148  const NekDouble y0 = 0.0;
149  const NekDouble beta = 5.0;
150  const NekDouble u0 = 1.0;
151  const NekDouble v0 = 0.5;
152  const NekDouble gamma = m_gamma;
153  NekDouble r, xbar, ybar, tmp;
154  NekDouble fac = 1.0 / (16.0 * gamma * M_PI * M_PI);
155 
156  // In 3D zero rhow field.
157  if (m_spacedim == 3)
158  {
159  Vmath::Zero(nq, &u[3][o], 1);
160  }
161 
162  // Fill storage
163  for (int i = 0; i < nq; ++i)
164  {
165  xbar = x[i] - u0 * time - x0;
166  ybar = y[i] - v0 * time - y0;
167  r = sqrt(xbar * xbar + ybar * ybar);
168  tmp = beta * exp(1 - r * r);
169  u[0][i + o] =
170  pow(1.0 - (gamma - 1.0) * tmp * tmp * fac, 1.0 / (gamma - 1.0));
171  u[1][i + o] = u[0][i + o] * (u0 - tmp * ybar / (2 * M_PI));
172  u[2][i + o] = u[0][i + o] * (v0 + tmp * xbar / (2 * M_PI));
173  u[m_spacedim + 1][i + o] =
174  pow(u[0][i + o], gamma) / (gamma - 1.0) +
175  0.5 * (u[1][i + o] * u[1][i + o] + u[2][i + o] * u[2][i + o]) /
176  u[0][i + o];
177  }
178 }
179 
180 } // namespace Nektar
virtual ~IsentropicVortex()
Destructor for EulerCFE class.
void EvaluateIsentropicVortex(const Array< OneD, NekDouble > &x, const Array< OneD, NekDouble > &y, const Array< OneD, NekDouble > &z, Array< OneD, Array< OneD, NekDouble >> &u, NekDouble time, const int o=0)
Isentropic Vortex Test Case.
virtual void v_EvaluateExactSolution(unsigned int field, Array< OneD, NekDouble > &outfield, const NekDouble time=0.0)
Get the exact solutions for isentropic vortex and Ringleb flow problems.
virtual void v_SetInitialConditions(NekDouble initialtime=0.0, bool dumpInitialConditions=true, const int domain=0)
Set the initial conditions.
virtual void v_GenerateSummary(SolverUtils::SummaryList &s)
Print a summary of time stepping parameters.
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:198
int m_spacedim
Spatial dimension (>= expansion dim).
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
SOLVER_UTILS_EXPORT void Checkpoint_Output(const int n)
Write checkpoint file of m_fields.
int m_checksteps
Number of steps between checkpoints.
SOLVER_UTILS_EXPORT int GetTotPoints()
Base class for unsteady solvers.
virtual SOLVER_UTILS_EXPORT void v_GenerateSummary(SummaryList &s)
Print a summary of time stepping parameters.
std::shared_ptr< SessionReader > SessionReaderSharedPtr
@ beta
Gauss Radau pinned at x=-1,.
Definition: PointsType.h:61
std::vector< std::pair< std::string, std::string > > SummaryList
Definition: Misc.h:48
EquationSystemFactory & GetEquationSystemFactory()
void AddSummaryItem(SummaryList &l, const std::string &name, const std::string &value)
Adds a summary item to the summary info list.
Definition: Misc.cpp:49
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:172
The above copyright notice and this permission notice shall be included.
Definition: CoupledSolver.h:1
double NekDouble
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:492
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1255
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:291