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 
227  std::string userDefStr;
228  int nreg = m_fields[0]->GetBndConditions().num_elements();
229  // Loop over Boundary Regions
230  for (int n = 0; n < nreg; ++n)
231  {
232  userDefStr = m_fields[0]->GetBndConditions()[n]->GetUserDefined();
233  if(!userDefStr.empty())
234  {
235  if(boost::iequals(userDefStr,"WallViscous"))
236  {
237  ASSERTL0(false, "WallViscous is a wrong bc for the "
238  "Euler equations");
239  }
240  else if(boost::iequals(userDefStr, "IsentropicVortex"))
241  {
242  // Isentropic Vortex Boundary Condition
243  SetBoundaryIsentropicVortex(n, time, cnt, inarray);
244  }
245  else if (boost::iequals(userDefStr,"RinglebFlow"))
246  {
247  SetBoundaryRinglebFlow(n, time, cnt, inarray);
248  }
249  else
250  {
251  // set up userdefined BC common to all solvers
252  SetCommonBC(userDefStr,n,time, cnt,inarray);
253  }
254  }
255 
256  // no User Defined conditions provided so skip cnt
257  cnt += m_fields[0]->GetBndCondExpansions()[n]->GetExpSize();
258  }
259  }
260 
261  /**
262  * @brief Get the exact solutions for isentropic vortex and Ringleb
263  * flow problems.
264  */
265  void EulerCFE::v_EvaluateExactSolution(
266  unsigned int field,
267  Array<OneD, NekDouble> &outfield,
268  const NekDouble time)
269  {
270  switch(m_problemType)
271  {
272  case eIsentropicVortex:
273  {
274  GetExactIsentropicVortex(field, outfield, time);
275  break;
276  }
277  case eRinglebFlow:
278  {
279  GetExactRinglebFlow(field, outfield);
280  break;
281  }
282  default:
283  {
284  EquationSystem::v_EvaluateExactSolution(field, outfield, time);
285  break;
286  }
287  }
288  }
289 
290  void EulerCFE::EvaluateIsentropicVortex(
291  const Array<OneD, NekDouble> &x,
292  const Array<OneD, NekDouble> &y,
293  const Array<OneD, NekDouble> &z,
294  Array<OneD, Array<OneD, NekDouble> > &u,
295  NekDouble time,
296  const int o)
297  {
298  int nq = x.num_elements();
299 
300  // Flow parameters
301  const NekDouble x0 = 5.0;
302  const NekDouble y0 = 0.0;
303  const NekDouble beta = 5.0;
304  const NekDouble u0 = 1.0;
305  const NekDouble v0 = 0.5;
306  const NekDouble gamma = m_gamma;
307  NekDouble r, xbar, ybar, tmp;
308  NekDouble fac = 1.0/(16.0*gamma*M_PI*M_PI);
309 
310  // In 3D zero rhow field.
311  if (m_spacedim == 3)
312  {
313  Vmath::Zero(nq, &u[3][o], 1);
314  }
315 
316  // Fill storage
317  for (int i = 0; i < nq; ++i)
318  {
319  xbar = x[i] - u0*time - x0;
320  ybar = y[i] - v0*time - y0;
321  r = sqrt(xbar*xbar + ybar*ybar);
322  tmp = beta*exp(1-r*r);
323  u[0][i+o] = pow(1.0 - (gamma-1.0)*tmp*tmp*fac, 1.0/(gamma-1.0));
324  u[1][i+o] = u[0][i+o]*(u0 - tmp*ybar/(2*M_PI));
325  u[2][i+o] = u[0][i+o]*(v0 + tmp*xbar/(2*M_PI));
326  u[m_spacedim+1][i+o] = pow(u[0][i+o], gamma)/(gamma-1.0) +
327  0.5*(u[1][i+o]*u[1][i+o] + u[2][i+o]*u[2][i+o]) / u[0][i+o];
328  }
329  }
330 
331  /**
332  * @brief Compute the exact solution for the isentropic vortex problem.
333  */
334  void EulerCFE::GetExactIsentropicVortex(
335  int field,
336  Array<OneD, NekDouble> &outarray,
337  NekDouble time)
338  {
339  int nTotQuadPoints = GetTotPoints();
340  Array<OneD, NekDouble> x(nTotQuadPoints);
341  Array<OneD, NekDouble> y(nTotQuadPoints);
342  Array<OneD, NekDouble> z(nTotQuadPoints);
343  Array<OneD, Array<OneD, NekDouble> > u(m_spacedim+2);
344 
345  m_fields[0]->GetCoords(x, y, z);
346 
347  for (int i = 0; i < m_spacedim + 2; ++i)
348  {
349  u[i] = Array<OneD, NekDouble>(nTotQuadPoints);
350  }
351 
352  EvaluateIsentropicVortex(x, y, z, u, time);
353 
354  Vmath::Vcopy(nTotQuadPoints, u[field], 1, outarray, 1);
355  }
356 
357  /**
358  * @brief Set the initial condition for the isentropic vortex problem.
359  */
360  void EulerCFE::SetInitialIsentropicVortex(NekDouble initialtime)
361  {
362  int nTotQuadPoints = GetTotPoints();
363  Array<OneD, NekDouble> x(nTotQuadPoints);
364  Array<OneD, NekDouble> y(nTotQuadPoints);
365  Array<OneD, NekDouble> z(nTotQuadPoints);
366  Array<OneD, Array<OneD, NekDouble> > u(m_spacedim+2);
367 
368  m_fields[0]->GetCoords(x, y, z);
369 
370  for (int i = 0; i < m_spacedim + 2; ++i)
371  {
372  u[i] = Array<OneD, NekDouble>(nTotQuadPoints);
373  }
374 
375  EvaluateIsentropicVortex(x, y, z, u, initialtime);
376 
377  // Forward transform to fill the coefficient space
378  for(int i = 0; i < m_fields.num_elements(); ++i)
379  {
380  Vmath::Vcopy(nTotQuadPoints, u[i], 1, m_fields[i]->UpdatePhys(), 1);
381  m_fields[i]->SetPhysState(true);
382  m_fields[i]->FwdTrans(m_fields[i]->GetPhys(),
383  m_fields[i]->UpdateCoeffs());
384  }
385  }
386 
387  /**
388  * @brief Set the boundary conditions for the isentropic vortex problem.
389  */
390  void EulerCFE::SetBoundaryIsentropicVortex(
391  int bcRegion,
392  NekDouble time,
393  int cnt,
394  Array<OneD, Array<OneD, NekDouble> > &physarray)
395  {
396  int nvariables = physarray.num_elements();
397  int nTraceNumPoints = GetTraceTotPoints();
398  Array<OneD, Array<OneD, NekDouble> > Fwd(nvariables);
399  const Array<OneD, const int> &bndTraceMap = m_fields[0]->GetTraceBndMap();
400 
401  // Get physical values of the forward trace (from exp to phys)
402  for (int i = 0; i < nvariables; ++i)
403  {
404  Fwd[i] = Array<OneD, NekDouble>(nTraceNumPoints);
405  m_fields[i]->ExtractTracePhys(physarray[i], Fwd[i]);
406  }
407 
408  int id2, e_max;
409  e_max = m_fields[0]->GetBndCondExpansions()[bcRegion]->GetExpSize();
410 
411  for(int e = 0; e < e_max; ++e)
412  {
413  int npoints = m_fields[0]->
414  GetBndCondExpansions()[bcRegion]->GetExp(e)->GetTotPoints();
415  int id1 = m_fields[0]->
416  GetBndCondExpansions()[bcRegion]->GetPhys_Offset(e);
417  id2 = m_fields[0]->GetTrace()->GetPhys_Offset(bndTraceMap[cnt++]);
418 
419  Array<OneD,NekDouble> x(npoints, 0.0);
420  Array<OneD,NekDouble> y(npoints, 0.0);
421  Array<OneD,NekDouble> z(npoints, 0.0);
422 
423  m_fields[0]->GetBndCondExpansions()[bcRegion]->
424  GetExp(e)->GetCoords(x, y, z);
425 
426  EvaluateIsentropicVortex(x, y, z, Fwd, time, id2);
427 
428  for (int i = 0; i < nvariables; ++i)
429  {
430  Vmath::Vcopy(npoints, &Fwd[i][id2], 1,
431  &(m_fields[i]->GetBndCondExpansions()[bcRegion]->
432  UpdatePhys())[id1], 1);
433  }
434  }
435  }
436 
437  /**
438  * @brief Compute the exact solution for the Ringleb flow problem.
439  */
440  void EulerCFE::GetExactRinglebFlow(
441  int field,
442  Array<OneD, NekDouble> &outarray)
443  {
444  int nTotQuadPoints = GetTotPoints();
445 
446  Array<OneD, NekDouble> rho(nTotQuadPoints, 100.0);
447  Array<OneD, NekDouble> rhou(nTotQuadPoints);
448  Array<OneD, NekDouble> rhov(nTotQuadPoints);
449  Array<OneD, NekDouble> E(nTotQuadPoints);
450  Array<OneD, NekDouble> x(nTotQuadPoints);
451  Array<OneD, NekDouble> y(nTotQuadPoints);
452  Array<OneD, NekDouble> z(nTotQuadPoints);
453 
454  m_fields[0]->GetCoords(x, y, z);
455 
456  // Flow parameters
457  NekDouble c, k, phi, r, J, VV, pp, sint, P, ss;
458  NekDouble J11, J12, J21, J22, det;
459  NekDouble Fx, Fy;
460  NekDouble xi, yi;
461  NekDouble dV;
462  NekDouble dtheta;
463  NekDouble par1;
464  NekDouble theta = M_PI / 4.0;
465  NekDouble kExt = 0.7;
466  NekDouble V = kExt * sin(theta);
467  NekDouble toll = 1.0e-8;
468  NekDouble errV = 1.0;
469  NekDouble errTheta = 1.0;
470  NekDouble gamma = m_gamma;
471  NekDouble gamma_1_2 = (gamma - 1.0) / 2.0;
472 
473  for (int i = 0; i < nTotQuadPoints; ++i)
474  {
475  while ((abs(errV) > toll) || (abs(errTheta) > toll))
476  {
477  VV = V * V;
478  sint = sin(theta);
479  c = sqrt(1.0 - gamma_1_2 * VV);
480  k = V / sint;
481  phi = 1.0 / k;
482  pp = phi * phi;
483  J = 1.0 / c + 1.0 / (3.0 * c * c * c) +
484  1.0 / (5.0 * c * c * c * c * c) -
485  0.5 * log((1.0 + c) / (1.0 - c));
486 
487  r = pow(c, 1.0 / gamma_1_2);
488  xi = 1.0 / (2.0 * r) * (1.0 / VV - 2.0 * pp) + J / 2.0;
489  yi = phi / (r * V) * sqrt(1.0 - VV * pp);
490  par1 = 25.0 - 5.0 * VV;
491  ss = sint * sint;
492 
493  Fx = xi - x[i];
494  Fy = yi - y[i];
495 
496  J11 = 39062.5 / pow(par1, 3.5) * (1.0 / VV - 2.0 / VV * ss) *
497  V + 1562.5 / pow(par1, 2.5) * (-2.0 / (VV * V) + 4.0 /
498  (VV * V) * ss) + 12.5 / pow(par1, 1.5) * V + 312.5 /
499  pow(par1, 2.5) * V + 7812.5 / pow(par1, 3.5) * V -
500  0.25 * (-1.0 / pow(par1, 0.5) * V/(1.0 - 0.2 *
501  pow(par1, 0.5)) - (1.0 + 0.2 * pow(par1, 0.5)) /
502  pow((1.0 - 0.2 * pow(par1, 0.5)), 2.0) /
503  pow(par1, 0.5) * V) / (1.0 + 0.2 * pow(par1, 0.5)) *
504  (1.0 - 0.2 * pow(par1, 0.5));
505 
506  J12 = -6250.0 / pow(par1, 2.5) / VV * sint * cos(theta);
507  J21 = -6250.0 / (VV * V) * sint /
508  pow(par1, 2.5) * pow((1.0 - ss), 0.5) +
509  78125.0 / V * sint / pow(par1, 3.5) *
510  pow((1.0 - ss), 0.5);
511 
512  // the matrix is singular when theta = pi/2
513  if(abs(y[i])<toll && abs(cos(theta))<toll)
514  {
515  J22 = -39062.5 / pow(par1, 3.5) / V + 3125 /
516  pow(par1, 2.5) / (VV * V) + 12.5 / pow(par1, 1.5) *
517  V + 312.5 / pow(par1, 2.5) * V + 7812.5 /
518  pow(par1, 3.5) * V - 0.25 * (-1.0 / pow(par1, 0.5) *
519  V / (1.0 - 0.2 * pow(par1, 0.5)) - (1.0 + 0.2 *
520  pow(par1, 0.5)) / pow((1.0 - 0.2 *
521  pow(par1, 0.5)), 2.0) / pow(par1, 0.5) * V) /
522  (1.0 + 0.2 * pow(par1, 0.5)) * (1.0 - 0.2 *
523  pow(par1,0.5));
524 
525  // dV = -dV/dx * Fx
526  dV = -1.0 / J22 * Fx;
527  dtheta = 0.0;
528  theta = M_PI / 2.0;
529  }
530  else
531  {
532  J22 = 3125.0 / VV * cos(theta) / pow(par1, 2.5) *
533  pow((1.0 - ss), 0.5) - 3125.0 / VV * ss /
534  pow(par1, 2.5) / pow((1.0 - ss), 0.5) * cos(theta);
535 
536  det = -1.0 / (J11 * J22 - J12 * J21);
537 
538  // [dV dtheta]' = -[invJ]*[Fx Fy]'
539  dV = det * ( J22 * Fx - J12 * Fy);
540  dtheta = det * (-J21 * Fx + J11 * Fy);
541  }
542 
543  V = V + dV;
544  theta = theta + dtheta;
545 
546  errV = abs(dV);
547  errTheta = abs(dtheta);
548 
549  }
550 
551  c = sqrt(1.0 - gamma_1_2 * VV);
552  r = pow(c, 1.0 / gamma_1_2);
553 
554  rho[i] = r;
555  rhou[i] = rho[i] * V * cos(theta);
556  rhov[i] = rho[i] * V * sin(theta);
557  P = (c * c) * rho[i] / gamma;
558  E[i] = P / (gamma - 1.0) + 0.5 *
559  (rhou[i] * rhou[i] / rho[i] + rhov[i] * rhov[i] / rho[i]);
560 
561  // Resetting the guess value
562  errV = 1.0;
563  errTheta = 1.0;
564  theta = M_PI/4.0;
565  V = kExt*sin(theta);
566  }
567 
568  switch (field)
569  {
570  case 0:
571  outarray = rho;
572  break;
573  case 1:
574  outarray = rhou;
575  break;
576  case 2:
577  outarray = rhov;
578  break;
579  case 3:
580  outarray = E;
581  break;
582  default:
583  ASSERTL0(false, "Error in variable number!");
584  break;
585  }
586  }
587 
588  /**
589  * @brief Set the initial condition for the Ringleb flow problem.
590  */
591  void EulerCFE::SetInitialRinglebFlow(void)
592  {
593  // Get number of different boundaries in the input file
594  int nbnd = m_fields[0]->GetBndConditions().num_elements();
595 
596  // Loop on all the edges of the input file
597  for(int bcRegion=0; bcRegion < nbnd; ++bcRegion)
598  {
599 
600  int npoints = m_fields[0]->
601  GetBndCondExpansions()[bcRegion]->GetNpoints();
602 
603  Array<OneD,NekDouble> x0(npoints, 0.0);
604  Array<OneD,NekDouble> x1(npoints, 0.0);
605  Array<OneD,NekDouble> x2(npoints, 0.0);
606 
607  Array<OneD, NekDouble> rho(npoints, 0.0);
608  Array<OneD, NekDouble> rhou(npoints, 0.0);
609  Array<OneD, NekDouble> rhov(npoints, 0.0);
610  Array<OneD, NekDouble> E(npoints, 0.0);
611 
612  m_fields[0]->GetBndCondExpansions()[bcRegion]->
613  GetCoords(x0, x1, x2);
614 
615  // Flow parameters
616  NekDouble c, k, phi, r, J, VV, pp, sint, P, ss;
617  NekDouble J11, J12, J21, J22, det;
618  NekDouble Fx, Fy;
619  NekDouble xi, yi;
620  NekDouble dV;
621  NekDouble dtheta;
622  NekDouble par1;
623  NekDouble theta = M_PI / 4.0;
624  NekDouble kExt = 0.7;
625  NekDouble V = kExt * sin(theta);
626  NekDouble toll = 1.0e-8;
627  NekDouble errV = 1.0;
628  NekDouble errTheta = 1.0;
629  NekDouble gamma = m_gamma;
630  NekDouble gamma_1_2 = (gamma - 1.0) / 2.0;
631 
632  // Loop on all the points of that edge
633  for (int j = 0; j < npoints; j++)
634  {
635  while ((abs(errV) > toll) || (abs(errTheta) > toll))
636  {
637 
638  VV = V * V;
639  sint = sin(theta);
640  c = sqrt(1.0 - gamma_1_2 * VV);
641  k = V / sint;
642  phi = 1.0 / k;
643  pp = phi * phi;
644  J = 1.0 / c + 1.0 / (3.0 * c * c * c) +
645  1.0 / (5.0 * c * c * c * c * c) -
646  0.5 * log((1.0 + c) / (1.0 - c));
647 
648  r = pow(c, 1.0 / gamma_1_2);
649  xi = 1.0 / (2.0 * r) * (1.0 / VV - 2.0 * pp) + J / 2.0;
650  yi = phi / (r * V) * sqrt(1.0 - VV * pp);
651  par1 = 25.0 - 5.0 * VV;
652  ss = sint * sint;
653 
654  Fx = xi - x0[j];
655  Fy = yi - x1[j];
656 
657  J11 = 39062.5 / pow(par1, 3.5) *
658  (1.0 / VV - 2.0 / VV * ss) * V +
659  1562.5 / pow(par1, 2.5) * (-2.0 /
660  (VV * V) + 4.0 / (VV * V) * ss) +
661  12.5 / pow(par1, 1.5) * V +
662  312.5 / pow(par1, 2.5) * V +
663  7812.5 / pow(par1, 3.5) * V -
664  0.25 * (-1.0 / pow(par1, 0.5) * V /
665  (1.0 - 0.2 * pow(par1, 0.5)) - (1.0 + 0.2 *
666  pow(par1, 0.5)) / pow((1.0 - 0.2 *
667  pow(par1, 0.5)), 2.0) /
668  pow(par1, 0.5) * V) /
669  (1.0 + 0.2 * pow(par1, 0.5)) *
670  (1.0 - 0.2 * pow(par1, 0.5));
671 
672  J12 = -6250.0 / pow(par1, 2.5) / VV * sint * cos(theta);
673  J21 = -6250.0 / (VV * V) * sint / pow(par1, 2.5) *
674  pow((1.0 - ss), 0.5) + 78125.0 / V * sint /
675  pow(par1, 3.5) * pow((1.0 - ss), 0.5);
676 
677  // the matrix is singular when theta = pi/2
678  if (abs(x1[j]) < toll && abs(cos(theta)) < toll)
679  {
680 
681  J22 = -39062.5 / pow(par1, 3.5) / V +
682  3125 / pow(par1, 2.5) / (VV * V) + 12.5 /
683  pow(par1, 1.5) * V + 312.5 / pow(par1, 2.5) * V +
684  7812.5 / pow(par1, 3.5) * V -
685  0.25 * (-1.0 / pow(par1, 0.5) * V /
686  (1.0 - 0.2 * pow(par1, 0.5)) -
687  (1.0 + 0.2 * pow(par1, 0.5)) /
688  pow((1.0 - 0.2* pow(par1, 0.5)), 2.0) /
689  pow(par1, 0.5) *V) /
690  (1.0 + 0.2 * pow(par1, 0.5)) *
691  (1.0 - 0.2 * pow(par1, 0.5));
692 
693  // dV = -dV/dx * Fx
694  dV = -1.0 / J22 * Fx;
695  dtheta = 0.0;
696  theta = M_PI / 2.0;
697  }
698  else
699  {
700 
701  J22 = 3125.0 / VV * cos(theta) / pow(par1, 2.5) *
702  pow((1.0 - ss), 0.5) - 3125.0 / VV * ss /
703  pow(par1, 2.5) / pow((1.0 - ss), 0.5) *
704  cos(theta);
705 
706  det = -1.0 / (J11 * J22 - J12 * J21);
707 
708  // [dV dtheta]' = -[invJ]*[Fx Fy]'
709  dV = det * ( J22 * Fx - J12 * Fy);
710  dtheta = det * (-J21 * Fx + J11 * Fy);
711  }
712 
713  V = V + dV;
714  theta = theta + dtheta;
715 
716  errV = abs(dV);
717  errTheta = abs(dtheta);
718 
719  }
720 
721  c = sqrt(1.0 - gamma_1_2 * VV);
722  rho[j] = pow(c, 1.0 / gamma_1_2) * exp(-1.0);
723  rhou[j] = rho[j] * V * cos(theta) * exp(-1.0);
724  rhov[j] = rho[j] * V * sin(theta) * exp(-1.0);
725  P = (c * c) * rho[j] / gamma;
726  E[j] = P / (gamma - 1.0) + 0.5 *
727  (rhou[j] * rhou[j] /
728  rho[j] + rhov[j] * rhov[j] / rho[j]);
729  errV = 1.0;
730  errTheta = 1.0;
731  theta = M_PI / 4.0;
732  V = kExt * sin(theta);
733 
734  }
735 
736  // Fill the physical space
737  m_fields[0]->GetBndCondExpansions()[bcRegion]->SetPhys(rho);
738  m_fields[1]->GetBndCondExpansions()[bcRegion]->SetPhys(rhou);
739  m_fields[2]->GetBndCondExpansions()[bcRegion]->SetPhys(rhov);
740  m_fields[3]->GetBndCondExpansions()[bcRegion]->SetPhys(E);
741 
742  // Forward transform to fill the coefficients space
743  for(int i = 0; i < m_fields.num_elements(); ++i)
744  {
745  m_fields[i]->GetBndCondExpansions()[bcRegion]->
746  FwdTrans_BndConstrained(
747  m_fields[i]->GetBndCondExpansions()[bcRegion]->
748  GetPhys(),
749  m_fields[i]->GetBndCondExpansions()[bcRegion]->
750  UpdateCoeffs());
751  }
752 
753  }
754 
755  // Calculation of the initial internal values as a weighted average
756  // over the distance from the boundaries
757  int nq = m_fields[0]->GetNpoints();
758 
759  Array<OneD,NekDouble> x0(nq);
760  Array<OneD,NekDouble> x1(nq);
761  Array<OneD,NekDouble> x2(nq);
762 
763  // Get the coordinates
764  // (assuming all fields have the same discretisation)
765  m_fields[0]->GetCoords(x0, x1, x2);
766 
767  for(int j = 0; j < nq; j++)
768  {
769  NekDouble Dist = 0.0;
770  NekDouble rho = 0.0;
771  NekDouble rhou = 0.0;
772  NekDouble rhov = 0.0;
773  NekDouble E = 0.0;
774  NekDouble SumDist = 0.0;
775 
776  // Calculation of all the distances
777  // Loop on all the edges of the input file
778  for (int bcRegion = 0; bcRegion < nbnd; ++bcRegion)
779  {
780  // Get quadrature points on the edge
781  int npoints = m_fields[0]->
782  GetBndCondExpansions()[bcRegion]->GetNpoints();
783 
784  Array<OneD,NekDouble> xb0(npoints, 0.0);
785  Array<OneD,NekDouble> xb1(npoints, 0.0);
786  Array<OneD,NekDouble> xb2(npoints, 0.0);
787 
788  m_fields[0]->GetBndCondExpansions()[bcRegion]->
789  GetCoords(xb0, xb1, xb2);
790 
791  for (int k = 0; k < npoints; k++)
792  {
793  Dist = sqrt((xb0[k] - x0[j]) * (xb0[k] - x0[j]) +
794  (xb1[k] - x1[j]) * (xb1[k] - x1[j]));
795 
796  SumDist += Dist;
797  rho += Dist * (m_fields[0]->
798  GetBndCondExpansions()[bcRegion]->GetPhys())[k];
799  rhou += Dist * (m_fields[1]->
800  GetBndCondExpansions()[bcRegion]->GetPhys())[k];
801  rhov += Dist * (m_fields[2]->
802  GetBndCondExpansions()[bcRegion]->GetPhys())[k];
803  E += Dist * (m_fields[3]->
804  GetBndCondExpansions()[bcRegion]->GetPhys())[k];
805  }
806  }
807 
808  rho = rho / SumDist;
809  rhou = rhou / SumDist;
810  rhov = rhov / SumDist;
811  E = E / SumDist;
812 
813  (m_fields[0]->UpdatePhys())[j] = rho;
814  (m_fields[1]->UpdatePhys())[j] = rhou;
815  (m_fields[2]->UpdatePhys())[j] = rhov;
816  (m_fields[3]->UpdatePhys())[j] = E;
817 
818  }
819 
820  for (int i = 0 ; i < m_fields.num_elements(); i++)
821  {
822  m_fields[i]->SetPhysState(true);
823  m_fields[i]->FwdTrans(m_fields[i]->GetPhys(),
824  m_fields[i]->UpdateCoeffs());
825  }
826 
827  // Dump initial conditions to file
828  std::string outname = m_sessionName + "_initialRingleb.chk";
829  WriteFld(outname);
830  }
831 
832  /**
833  * @brief Set the boundary conditions for the Ringleb flow problem.
834  */
835  void EulerCFE::SetBoundaryRinglebFlow(
836  int bcRegion,
837  NekDouble time,
838  int cnt,
839  Array<OneD, Array<OneD, NekDouble> > &physarray)
840  {
841  int nvariables = physarray.num_elements();
842  int nTraceNumPoints = GetTraceTotPoints();
843  Array<OneD, Array<OneD, NekDouble> > Fwd(nvariables);
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  nTraceNumPoints = nTraceNumPoints * n_planes;
853  }
854 
855  // Get physical values of the forward trace (from exp to phys)
856  for (int i = 0; i < nvariables; ++i)
857  {
858  Fwd[i] = Array<OneD, NekDouble>(nTraceNumPoints);
859  m_fields[i]->ExtractTracePhys(physarray[i], Fwd[i]);
860  }
861 
862  int id2, id2_plane, e_max;
863 
864  e_max = m_fields[0]->GetBndCondExpansions()[bcRegion]->GetExpSize();
865 
866  for(int e = 0; e < e_max; ++e)
867  {
868  int npoints = m_fields[0]->
869  GetBndCondExpansions()[bcRegion]->GetExp(e)->GetTotPoints();
870  int id1 = m_fields[0]->
871  GetBndCondExpansions()[bcRegion]->GetPhys_Offset(e);
872 
873  // For 3DHomogenoeus1D
874  if (m_expdim == 2 && m_HomogeneousType == eHomogeneous1D)
875  {
876  int cnt_plane = cnt/n_planes;
877  int e_plane;
878  int e_max_plane = e_max/n_planes;
879  int nTracePts_plane = GetTraceTotPoints();
880 
881  int planeID = floor((e + 0.5 )/ e_max_plane );
882  e_plane = e - e_max_plane*planeID;
883 
884  id2_plane = m_fields[0]->GetTrace()->GetPhys_Offset(
885  m_fields[0]->GetTraceMap()->
886  GetBndCondCoeffsToGlobalCoeffsMap(
887  cnt_plane + e_plane));
888  id2 = id2_plane + planeID*nTracePts_plane;
889  }
890  else // For general case
891  {
892  id2 = m_fields[0]->
893  GetTrace()->GetPhys_Offset(m_fields[0]->GetTraceMap()->
894  GetBndCondTraceToGlobalTraceMap(cnt++));
895  }
896 
897  Array<OneD,NekDouble> x0(npoints, 0.0);
898  Array<OneD,NekDouble> x1(npoints, 0.0);
899  Array<OneD,NekDouble> x2(npoints, 0.0);
900 
901  m_fields[0]->GetBndCondExpansions()[bcRegion]->
902  GetExp(e)->GetCoords(x0, x1, x2);
903 
904  // Flow parameters
905  NekDouble c, k, phi, r, J, VV, pp, sint, P, ss;
906  NekDouble J11, J12, J21, J22, det;
907  NekDouble Fx, Fy;
908  NekDouble xi, yi;
909  NekDouble dV;
910  NekDouble dtheta;
911  NekDouble par1;
912  NekDouble theta = M_PI / 4.0;
913  NekDouble kExt = 0.7;
914  NekDouble V = kExt * sin(theta);
915  NekDouble toll = 1.0e-8;
916  NekDouble errV = 1.0;
917  NekDouble errTheta = 1.0;
918  NekDouble gamma = m_gamma;
919  NekDouble gamma_1_2 = (gamma - 1.0) / 2.0;
920 
921  // Loop on all the points of that edge
922  for (int j = 0; j < npoints; j++)
923  {
924 
925  while ((abs(errV) > toll) || (abs(errTheta) > toll))
926  {
927  VV = V * V;
928  sint = sin(theta);
929  c = sqrt(1.0 - gamma_1_2 * VV);
930  k = V / sint;
931  phi = 1.0 / k;
932  pp = phi * phi;
933  J = 1.0 / c + 1.0 / (3.0 * c * c * c) +
934  1.0 / (5.0 * c * c * c * c * c) -
935  0.5 * log((1.0 + c) / (1.0 - c));
936 
937  r = pow(c, 1.0 / gamma_1_2);
938  xi = 1.0 / (2.0 * r) * (1.0 / VV - 2.0 * pp) + J / 2.0;
939  yi = phi / (r * V) * sqrt(1.0 - VV * pp);
940  par1 = 25.0 - 5.0 * VV;
941  ss = sint * sint;
942 
943  Fx = xi - x0[j];
944  Fy = yi - x1[j];
945 
946  J11 = 39062.5 / pow(par1, 3.5) *
947  (1.0 / VV - 2.0 / VV * ss) * V + 1562.5 /
948  pow(par1, 2.5) * (-2.0 / (VV * V) + 4.0 /
949  (VV * V) * ss) + 12.5 / pow(par1, 1.5) * V +
950  312.5 / pow(par1, 2.5) * V + 7812.5 /
951  pow(par1, 3.5) * V - 0.25 *
952  (-1.0 / pow(par1, 0.5) * V / (1.0 - 0.2 *
953  pow(par1, 0.5)) - (1.0 + 0.2 * pow(par1, 0.5)) /
954  pow((1.0 - 0.2 * pow(par1, 0.5)), 2.0) /
955  pow(par1, 0.5) * V) / (1.0 + 0.2 * pow(par1, 0.5)) *
956  (1.0 - 0.2 * pow(par1, 0.5));
957 
958  J12 = -6250.0 / pow(par1, 2.5) / VV * sint * cos(theta);
959  J21 = -6250.0 / (VV * V) * sint / pow(par1, 2.5) *
960  pow((1.0 - ss), 0.5) + 78125.0 / V * sint /
961  pow(par1, 3.5) * pow((1.0 - ss), 0.5);
962 
963  // the matrix is singular when theta = pi/2
964  if (abs(x1[j]) < toll && abs(cos(theta)) < toll)
965  {
966  J22 = -39062.5 / pow(par1, 3.5) / V + 3125 /
967  pow(par1, 2.5) / (VV * V) + 12.5 /
968  pow(par1, 1.5) * V + 312.5 / pow(par1, 2.5) *
969  V + 7812.5 / pow(par1, 3.5) * V - 0.25 *
970  (-1.0 / pow(par1, 0.5) * V / (1.0 - 0.2 *
971  pow(par1, 0.5)) - (1.0 + 0.2 * pow(par1, 0.5)) /
972  pow((1.0 - 0.2 * pow(par1, 0.5)), 2.0) /
973  pow(par1, 0.5) * V) / (1.0 + 0.2 *
974  pow(par1, 0.5)) * (1.0 - 0.2 * pow(par1, 0.5));
975 
976  // dV = -dV/dx * Fx
977  dV = -1.0 / J22 * Fx;
978  dtheta = 0.0;
979  theta = M_PI / 2.0;
980  }
981  else
982  {
983  J22 = 3125.0 / VV * cos(theta) / pow(par1, 2.5) *
984  pow((1.0 - ss), 0.5) - 3125.0 / VV * ss /
985  pow(par1, 2.5) / pow((1.0 - ss), 0.5) *
986  cos(theta);
987 
988  det = -1.0 / (J11 * J22 - J12 * J21);
989 
990  // [dV dtheta]' = -[invJ]*[Fx Fy]'
991  dV = det * ( J22 * Fx - J12 * Fy);
992  dtheta = det * (-J21 * Fx + J11 * Fy);
993  }
994 
995  V = V + dV;
996  theta = theta + dtheta;
997 
998  errV = abs(dV);
999  errTheta = abs(dtheta);
1000  }
1001 
1002  c = sqrt(1.0 - gamma_1_2 * VV);
1003  int kk = id2 + j;
1004  NekDouble timeramp = 200.0;
1005  std::string restartstr = "RESTART";
1006  if (time<timeramp &&
1007  !(m_session->DefinesFunction("InitialConditions") &&
1008  m_session->GetFunctionType("InitialConditions", 0) ==
1009  LibUtilities::eFunctionTypeFile))
1010  {
1011  Fwd[0][kk] = pow(c, 1.0 / gamma_1_2) *
1012  exp(-1.0 + time /timeramp);
1013 
1014  Fwd[1][kk] = Fwd[0][kk] * V * cos(theta) *
1015  exp(-1 + time / timeramp);
1016 
1017  Fwd[2][kk] = Fwd[0][kk] * V * sin(theta) *
1018  exp(-1 + time / timeramp);
1019  }
1020  else
1021  {
1022  Fwd[0][kk] = pow(c, 1.0 / gamma_1_2);
1023  Fwd[1][kk] = Fwd[0][kk] * V * cos(theta);
1024  Fwd[2][kk] = Fwd[0][kk] * V * sin(theta);
1025  }
1026 
1027  P = (c * c) * Fwd[0][kk] / gamma;
1028  Fwd[3][kk] = P / (gamma - 1.0) + 0.5 *
1029  (Fwd[1][kk] * Fwd[1][kk] / Fwd[0][kk] +
1030  Fwd[2][kk] * Fwd[2][kk] / Fwd[0][kk]);
1031 
1032  errV = 1.0;
1033  errTheta = 1.0;
1034  theta = M_PI / 4.0;
1035  V = kExt * sin(theta);
1036  }
1037 
1038  for (int i = 0; i < nvariables; ++i)
1039  {
1040  Vmath::Vcopy(npoints, &Fwd[i][id2], 1,
1041  &(m_fields[i]->GetBndCondExpansions()[bcRegion]->
1042  UpdatePhys())[id1],1);
1043  }
1044  }
1045  }
1046 }
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:360
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:290
Nektar::SolverUtils::EquationSystem::m_projectionType
enum MultiRegions::ProjectionType m_projectionType
Type of projection; e.g continuous or discontinuous.
Definition: EquationSystem.h:514
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:591
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: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::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:324
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::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:835
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:334
Nektar::EulerCFE::GetExactRinglebFlow
void GetExactRinglebFlow(int field, Array< OneD, NekDouble > &outarray)
Ringleb Flow Test Case.
Definition: EulerCFE.cpp:440
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:265
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::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::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:390
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 Sat Jun 4 2016 00:02:20 for Nektar++ by   doxygen 1.8.8