Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
APE.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File APE.cpp
4 //
5 // For more information, please see: http://www.nektar.info
6 //
7 // The MIT License
8 //
9 // Copyright (c) 2015 Kilian Lackhove
10 // Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
11 // Department of Aeronautics, Imperial College London (UK), and Scientific
12 // Computing and Imaging Institute, University of Utah (USA).
13 //
14 // License for the specific language governing rights and limitations under
15 // Permission is hereby granted, free of charge, to any person obtaining a
16 // copy of this software and associated documentation files (the "Software"),
17 // to deal in the Software without restriction, including without limitation
18 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
19 // and/or sell copies of the Software, and to permit persons to whom the
20 // Software is furnished to do so, subject to the following conditions:
21 //
22 // The above copyright notice and this permission notice shall be included
23 // in all copies or substantial portions of the Software.
24 //
25 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
26 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
27 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
28 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
29 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
30 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
31 // DEALINGS IN THE SOFTWARE.
32 //
33 // Description: APE1/APE4 (Acoustic Perturbation Equations)
34 //
35 ///////////////////////////////////////////////////////////////////////////////
36 
37 #include <iostream>
38 
41 #include <LocalRegions/TriExp.h>
42 #include <LocalRegions/QuadExp.h>
43 #include <LocalRegions/HexExp.h>
44 
48 
49 using namespace std;
50 
51 namespace Nektar
52 {
53 string APE::className = GetEquationSystemFactory().RegisterCreatorFunction(
54  "APE", APE::create,
55  "APE1/APE4 (Acoustic Perturbation Equations)");
56 
57 
58 APE::APE(
60  : UnsteadySystem(pSession)
61 {
62 }
63 
64 
65 /**
66  * @brief Initialization object for the APE class.
67  */
69 {
71 
72  // TODO: We have a bug somewhere in the 1D boundary conditions. Therefore 1D
73  // problems are currently disabled. This should get fixed in the future.
74  ASSERTL0(m_spacedim > 1, "1D problems currently not supported by the APE class.");
75 
77  "Only Projection=DisContinuous supported by the APE class.");
78 
79  // Load isentropic coefficient, Ratio of specific heats
80  m_session->LoadParameter("Gamma", m_gamma, 1.4);
81 
82  m_session->LoadParameter("IO_CFLSteps", m_cflsteps, 0);
83 
84  // Define Baseflow and source term fields
85  switch (m_spacedim)
86  {
87  case 1:
88  {
89  for (int i = 0; i < m_spacedim + 2; ++i)
90  {
93  }
94  break;
95  }
96 
97  case 2:
98  {
99  for (int i = 0; i < m_spacedim + 2; ++i)
100  {
103  }
104  break;
105  }
106 
107  case 3:
108  {
109  for (int i = 0; i < m_spacedim + 2; ++i)
110  {
112  AllocateSharedPtr(m_session, m_graph);
113  }
114  break;
115  }
116 
117  default:
118  {
119 
120  ASSERTL0(false, "Expansion dimension not recognised");
121  break;
122  }
123  }
124 
125  m_bfNames.push_back("p0");
126  m_bfNames.push_back("rho0");
127  m_bfNames.push_back("u0");
128  m_bfNames.push_back("v0");
129  m_bfNames.push_back("w0");
130 
131  // Resize the advection velocities vector to dimension of the problem
132  m_bfNames.resize(m_spacedim + 2);
133 
134  // Initialize basefield
136  for (int i = 0; i < m_bf.num_elements(); ++i)
137  {
139  }
140  EvaluateFunction(m_bfNames, m_bf, "Baseflow", m_time);
141 
143 
144  // Do not forwards transform initial condition
145  m_homoInitialFwd = false;
146 
147  // Define the normal velocity fields
148  if (m_fields[0]->GetTrace())
149  {
151  for (int i = 0; i < m_spacedim + 2; i++)
152  {
154  }
155  }
156 
157  // Set up locations of velocity and base velocity vectors.
160  for (int i = 0; i < m_spacedim; ++i)
161  {
162  // u', v', w'
163  m_vecLocs[0][i] = 1 + i;
164  }
165 
166  string riemName;
167  m_session->LoadSolverInfo("UpwindType", riemName, "APEUpwind");
169  riemName);
170  m_riemannSolver->SetVector("N", &APE::GetNormals, this);
171  m_riemannSolver->SetVector("basefield", &APE::GetBasefield, this);
172  m_riemannSolver->SetAuxVec("vecLocs", &APE::GetVecLocs, this);
173  m_riemannSolver->SetParam("Gamma", &APE::GetGamma, this);
174 
175  // Set up advection operator
176  string advName;
177  m_session->LoadSolverInfo("AdvectionType", advName, "WeakDG");
179  .CreateInstance(advName, advName);
180  m_advection->SetFluxVector(&APE::GetFluxVector, this);
181  m_advection->SetRiemannSolver(m_riemannSolver);
182  m_advection->InitObject(m_session, m_fields);
183 
185  {
188  }
189  else
190  {
191  ASSERTL0(false, "Implicit APE not set up.");
192  }
193 }
194 
195 
196 /**
197  * @brief Destructor for APE class.
198  */
200 {
201 
202 }
203 
204 
206 {
207  int nElm = m_fields[0]->GetExpSize();
208  const Array<OneD, int> expOrder = GetNumExpModesPerExp();
209 
210  Array<OneD, NekDouble> cfl(nElm, 0.0);
211  Array<OneD, NekDouble> stdVelocity(nElm, 0.0);
212 
213  // Get standard velocity to compute the time-step limit
214  GetStdVelocity(stdVelocity);
215 
216  // Factors to compute the time-step limit
218  NekDouble cLambda = 0.2; // Spencer book-317
219 
220  // Loop over elements to compute the time-step limit for each element
221  for (int el = 0; el < nElm; ++el)
222  {
223  NekDouble lambdaMax = stdVelocity[el] * cLambda
224  * (expOrder[el] - 1) * (expOrder[el] - 1);
225  cfl[el] = m_timestep * lambdaMax / alpha;
226  }
227 
228  // Get the minimum time-step limit and return the time-step
229  NekDouble maxCFL = Vmath::Vmax(nElm, cfl, 1);
230  m_comm->AllReduce(maxCFL, LibUtilities::ReduceMax);
231 
232  return maxCFL;
233 }
234 
235 
236 /**
237  * @brief Return the flux vector for the APE equations.
238  *
239  * @param physfield Fields.
240  * @param flux Resulting flux. flux[eq][dir][pt]
241  */
243  const Array<OneD, Array<OneD, NekDouble> > &physfield,
245 {
246  int nq = physfield[0].num_elements();
247  Array<OneD, NekDouble> tmp1(nq);
248  Array<OneD, NekDouble> tmp2(nq);
249 
250  ASSERTL1(flux[0].num_elements() == m_spacedim,
251  "Dimension of flux array and velocity array do not match");
252 
253  // F_{adv,p',j} = \gamma p_0 u'_j + p' \bar{u}_j
254  for (int j = 0; j < m_spacedim; ++j)
255  {
256  Vmath::Zero(nq, flux[0][j], 1);
257 
258  // construct \gamma p_0 u'_j term
259  Vmath::Smul(nq, m_gamma, m_bf[0], 1, tmp1, 1);
260  Vmath::Vmul(nq, tmp1, 1, physfield[j+1], 1, tmp1, 1);
261 
262  // construct p' \bar{u}_j term
263  Vmath::Vmul(nq, physfield[0], 1, m_bf[j+2], 1, tmp2, 1);
264 
265  // add both terms
266  Vmath::Vadd(nq, tmp1, 1, tmp2, 1, flux[0][j], 1);
267  }
268 
269  for (int i = 1; i < flux.num_elements(); ++i)
270  {
271  ASSERTL1(flux[i].num_elements() == m_spacedim,
272  "Dimension of flux array and velocity array do not match");
273 
274  // F_{adv,u'_i,j} = (p'/ \bar{rho} + \bar{u}_k u'_k) \delta_{ij}
275  for (int j = 0; j < m_spacedim; ++j)
276  {
277  Vmath::Zero(nq, flux[i][j], 1);
278 
279  if (i - 1 == j)
280  {
281  // contruct p'/ \bar{rho} term
282  Vmath::Vdiv(nq, physfield[0], 1, m_bf[1], 1, flux[i][j], 1);
283 
284  // construct \bar{u}_k u'_k term
285  Vmath::Zero(nq, tmp1, 1);
286  for (int k = 0; k < m_spacedim; ++k)
287  {
288  Vmath::Vvtvp(nq, physfield[k + 1], 1, m_bf[k + 2 ], 1, tmp1, 1, tmp1, 1);
289  }
290 
291  // add terms
292  Vmath::Vadd(nq, flux[i][j], 1, tmp1, 1, flux[i][j], 1);
293  }
294  }
295  }
296 }
297 
298 
299 /**
300  * @brief v_PostIntegrate
301  */
302 bool APE::v_PreIntegrate(int step)
303 {
304  EvaluateFunction(m_bfNames, m_bf, "Baseflow", m_time);
305 
307  std::vector<SolverUtils::ForcingSharedPtr>::const_iterator x;
308  for (x = m_forcing.begin(); x != m_forcing.end(); ++x)
309  {
310  for (int i = 0; i < (*x)->GetForces().num_elements(); ++i)
311  {
312  m_bfField->IProductWRTBase((*x)->GetForces()[i], tmpC);
313  m_bfField->MultiplyByElmtInvMass(tmpC, tmpC);
314  m_bfField->LocalToGlobal(tmpC, tmpC);
315  m_bfField->GlobalToLocal(tmpC, tmpC);
316  m_bfField->BwdTrans(tmpC, (*x)->UpdateForces()[i]);
317  }
318  }
319 
320  for (int i = 0; i < m_spacedim + 2; ++i)
321  {
322  // ensure the field is C0-continuous
323  m_bfField->IProductWRTBase(m_bf[i], tmpC);
324  m_bfField->MultiplyByElmtInvMass(tmpC, tmpC);
325  m_bfField->LocalToGlobal(tmpC, tmpC);
326  m_bfField->GlobalToLocal(tmpC, tmpC);
327  m_bfField->BwdTrans(tmpC, m_bf[i]);
328  }
329 
330  return UnsteadySystem::v_PreIntegrate(step);
331 }
332 
333 
334 /**
335  * @brief v_PostIntegrate
336  */
337 bool APE::v_PostIntegrate(int step)
338 {
339  if (m_cflsteps && !((step + 1) % m_cflsteps))
340  {
341  NekDouble cfl = GetCFLEstimate();
342  if (m_comm->GetRank() == 0)
343  {
344  cout << "CFL: " << cfl << endl;
345  }
346  }
347 
348  return UnsteadySystem::v_PostIntegrate(step);
349 }
350 
351 
352 /**
353  * @brief Compute the right-hand side.
354  */
355 void APE::DoOdeRhs(const Array<OneD, const Array<OneD, NekDouble> >&inarray,
356  Array<OneD, Array<OneD, NekDouble> >&outarray,
357  const NekDouble time)
358 {
359  int nVariables = inarray.num_elements();
360  int nq = GetTotPoints();
361 
362  // WeakDG does not use advVel, so we only provide a dummy array
364  m_advection->Advect(nVariables, m_fields, advVel, inarray, outarray, time);
365 
366  // Negate the LHS terms
367  for (int i = 0; i < nVariables; ++i)
368  {
369  Vmath::Neg(nq, outarray[i], 1);
370  }
371 
372  std::vector<SolverUtils::ForcingSharedPtr>::const_iterator x;
373  for (x = m_forcing.begin(); x != m_forcing.end(); ++x)
374  {
375  (*x)->Apply(m_fields, outarray, outarray, m_time);
376  }
377 }
378 
379 
380 /**
381  * @brief Compute the projection and call the method for imposing the
382  * boundary conditions in case of discontinuous projection.
383  */
385  Array<OneD, Array<OneD, NekDouble> >&outarray,
386  const NekDouble time)
387 {
388  int nvariables = inarray.num_elements();
389  int nq = m_fields[0]->GetNpoints();
390 
391  // deep copy
392  for (int i = 0; i < nvariables; ++i)
393  {
394  Vmath::Vcopy(nq, inarray[i], 1, outarray[i], 1);
395  }
396 
397  SetBoundaryConditions(outarray, time);
398 }
399 
400 
401 /**
402  * @brief Apply the Boundary Conditions to the APE equations.
403  */
405  NekDouble time)
406 {
407  std::string varName;
408  int nvariables = m_fields.num_elements();
409  int cnt = 0;
410  int nTracePts = GetTraceTotPoints();
411 
412  // Extract trace for boundaries. Needs to be done on all processors to avoid
413  // deadlock.
414  Array<OneD, Array<OneD, NekDouble> > Fwd(nvariables);
415  for (int i = 0; i < nvariables; ++i)
416  {
417  Fwd[i] = Array<OneD, NekDouble>(nTracePts);
418  m_fields[i]->ExtractTracePhys(inarray[i], Fwd[i]);
419  }
420 
421  // loop over Boundary Regions
422  for(int n = 0; n < m_fields[0]->GetBndConditions().num_elements(); ++n)
423  {
424  // Wall Boundary Condition
425  if (boost::iequals(m_fields[0]->GetBndConditions()[n]->GetUserDefined(),"Wall"))
426  {
427  WallBC(n, cnt, Fwd, inarray);
428  }
429 
430  // Time Dependent Boundary Condition (specified in meshfile)
431  if (m_fields[0]->GetBndConditions()[n]->IsTimeDependent())
432  {
433  for (int i = 0; i < nvariables; ++i)
434  {
435  varName = m_session->GetVariable(i);
436  m_fields[i]->EvaluateBoundaryConditions(time, varName);
437  }
438  }
439  cnt +=m_fields[0]->GetBndCondExpansions()[n]->GetExpSize();
440  }
441 }
442 
443 
444 /**
445  * @brief Wall boundary conditions for the APE equations.
446  */
447 void APE::WallBC(int bcRegion, int cnt,
449  Array<OneD, Array<OneD, NekDouble> > &physarray)
450 {
451  int nVariables = physarray.num_elements();
452 
453  const Array<OneD, const int> &traceBndMap = m_fields[0]->GetTraceBndMap();
454 
455  // Adjust the physical values of the trace to take
456  // user defined boundaries into account
457  int id1, id2, nBCEdgePts;
458  int eMax = m_fields[0]->GetBndCondExpansions()[bcRegion]->GetExpSize();
459 
460  for (int e = 0; e < eMax; ++e)
461  {
462  nBCEdgePts = m_fields[0]->GetBndCondExpansions()[bcRegion]->
463  GetExp(e)->GetTotPoints();
464  id1 = m_fields[0]->GetBndCondExpansions()[bcRegion]->GetPhys_Offset(e);
465  id2 = m_fields[0]->GetTrace()->GetPhys_Offset(traceBndMap[cnt+e]);
466 
467  // For 2D/3D, define: v* = v - 2(v.n)n
468  Array<OneD, NekDouble> tmp(nBCEdgePts, 0.0);
469 
470  // Calculate (v.n)
471  for (int i = 0; i < m_spacedim; ++i)
472  {
473  Vmath::Vvtvp(nBCEdgePts,
474  &Fwd[1+i][id2], 1,
475  &m_traceNormals[i][id2], 1,
476  &tmp[0], 1,
477  &tmp[0], 1);
478  }
479 
480  // Calculate 2.0(v.n)
481  Vmath::Smul(nBCEdgePts, -2.0, &tmp[0], 1, &tmp[0], 1);
482 
483  // Calculate v* = v - 2.0(v.n)n
484  for (int i = 0; i < m_spacedim; ++i)
485  {
486  Vmath::Vvtvp(nBCEdgePts,
487  &tmp[0], 1,
488  &m_traceNormals[i][id2], 1,
489  &Fwd[1+i][id2], 1,
490  &Fwd[1+i][id2], 1);
491  }
492 
493  // Copy boundary adjusted values into the boundary expansion
494  for (int i = 0; i < nVariables; ++i)
495  {
496  Vmath::Vcopy(nBCEdgePts,
497  &Fwd[i][id2], 1,
498  &(m_fields[i]->GetBndCondExpansions()[bcRegion]->UpdatePhys())[id1], 1);
499  }
500  }
501 }
502 
503 
504 /**
505  * @brief Compute the advection velocity in the standard space
506  * for each element of the expansion.
507  *
508  * @param stdV Standard velocity field.
509  */
511 {
512  int nElm = m_fields[0]->GetExpSize();
513 
514  ASSERTL1(stdV.num_elements() == nElm, "stdV malformed");
515 
519 
520  // Zero output array
521  Vmath::Zero(stdV.num_elements(), stdV, 1);
522 
523  int cnt = 0;
524 
525  for (int el = 0; el < nElm; ++el)
526  {
527  ptsKeys = m_fields[0]->GetExp(el)->GetPointsKeys();
528 
529  // Possible bug: not multiply by jacobian??
530  const SpatialDomains::GeomFactorsSharedPtr metricInfo =
531  m_fields[0]->GetExp(el)->GetGeom()->GetMetricInfo();
532  const Array<TwoD, const NekDouble> &gmat =
533  m_fields[0]->GetExp(el)->GetGeom()->GetMetricInfo()
534  ->GetDerivFactors(ptsKeys);
535 
536  int nq = m_fields[0]->GetExp(el)->GetTotPoints();
537 
538  for (int i = 0; i < m_spacedim; ++i)
539  {
540  stdVelocity[i] = Array<OneD, NekDouble>(nq, 0.0);
541 
542  velocity[i] = Array<OneD, NekDouble>(nq, 0.0);
543  for (int j = 0; j < nq; ++j)
544  {
545  // The total advection velocity is v+c, so we need to scale c by
546  // adding it before we do the transformation.
547  NekDouble c = sqrt(m_gamma * m_bf[0][cnt+j] / m_bf[1][cnt+j]);
548  velocity[i][j] = m_bf[i+2][cnt+j] + c;
549  }
550  }
551 
552  // scale the velocity components
553  if (metricInfo->GetGtype() == SpatialDomains::eDeformed)
554  {
555  // d xi/ dx = gmat = 1/J * d x/d xi
556  for (int i = 0; i < m_spacedim; ++i)
557  {
558  Vmath::Vmul(nq, gmat[i], 1, velocity[0], 1, stdVelocity[i], 1);
559  for (int j = 1; j < m_spacedim; ++j)
560  {
561  Vmath::Vvtvp(nq, gmat[m_spacedim * j + i], 1, velocity[j],
562  1, stdVelocity[i], 1, stdVelocity[i], 1);
563  }
564  }
565  }
566  else
567  {
568  for (int i = 0; i < m_spacedim; ++i)
569  {
570  Vmath::Smul(nq, gmat[i][0], velocity[0], 1, stdVelocity[i], 1);
571  for (int j = 1; j < m_spacedim; ++j)
572  {
573  Vmath::Svtvp(nq, gmat[m_spacedim * j + i][0], velocity[j],
574  1, stdVelocity[i], 1, stdVelocity[i], 1);
575  }
576  }
577  }
578 
579  // compute the max absolute velocity of the element
580  for (int i = 0; i < nq; ++i)
581  {
582  NekDouble pntVelocity = 0.0;
583  for (int j = 0; j < m_spacedim; ++j)
584  {
585  pntVelocity += stdVelocity[j][i] * stdVelocity[j][i];
586  }
587  pntVelocity = sqrt(pntVelocity);
588 
589  if (pntVelocity > stdV[el])
590  {
591  stdV[el] = pntVelocity;
592  }
593  }
594 
595  cnt += nq;
596  }
597 }
598 
599 
601  std::vector<Array<OneD, NekDouble> > &fieldcoeffs,
602  std::vector<std::string> &variables)
603 {
604  for (int i = 0; i < m_spacedim + 2; i++)
605  {
607 
608  // ensure the field is C0-continuous
609  m_bfField->IProductWRTBase(m_bf[i], tmpC);
610  m_bfField->MultiplyByElmtInvMass(tmpC, tmpC);
611  m_bfField->LocalToGlobal(tmpC, tmpC);
612  m_bfField->GlobalToLocal(tmpC, tmpC);
613 
614  variables.push_back(m_bfNames[i]);
615  fieldcoeffs.push_back(tmpC);
616  }
617 
618  int f = 0;
619  std::vector<SolverUtils::ForcingSharedPtr>::const_iterator x;
620  for (x = m_forcing.begin(); x != m_forcing.end(); ++x)
621  {
622  for (int i = 0; i < (*x)->GetForces().num_elements(); ++i)
623  {
625 
626  m_bfField->IProductWRTBase((*x)->GetForces()[i], tmpC);
627  m_bfField->MultiplyByElmtInvMass(tmpC, tmpC);
628  m_bfField->LocalToGlobal(tmpC, tmpC);
629  m_bfField->GlobalToLocal(tmpC, tmpC);
630 
631  variables.push_back("F_" + boost::lexical_cast<string>(f) +
632  "_" + m_session->GetVariable(i));
633  fieldcoeffs.push_back(tmpC);
634  }
635  f++;
636  }
637 }
638 
639 
640 /**
641  * @brief Get the normal vectors.
642  */
644 {
645  return m_traceNormals;
646 }
647 
648 
649 /**
650  * @brief Get the locations of the components of the directed fields within the fields array.
651  */
653 {
654  return m_vecLocs;
655 }
656 
657 
658 /**
659  * @brief Get the baseflow field.
660  */
662 {
663  for (int i = 0; i < m_spacedim + 2; i++)
664  {
665  m_fields[0]->ExtractTracePhys(m_bf[i], m_traceBasefield[i]);
666  }
667  return m_traceBasefield;
668 }
669 
670 
671 /**
672  * @brief Get the heat capacity ratio.
673  */
675 {
676  return m_gamma;
677 }
678 
679 } //end of namespace
680 
void WallBC(int bcRegion, int cnt, Array< OneD, Array< OneD, NekDouble > > &Fwd, Array< OneD, Array< OneD, NekDouble > > &physarray)
Wall boundary conditions for the APE equations.
Definition: APE.cpp:447
const Array< OneD, const Array< OneD, NekDouble > > & GetNormals()
Get the normal vectors.
Definition: APE.cpp:643
virtual bool v_PostIntegrate(int step)
v_PostIntegrate
Definition: APE.cpp:337
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
virtual SOLVER_UTILS_EXPORT bool v_PreIntegrate(int step)
SolverUtils::AdvectionSharedPtr m_advection
Definition: APE.h:75
std::vector< PointsKey > PointsKeyVector
Definition: Points.h:242
bool m_homoInitialFwd
Flag to determine if simulation should start in homogeneous forward transformed state.
tBaseSharedPtr CreateInstance(tKey idKey BOOST_PP_COMMA_IF(MAX_PARAM) BOOST_PP_ENUM_BINARY_PARAMS(MAX_PARAM, tParam, x))
Create an instance of the class referred to by idKey.
Definition: NekFactory.hpp:162
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
MultiRegions::ExpListSharedPtr m_bfField
Definition: APE.h:83
NekDouble m_time
Current time of simulation.
T Vmax(int n, const T *x, const int incx)
Return the maximum element in x – called vmax to avoid conflict with max.
Definition: Vmath.cpp:779
LibUtilities::TimeIntegrationSchemeOperators m_ode
The time integration scheme operators to use.
NekDouble m_timestep
Time step size.
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
SOLVER_UTILS_EXPORT const Array< OneD, int > GetNumExpModesPerExp()
Array< OneD, Array< OneD, NekDouble > > m_vecLocs
Definition: APE.h:79
void Svtvp(int n, const T alpha, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
svtvp (scalar times vector plus vector): z = alpha*x + y
Definition: Vmath.cpp:485
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
Definition: Vmath.cpp:442
enum MultiRegions::ProjectionType m_projectionType
Type of projection; e.g continuous or discontinuous.
STL namespace.
virtual bool v_PreIntegrate(int step)
v_PostIntegrate
Definition: APE.cpp:302
std::vector< SolverUtils::ForcingSharedPtr > m_forcing
Definition: APE.h:76
static SOLVER_UTILS_EXPORT std::vector< ForcingSharedPtr > Load(const LibUtilities::SessionReaderSharedPtr &pSession, const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const unsigned int &pNumForcingFields=0)
Definition: Forcing.cpp:86
void Vdiv(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x/y.
Definition: Vmath.cpp:241
SOLVER_UTILS_EXPORT NekDouble MaxTimeStepEstimator()
Get the maximum timestep estimator for cfl control.
void GetFluxVector(const Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &flux)
Return the flux vector for the APE equations.
Definition: APE.cpp:242
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:51
LibUtilities::CommSharedPtr m_comm
Communicator.
SOLVER_UTILS_EXPORT int GetTotPoints()
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Array holding trace normals for DG simulations in the forwards direction.
void DefineProjection(FuncPointerT func, ObjectPointerT obj)
NekDouble GetGamma()
Get the heat capacity ratio.
Definition: APE.cpp:674
virtual void v_ExtraFldOutput(std::vector< Array< OneD, NekDouble > > &fieldcoeffs, std::vector< std::string > &variables)
Definition: APE.cpp:600
void GetStdVelocity(Array< OneD, NekDouble > &stdV)
Compute the advection velocity in the standard space for each element of the expansion.
Definition: APE.cpp:510
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
Definition: Vmath.cpp:213
bool m_explicitAdvection
Indicates if explicit or implicit treatment of advection is used.
const Array< OneD, const Array< OneD, NekDouble > > & GetVecLocs()
Get the locations of the components of the directed fields within the fields array.
Definition: APE.cpp:652
Array< OneD, Array< OneD, NekDouble > > m_traceBasefield
Definition: APE.h:78
virtual SOLVER_UTILS_EXPORT bool v_PostIntegrate(int step)
void DefineOdeRhs(FuncPointerT func, ObjectPointerT obj)
virtual ~APE()
Destructor.
Definition: APE.cpp:199
Base class for unsteady solvers.
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: APE.cpp:355
RiemannSolverFactory & GetRiemannSolverFactory()
int m_cflsteps
dump cfl estimate
Definition: APE.h:86
int m_spacedim
Spatial dimension (>= expansion dim).
AdvectionFactory & GetAdvectionFactory()
Gets the factory for initialising advection objects.
Definition: Advection.cpp:46
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:396
double NekDouble
virtual void v_InitObject()
Initialization object for the APE class.
Definition: APE.cpp:68
virtual SOLVER_UTILS_EXPORT void v_InitObject()
Init object for UnsteadySystem class.
void SetBoundaryConditions(Array< OneD, Array< OneD, NekDouble > > &physarray, NekDouble time)
Apply the Boundary Conditions to the APE equations.
Definition: APE.cpp:404
SOLVER_UTILS_EXPORT void EvaluateFunction(Array< OneD, Array< OneD, NekDouble > > &pArray, std::string pFunctionName, const NekDouble pTime=0.0, const int domain=0)
Evaluates a function as specified in the session file.
NekDouble m_gamma
Isentropic coefficient, Ratio of specific heats (APE)
Definition: APE.h:81
const Array< OneD, const Array< OneD, NekDouble > > & GetBasefield()
Get the baseflow field.
Definition: APE.cpp:661
EquationSystemFactory & GetEquationSystemFactory()
SOLVER_UTILS_EXPORT int GetTraceTotPoints()
SolverUtils::RiemannSolverSharedPtr m_riemannSolver
Definition: APE.h:77
NekDouble GetCFLEstimate()
Definition: APE.cpp:205
boost::shared_ptr< GeomFactors > GeomFactorsSharedPtr
Pointer to a GeomFactors object.
Definition: GeomFactors.h:62
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
SOLVER_UTILS_EXPORT int GetTraceNpoints()
SOLVER_UTILS_EXPORT int GetNcoeffs()
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: APE.cpp:384
SpatialDomains::MeshGraphSharedPtr m_graph
Pointer to graph defining mesh.
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:373
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:228
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1061
Geometry is curved or has non-constant factors.
Array< OneD, Array< OneD, NekDouble > > m_bf
Definition: APE.h:82
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.cpp:299
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.cpp:183
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215
std::vector< std::string > m_bfNames
Definition: APE.h:84