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