Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
RinglebFlow.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File RinglebFlow.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 // License for the specific language governing rights and limitations under
14 // Permission is hereby granted, free of charge, to any person obtaining a
15 // copy of this software and associated documentation files (the "Software"),
16 // to deal in the Software without restriction, including without limitation
17 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
18 // and/or sell copies of the Software, and to permit persons to whom the
19 // Software is furnished to do so, subject to the following conditions:
20 //
21 // The above copyright notice and this permission notice shall be included
22 // in all copies or substantial portions of the Software.
23 //
24 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
25 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
26 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
27 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
28 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
29 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
30 // DEALINGS IN THE SOFTWARE.
31 //
32 // Description: Euler equations for Ringleb flow
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
37 
38 using namespace std;
39 
40 namespace Nektar
41 {
42  string RinglebFlow::className =
44  "RinglebFlow", RinglebFlow::create,
45  "Euler equations for Ringleb flow.");
46 
47  RinglebFlow::RinglebFlow(
49  : EulerCFE(pSession)
50  {
51  }
52 
53  /**
54  * @brief Destructor for EulerCFE class.
55  */
57  {
58  }
59 
60  /**
61  * @brief Print out a summary with some relevant information.
62  */
64  {
66  SolverUtils::AddSummaryItem(s, "Problem Type", "RinglebFlow");
67  }
68 
69  /**
70  * @brief Get the exact solutions for isentropic vortex and Ringleb
71  * flow problems.
72  */
74  unsigned int field,
75  Array<OneD, NekDouble> &outfield,
76  const NekDouble time)
77  {
78  GetExactRinglebFlow( field, outfield);
79  }
80 
81  /**
82  * @brief Compute the exact solution for the Ringleb flow problem.
83  */
85  int field,
86  Array<OneD, NekDouble> &outarray)
87  {
88  int nTotQuadPoints = GetTotPoints();
89 
90  Array<OneD, NekDouble> rho(nTotQuadPoints, 100.0);
91  Array<OneD, NekDouble> rhou(nTotQuadPoints);
92  Array<OneD, NekDouble> rhov(nTotQuadPoints);
93  Array<OneD, NekDouble> E(nTotQuadPoints);
94  Array<OneD, NekDouble> x(nTotQuadPoints);
95  Array<OneD, NekDouble> y(nTotQuadPoints);
96  Array<OneD, NekDouble> z(nTotQuadPoints);
97 
98  m_fields[0]->GetCoords(x, y, z);
99 
100  // Flow parameters
101  NekDouble c, k, phi, r, J, VV, pp, sint, P, ss;
102  NekDouble J11, J12, J21, J22, det;
103  NekDouble Fx, Fy;
104  NekDouble xi, yi;
105  NekDouble dV;
106  NekDouble dtheta;
107  NekDouble par1;
108  NekDouble theta = M_PI / 4.0;
109  NekDouble kExt = 0.7;
110  NekDouble V = kExt * sin(theta);
111  NekDouble toll = 1.0e-8;
112  NekDouble errV = 1.0;
113  NekDouble errTheta = 1.0;
114  NekDouble gamma = m_gamma;
115  NekDouble gamma_1_2 = (gamma - 1.0) / 2.0;
116 
117  for (int i = 0; i < nTotQuadPoints; ++i)
118  {
119  while ((abs(errV) > toll) || (abs(errTheta) > toll))
120  {
121  VV = V * V;
122  sint = sin(theta);
123  c = sqrt(1.0 - gamma_1_2 * VV);
124  k = V / sint;
125  phi = 1.0 / k;
126  pp = phi * phi;
127  J = 1.0 / c + 1.0 / (3.0 * c * c * c) +
128  1.0 / (5.0 * c * c * c * c * c) -
129  0.5 * log((1.0 + c) / (1.0 - c));
130 
131  r = pow(c, 1.0 / gamma_1_2);
132  xi = 1.0 / (2.0 * r) * (1.0 / VV - 2.0 * pp) + J / 2.0;
133  yi = phi / (r * V) * sqrt(1.0 - VV * pp);
134  par1 = 25.0 - 5.0 * VV;
135  ss = sint * sint;
136 
137  Fx = xi - x[i];
138  Fy = yi - y[i];
139 
140  J11 = 39062.5 / pow(par1, 3.5) * (1.0 / VV - 2.0 / VV * ss) *
141  V + 1562.5 / pow(par1, 2.5) * (-2.0 / (VV * V) + 4.0 /
142  (VV * V) * ss) + 12.5 / pow(par1, 1.5) * V + 312.5 /
143  pow(par1, 2.5) * V + 7812.5 / pow(par1, 3.5) * V -
144  0.25 * (-1.0 / pow(par1, 0.5) * V/(1.0 - 0.2 *
145  pow(par1, 0.5)) - (1.0 + 0.2 * pow(par1, 0.5)) /
146  pow((1.0 - 0.2 * pow(par1, 0.5)), 2.0) /
147  pow(par1, 0.5) * V) / (1.0 + 0.2 * pow(par1, 0.5)) *
148  (1.0 - 0.2 * pow(par1, 0.5));
149 
150  J12 = -6250.0 / pow(par1, 2.5) / VV * sint * cos(theta);
151  J21 = -6250.0 / (VV * V) * sint /
152  pow(par1, 2.5) * pow((1.0 - ss), 0.5) +
153  78125.0 / V * sint / pow(par1, 3.5) *
154  pow((1.0 - ss), 0.5);
155 
156  // the matrix is singular when theta = pi/2
157  if(abs(y[i])<toll && abs(cos(theta))<toll)
158  {
159  J22 = -39062.5 / pow(par1, 3.5) / V + 3125 /
160  pow(par1, 2.5) / (VV * V) + 12.5 / pow(par1, 1.5) *
161  V + 312.5 / pow(par1, 2.5) * V + 7812.5 /
162  pow(par1, 3.5) * V - 0.25 * (-1.0 / pow(par1, 0.5) *
163  V / (1.0 - 0.2 * pow(par1, 0.5)) - (1.0 + 0.2 *
164  pow(par1, 0.5)) / pow((1.0 - 0.2 *
165  pow(par1, 0.5)), 2.0) / pow(par1, 0.5) * V) /
166  (1.0 + 0.2 * pow(par1, 0.5)) * (1.0 - 0.2 *
167  pow(par1,0.5));
168 
169  // dV = -dV/dx * Fx
170  dV = -1.0 / J22 * Fx;
171  dtheta = 0.0;
172  theta = M_PI / 2.0;
173  }
174  else
175  {
176  J22 = 3125.0 / VV * cos(theta) / pow(par1, 2.5) *
177  pow((1.0 - ss), 0.5) - 3125.0 / VV * ss /
178  pow(par1, 2.5) / pow((1.0 - ss), 0.5) * cos(theta);
179 
180  det = -1.0 / (J11 * J22 - J12 * J21);
181 
182  // [dV dtheta]' = -[invJ]*[Fx Fy]'
183  dV = det * ( J22 * Fx - J12 * Fy);
184  dtheta = det * (-J21 * Fx + J11 * Fy);
185  }
186 
187  V = V + dV;
188  theta = theta + dtheta;
189 
190  errV = abs(dV);
191  errTheta = abs(dtheta);
192 
193  }
194 
195  c = sqrt(1.0 - gamma_1_2 * VV);
196  r = pow(c, 1.0 / gamma_1_2);
197 
198  rho[i] = r;
199  rhou[i] = rho[i] * V * cos(theta);
200  rhov[i] = rho[i] * V * sin(theta);
201  P = (c * c) * rho[i] / gamma;
202  E[i] = P / (gamma - 1.0) + 0.5 *
203  (rhou[i] * rhou[i] / rho[i] + rhov[i] * rhov[i] / rho[i]);
204 
205  // Resetting the guess value
206  errV = 1.0;
207  errTheta = 1.0;
208  theta = M_PI/4.0;
209  V = kExt*sin(theta);
210  }
211 
212  switch (field)
213  {
214  case 0:
215  outarray = rho;
216  break;
217  case 1:
218  outarray = rhou;
219  break;
220  case 2:
221  outarray = rhov;
222  break;
223  case 3:
224  outarray = E;
225  break;
226  default:
227  ASSERTL0(false, "Error in variable number!");
228  break;
229  }
230  }
231 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
std::vector< std::pair< std::string, std::string > > SummaryList
Definition: Misc.h:47
STL namespace.
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:51
SOLVER_UTILS_EXPORT int GetTotPoints()
virtual SOLVER_UTILS_EXPORT void v_GenerateSummary(SummaryList &s)
Print a summary of time stepping parameters.
void AddSummaryItem(SummaryList &l, const std::string &name, const std::string &value)
Adds a summary item to the summary info list.
Definition: Misc.cpp:50
void GetExactRinglebFlow(int field, Array< OneD, NekDouble > &outarray)
Ringleb Flow Test Case.
Definition: RinglebFlow.cpp:84
double NekDouble
EquationSystemFactory & GetEquationSystemFactory()
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
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.
Definition: RinglebFlow.cpp:73
virtual void v_GenerateSummary(SolverUtils::SummaryList &s)
Print a summary of time stepping parameters.
Definition: RinglebFlow.cpp:63
virtual ~RinglebFlow()
Destructor for EulerCFE class.
Definition: RinglebFlow.cpp:56
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215