Nektar++
CompressibleFlowSystem.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File CompressibleFlowSystem.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 // Permission is hereby granted, free of charge, to any person obtaining a
14 // copy of this software and associated documentation files (the "Software"),
15 // to deal in the Software without restriction, including without limitation
16 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
17 // and/or sell copies of the Software, and to permit persons to whom the
18 // Software is furnished to do so, subject to the following conditions:
19 //
20 // The above copyright notice and this permission notice shall be included
21 // in all copies or substantial portions of the Software.
22 //
23 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29 // DEALINGS IN THE SOFTWARE.
30 //
31 // Description: Compressible flow system base class with auxiliary functions
32 //
33 ///////////////////////////////////////////////////////////////////////////////
34 
35 #include <boost/core/ignore_unused.hpp>
36 
38 
40 
41 using namespace std;
42 
43 namespace Nektar
44 {
45 CompressibleFlowSystem::CompressibleFlowSystem(
48  : UnsteadySystem(pSession, pGraph), AdvectionSystem(pSession, pGraph)
49 {
50 }
51 
52 /**
53  * @brief Initialization object for CompressibleFlowSystem class.
54  */
55 void CompressibleFlowSystem::v_InitObject(bool DeclareFields)
56 {
57  AdvectionSystem::v_InitObject(DeclareFields);
58 
59  for (int i = 0; i < m_fields.size(); i++)
60  {
61  // Use BwdTrans to make sure initial condition is in solution space
62  m_fields[i]->BwdTrans(m_fields[i]->GetCoeffs(),
63  m_fields[i]->UpdatePhys());
64  }
65 
68 
69  ASSERTL0(m_session->DefinesSolverInfo("UPWINDTYPE"),
70  "No UPWINDTYPE defined in session.");
71 
72  // Do not forwards transform initial condition
73  m_homoInitialFwd = false;
74 
75  // Set up locations of velocity vector.
78  for (int i = 0; i < m_spacedim; ++i)
79  {
80  m_vecLocs[0][i] = 1 + i;
81  }
82 
83  // Loading parameters from session file
85 
86  // Setting up advection and diffusion operators
87  InitAdvection();
88 
89  // Create artificial diffusion with laplacian operator
90  if (m_shockCaptureType == "NonSmooth")
91  {
94  }
95 
96  // Forcing terms for the sponge region
97  m_forcing = SolverUtils::Forcing::Load(m_session, shared_from_this(),
98  m_fields, m_fields.size());
99 
100  // User-defined boundary conditions
101  int cnt = 0;
102  for (int n = 0; n < m_fields[0]->GetBndConditions().size(); ++n)
103  {
104  std::string type = m_fields[0]->GetBndConditions()[n]->GetUserDefined();
105 
106  if (m_fields[0]->GetBndConditions()[n]->GetBoundaryConditionType() ==
108  {
109  continue;
110  }
111 
112  if (!type.empty())
113  {
114  m_bndConds.push_back(GetCFSBndCondFactory().CreateInstance(
115  type, m_session, m_fields, m_traceNormals, m_spacedim, n, cnt));
116  }
117  cnt += m_fields[0]->GetBndCondExpansions()[n]->GetExpSize();
118  }
119 
122 
124 }
125 
126 /**
127  * @brief Destructor for CompressibleFlowSystem class.
128  */
130 {
131 }
132 
133 /**
134  * @brief Load CFS parameters from the session file
135  */
137 {
138  // Get gamma parameter from session file.
139  m_session->LoadParameter("Gamma", m_gamma, 1.4);
140 
141  // Shock capture
142  m_session->LoadSolverInfo("ShockCaptureType", m_shockCaptureType, "Off");
143 
144  // Load parameters for exponential filtering
145  m_session->MatchSolverInfo("ExponentialFiltering", "True", m_useFiltering,
146  false);
147  if (m_useFiltering)
148  {
149  m_session->LoadParameter("FilterAlpha", m_filterAlpha, 36);
150  m_session->LoadParameter("FilterExponent", m_filterExponent, 16);
151  m_session->LoadParameter("FilterCutoff", m_filterCutoff, 0);
152  }
153 
154  // Load CFL for local time-stepping (for steady state)
155  m_session->MatchSolverInfo("LocalTimeStep", "True", m_useLocalTimeStep,
156  false);
157  if (m_useLocalTimeStep)
158  {
160  "Local time stepping requires CFL parameter.");
161  }
162 }
163 
164 /**
165  * @brief Create advection and diffusion objects for CFS
166  */
168 {
169  // Check if projection type is correct
171  "Unsupported projection type.");
172 
173  string advName, riemName;
174  m_session->LoadSolverInfo("AdvectionType", advName, "WeakDG");
175 
176  m_advObject =
178 
180  {
181  m_advObject->SetFluxVector(
183  }
184  else
185  {
187  this);
188  }
189 
190  // Setting up Riemann solver for advection operator
191  m_session->LoadSolverInfo("UpwindType", riemName, "Average");
192 
195  riemName, m_session);
196 
197  // Setting up parameters for advection operator Riemann solver
198  riemannSolver->SetParam("gamma", &CompressibleFlowSystem::GetGamma, this);
199  riemannSolver->SetAuxVec("vecLocs", &CompressibleFlowSystem::GetVecLocs,
200  this);
201  riemannSolver->SetVector("N", &CompressibleFlowSystem::GetNormals, this);
202 
203  // Concluding initialisation of advection / diffusion operators
204  m_advObject->SetRiemannSolver(riemannSolver);
205  m_advObject->InitObject(m_session, m_fields);
206 }
207 
208 /**
209  * @brief Compute the right-hand side.
210  */
212  const Array<OneD, const Array<OneD, NekDouble>> &inarray,
213  Array<OneD, Array<OneD, NekDouble>> &outarray, const NekDouble time)
214 {
215  int nvariables = inarray.size();
216  int npoints = GetNpoints();
217  int nTracePts = GetTraceTotPoints();
218 
219  m_bndEvaluateTime = time;
220 
221  // Store forwards/backwards space along trace space
222  Array<OneD, Array<OneD, NekDouble>> Fwd(nvariables);
223  Array<OneD, Array<OneD, NekDouble>> Bwd(nvariables);
224 
226  {
229  }
230  else
231  {
232  for (int i = 0; i < nvariables; ++i)
233  {
234  Fwd[i] = Array<OneD, NekDouble>(nTracePts, 0.0);
235  Bwd[i] = Array<OneD, NekDouble>(nTracePts, 0.0);
236  m_fields[i]->GetFwdBwdTracePhys(inarray[i], Fwd[i], Bwd[i]);
237  }
238  }
239 
240  // Calculate advection
241  LibUtilities::Timer timer;
242  timer.Start();
243  DoAdvection(inarray, outarray, time, Fwd, Bwd);
244  timer.Stop();
245  timer.AccumulateRegion("DoAdvection");
246 
247  // Negate results
248  for (int i = 0; i < nvariables; ++i)
249  {
250  Vmath::Neg(npoints, outarray[i], 1);
251  }
252 
253  // Add diffusion terms
254  timer.Start();
255  DoDiffusion(inarray, outarray, Fwd, Bwd);
256  timer.Stop();
257  timer.AccumulateRegion("DoDiffusion");
258 
259  // Add forcing terms
260  for (auto &x : m_forcing)
261  {
262  x->Apply(m_fields, inarray, outarray, time);
263  }
264 
265  if (m_useLocalTimeStep)
266  {
267  int nElements = m_fields[0]->GetExpSize();
268  int nq, offset;
269  NekDouble fac;
271 
272  Array<OneD, NekDouble> tstep(nElements, 0.0);
273  GetElmtTimeStep(inarray, tstep);
274 
275  // Loop over elements
276  for (int n = 0; n < nElements; ++n)
277  {
278  nq = m_fields[0]->GetExp(n)->GetTotPoints();
279  offset = m_fields[0]->GetPhys_Offset(n);
280  fac = tstep[n] / m_timestep;
281  for (int i = 0; i < nvariables; ++i)
282  {
283  Vmath::Smul(nq, fac, outarray[i] + offset, 1,
284  tmp = outarray[i] + offset, 1);
285  }
286  }
287  }
288 }
289 
290 /**
291  * @brief Compute the projection and call the method for imposing the
292  * boundary conditions in case of discontinuous projection.
293  */
295  const Array<OneD, const Array<OneD, NekDouble>> &inarray,
296  Array<OneD, Array<OneD, NekDouble>> &outarray, const NekDouble time)
297 {
298  int nvariables = inarray.size();
299 
300  switch (m_projectionType)
301  {
303  {
304  // Just copy over array
305  int npoints = GetNpoints();
306 
307  for (int i = 0; i < nvariables; ++i)
308  {
309  Vmath::Vcopy(npoints, inarray[i], 1, outarray[i], 1);
310  if (m_useFiltering)
311  {
312  m_fields[i]->ExponentialFilter(outarray[i], m_filterAlpha,
315  }
316  }
317  SetBoundaryConditions(outarray, time);
318  break;
319  }
322  {
323  NEKERROR(ErrorUtil::efatal, "No Continuous Galerkin for full "
324  "compressible Navier-Stokes equations");
325  break;
326  }
327  default:
328  NEKERROR(ErrorUtil::efatal, "Unknown projection scheme");
329  break;
330  }
331 }
332 
333 /**
334  * @brief Compute the advection terms for the right-hand side
335  */
337  const Array<OneD, Array<OneD, NekDouble>> &inarray,
338  Array<OneD, Array<OneD, NekDouble>> &outarray, const NekDouble time,
339  const Array<OneD, Array<OneD, NekDouble>> &pFwd,
340  const Array<OneD, Array<OneD, NekDouble>> &pBwd)
341 {
342  int nvariables = inarray.size();
344 
345  m_advObject->Advect(nvariables, m_fields, advVel, inarray, outarray, time,
346  pFwd, pBwd);
347 }
348 
349 /**
350  * @brief Add the diffusions terms to the right-hand side
351  */
353  const Array<OneD, Array<OneD, NekDouble>> &inarray,
354  Array<OneD, Array<OneD, NekDouble>> &outarray,
355  const Array<OneD, Array<OneD, NekDouble>> &pFwd,
356  const Array<OneD, Array<OneD, NekDouble>> &pBwd)
357 {
358  v_DoDiffusion(inarray, outarray, pFwd, pBwd);
359 }
360 
362  Array<OneD, Array<OneD, NekDouble>> &physarray, NekDouble time)
363 {
364  int nTracePts = GetTraceTotPoints();
365  int nvariables = physarray.size();
366 
367  Array<OneD, Array<OneD, NekDouble>> Fwd(nvariables);
368  for (int i = 0; i < nvariables; ++i)
369  {
370  Fwd[i] = Array<OneD, NekDouble>(nTracePts);
371  m_fields[i]->ExtractTracePhys(physarray[i], Fwd[i]);
372  }
373 
374  if (m_bndConds.size())
375  {
376  // Loop over user-defined boundary conditions
377  for (auto &x : m_bndConds)
378  {
379  x->Apply(Fwd, physarray, time);
380  }
381  }
382 }
383 /**
384  * @brief Set up a weight on physical boundaries for boundary condition
385  * applications
386  */
388 {
389  if (m_bndConds.size())
390  {
391  // Loop over user-defined boundary conditions
392  for (auto &x : m_bndConds)
393  {
394  x->ApplyBwdWeight();
395  }
396  }
397 }
398 
399 /**
400  * @brief Return the flux vector for the compressible Euler equations.
401  *
402  * @param physfield Fields.
403  * @param flux Resulting flux.
404  */
406  const Array<OneD, const Array<OneD, NekDouble>> &physfield,
408 {
409  auto nVariables = physfield.size();
410  auto nPts = physfield[0].size();
411 
412  constexpr unsigned short maxVel = 3;
413  constexpr unsigned short maxFld = 5;
414 
415  // hardcoding done for performance reasons
416  ASSERTL1(nVariables <= maxFld, "GetFluxVector, hard coded max fields");
417 
418  for (size_t p = 0; p < nPts; ++p)
419  {
420  // local storage
421  std::array<NekDouble, maxFld> fieldTmp;
422  std::array<NekDouble, maxVel> velocity;
423 
424  // rearrenge and load data
425  for (size_t f = 0; f < nVariables; ++f)
426  {
427  fieldTmp[f] = physfield[f][p]; // load
428  }
429 
430  // 1 / rho
431  NekDouble oneOrho = 1.0 / fieldTmp[0];
432 
433  for (size_t d = 0; d < m_spacedim; ++d)
434  {
435  // Flux vector for the rho equation
436  flux[0][d][p] = fieldTmp[d + 1]; // store
437  // compute velocity
438  velocity[d] = fieldTmp[d + 1] * oneOrho;
439  }
440 
441  NekDouble pressure = m_varConv->GetPressure(fieldTmp.data());
442  NekDouble ePlusP = fieldTmp[m_spacedim + 1] + pressure;
443  for (size_t f = 0; f < m_spacedim; ++f)
444  {
445  // Flux vector for the velocity fields
446  for (size_t d = 0; d < m_spacedim; ++d)
447  {
448  flux[f + 1][d][p] = velocity[d] * fieldTmp[f + 1]; // store
449  }
450 
451  // Add pressure to appropriate field
452  flux[f + 1][f][p] += pressure;
453 
454  // Flux vector for energy
455  flux[m_spacedim + 1][f][p] = ePlusP * velocity[f]; // store
456  }
457  }
458 }
459 
460 /**
461  * @brief Return the flux vector for the compressible Euler equations
462  * by using the de-aliasing technique.
463  *
464  * @param physfield Fields.
465  * @param flux Resulting flux.
466  */
468  const Array<OneD, const Array<OneD, NekDouble>> &physfield,
470 {
471  int i, j;
472  int nq = physfield[0].size();
473  int nVariables = m_fields.size();
474 
475  // Factor to rescale 1d points in dealiasing
476  NekDouble OneDptscale = 2;
477  nq = m_fields[0]->Get1DScaledTotPoints(OneDptscale);
478 
481 
482  Array<OneD, Array<OneD, NekDouble>> physfield_interp(nVariables);
483  TensorOfArray3D<NekDouble> flux_interp(nVariables);
484 
485  for (i = 0; i < nVariables; ++i)
486  {
487  physfield_interp[i] = Array<OneD, NekDouble>(nq);
489  m_fields[0]->PhysInterp1DScaled(OneDptscale, physfield[i],
490  physfield_interp[i]);
491 
492  for (j = 0; j < m_spacedim; ++j)
493  {
494  flux_interp[i][j] = Array<OneD, NekDouble>(nq);
495  }
496  }
497 
498  // Flux vector for the rho equation
499  for (i = 0; i < m_spacedim; ++i)
500  {
501  velocity[i] = Array<OneD, NekDouble>(nq);
502 
503  // Galerkin project solution back to original space
504  m_fields[0]->PhysGalerkinProjection1DScaled(
505  OneDptscale, physfield_interp[i + 1], flux[0][i]);
506  }
507 
508  m_varConv->GetVelocityVector(physfield_interp, velocity);
509  m_varConv->GetPressure(physfield_interp, pressure);
510 
511  // Evaluation of flux vector for the velocity fields
512  for (i = 0; i < m_spacedim; ++i)
513  {
514  for (j = 0; j < m_spacedim; ++j)
515  {
516  Vmath::Vmul(nq, velocity[j], 1, physfield_interp[i + 1], 1,
517  flux_interp[i + 1][j], 1);
518  }
519 
520  // Add pressure to appropriate field
521  Vmath::Vadd(nq, flux_interp[i + 1][i], 1, pressure, 1,
522  flux_interp[i + 1][i], 1);
523  }
524 
525  // Galerkin project solution back to original space
526  for (i = 0; i < m_spacedim; ++i)
527  {
528  for (j = 0; j < m_spacedim; ++j)
529  {
530  m_fields[0]->PhysGalerkinProjection1DScaled(
531  OneDptscale, flux_interp[i + 1][j], flux[i + 1][j]);
532  }
533  }
534 
535  // Evaluation of flux vector for energy
536  Vmath::Vadd(nq, physfield_interp[m_spacedim + 1], 1, pressure, 1, pressure,
537  1);
538 
539  for (j = 0; j < m_spacedim; ++j)
540  {
541  Vmath::Vmul(nq, velocity[j], 1, pressure, 1,
542  flux_interp[m_spacedim + 1][j], 1);
543 
544  // Galerkin project solution back to original space
545  m_fields[0]->PhysGalerkinProjection1DScaled(
546  OneDptscale, flux_interp[m_spacedim + 1][j],
547  flux[m_spacedim + 1][j]);
548  }
549 }
550 
551 /**
552  * @brief Calculate the maximum timestep on each element
553  * subject to CFL restrictions.
554  */
556  const Array<OneD, const Array<OneD, NekDouble>> &inarray,
557  Array<OneD, NekDouble> &tstep)
558 {
559  boost::ignore_unused(inarray);
560 
561  int nElements = m_fields[0]->GetExpSize();
562 
563  // Change value of m_timestep (in case it is set to zero)
564  NekDouble tmp = m_timestep;
565  m_timestep = 1.0;
566 
567  Array<OneD, NekDouble> cfl(nElements);
568  cfl = GetElmtCFLVals();
569 
570  // Factors to compute the time-step limit
572 
573  // Loop over elements to compute the time-step limit for each element
574  for (int n = 0; n < nElements; ++n)
575  {
576  tstep[n] = m_cflSafetyFactor * alpha / cfl[n];
577  }
578 
579  // Restore value of m_timestep
580  m_timestep = tmp;
581 }
582 
583 /**
584  * @brief Calculate the maximum timestep subject to CFL restrictions.
585  */
587  const Array<OneD, const Array<OneD, NekDouble>> &inarray)
588 {
589  int nElements = m_fields[0]->GetExpSize();
590  Array<OneD, NekDouble> tstep(nElements, 0.0);
591 
592  GetElmtTimeStep(inarray, tstep);
593 
594  // Get the minimum time-step limit and return the time-step
595  NekDouble TimeStep = Vmath::Vmin(nElements, tstep, 1);
596  m_comm->AllReduce(TimeStep, LibUtilities::ReduceMin);
597 
598  NekDouble tmp = m_timestep;
599  m_timestep = TimeStep;
600 
601  Array<OneD, NekDouble> cflNonAcoustic(nElements, 0.0);
602  cflNonAcoustic = GetElmtCFLVals(false);
603 
604  // Get the minimum time-step limit and return the time-step
605  NekDouble MaxcflNonAcoustic = Vmath::Vmax(nElements, cflNonAcoustic, 1);
606  m_comm->AllReduce(MaxcflNonAcoustic, LibUtilities::ReduceMax);
607 
608  m_cflNonAcoustic = MaxcflNonAcoustic;
609  m_timestep = tmp;
610 
611  return TimeStep;
612 }
613 
614 /**
615  * @brief Apply artificial diffusion (Laplacian operator)
616  */
618  const Array<OneD, Array<OneD, NekDouble>> &inarray,
619  Array<OneD, Array<OneD, NekDouble>> &outarray,
620  const Array<OneD, Array<OneD, NekDouble>> &pFwd,
621  const Array<OneD, Array<OneD, NekDouble>> &pBwd)
622 {
623  boost::ignore_unused(pFwd, pBwd);
625  {
626  m_artificialDiffusion->DoArtificialDiffusion(inarray, outarray);
627  }
628 }
629 
630 /**
631  * @brief Set up logic for residual calculation.
632  */
634  bool dumpInitialConditions,
635  const int domain)
636 {
637  boost::ignore_unused(domain);
638 
639  EquationSystem::v_SetInitialConditions(initialtime, false);
640 
641  // insert white noise in initial condition
642  NekDouble Noise;
643  int phystot = m_fields[0]->GetTotPoints();
644  Array<OneD, NekDouble> noise(phystot);
645 
646  m_session->LoadParameter("Noise", Noise, 0.0);
647  int m_nConvectiveFields = m_fields.size();
648 
649  if (Noise > 0.0)
650  {
651  int seed = -m_comm->GetRank() * m_nConvectiveFields;
652  for (int i = 0; i < m_nConvectiveFields; i++)
653  {
654  Vmath::FillWhiteNoise(phystot, Noise, noise, 1, seed);
655  --seed;
656  Vmath::Vadd(phystot, m_fields[i]->GetPhys(), 1, noise, 1,
657  m_fields[i]->UpdatePhys(), 1);
658  m_fields[i]->FwdTransLocalElmt(m_fields[i]->GetPhys(),
659  m_fields[i]->UpdateCoeffs());
660  }
661  }
662 
663  if (dumpInitialConditions && m_checksteps)
664  {
666  m_nchk++;
667  }
668 }
669 
670 /**
671  * @brief Compute the advection velocity in the standard space
672  * for each element of the expansion.
673  */
675  const NekDouble SpeedSoundFactor)
676 {
677  int nTotQuadPoints = GetTotPoints();
678  int n_element = m_fields[0]->GetExpSize();
679  int expdim = m_fields[0]->GetGraph()->GetMeshDimension();
680  int nfields = m_fields.size();
681  int offset;
683 
684  Array<OneD, Array<OneD, NekDouble>> physfields(nfields);
685  for (int i = 0; i < nfields; ++i)
686  {
687  physfields[i] = m_fields[i]->GetPhys();
688  }
689 
690  Array<OneD, NekDouble> stdV(n_element, 0.0);
691 
692  // Getting the velocity vector on the 2D normal space
696  Array<OneD, NekDouble> soundspeed(nTotQuadPoints);
698 
699  for (int i = 0; i < m_spacedim; ++i)
700  {
701  velocity[i] = Array<OneD, NekDouble>(nTotQuadPoints);
702  stdVelocity[i] = Array<OneD, NekDouble>(nTotQuadPoints, 0.0);
703  stdSoundSpeed[i] = Array<OneD, NekDouble>(nTotQuadPoints, 0.0);
704  }
705 
706  m_varConv->GetVelocityVector(physfields, velocity);
707  m_varConv->GetSoundSpeed(physfields, soundspeed);
708 
709  for (int el = 0; el < n_element; ++el)
710  {
711  ptsKeys = m_fields[0]->GetExp(el)->GetPointsKeys();
712  offset = m_fields[0]->GetPhys_Offset(el);
713  int nq = m_fields[0]->GetExp(el)->GetTotPoints();
714 
715  const SpatialDomains::GeomFactorsSharedPtr metricInfo =
716  m_fields[0]->GetExp(el)->GetGeom()->GetMetricInfo();
717  const Array<TwoD, const NekDouble> &gmat =
718  m_fields[0]
719  ->GetExp(el)
720  ->GetGeom()
721  ->GetMetricInfo()
722  ->GetDerivFactors(ptsKeys);
723 
724  // Convert to standard element
725  // consider soundspeed in all directions
726  // (this might overestimate the cfl)
727  if (metricInfo->GetGtype() == SpatialDomains::eDeformed)
728  {
729  // d xi/ dx = gmat = 1/J * d x/d xi
730  for (int i = 0; i < expdim; ++i)
731  {
732  Vmath::Vmul(nq, gmat[i], 1, velocity[0] + offset, 1,
733  tmp = stdVelocity[i] + offset, 1);
734  Vmath::Vmul(nq, gmat[i], 1, soundspeed + offset, 1,
735  tmp = stdSoundSpeed[i] + offset, 1);
736  for (int j = 1; j < expdim; ++j)
737  {
738  Vmath::Vvtvp(nq, gmat[expdim * j + i], 1,
739  velocity[j] + offset, 1,
740  stdVelocity[i] + offset, 1,
741  tmp = stdVelocity[i] + offset, 1);
742  Vmath::Vvtvp(nq, gmat[expdim * j + i], 1,
743  soundspeed + offset, 1,
744  stdSoundSpeed[i] + offset, 1,
745  tmp = stdSoundSpeed[i] + offset, 1);
746  }
747  }
748  }
749  else
750  {
751  for (int i = 0; i < expdim; ++i)
752  {
753  Vmath::Smul(nq, gmat[i][0], velocity[0] + offset, 1,
754  tmp = stdVelocity[i] + offset, 1);
755  Vmath::Smul(nq, gmat[i][0], soundspeed + offset, 1,
756  tmp = stdSoundSpeed[i] + offset, 1);
757  for (int j = 1; j < expdim; ++j)
758  {
759  Vmath::Svtvp(nq, gmat[expdim * j + i][0],
760  velocity[j] + offset, 1,
761  stdVelocity[i] + offset, 1,
762  tmp = stdVelocity[i] + offset, 1);
763  Vmath::Svtvp(nq, gmat[expdim * j + i][0],
764  soundspeed + offset, 1,
765  stdSoundSpeed[i] + offset, 1,
766  tmp = stdSoundSpeed[i] + offset, 1);
767  }
768  }
769  }
770 
771  NekDouble vel;
772  for (int i = 0; i < nq; ++i)
773  {
774  NekDouble pntVelocity = 0.0;
775  for (int j = 0; j < expdim; ++j)
776  {
777  // Add sound speed
778  vel = std::abs(stdVelocity[j][offset + i]) +
779  SpeedSoundFactor * std::abs(stdSoundSpeed[j][offset + i]);
780  pntVelocity += vel * vel;
781  }
782  pntVelocity = sqrt(pntVelocity);
783  if (pntVelocity > stdV[el])
784  {
785  stdV[el] = pntVelocity;
786  }
787  }
788  }
789 
790  return stdV;
791 }
792 
793 /**
794  * @brief Set the denominator to compute the time step when a cfl
795  * control is employed. This function is no longer used but is still
796  * here for being utilised in the future.
797  *
798  * @param n Order of expansion element by element.
799  */
801 {
802  ASSERTL0(n <= 20, "Illegal modes dimension for CFL calculation "
803  "(P has to be less then 20)");
804 
805  NekDouble CFLDG[21] = {2.0000, 6.0000, 11.8424, 19.1569, 27.8419,
806  37.8247, 49.0518, 61.4815, 75.0797, 89.8181,
807  105.6700, 122.6200, 140.6400, 159.7300, 179.8500,
808  201.0100, 223.1800, 246.3600, 270.5300, 295.6900,
809  321.8300}; // CFLDG 1D [0-20]
810  NekDouble CFL = 0.0;
811 
813  {
814  CFL = CFLDG[n];
815  }
816  else
817  {
818  NEKERROR(ErrorUtil::efatal, "Continuous Galerkin stability "
819  "coefficients not introduced yet.");
820  }
821 
822  return CFL;
823 }
824 
825 /**
826  * @brief Compute the vector of denominators to compute the time step
827  * when a cfl control is employed. This function is no longer used but
828  * is still here for being utilised in the future.
829  *
830  * @param ExpOrder Order of expansion element by element.
831  */
833  const Array<OneD, int> &ExpOrder)
834 {
835  int i;
836  Array<OneD, NekDouble> returnval(m_fields[0]->GetExpSize(), 0.0);
837  for (i = 0; i < m_fields[0]->GetExpSize(); i++)
838  {
839  returnval[i] = GetStabilityLimit(ExpOrder[i]);
840  }
841  return returnval;
842 }
843 
845  std::vector<Array<OneD, NekDouble>> &fieldcoeffs,
846  std::vector<std::string> &variables)
847 {
848  bool extraFields;
849  m_session->MatchSolverInfo("OutputExtraFields", "True", extraFields, true);
850  if (extraFields)
851  {
852  const int nPhys = m_fields[0]->GetNpoints();
853  const int nCoeffs = m_fields[0]->GetNcoeffs();
855 
856  for (int i = 0; i < m_fields.size(); ++i)
857  {
858  tmp[i] = m_fields[i]->GetPhys();
859  }
860 
863  for (int i = 0; i < m_spacedim; ++i)
864  {
865  velocity[i] = Array<OneD, NekDouble>(nPhys);
866  velFwd[i] = Array<OneD, NekDouble>(nCoeffs);
867  }
868 
869  Array<OneD, NekDouble> pressure(nPhys), temperature(nPhys);
870  Array<OneD, NekDouble> entropy(nPhys);
871  Array<OneD, NekDouble> soundspeed(nPhys), mach(nPhys);
872  Array<OneD, NekDouble> sensor(nPhys), SensorKappa(nPhys);
873 
874  m_varConv->GetVelocityVector(tmp, velocity);
875  m_varConv->GetPressure(tmp, pressure);
876  m_varConv->GetTemperature(tmp, temperature);
877  m_varConv->GetEntropy(tmp, entropy);
878  m_varConv->GetSoundSpeed(tmp, soundspeed);
879  m_varConv->GetMach(tmp, soundspeed, mach);
880 
881  int sensorOffset;
882  m_session->LoadParameter("SensorOffset", sensorOffset, 1);
883  m_varConv->GetSensor(m_fields[0], tmp, sensor, SensorKappa,
884  sensorOffset);
885 
886  Array<OneD, NekDouble> pFwd(nCoeffs), TFwd(nCoeffs);
887  Array<OneD, NekDouble> sFwd(nCoeffs);
888  Array<OneD, NekDouble> aFwd(nCoeffs), mFwd(nCoeffs);
889  Array<OneD, NekDouble> sensFwd(nCoeffs);
890 
891  string velNames[3] = {"u", "v", "w"};
892  for (int i = 0; i < m_spacedim; ++i)
893  {
894  m_fields[0]->FwdTransLocalElmt(velocity[i], velFwd[i]);
895  variables.push_back(velNames[i]);
896  fieldcoeffs.push_back(velFwd[i]);
897  }
898 
899  m_fields[0]->FwdTransLocalElmt(pressure, pFwd);
900  m_fields[0]->FwdTransLocalElmt(temperature, TFwd);
901  m_fields[0]->FwdTransLocalElmt(entropy, sFwd);
902  m_fields[0]->FwdTransLocalElmt(soundspeed, aFwd);
903  m_fields[0]->FwdTransLocalElmt(mach, mFwd);
904  m_fields[0]->FwdTransLocalElmt(sensor, sensFwd);
905 
906  variables.push_back("p");
907  variables.push_back("T");
908  variables.push_back("s");
909  variables.push_back("a");
910  variables.push_back("Mach");
911  variables.push_back("Sensor");
912  fieldcoeffs.push_back(pFwd);
913  fieldcoeffs.push_back(TFwd);
914  fieldcoeffs.push_back(sFwd);
915  fieldcoeffs.push_back(aFwd);
916  fieldcoeffs.push_back(mFwd);
917  fieldcoeffs.push_back(sensFwd);
918 
920  {
921  // reuse pressure
922  Array<OneD, NekDouble> sensorFwd(nCoeffs);
923  m_artificialDiffusion->GetArtificialViscosity(tmp, pressure);
924  m_fields[0]->FwdTransLocalElmt(pressure, sensorFwd);
925 
926  variables.push_back("ArtificialVisc");
927  fieldcoeffs.push_back(sensorFwd);
928  }
929  }
930 }
931 
932 /**
933  *
934  */
936  const Array<OneD, const Array<OneD, NekDouble>> &physfield,
938 {
939  m_varConv->GetPressure(physfield, pressure);
940 }
941 
942 /**
943  *
944  */
946  const Array<OneD, const Array<OneD, NekDouble>> &physfield,
947  Array<OneD, NekDouble> &density)
948 {
949  density = physfield[0];
950 }
951 
952 /**
953  *
954  */
956  const Array<OneD, const Array<OneD, NekDouble>> &physfield,
957  Array<OneD, Array<OneD, NekDouble>> &velocity)
958 {
959  m_varConv->GetVelocityVector(physfield, velocity);
960 }
961 
964 {
965  boost::ignore_unused(step);
966  const int nPoints = GetTotPoints();
967  const int nFields = m_fields.size();
969  Array<OneD, Array<OneD, NekDouble>> inarray(nFields);
970  for (int i = 0; i < nFields; ++i)
971  {
972  rhs[i] = Array<OneD, NekDouble>(nPoints, 0.0);
973  inarray[i] = m_fields[i]->UpdatePhys();
974  }
975 
976  DoOdeRhs(inarray, rhs, m_time);
977 
978  // Holds L2 errors.
980  Array<OneD, NekDouble> RHSL2(nFields);
981  Array<OneD, NekDouble> residual(nFields);
982 
983  for (int i = 0; i < nFields; ++i)
984  {
985  tmp = rhs[i];
986 
987  Vmath::Vmul(nPoints, tmp, 1, tmp, 1, tmp, 1);
988  residual[i] = Vmath::Vsum(nPoints, tmp, 1);
989  }
990 
991  m_comm->AllReduce(residual, LibUtilities::ReduceSum);
992 
993  NekDouble onPoints = 1.0 / NekDouble(nPoints);
994  for (int i = 0; i < nFields; ++i)
995  {
996  L2[i] = sqrt(residual[i] * onPoints);
997  }
998 }
999 } // namespace Nektar
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:215
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mode...
Definition: ErrorUtil.hpp:209
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
Definition: ErrorUtil.hpp:249
virtual void v_InitObject(bool DeclareFields=true)
Initialization object for CompressibleFlowSystem class.
void DoAdvection(const Array< OneD, Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray, const NekDouble time, const Array< OneD, Array< OneD, NekDouble >> &pFwd, const Array< OneD, Array< OneD, NekDouble >> &pBwd)
Compute the advection terms for the right-hand side.
virtual void GetVelocity(const Array< OneD, const Array< OneD, NekDouble >> &physfield, Array< OneD, Array< OneD, NekDouble >> &velocity)
Extract array with velocity from physfield.
void DoOdeRhs(const Array< OneD, const Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray, const NekDouble time)
Compute the right-hand side.
void DoDiffusion(const Array< OneD, Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray, const Array< OneD, Array< OneD, NekDouble >> &pFwd, const Array< OneD, Array< OneD, NekDouble >> &pBwd)
Add the diffusions terms to the right-hand side.
void GetFluxVector(const Array< OneD, const Array< OneD, NekDouble >> &physfield, TensorOfArray3D< NekDouble > &flux)
Return the flux vector for the compressible Euler equations.
Array< OneD, NekDouble > GetStabilityLimitVector(const Array< OneD, int > &ExpOrder)
Function to calculate the stability limit for DG/CG (a vector of them).
void InitAdvection()
Create advection and diffusion objects for CFS.
virtual ~CompressibleFlowSystem()
Destructor for CompressibleFlowSystem class.
virtual void v_DoDiffusion(const Array< OneD, Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray, const Array< OneD, Array< OneD, NekDouble >> &pFwd, const Array< OneD, Array< OneD, NekDouble >> &pBwd)
Apply artificial diffusion (Laplacian operator)
virtual NekDouble v_GetTimeStep(const Array< OneD, const Array< OneD, NekDouble >> &inarray)
Calculate the maximum timestep subject to CFL restrictions.
virtual Array< OneD, NekDouble > v_GetMaxStdVelocity(const NekDouble SpeedSoundFactor)
Compute the advection velocity in the standard space for each element of the expansion.
void GetFluxVectorDeAlias(const Array< OneD, const Array< OneD, NekDouble >> &physfield, TensorOfArray3D< NekDouble > &flux)
Return the flux vector for the compressible Euler equations by using the de-aliasing technique.
void SetBoundaryConditions(Array< OneD, Array< OneD, NekDouble >> &physarray, NekDouble time)
virtual void v_SetInitialConditions(NekDouble initialtime=0.0, bool dumpInitialConditions=true, const int domain=0)
Set up logic for residual calculation.
NekDouble GetStabilityLimit(int n)
Function to calculate the stability limit for DG/CG.
std::vector< CFSBndCondSharedPtr > m_bndConds
void GetElmtTimeStep(const Array< OneD, const Array< OneD, NekDouble >> &inarray, Array< OneD, NekDouble > &tstep)
Calculate the maximum timestep on each element subject to CFL restrictions.
Array< OneD, Array< OneD, NekDouble > > m_vecLocs
virtual void v_SteadyStateResidual(int step, Array< OneD, NekDouble > &L2)
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...
virtual void v_ExtraFldOutput(std::vector< Array< OneD, NekDouble >> &fieldcoeffs, std::vector< std::string > &variables)
VariableConverterSharedPtr m_varConv
void SetBoundaryConditionsBwdWeight()
Set up a weight on physical boundaries for boundary condition applications.
std::vector< SolverUtils::ForcingSharedPtr > m_forcing
ArtificialDiffusionSharedPtr m_artificialDiffusion
virtual void GetDensity(const Array< OneD, const Array< OneD, NekDouble >> &physfield, Array< OneD, NekDouble > &density)
Extract array with density from physfield.
const Array< OneD, const Array< OneD, NekDouble > > & GetNormals()
void InitialiseParameters()
Load CFS parameters from the session file.
const Array< OneD, const Array< OneD, NekDouble > > & GetVecLocs()
tBaseSharedPtr CreateInstance(tKey idKey, tParam... args)
Create an instance of the class referred to by idKey.
Definition: NekFactory.hpp:144
void DefineProjection(FuncPointerT func, ObjectPointerT obj)
void DefineOdeRhs(FuncPointerT func, ObjectPointerT obj)
void AccumulateRegion(std::string, int iolevel=0)
Accumulate elapsed time for a region.
Definition: Timer.cpp:73
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
A base class for PDEs which include an advection component.
SolverUtils::AdvectionSharedPtr m_advObject
Advection term.
virtual SOLVER_UTILS_EXPORT void v_InitObject(bool DeclareField=true)
Init object for UnsteadySystem class.
SOLVER_UTILS_EXPORT Array< OneD, NekDouble > GetElmtCFLVals(const bool FlagAcousticCFL=true)
int m_spacedim
Spatial dimension (>= expansion dim).
virtual SOLVER_UTILS_EXPORT void v_SetInitialConditions(NekDouble initialtime=0.0, bool dumpInitialConditions=true, const int domain=0)
SpatialDomains::MeshGraphSharedPtr m_graph
Pointer to graph defining mesh.
LibUtilities::CommSharedPtr m_comm
Communicator.
NekDouble m_timestep
Time step size.
NekDouble m_time
Current time of simulation.
SOLVER_UTILS_EXPORT int GetTraceTotPoints()
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
SOLVER_UTILS_EXPORT void Checkpoint_Output(const int n)
Write checkpoint file of m_fields.
SOLVER_UTILS_EXPORT int GetExpSize()
bool m_specHP_dealiasing
Flag to determine if dealisising is usde for the Spectral/hp element discretisation.
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Array holding trace normals for DG simulations in the forwards direction.
enum HomogeneousType m_HomogeneousType
SOLVER_UTILS_EXPORT int GetNpoints()
int m_nchk
Number of checkpoints written so far.
enum MultiRegions::ProjectionType m_projectionType
Type of projection; e.g continuous or discontinuous.
SOLVER_UTILS_EXPORT MultiRegions::ExpListSharedPtr GetPressure()
Get pressure field if available.
int m_checksteps
Number of steps between checkpoints.
SOLVER_UTILS_EXPORT int GetTotPoints()
static SOLVER_UTILS_EXPORT std::vector< ForcingSharedPtr > Load(const LibUtilities::SessionReaderSharedPtr &pSession, const std::weak_ptr< EquationSystem > &pEquation, const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const unsigned int &pNumForcingFields=0)
Definition: Forcing.cpp:120
Base class for unsteady solvers.
LibUtilities::TimeIntegrationSchemeOperators m_ode
The time integration scheme operators to use.
NekDouble m_cflSafetyFactor
CFL safety factor (comprise between 0 to 1).
SOLVER_UTILS_EXPORT NekDouble MaxTimeStepEstimator()
Get the maximum timestep estimator for cfl control.
bool m_homoInitialFwd
Flag to determine if simulation should start in homogeneous forward transformed state.
std::shared_ptr< SessionReader > SessionReaderSharedPtr
std::vector< PointsKey > PointsKeyVector
Definition: Points.h:250
std::shared_ptr< RiemannSolver > RiemannSolverSharedPtr
A shared pointer to an EquationSystem object.
AdvectionFactory & GetAdvectionFactory()
Gets the factory for initialising advection objects.
Definition: Advection.cpp:47
RiemannSolverFactory & GetRiemannSolverFactory()
std::shared_ptr< GeomFactors > GeomFactorsSharedPtr
Pointer to a GeomFactors object.
Definition: GeomFactors.h:62
@ eDeformed
Geometry is curved or has non-constant factors.
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:172
The above copyright notice and this permission notice shall be included.
Definition: CoupledSolver.h:1
CFSBndCondFactory & GetCFSBndCondFactory()
Declaration of the boundary condition factory singleton.
Definition: CFSBndCond.cpp:41
static Array< OneD, Array< OneD, NekDouble > > NullNekDoubleArrayOfArray
ArtificialDiffusionFactory & GetArtificialDiffusionFactory()
Declaration of the artificial diffusion factory singleton.
double NekDouble
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:209
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:622
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:518
T Vmin(int n, const T *x, const int incx)
Return the minimum element in x - called vmin to avoid conflict with min.
Definition: Vmath.cpp:1050
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:574
T Vsum(int n, const T *x, const int incx)
Subtract return sum(x)
Definition: Vmath.cpp:895
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:359
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
Definition: Vmath.cpp:248
void FillWhiteNoise(int n, const T eps, T *x, const int incx, int outseed)
Fills a vector with white noise.
Definition: Vmath.cpp:155
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:945
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1255
scalarT< T > abs(scalarT< T > in)
Definition: scalar.hpp:295
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:291