Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
EulerCFE.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File EulerCFE.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 in consƒervative variables without artificial
33 // diffusion
34 //
35 ///////////////////////////////////////////////////////////////////////////////
36 
37 #include <iostream>
38 #include <iomanip>
39 #include <boost/algorithm/string.hpp>
40 
41 #include <MultiRegions/AssemblyMap/AssemblyMapDG.h>
42 #include <CompressibleFlowSolver/EquationSystems/EulerCFE.h>
43 
44 using namespace std;
45 
46 namespace Nektar
47 {
48  string EulerCFE::className =
49  SolverUtils::GetEquationSystemFactory().RegisterCreatorFunction(
50  "EulerCFE", EulerCFE::create,
51  "Euler equations in conservative variables without "
52  "artificial diffusion.");
53 
54  EulerCFE::EulerCFE(
55  const LibUtilities::SessionReaderSharedPtr& pSession)
56  : CompressibleFlowSystem(pSession)
57  {
58  }
59 
60  void EulerCFE::v_InitObject()
61  {
62  CompressibleFlowSystem::v_InitObject();
63 
64  if(m_session->DefinesSolverInfo("PROBLEMTYPE"))
65  {
66  int i;
67  std::string ProblemTypeStr =
68  m_session->GetSolverInfo("PROBLEMTYPE");
69  for (i = 0; i < (int) SIZE_ProblemType; ++i)
70  {
71  if (boost::iequals(ProblemTypeMap[i], ProblemTypeStr))
72  {
73  m_problemType = (ProblemType)i;
74  break;
75  }
76  }
77  }
78  else
79  {
80  m_problemType = (ProblemType)0;
81  }
82 
83  if (m_explicitAdvection)
84  {
85  m_ode.DefineOdeRhs (&EulerCFE::DoOdeRhs, this);
86  m_ode.DefineProjection (&EulerCFE::DoOdeProjection, this);
87  }
88  else
89  {
90  ASSERTL0(false, "Implicit CFE not set up.");
91  }
92  }
93 
94  /**
95  * @brief Destructor for EulerCFE class.
96  */
97  EulerCFE::~EulerCFE()
98  {
99  }
100 
101  /**
102  * @brief Print out a summary with some relevant information.
103  */
104  void EulerCFE::v_GenerateSummary(SolverUtils::SummaryList& s)
105  {
106  CompressibleFlowSystem::v_GenerateSummary(s);
107  SolverUtils::AddSummaryItem(s, "Problem Type", ProblemTypeMap[m_problemType]);
108  }
109 
110  /**
111  * @brief Set the initial conditions.
112  */
113  void EulerCFE::v_SetInitialConditions(
114  NekDouble initialtime,
115  bool dumpInitialConditions,
116  const int domain)
117  {
118  switch (m_problemType)
119  {
120  case eIsentropicVortex:
121  {
122  SetInitialIsentropicVortex(initialtime);
123  break;
124  }
125  default:
126  {
127  EquationSystem::v_SetInitialConditions(initialtime, false);
128  break;
129  }
130  }
131 
132  CompressibleFlowSystem::v_SetInitialConditions();
133 
134  if (dumpInitialConditions && m_checksteps)
135  {
136  // Dump initial conditions to file
137  Checkpoint_Output(0);
138  }
139  }
140 
141  /**
142  * @brief Compute the right-hand side.
143  */
144  void EulerCFE::DoOdeRhs(
145  const Array<OneD, const Array<OneD, NekDouble> > &inarray,
146  Array<OneD, Array<OneD, NekDouble> > &outarray,
147  const NekDouble time)
148  {
149  int i;
150  int nvariables = inarray.num_elements();
151  int npoints = GetNpoints();
152 
153  Array<OneD, Array<OneD, NekDouble> > advVel(m_spacedim);
154 
155  m_advection->Advect(nvariables, m_fields, advVel, inarray,
156  outarray, time);
157 
158  for (i = 0; i < nvariables; ++i)
159  {
160  Vmath::Neg(npoints, outarray[i], 1);
161  }
162 
163  // Add sponge layer if defined in the session file
164  std::vector<SolverUtils::ForcingSharedPtr>::const_iterator x;
165  for (x = m_forcing.begin(); x != m_forcing.end(); ++x)
166  {
167  (*x)->Apply(m_fields, inarray, outarray, time);
168  }
169  }
170 
171  /**
172  * @brief Compute the projection and call the method for imposing the
173  * boundary conditions in case of discontinuous projection.
174  */
175  void EulerCFE::DoOdeProjection(
176  const Array<OneD, const Array<OneD, NekDouble> > &inarray,
177  Array<OneD, Array<OneD, NekDouble> > &outarray,
178  const NekDouble time)
179  {
180  int i;
181  int nvariables = inarray.num_elements();
182 
183  switch (m_projectionType)
184  {
185  case MultiRegions::eDiscontinuous:
186  {
187  // Just copy over array
188  int npoints = GetNpoints();
189 
190  for(i = 0; i < nvariables; ++i)
191  {
192  Vmath::Vcopy(npoints, inarray[i], 1, outarray[i], 1);
193  }
194  SetBoundaryConditions(outarray, time);
195  break;
196  }
197  case MultiRegions::eGalerkin:
198  case MultiRegions::eMixed_CG_Discontinuous:
199  {
200  ASSERTL0(false, "No Continuous Galerkin for Euler equations");
201  break;
202  }
203  default:
204  ASSERTL0(false, "Unknown projection scheme");
205  break;
206  }
207  }
208 
209  /**
210  * @brief Set boundary conditions which can be:
211  * a) Wall and Symmerty BCs implemented at CompressibleFlowSystem level
212  * since they are compressible solver specific;
213  * b) Isentropic vortex and Ringleb flow BCs implemented at EulerCFE level
214  * since they are Euler solver specific;
215  * c) Time dependent BCs.
216  *
217  * @param inarray: fields array.
218  * @param time: time.
219  */
220  void EulerCFE::SetBoundaryConditions(
221  Array<OneD, Array<OneD, NekDouble> > &inarray,
222  NekDouble time)
223  {
224  std::string varName;
225  int cnt = 0;
226  int nvariables = inarray.num_elements();
227  int nTracePts = GetTraceTotPoints();
228 
229  Array<OneD, Array<OneD, NekDouble> > Fwd(nvariables);
230  for (int i = 0; i < nvariables; ++i)
231  {
232  Fwd[i] = Array<OneD, NekDouble>(nTracePts);
233  m_fields[i]->ExtractTracePhys(inarray[i], Fwd[i]);
234  }
235 
236  std::string userDefStr;
237  int nreg = m_fields[0]->GetBndConditions().num_elements();
238  // Loop over Boundary Regions
239  for (int n = 0; n < nreg; ++n)
240  {
241  userDefStr = m_fields[0]->GetBndConditions()[n]->GetUserDefined();
242  if(!userDefStr.empty())
243  {
244  if(boost::iequals(userDefStr,"WallViscous"))
245  {
246  ASSERTL0(false, "WallViscous is a wrong bc for the "
247  "Euler equations");
248  }
249  else if(boost::iequals(userDefStr, "IsentropicVortex"))
250  {
251  // Isentropic Vortex Boundary Condition
252  SetBoundaryIsentropicVortex(n, time, cnt, Fwd, inarray);
253  }
254  else if (boost::iequals(userDefStr,"RinglebFlow"))
255  {
256  SetBoundaryRinglebFlow(n, time, cnt, Fwd, inarray);
257  }
258  else
259  {
260  // set up userdefined BC common to all solvers
261  SetCommonBC(userDefStr, n, time, cnt, Fwd, inarray);
262  }
263  }
264 
265  // no User Defined conditions provided so skip cnt
266  cnt += m_fields[0]->GetBndCondExpansions()[n]->GetExpSize();
267  }
268  }
269 
270  /**
271  * @brief Get the exact solutions for isentropic vortex and Ringleb
272  * flow problems.
273  */
274  void EulerCFE::v_EvaluateExactSolution(
275  unsigned int field,
276  Array<OneD, NekDouble> &outfield,
277  const NekDouble time)
278  {
279  switch(m_problemType)
280  {
281  case eIsentropicVortex:
282  {
283  GetExactIsentropicVortex(field, outfield, time);
284  break;
285  }
286  case eRinglebFlow:
287  {
288  GetExactRinglebFlow(field, outfield);
289  break;
290  }
291  default:
292  {
293  EquationSystem::v_EvaluateExactSolution(field, outfield, time);
294  break;
295  }
296  }
297  }
298 
299  void EulerCFE::EvaluateIsentropicVortex(
300  const Array<OneD, NekDouble> &x,
301  const Array<OneD, NekDouble> &y,
302  const Array<OneD, NekDouble> &z,
303  Array<OneD, Array<OneD, NekDouble> > &u,
304  NekDouble time,
305  const int o)
306  {
307  int nq = x.num_elements();
308 
309  // Flow parameters
310  const NekDouble x0 = 5.0;
311  const NekDouble y0 = 0.0;
312  const NekDouble beta = 5.0;
313  const NekDouble u0 = 1.0;
314  const NekDouble v0 = 0.5;
315  const NekDouble gamma = m_gamma;
316  NekDouble r, xbar, ybar, tmp;
317  NekDouble fac = 1.0/(16.0*gamma*M_PI*M_PI);
318 
319  // In 3D zero rhow field.
320  if (m_spacedim == 3)
321  {
322  Vmath::Zero(nq, &u[3][o], 1);
323  }
324 
325  // Fill storage
326  for (int i = 0; i < nq; ++i)
327  {
328  xbar = x[i] - u0*time - x0;
329  ybar = y[i] - v0*time - y0;
330  r = sqrt(xbar*xbar + ybar*ybar);
331  tmp = beta*exp(1-r*r);
332  u[0][i+o] = pow(1.0 - (gamma-1.0)*tmp*tmp*fac, 1.0/(gamma-1.0));
333  u[1][i+o] = u[0][i+o]*(u0 - tmp*ybar/(2*M_PI));
334  u[2][i+o] = u[0][i+o]*(v0 + tmp*xbar/(2*M_PI));
335  u[m_spacedim+1][i+o] = pow(u[0][i+o], gamma)/(gamma-1.0) +
336  0.5*(u[1][i+o]*u[1][i+o] + u[2][i+o]*u[2][i+o]) / u[0][i+o];
337  }
338  }
339 
340  /**
341  * @brief Compute the exact solution for the isentropic vortex problem.
342  */
343  void EulerCFE::GetExactIsentropicVortex(
344  int field,
345  Array<OneD, NekDouble> &outarray,
346  NekDouble time)
347  {
348  int nTotQuadPoints = GetTotPoints();
349  Array<OneD, NekDouble> x(nTotQuadPoints);
350  Array<OneD, NekDouble> y(nTotQuadPoints);
351  Array<OneD, NekDouble> z(nTotQuadPoints);
352  Array<OneD, Array<OneD, NekDouble> > u(m_spacedim+2);
353 
354  m_fields[0]->GetCoords(x, y, z);
355 
356  for (int i = 0; i < m_spacedim + 2; ++i)
357  {
358  u[i] = Array<OneD, NekDouble>(nTotQuadPoints);
359  }
360 
361  EvaluateIsentropicVortex(x, y, z, u, time);
362 
363  Vmath::Vcopy(nTotQuadPoints, u[field], 1, outarray, 1);
364  }
365 
366  /**
367  * @brief Set the initial condition for the isentropic vortex problem.
368  */
369  void EulerCFE::SetInitialIsentropicVortex(NekDouble initialtime)
370  {
371  int nTotQuadPoints = GetTotPoints();
372  Array<OneD, NekDouble> x(nTotQuadPoints);
373  Array<OneD, NekDouble> y(nTotQuadPoints);
374  Array<OneD, NekDouble> z(nTotQuadPoints);
375  Array<OneD, Array<OneD, NekDouble> > u(m_spacedim+2);
376 
377  m_fields[0]->GetCoords(x, y, z);
378 
379  for (int i = 0; i < m_spacedim + 2; ++i)
380  {
381  u[i] = Array<OneD, NekDouble>(nTotQuadPoints);
382  }
383 
384  EvaluateIsentropicVortex(x, y, z, u, initialtime);
385 
386  // Forward transform to fill the coefficient space
387  for(int i = 0; i < m_fields.num_elements(); ++i)
388  {
389  Vmath::Vcopy(nTotQuadPoints, u[i], 1, m_fields[i]->UpdatePhys(), 1);
390  m_fields[i]->SetPhysState(true);
391  m_fields[i]->FwdTrans(m_fields[i]->GetPhys(),
392  m_fields[i]->UpdateCoeffs());
393  }
394  }
395 
396  /**
397  * @brief Set the boundary conditions for the isentropic vortex problem.
398  */
399  void EulerCFE::SetBoundaryIsentropicVortex(
400  int bcRegion,
401  NekDouble time,
402  int cnt,
403  Array<OneD, Array<OneD, NekDouble> > &Fwd,
404  Array<OneD, Array<OneD, NekDouble> > &physarray)
405  {
406  int nvariables = physarray.num_elements();
407  const Array<OneD, const int> &bndTraceMap = m_fields[0]->GetTraceBndMap();
408 
409  int id2, e_max;
410  e_max = m_fields[0]->GetBndCondExpansions()[bcRegion]->GetExpSize();
411 
412  for(int e = 0; e < e_max; ++e)
413  {
414  int npoints = m_fields[0]->
415  GetBndCondExpansions()[bcRegion]->GetExp(e)->GetTotPoints();
416  int id1 = m_fields[0]->
417  GetBndCondExpansions()[bcRegion]->GetPhys_Offset(e);
418  id2 = m_fields[0]->GetTrace()->GetPhys_Offset(bndTraceMap[cnt++]);
419 
420  Array<OneD,NekDouble> x(npoints, 0.0);
421  Array<OneD,NekDouble> y(npoints, 0.0);
422  Array<OneD,NekDouble> z(npoints, 0.0);
423 
424  m_fields[0]->GetBndCondExpansions()[bcRegion]->
425  GetExp(e)->GetCoords(x, y, z);
426 
427  EvaluateIsentropicVortex(x, y, z, Fwd, time, id2);
428 
429  for (int i = 0; i < nvariables; ++i)
430  {
431  Vmath::Vcopy(npoints, &Fwd[i][id2], 1,
432  &(m_fields[i]->GetBndCondExpansions()[bcRegion]->
433  UpdatePhys())[id1], 1);
434  }
435  }
436  }
437 
438  /**
439  * @brief Compute the exact solution for the Ringleb flow problem.
440  */
441  void EulerCFE::GetExactRinglebFlow(
442  int field,
443  Array<OneD, NekDouble> &outarray)
444  {
445  int nTotQuadPoints = GetTotPoints();
446 
447  Array<OneD, NekDouble> rho(nTotQuadPoints, 100.0);
448  Array<OneD, NekDouble> rhou(nTotQuadPoints);
449  Array<OneD, NekDouble> rhov(nTotQuadPoints);
450  Array<OneD, NekDouble> E(nTotQuadPoints);
451  Array<OneD, NekDouble> x(nTotQuadPoints);
452  Array<OneD, NekDouble> y(nTotQuadPoints);
453  Array<OneD, NekDouble> z(nTotQuadPoints);
454 
455  m_fields[0]->GetCoords(x, y, z);
456 
457  // Flow parameters
458  NekDouble c, k, phi, r, J, VV, pp, sint, P, ss;
459  NekDouble J11, J12, J21, J22, det;
460  NekDouble Fx, Fy;
461  NekDouble xi, yi;
462  NekDouble dV;
463  NekDouble dtheta;
464  NekDouble par1;
465  NekDouble theta = M_PI / 4.0;
466  NekDouble kExt = 0.7;
467  NekDouble V = kExt * sin(theta);
468  NekDouble toll = 1.0e-8;
469  NekDouble errV = 1.0;
470  NekDouble errTheta = 1.0;
471  NekDouble gamma = m_gamma;
472  NekDouble gamma_1_2 = (gamma - 1.0) / 2.0;
473 
474  for (int i = 0; i < nTotQuadPoints; ++i)
475  {
476  while ((abs(errV) > toll) || (abs(errTheta) > toll))
477  {
478  VV = V * V;
479  sint = sin(theta);
480  c = sqrt(1.0 - gamma_1_2 * VV);
481  k = V / sint;
482  phi = 1.0 / k;
483  pp = phi * phi;
484  J = 1.0 / c + 1.0 / (3.0 * c * c * c) +
485  1.0 / (5.0 * c * c * c * c * c) -
486  0.5 * log((1.0 + c) / (1.0 - c));
487 
488  r = pow(c, 1.0 / gamma_1_2);
489  xi = 1.0 / (2.0 * r) * (1.0 / VV - 2.0 * pp) + J / 2.0;
490  yi = phi / (r * V) * sqrt(1.0 - VV * pp);
491  par1 = 25.0 - 5.0 * VV;
492  ss = sint * sint;
493 
494  Fx = xi - x[i];
495  Fy = yi - y[i];
496 
497  J11 = 39062.5 / pow(par1, 3.5) * (1.0 / VV - 2.0 / VV * ss) *
498  V + 1562.5 / pow(par1, 2.5) * (-2.0 / (VV * V) + 4.0 /
499  (VV * V) * ss) + 12.5 / pow(par1, 1.5) * V + 312.5 /
500  pow(par1, 2.5) * V + 7812.5 / pow(par1, 3.5) * V -
501  0.25 * (-1.0 / pow(par1, 0.5) * V/(1.0 - 0.2 *
502  pow(par1, 0.5)) - (1.0 + 0.2 * pow(par1, 0.5)) /
503  pow((1.0 - 0.2 * pow(par1, 0.5)), 2.0) /
504  pow(par1, 0.5) * V) / (1.0 + 0.2 * pow(par1, 0.5)) *
505  (1.0 - 0.2 * pow(par1, 0.5));
506 
507  J12 = -6250.0 / pow(par1, 2.5) / VV * sint * cos(theta);
508  J21 = -6250.0 / (VV * V) * sint /
509  pow(par1, 2.5) * pow((1.0 - ss), 0.5) +
510  78125.0 / V * sint / pow(par1, 3.5) *
511  pow((1.0 - ss), 0.5);
512 
513  // the matrix is singular when theta = pi/2
514  if(abs(y[i])<toll && abs(cos(theta))<toll)
515  {
516  J22 = -39062.5 / pow(par1, 3.5) / V + 3125 /
517  pow(par1, 2.5) / (VV * V) + 12.5 / pow(par1, 1.5) *
518  V + 312.5 / pow(par1, 2.5) * V + 7812.5 /
519  pow(par1, 3.5) * V - 0.25 * (-1.0 / pow(par1, 0.5) *
520  V / (1.0 - 0.2 * pow(par1, 0.5)) - (1.0 + 0.2 *
521  pow(par1, 0.5)) / pow((1.0 - 0.2 *
522  pow(par1, 0.5)), 2.0) / pow(par1, 0.5) * V) /
523  (1.0 + 0.2 * pow(par1, 0.5)) * (1.0 - 0.2 *
524  pow(par1,0.5));
525 
526  // dV = -dV/dx * Fx
527  dV = -1.0 / J22 * Fx;
528  dtheta = 0.0;
529  theta = M_PI / 2.0;
530  }
531  else
532  {
533  J22 = 3125.0 / VV * cos(theta) / pow(par1, 2.5) *
534  pow((1.0 - ss), 0.5) - 3125.0 / VV * ss /
535  pow(par1, 2.5) / pow((1.0 - ss), 0.5) * cos(theta);
536 
537  det = -1.0 / (J11 * J22 - J12 * J21);
538 
539  // [dV dtheta]' = -[invJ]*[Fx Fy]'
540  dV = det * ( J22 * Fx - J12 * Fy);
541  dtheta = det * (-J21 * Fx + J11 * Fy);
542  }
543 
544  V = V + dV;
545  theta = theta + dtheta;
546 
547  errV = abs(dV);
548  errTheta = abs(dtheta);
549 
550  }
551 
552  c = sqrt(1.0 - gamma_1_2 * VV);
553  r = pow(c, 1.0 / gamma_1_2);
554 
555  rho[i] = r;
556  rhou[i] = rho[i] * V * cos(theta);
557  rhov[i] = rho[i] * V * sin(theta);
558  P = (c * c) * rho[i] / gamma;
559  E[i] = P / (gamma - 1.0) + 0.5 *
560  (rhou[i] * rhou[i] / rho[i] + rhov[i] * rhov[i] / rho[i]);
561 
562  // Resetting the guess value
563  errV = 1.0;
564  errTheta = 1.0;
565  theta = M_PI/4.0;
566  V = kExt*sin(theta);
567  }
568 
569  switch (field)
570  {
571  case 0:
572  outarray = rho;
573  break;
574  case 1:
575  outarray = rhou;
576  break;
577  case 2:
578  outarray = rhov;
579  break;
580  case 3:
581  outarray = E;
582  break;
583  default:
584  ASSERTL0(false, "Error in variable number!");
585  break;
586  }
587  }
588 
589  /**
590  * @brief Set the initial condition for the Ringleb flow problem.
591  */
592  void EulerCFE::SetInitialRinglebFlow(void)
593  {
594  // Get number of different boundaries in the input file
595  int nbnd = m_fields[0]->GetBndConditions().num_elements();
596 
597  // Loop on all the edges of the input file
598  for(int bcRegion=0; bcRegion < nbnd; ++bcRegion)
599  {
600 
601  int npoints = m_fields[0]->
602  GetBndCondExpansions()[bcRegion]->GetNpoints();
603 
604  Array<OneD,NekDouble> x0(npoints, 0.0);
605  Array<OneD,NekDouble> x1(npoints, 0.0);
606  Array<OneD,NekDouble> x2(npoints, 0.0);
607 
608  Array<OneD, NekDouble> rho(npoints, 0.0);
609  Array<OneD, NekDouble> rhou(npoints, 0.0);
610  Array<OneD, NekDouble> rhov(npoints, 0.0);
611  Array<OneD, NekDouble> E(npoints, 0.0);
612 
613  m_fields[0]->GetBndCondExpansions()[bcRegion]->
614  GetCoords(x0, x1, x2);
615 
616  // Flow parameters
617  NekDouble c, k, phi, r, J, VV, pp, sint, P, ss;
618  NekDouble J11, J12, J21, J22, det;
619  NekDouble Fx, Fy;
620  NekDouble xi, yi;
621  NekDouble dV;
622  NekDouble dtheta;
623  NekDouble par1;
624  NekDouble theta = M_PI / 4.0;
625  NekDouble kExt = 0.7;
626  NekDouble V = kExt * sin(theta);
627  NekDouble toll = 1.0e-8;
628  NekDouble errV = 1.0;
629  NekDouble errTheta = 1.0;
630  NekDouble gamma = m_gamma;
631  NekDouble gamma_1_2 = (gamma - 1.0) / 2.0;
632 
633  // Loop on all the points of that edge
634  for (int j = 0; j < npoints; j++)
635  {
636  while ((abs(errV) > toll) || (abs(errTheta) > toll))
637  {
638 
639  VV = V * V;
640  sint = sin(theta);
641  c = sqrt(1.0 - gamma_1_2 * VV);
642  k = V / sint;
643  phi = 1.0 / k;
644  pp = phi * phi;
645  J = 1.0 / c + 1.0 / (3.0 * c * c * c) +
646  1.0 / (5.0 * c * c * c * c * c) -
647  0.5 * log((1.0 + c) / (1.0 - c));
648 
649  r = pow(c, 1.0 / gamma_1_2);
650  xi = 1.0 / (2.0 * r) * (1.0 / VV - 2.0 * pp) + J / 2.0;
651  yi = phi / (r * V) * sqrt(1.0 - VV * pp);
652  par1 = 25.0 - 5.0 * VV;
653  ss = sint * sint;
654 
655  Fx = xi - x0[j];
656  Fy = yi - x1[j];
657 
658  J11 = 39062.5 / pow(par1, 3.5) *
659  (1.0 / VV - 2.0 / VV * ss) * V +
660  1562.5 / pow(par1, 2.5) * (-2.0 /
661  (VV * V) + 4.0 / (VV * V) * ss) +
662  12.5 / pow(par1, 1.5) * V +
663  312.5 / pow(par1, 2.5) * V +
664  7812.5 / pow(par1, 3.5) * V -
665  0.25 * (-1.0 / pow(par1, 0.5) * V /
666  (1.0 - 0.2 * pow(par1, 0.5)) - (1.0 + 0.2 *
667  pow(par1, 0.5)) / pow((1.0 - 0.2 *
668  pow(par1, 0.5)), 2.0) /
669  pow(par1, 0.5) * V) /
670  (1.0 + 0.2 * pow(par1, 0.5)) *
671  (1.0 - 0.2 * pow(par1, 0.5));
672 
673  J12 = -6250.0 / pow(par1, 2.5) / VV * sint * cos(theta);
674  J21 = -6250.0 / (VV * V) * sint / pow(par1, 2.5) *
675  pow((1.0 - ss), 0.5) + 78125.0 / V * sint /
676  pow(par1, 3.5) * pow((1.0 - ss), 0.5);
677 
678  // the matrix is singular when theta = pi/2
679  if (abs(x1[j]) < toll && abs(cos(theta)) < toll)
680  {
681 
682  J22 = -39062.5 / pow(par1, 3.5) / V +
683  3125 / pow(par1, 2.5) / (VV * V) + 12.5 /
684  pow(par1, 1.5) * V + 312.5 / pow(par1, 2.5) * V +
685  7812.5 / pow(par1, 3.5) * V -
686  0.25 * (-1.0 / pow(par1, 0.5) * V /
687  (1.0 - 0.2 * pow(par1, 0.5)) -
688  (1.0 + 0.2 * pow(par1, 0.5)) /
689  pow((1.0 - 0.2* pow(par1, 0.5)), 2.0) /
690  pow(par1, 0.5) *V) /
691  (1.0 + 0.2 * pow(par1, 0.5)) *
692  (1.0 - 0.2 * pow(par1, 0.5));
693 
694  // dV = -dV/dx * Fx
695  dV = -1.0 / J22 * Fx;
696  dtheta = 0.0;
697  theta = M_PI / 2.0;
698  }
699  else
700  {
701 
702  J22 = 3125.0 / VV * cos(theta) / pow(par1, 2.5) *
703  pow((1.0 - ss), 0.5) - 3125.0 / VV * ss /
704  pow(par1, 2.5) / pow((1.0 - ss), 0.5) *
705  cos(theta);
706 
707  det = -1.0 / (J11 * J22 - J12 * J21);
708 
709  // [dV dtheta]' = -[invJ]*[Fx Fy]'
710  dV = det * ( J22 * Fx - J12 * Fy);
711  dtheta = det * (-J21 * Fx + J11 * Fy);
712  }
713 
714  V = V + dV;
715  theta = theta + dtheta;
716 
717  errV = abs(dV);
718  errTheta = abs(dtheta);
719 
720  }
721 
722  c = sqrt(1.0 - gamma_1_2 * VV);
723  rho[j] = pow(c, 1.0 / gamma_1_2) * exp(-1.0);
724  rhou[j] = rho[j] * V * cos(theta) * exp(-1.0);
725  rhov[j] = rho[j] * V * sin(theta) * exp(-1.0);
726  P = (c * c) * rho[j] / gamma;
727  E[j] = P / (gamma - 1.0) + 0.5 *
728  (rhou[j] * rhou[j] /
729  rho[j] + rhov[j] * rhov[j] / rho[j]);
730  errV = 1.0;
731  errTheta = 1.0;
732  theta = M_PI / 4.0;
733  V = kExt * sin(theta);
734 
735  }
736 
737  // Fill the physical space
738  m_fields[0]->GetBndCondExpansions()[bcRegion]->SetPhys(rho);
739  m_fields[1]->GetBndCondExpansions()[bcRegion]->SetPhys(rhou);
740  m_fields[2]->GetBndCondExpansions()[bcRegion]->SetPhys(rhov);
741  m_fields[3]->GetBndCondExpansions()[bcRegion]->SetPhys(E);
742 
743  // Forward transform to fill the coefficients space
744  for(int i = 0; i < m_fields.num_elements(); ++i)
745  {
746  m_fields[i]->GetBndCondExpansions()[bcRegion]->
747  FwdTrans_BndConstrained(
748  m_fields[i]->GetBndCondExpansions()[bcRegion]->
749  GetPhys(),
750  m_fields[i]->GetBndCondExpansions()[bcRegion]->
751  UpdateCoeffs());
752  }
753 
754  }
755 
756  // Calculation of the initial internal values as a weighted average
757  // over the distance from the boundaries
758  int nq = m_fields[0]->GetNpoints();
759 
760  Array<OneD,NekDouble> x0(nq);
761  Array<OneD,NekDouble> x1(nq);
762  Array<OneD,NekDouble> x2(nq);
763 
764  // Get the coordinates
765  // (assuming all fields have the same discretisation)
766  m_fields[0]->GetCoords(x0, x1, x2);
767 
768  for(int j = 0; j < nq; j++)
769  {
770  NekDouble Dist = 0.0;
771  NekDouble rho = 0.0;
772  NekDouble rhou = 0.0;
773  NekDouble rhov = 0.0;
774  NekDouble E = 0.0;
775  NekDouble SumDist = 0.0;
776 
777  // Calculation of all the distances
778  // Loop on all the edges of the input file
779  for (int bcRegion = 0; bcRegion < nbnd; ++bcRegion)
780  {
781  // Get quadrature points on the edge
782  int npoints = m_fields[0]->
783  GetBndCondExpansions()[bcRegion]->GetNpoints();
784 
785  Array<OneD,NekDouble> xb0(npoints, 0.0);
786  Array<OneD,NekDouble> xb1(npoints, 0.0);
787  Array<OneD,NekDouble> xb2(npoints, 0.0);
788 
789  m_fields[0]->GetBndCondExpansions()[bcRegion]->
790  GetCoords(xb0, xb1, xb2);
791 
792  for (int k = 0; k < npoints; k++)
793  {
794  Dist = sqrt((xb0[k] - x0[j]) * (xb0[k] - x0[j]) +
795  (xb1[k] - x1[j]) * (xb1[k] - x1[j]));
796 
797  SumDist += Dist;
798  rho += Dist * (m_fields[0]->
799  GetBndCondExpansions()[bcRegion]->GetPhys())[k];
800  rhou += Dist * (m_fields[1]->
801  GetBndCondExpansions()[bcRegion]->GetPhys())[k];
802  rhov += Dist * (m_fields[2]->
803  GetBndCondExpansions()[bcRegion]->GetPhys())[k];
804  E += Dist * (m_fields[3]->
805  GetBndCondExpansions()[bcRegion]->GetPhys())[k];
806  }
807  }
808 
809  rho = rho / SumDist;
810  rhou = rhou / SumDist;
811  rhov = rhov / SumDist;
812  E = E / SumDist;
813 
814  (m_fields[0]->UpdatePhys())[j] = rho;
815  (m_fields[1]->UpdatePhys())[j] = rhou;
816  (m_fields[2]->UpdatePhys())[j] = rhov;
817  (m_fields[3]->UpdatePhys())[j] = E;
818 
819  }
820 
821  for (int i = 0 ; i < m_fields.num_elements(); i++)
822  {
823  m_fields[i]->SetPhysState(true);
824  m_fields[i]->FwdTrans(m_fields[i]->GetPhys(),
825  m_fields[i]->UpdateCoeffs());
826  }
827 
828  // Dump initial conditions to file
829  std::string outname = m_sessionName + "_initialRingleb.chk";
830  WriteFld(outname);
831  }
832 
833  /**
834  * @brief Set the boundary conditions for the Ringleb flow problem.
835  */
836  void EulerCFE::SetBoundaryRinglebFlow(
837  int bcRegion,
838  NekDouble time,
839  int cnt,
840  Array<OneD, Array<OneD, NekDouble> > &Fwd,
841  Array<OneD, Array<OneD, NekDouble> > &physarray)
842  {
843  int nvariables = physarray.num_elements();
844 
845  // For 3DHomogenoeus1D
846  int n_planes = 1;
847  if (m_expdim == 2 && m_HomogeneousType == eHomogeneous1D)
848  {
849  int nPointsTot = m_fields[0]->GetTotPoints();
850  int nPointsTot_plane = m_fields[0]->GetPlane(0)->GetTotPoints();
851  n_planes = nPointsTot/nPointsTot_plane;
852  }
853 
854  int id2, id2_plane, e_max;
855 
856  e_max = m_fields[0]->GetBndCondExpansions()[bcRegion]->GetExpSize();
857 
858  for(int e = 0; e < e_max; ++e)
859  {
860  int npoints = m_fields[0]->
861  GetBndCondExpansions()[bcRegion]->GetExp(e)->GetTotPoints();
862  int id1 = m_fields[0]->
863  GetBndCondExpansions()[bcRegion]->GetPhys_Offset(e);
864 
865  // For 3DHomogenoeus1D
866  if (m_expdim == 2 && m_HomogeneousType == eHomogeneous1D)
867  {
868  int cnt_plane = cnt/n_planes;
869  int e_plane;
870  int e_max_plane = e_max/n_planes;
871  int nTracePts_plane = GetTraceTotPoints();
872 
873  int planeID = floor((e + 0.5 )/ e_max_plane );
874  e_plane = e - e_max_plane*planeID;
875 
876  id2_plane = m_fields[0]->GetTrace()->GetPhys_Offset(
877  m_fields[0]->GetTraceMap()->
878  GetBndCondCoeffsToGlobalCoeffsMap(
879  cnt_plane + e_plane));
880  id2 = id2_plane + planeID*nTracePts_plane;
881  }
882  else // For general case
883  {
884  id2 = m_fields[0]->
885  GetTrace()->GetPhys_Offset(m_fields[0]->GetTraceMap()->
886  GetBndCondTraceToGlobalTraceMap(cnt++));
887  }
888 
889  Array<OneD,NekDouble> x0(npoints, 0.0);
890  Array<OneD,NekDouble> x1(npoints, 0.0);
891  Array<OneD,NekDouble> x2(npoints, 0.0);
892 
893  m_fields[0]->GetBndCondExpansions()[bcRegion]->
894  GetExp(e)->GetCoords(x0, x1, x2);
895 
896  // Flow parameters
897  NekDouble c, k, phi, r, J, VV, pp, sint, P, ss;
898  NekDouble J11, J12, J21, J22, det;
899  NekDouble Fx, Fy;
900  NekDouble xi, yi;
901  NekDouble dV;
902  NekDouble dtheta;
903  NekDouble par1;
904  NekDouble theta = M_PI / 4.0;
905  NekDouble kExt = 0.7;
906  NekDouble V = kExt * sin(theta);
907  NekDouble toll = 1.0e-8;
908  NekDouble errV = 1.0;
909  NekDouble errTheta = 1.0;
910  NekDouble gamma = m_gamma;
911  NekDouble gamma_1_2 = (gamma - 1.0) / 2.0;
912 
913  // Loop on all the points of that edge
914  for (int j = 0; j < npoints; j++)
915  {
916 
917  while ((abs(errV) > toll) || (abs(errTheta) > toll))
918  {
919  VV = V * V;
920  sint = sin(theta);
921  c = sqrt(1.0 - gamma_1_2 * VV);
922  k = V / sint;
923  phi = 1.0 / k;
924  pp = phi * phi;
925  J = 1.0 / c + 1.0 / (3.0 * c * c * c) +
926  1.0 / (5.0 * c * c * c * c * c) -
927  0.5 * log((1.0 + c) / (1.0 - c));
928 
929  r = pow(c, 1.0 / gamma_1_2);
930  xi = 1.0 / (2.0 * r) * (1.0 / VV - 2.0 * pp) + J / 2.0;
931  yi = phi / (r * V) * sqrt(1.0 - VV * pp);
932  par1 = 25.0 - 5.0 * VV;
933  ss = sint * sint;
934 
935  Fx = xi - x0[j];
936  Fy = yi - x1[j];
937 
938  J11 = 39062.5 / pow(par1, 3.5) *
939  (1.0 / VV - 2.0 / VV * ss) * V + 1562.5 /
940  pow(par1, 2.5) * (-2.0 / (VV * V) + 4.0 /
941  (VV * V) * ss) + 12.5 / pow(par1, 1.5) * V +
942  312.5 / pow(par1, 2.5) * V + 7812.5 /
943  pow(par1, 3.5) * V - 0.25 *
944  (-1.0 / pow(par1, 0.5) * V / (1.0 - 0.2 *
945  pow(par1, 0.5)) - (1.0 + 0.2 * pow(par1, 0.5)) /
946  pow((1.0 - 0.2 * pow(par1, 0.5)), 2.0) /
947  pow(par1, 0.5) * V) / (1.0 + 0.2 * pow(par1, 0.5)) *
948  (1.0 - 0.2 * pow(par1, 0.5));
949 
950  J12 = -6250.0 / pow(par1, 2.5) / VV * sint * cos(theta);
951  J21 = -6250.0 / (VV * V) * sint / pow(par1, 2.5) *
952  pow((1.0 - ss), 0.5) + 78125.0 / V * sint /
953  pow(par1, 3.5) * pow((1.0 - ss), 0.5);
954 
955  // the matrix is singular when theta = pi/2
956  if (abs(x1[j]) < toll && abs(cos(theta)) < toll)
957  {
958  J22 = -39062.5 / pow(par1, 3.5) / V + 3125 /
959  pow(par1, 2.5) / (VV * V) + 12.5 /
960  pow(par1, 1.5) * V + 312.5 / pow(par1, 2.5) *
961  V + 7812.5 / pow(par1, 3.5) * V - 0.25 *
962  (-1.0 / pow(par1, 0.5) * V / (1.0 - 0.2 *
963  pow(par1, 0.5)) - (1.0 + 0.2 * pow(par1, 0.5)) /
964  pow((1.0 - 0.2 * pow(par1, 0.5)), 2.0) /
965  pow(par1, 0.5) * V) / (1.0 + 0.2 *
966  pow(par1, 0.5)) * (1.0 - 0.2 * pow(par1, 0.5));
967 
968  // dV = -dV/dx * Fx
969  dV = -1.0 / J22 * Fx;
970  dtheta = 0.0;
971  theta = M_PI / 2.0;
972  }
973  else
974  {
975  J22 = 3125.0 / VV * cos(theta) / pow(par1, 2.5) *
976  pow((1.0 - ss), 0.5) - 3125.0 / VV * ss /
977  pow(par1, 2.5) / pow((1.0 - ss), 0.5) *
978  cos(theta);
979 
980  det = -1.0 / (J11 * J22 - J12 * J21);
981 
982  // [dV dtheta]' = -[invJ]*[Fx Fy]'
983  dV = det * ( J22 * Fx - J12 * Fy);
984  dtheta = det * (-J21 * Fx + J11 * Fy);
985  }
986 
987  V = V + dV;
988  theta = theta + dtheta;
989 
990  errV = abs(dV);
991  errTheta = abs(dtheta);
992  }
993 
994  c = sqrt(1.0 - gamma_1_2 * VV);
995  int kk = id2 + j;
996  NekDouble timeramp = 200.0;
997  std::string restartstr = "RESTART";
998  if (time<timeramp &&
999  !(m_session->DefinesFunction("InitialConditions") &&
1000  m_session->GetFunctionType("InitialConditions", 0) ==
1001  LibUtilities::eFunctionTypeFile))
1002  {
1003  Fwd[0][kk] = pow(c, 1.0 / gamma_1_2) *
1004  exp(-1.0 + time /timeramp);
1005 
1006  Fwd[1][kk] = Fwd[0][kk] * V * cos(theta) *
1007  exp(-1 + time / timeramp);
1008 
1009  Fwd[2][kk] = Fwd[0][kk] * V * sin(theta) *
1010  exp(-1 + time / timeramp);
1011  }
1012  else
1013  {
1014  Fwd[0][kk] = pow(c, 1.0 / gamma_1_2);
1015  Fwd[1][kk] = Fwd[0][kk] * V * cos(theta);
1016  Fwd[2][kk] = Fwd[0][kk] * V * sin(theta);
1017  }
1018 
1019  P = (c * c) * Fwd[0][kk] / gamma;
1020  Fwd[3][kk] = P / (gamma - 1.0) + 0.5 *
1021  (Fwd[1][kk] * Fwd[1][kk] / Fwd[0][kk] +
1022  Fwd[2][kk] * Fwd[2][kk] / Fwd[0][kk]);
1023 
1024  errV = 1.0;
1025  errTheta = 1.0;
1026  theta = M_PI / 4.0;
1027  V = kExt * sin(theta);
1028  }
1029 
1030  for (int i = 0; i < nvariables; ++i)
1031  {
1032  Vmath::Vcopy(npoints, &Fwd[i][id2], 1,
1033  &(m_fields[i]->GetBndCondExpansions()[bcRegion]->
1034  UpdatePhys())[id1],1);
1035  }
1036  }
1037  }
1038 }
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:188
Nektar::EulerCFE::v_SetInitialConditions
virtual void v_SetInitialConditions(NekDouble initialtime=0.0, bool dumpInitialConditions=true, const int domain=0)
Set the initial conditions.
Definition: EulerCFE.cpp:113
Nektar::CompressibleFlowSystem::v_InitObject
virtual void v_InitObject()
Initialization object for CompressibleFlowSystem class.
Definition: CompressibleFlowSystem.cpp:65
Nektar::EulerCFE::v_GenerateSummary
virtual void v_GenerateSummary(SolverUtils::SummaryList &s)
Print a summary of time stepping parameters.
Definition: EulerCFE.cpp:104
Nektar::SolverUtils::UnsteadySystem::m_ode
LibUtilities::TimeIntegrationSchemeOperators m_ode
The time integration scheme operators to use.
Definition: UnsteadySystem.h:67
Nektar::SolverUtils::SummaryList
std::vector< std::pair< std::string, std::string > > SummaryList
Definition: Misc.h:47
Nektar::EulerCFE::SetInitialIsentropicVortex
void SetInitialIsentropicVortex(NekDouble initialtime)
Set the initial condition for the isentropic vortex problem.
Definition: EulerCFE.cpp:369
Nektar::SolverUtils::EquationSystem::m_expdim
int m_expdim
Expansion dimension.
Definition: EquationSystem.h:494
Nektar::SolverUtils::EquationSystem::Checkpoint_Output
SOLVER_UTILS_EXPORT void Checkpoint_Output(const int n)
Write checkpoint file of m_fields.
Definition: EquationSystem.cpp:1969
Nektar::eRinglebFlow
Ringleb Flow.
Definition: EulerCFE.h:47
Nektar::EulerCFE::EvaluateIsentropicVortex
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.
Definition: EulerCFE.cpp:299
Nektar::SolverUtils::EquationSystem::m_projectionType
enum MultiRegions::ProjectionType m_projectionType
Type of projection; e.g continuous or discontinuous.
Definition: EquationSystem.h:514
Nektar::EulerCFE::SetBoundaryIsentropicVortex
void SetBoundaryIsentropicVortex(int bcRegion, NekDouble time, int cnt, Array< OneD, Array< OneD, NekDouble > > &Fwd, Array< OneD, Array< OneD, NekDouble > > &physarray)
Set the boundary conditions for the isentropic vortex problem.
Definition: EulerCFE.cpp:399
std
STL namespace.
Nektar::SolverUtils::EquationSystem::m_sessionName
std::string m_sessionName
Name of the session.
Definition: EquationSystem.h:470
Nektar::EulerCFE::SetInitialRinglebFlow
void SetInitialRinglebFlow(void)
Set the initial condition for the Ringleb flow problem.
Definition: EulerCFE.cpp:592
Nektar::LibUtilities::SessionReaderSharedPtr
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:51
Nektar::SolverUtils::EquationSystem::m_checksteps
int m_checksteps
Number of steps between checkpoints.
Definition: EquationSystem.h:490
Nektar::EulerCFE::SetBoundaryRinglebFlow
void SetBoundaryRinglebFlow(int bcRegion, NekDouble time, int cnt, Array< OneD, Array< OneD, NekDouble > > &Fwd, Array< OneD, Array< OneD, NekDouble > > &physarray)
Set the boundary conditions for the Ringleb flow problem.
Definition: EulerCFE.cpp:836
Nektar::EulerCFE::~EulerCFE
virtual ~EulerCFE()
problem type selector
Definition: EulerCFE.cpp:97
Nektar::SolverUtils::EquationSystem::GetTotPoints
SOLVER_UTILS_EXPORT int GetTotPoints()
Definition: EquationSystem.h:866
Nektar::CompressibleFlowSystem::v_GenerateSummary
virtual void v_GenerateSummary(SolverUtils::SummaryList &s)
Print a summary of time stepping parameters.
Definition: CompressibleFlowSystem.cpp:310
Nektar::LibUtilities::TimeIntegrationSchemeOperators::DefineProjection
void DefineProjection(FuncPointerT func, ObjectPointerT obj)
Definition: TimeIntegrationScheme.h:202
Nektar::EulerCFE::DoOdeProjection
void DoOdeProjection(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
Compute the projection and call the method for imposing the boundary conditions in case of discontinu...
Definition: EulerCFE.cpp:175
Nektar::SolverUtils::UnsteadySystem::m_explicitAdvection
bool m_explicitAdvection
Indicates if explicit or implicit treatment of advection is used.
Definition: UnsteadySystem.h:75
Nektar::LibUtilities::TimeIntegrationSchemeOperators::DefineOdeRhs
void DefineOdeRhs(FuncPointerT func, ObjectPointerT obj)
Definition: TimeIntegrationScheme.h:184
Nektar::SolverUtils::AddSummaryItem
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
Nektar::SolverUtils::EquationSystem::m_spacedim
int m_spacedim
Spatial dimension (>= expansion dim).
Definition: EquationSystem.h:492
Vmath::Neg
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:382
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:44
Nektar::SolverUtils::EquationSystem::v_SetInitialConditions
virtual SOLVER_UTILS_EXPORT void v_SetInitialConditions(NekDouble initialtime=0.0, bool dumpInitialConditions=true, const int domain=0)
Definition: EquationSystem.cpp:1180
Nektar::SIZE_ProblemType
Length of enum list.
Definition: EulerADCFE.h:47
Nektar::EulerCFE::DoOdeRhs
void DoOdeRhs(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
Compute the right-hand side.
Definition: EulerCFE.cpp:144
Nektar::EulerCFE::GetExactIsentropicVortex
void GetExactIsentropicVortex(int field, Array< OneD, NekDouble > &outarray, NekDouble time)
Compute the exact solution for the isentropic vortex problem.
Definition: EulerCFE.cpp:343
Nektar::EulerCFE::GetExactRinglebFlow
void GetExactRinglebFlow(int field, Array< OneD, NekDouble > &outarray)
Ringleb Flow Test Case.
Definition: EulerCFE.cpp:441
Nektar::SolverUtils::GetEquationSystemFactory
EquationSystemFactory & GetEquationSystemFactory()
Definition: EquationSystem.cpp:83
Nektar::CompressibleFlowSystem::m_advection
SolverUtils::AdvectionSharedPtr m_advection
Definition: CompressibleFlowSystem.h:81
Nektar::SolverUtils::EquationSystem::GetTraceTotPoints
SOLVER_UTILS_EXPORT int GetTraceTotPoints()
Definition: EquationSystem.h:841
Nektar::eIsentropicVortex
Isentropic Vortex.
Definition: EulerCFE.h:46
Nektar::SolverUtils::EquationSystem::WriteFld
SOLVER_UTILS_EXPORT void WriteFld(const std::string &outname)
Write field data to the given filename.
Definition: EquationSystem.cpp:2009
Nektar::EulerCFE::v_EvaluateExactSolution
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: EulerCFE.cpp:274
Nektar::SolverUtils::EquationSystem::GetNpoints
SOLVER_UTILS_EXPORT int GetNpoints()
Definition: EquationSystem.h:876
EulerCFE.h
Nektar::EulerCFE::m_problemType
ProblemType m_problemType
Definition: EulerCFE.h:77
Nektar::SolverUtils::EquationSystem::m_fields
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
Definition: EquationSystem.h:460
Nektar::CompressibleFlowSystem::v_SetInitialConditions
virtual void v_SetInitialConditions(NekDouble initialtime=0.0, bool dumpInitialConditions=true, const int domain=0)
Nektar::EulerCFE::v_InitObject
virtual void v_InitObject()
Initialization object for CompressibleFlowSystem class.
Definition: EulerCFE.cpp:60
Nektar::SolverUtils::EquationSystem::m_session
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
Definition: EquationSystem.h:452
Nektar::CompressibleFlowSystem::SetCommonBC
void SetCommonBC(const std::string &userDefStr, const int n, const NekDouble time, int &cnt, Array< OneD, Array< OneD, NekDouble > > &Fwd, Array< OneD, Array< OneD, NekDouble > > &inarray)
Set boundary conditions which can be: a) Wall and Symmerty BCs implemented at CompressibleFlowSystem ...
Definition: CompressibleFlowSystem.cpp:324
Nektar::ProblemType
ProblemType
Definition: EulerADCFE.h:44
AssemblyMapDG.h
Vmath::Zero
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:359
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1047
Nektar::CompressibleFlowSystem::m_forcing
std::vector< SolverUtils::ForcingSharedPtr > m_forcing
Definition: CompressibleFlowSystem.h:122
Nektar::SolverUtils::EquationSystem::v_EvaluateExactSolution
virtual SOLVER_UTILS_EXPORT void v_EvaluateExactSolution(unsigned int field, Array< OneD, NekDouble > &outfield, const NekDouble time)
Definition: EquationSystem.cpp:1231
Nektar::EulerCFE::SetBoundaryConditions
void SetBoundaryConditions(Array< OneD, Array< OneD, NekDouble > > &physarray, NekDouble time)
Set boundary conditions which can be: a) Wall and Symmerty BCs implemented at CompressibleFlowSystem ...
Definition: EulerCFE.cpp:220
Nektar::SolverUtils::EquationSystem::m_HomogeneousType
enum HomogeneousType m_HomogeneousType
Definition: EquationSystem.h:540
Nektar::LibUtilities::NekFactory::RegisterCreatorFunction
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215
Nektar::ProblemTypeMap
const char *const ProblemTypeMap[]
Definition: EulerADCFE.h:50
Generated on Fri Sep 9 2016 12:03:14 for Nektar++ by   doxygen 1.8.8