Nektar++
VelocityCorrectionScheme.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File VelocityCorrection.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: Velocity Correction Scheme for the Incompressible
32 // Navier Stokes equations
33 ///////////////////////////////////////////////////////////////////////////////
34 
39 #include <MultiRegions/ContField.h>
41 #include <SolverUtils/Core/Misc.h>
42 
43 #include <boost/algorithm/string.hpp>
44 
45 using namespace std;
46 
47 namespace Nektar
48 {
49 using namespace MultiRegions;
50 
51 string VelocityCorrectionScheme::className =
53  "VelocityCorrectionScheme", VelocityCorrectionScheme::create);
54 
55 /**
56  * Constructor. Creates ...
57  *
58  * \param
59  * \param
60  */
61 VelocityCorrectionScheme::VelocityCorrectionScheme(
64  : UnsteadySystem(pSession, pGraph), IncNavierStokes(pSession, pGraph),
65  m_varCoeffLap(StdRegions::NullVarCoeffMap)
66 {
67 }
68 
70 {
71  int n;
72 
73  IncNavierStokes::v_InitObject(DeclareField);
74  m_explicitDiffusion = false;
75 
76  // Set m_pressure to point to last field of m_fields;
77  if (boost::iequals(m_session->GetVariable(m_fields.size() - 1), "p"))
78  {
79  m_nConvectiveFields = m_fields.size() - 1;
81  }
82  else
83  {
84  ASSERTL0(false, "Need to set up pressure field definition");
85  }
86 
87  // Determine diffusion coefficients for each field
89  for (n = 0; n < m_nConvectiveFields; ++n)
90  {
91  std::string varName = m_session->GetVariable(n);
92  if (m_session->DefinesFunction("DiffusionCoefficient", varName))
93  {
95  m_session->GetFunction("DiffusionCoefficient", varName);
96  m_diffCoeff[n] = ffunc->Evaluate();
97  }
98  }
99 
100  // Integrate only the convective fields
101  for (n = 0; n < m_nConvectiveFields; ++n)
102  {
103  m_intVariables.push_back(n);
104  }
105 
107  SetUpSVV();
108 
109  // check to see if it is explicity turned off
110  m_session->MatchSolverInfo("GJPStabilisation", "False",
111  m_useGJPStabilisation, true);
112 
113  // if GJPStabilisation set to False bool will be true and
114  // if not false so negate/revese bool
116 
117  m_session->MatchSolverInfo("GJPNormalVelocity", "True", m_useGJPNormalVel,
118  false);
119 
120  if (m_useGJPNormalVel)
121  {
122  ASSERTL0(boost::iequals(m_session->GetSolverInfo("GJPStabilisation"),
123  "Explicit"),
124  "Can only specify GJPNormalVelocity with"
125  " GJPStabilisation set to Explicit currently");
126  }
127 
128  m_session->LoadParameter("GJPJumpScale", m_GJPJumpScale, 1.0);
129 
130  m_session->MatchSolverInfo("SmoothAdvection", "True", m_SmoothAdvection,
131  false);
132 
133  // set explicit time-intregration class operators
136 
137  // set implicit time-intregration class operators
140 
141  // Set up bits for flowrate.
142  m_session->LoadParameter("Flowrate", m_flowrate, 0.0);
143  m_session->LoadParameter("IO_FlowSteps", m_flowrateSteps, 0);
144 }
145 
147 {
148  // creation of the extrapolation object
151  {
152  std::string vExtrapolation = v_GetExtrapolateStr();
153  if (m_session->DefinesSolverInfo("Extrapolation"))
154  {
155  vExtrapolation = v_GetSubSteppingExtrapolateStr(
156  m_session->GetSolverInfo("Extrapolation"));
157  }
159  vExtrapolation, m_session, m_fields, m_pressure, m_velocity,
160  m_advObject);
161 
162  m_extrapolation->SetForcing(m_forcing);
163  m_extrapolation->SubSteppingTimeIntegration(m_intScheme);
164  m_extrapolation->GenerateHOPBCMap(m_session);
165  }
166 }
167 
168 /**
169  * @brief Set up the Stokes solution used to impose constant flowrate
170  * through a boundary.
171  *
172  * This routine solves a Stokes equation using a unit forcing direction,
173  * specified by the user to be in the desired flow direction. This field can
174  * then be used to correct the end of each timestep to impose a constant
175  * volumetric flow rate through a user-defined boundary.
176  *
177  * There are three modes of operation:
178  *
179  * - Standard two-dimensional or three-dimensional simulations (e.g. pipes
180  * or channels)
181  * - 3DH1D simulations where the forcing is not in the homogeneous
182  * direction (e.g. channel flow, where the y-direction of the 2D mesh
183  * is perpendicular to the wall);
184  * - 3DH1D simulations where the forcing is in the homogeneous direction
185  * (e.g. pipe flow in the z-direction).
186  *
187  * In the first two cases, the user should define:
188  * - the `Flowrate` parameter, which dictates the volumetric flux through
189  * the reference area
190  * - tag a boundary region with the `Flowrate` user-defined type to define
191  * the reference area
192  * - define a `FlowrateForce` function with components `ForceX`, `ForceY`
193  * and `ForceZ` that defines a unit forcing in the appropriate direction.
194  *
195  * In the latter case, the user should define only the `Flowrate`; the
196  * reference area is taken to be the homogeneous plane and the force is
197  * assumed to be the unit z-vector \f$ \hat{e}_z \f$.
198  *
199  * This routine solves a single timestep of the Stokes problem
200  * (premultiplied by the backwards difference coefficient):
201  *
202  * \f[ \frac{\partial\mathbf{u}}{\partial t} = -\nabla p +
203  * \nu\nabla^2\mathbf{u} + \mathbf{f} \f]
204  *
205  * with a zero initial condition to obtain a field \f$ \mathbf{u}_s \f$. The
206  * flowrate is then corrected at each timestep \f$ n \f$ by adding the
207  * correction \f$ \alpha\mathbf{u}_s \f$ where
208  *
209  * \f[ \alpha = \frac{\overline{Q} - Q(\mathbf{u^n})}{Q(\mathbf{u}_s)} \f]
210  *
211  * where \f$ Q(\cdot)\f$ is the volumetric flux through the appropriate
212  * surface or line, which is implemented in
213  * VelocityCorrectionScheme::MeasureFlowrate. For more details, see chapter
214  * 3.2 of the thesis of D. Moxey (University of Warwick, 2011).
215  */
217 {
218  m_flowrateBndID = -1;
219  m_flowrateArea = 0.0;
220 
222  m_fields[0]->GetBndConditions();
223 
224  std::string forces[] = {"X", "Y", "Z"};
225  Array<OneD, NekDouble> flowrateForce(m_spacedim, 0.0);
226 
227  // Set up flowrate forces.
228  bool defined = true;
229  for (int i = 0; i < m_spacedim; ++i)
230  {
231  std::string varName = std::string("Force") + forces[i];
232  defined = m_session->DefinesFunction("FlowrateForce", varName);
233 
234  if (!defined && m_HomogeneousType == eHomogeneous1D)
235  {
236  break;
237  }
238 
239  ASSERTL0(defined,
240  "A 'FlowrateForce' function must defined with components "
241  "[ForceX, ...] to define direction of flowrate forcing");
242 
244  m_session->GetFunction("FlowrateForce", varName);
245  flowrateForce[i] = ffunc->Evaluate();
246  }
247 
248  // Define flag for case with homogeneous expansion and forcing not in the
249  // z-direction
250  m_homd1DFlowinPlane = false;
251  if (defined && m_HomogeneousType == eHomogeneous1D)
252  {
253  m_homd1DFlowinPlane = true;
254  }
255 
256  // For 3DH1D simulations, if force isn't defined then assume in
257  // z-direction.
258  if (!defined)
259  {
260  flowrateForce[2] = 1.0;
261  }
262 
263  // Find the boundary condition that is tagged as the flowrate boundary.
264  for (int i = 0; i < bcs.size(); ++i)
265  {
266  if (boost::iequals(bcs[i]->GetUserDefined(), "Flowrate"))
267  {
268  m_flowrateBndID = i;
269  break;
270  }
271  }
272 
273  int tmpBr = m_flowrateBndID;
274  m_comm->AllReduce(tmpBr, LibUtilities::ReduceMax);
275  ASSERTL0(tmpBr >= 0 || m_HomogeneousType == eHomogeneous1D,
276  "One boundary region must be marked using the 'Flowrate' "
277  "user-defined type to monitor the volumetric flowrate.");
278 
279  // Extract an appropriate expansion list to represents the boundary.
280  if (m_flowrateBndID >= 0)
281  {
282  // For a boundary, extract the boundary itself.
283  m_flowrateBnd = m_fields[0]->GetBndCondExpansions()[m_flowrateBndID];
284  }
286  {
287  // For 3DH1D simulations with no force specified, find the mean
288  // (0th) plane.
289  Array<OneD, unsigned int> zIDs = m_fields[0]->GetZIDs();
290  int tmpId = -1;
291 
292  for (int i = 0; i < zIDs.size(); ++i)
293  {
294  if (zIDs[i] == 0)
295  {
296  tmpId = i;
297  break;
298  }
299  }
300 
301  ASSERTL1(tmpId <= 0, "Should be either at location 0 or -1 if not "
302  "found");
303 
304  if (tmpId != -1)
305  {
306  m_flowrateBnd = m_fields[0]->GetPlane(tmpId);
307  }
308  }
309 
310  // At this point, some processors may not have m_flowrateBnd
311  // set if they don't contain the appropriate boundary. To
312  // calculate the area, we integrate 1.0 over the boundary
313  // (which has been set up with the appropriate subcommunicator
314  // to avoid deadlock), and then communicate this to the other
315  // processors with an AllReduce.
316  if (m_flowrateBnd)
317  {
318  Array<OneD, NekDouble> inArea(m_flowrateBnd->GetNpoints(), 1.0);
319  m_flowrateArea = m_flowrateBnd->Integral(inArea);
320  }
322 
323  // In homogeneous case with forcing not aligned to the z-direction,
324  // redefine m_flowrateBnd so it is a 1D expansion
327  {
328  // For 3DH1D simulations with no force specified, find the mean
329  // (0th) plane.
330  Array<OneD, unsigned int> zIDs = m_fields[0]->GetZIDs();
331  m_planeID = -1;
332 
333  for (int i = 0; i < zIDs.size(); ++i)
334  {
335  if (zIDs[i] == 0)
336  {
337  m_planeID = i;
338  break;
339  }
340  }
341 
342  ASSERTL1(m_planeID <= 0, "Should be either at location 0 or -1 "
343  "if not found");
344 
345  if (m_planeID != -1)
346  {
347  m_flowrateBnd =
348  m_fields[0]->GetBndCondExpansions()[m_flowrateBndID]->GetPlane(
349  m_planeID);
350  }
351  }
352 
353  // Set up some storage for the Stokes solution (to be stored in
354  // m_flowrateStokes) and its initial condition (inTmp), which holds the
355  // unit forcing.
356  int nqTot = m_fields[0]->GetNpoints();
359 
360  for (int i = 0; i < m_spacedim; ++i)
361  {
362  inTmp[i] = Array<OneD, NekDouble>(nqTot, flowrateForce[i] * aii_dt);
363  m_flowrateStokes[i] = Array<OneD, NekDouble>(nqTot, 0.0);
364 
366  {
367  Array<OneD, NekDouble> inTmp2(nqTot);
368  m_fields[i]->HomogeneousFwdTrans(inTmp[i], inTmp2);
369  m_fields[i]->SetWaveSpace(true);
370  inTmp[i] = inTmp2;
371  }
372 
373  Vmath::Zero(m_fields[i]->GetNcoeffs(), m_fields[i]->UpdateCoeffs(), 1);
374  }
375 
376  // Create temporary extrapolation object to avoid issues with
377  // m_extrapolation for HOPBCs using higher order timestepping schemes.
378  // Zero pressure BCs in Neumann boundaries that may have been
379  // set in the advection step.
381  m_pressure->GetBndConditions();
383  m_pressure->GetBndCondExpansions();
384  for (int n = 0; n < PBndConds.size(); ++n)
385  {
386  if (PBndConds[n]->GetBoundaryConditionType() ==
388  {
389  Vmath::Zero(PBndExp[n]->GetNcoeffs(), PBndExp[n]->UpdateCoeffs(),
390  1);
391  }
392  }
393 
394  // Finally, calculate the solution and the flux of the Stokes
395  // solution. We set m_greenFlux to maximum numeric limit, which signals
396  // to SolveUnsteadyStokesSystem that we don't need to apply a flowrate
397  // force.
398  m_greenFlux = numeric_limits<NekDouble>::max();
399  m_flowrateAiidt = aii_dt;
400 
401  // Save the number of convective field in case it is not set
402  // to spacedim. Only need velocity components for stokes forcing
403  int SaveNConvectiveFields = m_nConvectiveFields;
405  SolveUnsteadyStokesSystem(inTmp, m_flowrateStokes, 0.0, aii_dt);
406  m_nConvectiveFields = SaveNConvectiveFields;
408 
409  // If the user specified IO_FlowSteps, open a handle to store output.
410  if (m_comm->GetRank() == 0 && m_flowrateSteps &&
411  !m_flowrateStream.is_open())
412  {
413  std::string filename = m_session->GetSessionName();
414  filename += ".prs";
415  m_flowrateStream.open(filename.c_str());
416  m_flowrateStream.setf(ios::scientific, ios::floatfield);
417  m_flowrateStream << "# step time dP" << endl
418  << "# -------------------------------------------"
419  << endl;
420  }
421 
422  // Replace pressure BCs with those evaluated from advection step
423  m_extrapolation->CopyPressureHBCsToPbndExp();
424 }
425 
426 /**
427  * @brief Measure the volumetric flow rate through the volumetric flow rate
428  * reference surface.
429  *
430  * This routine computes the volumetric flow rate
431  *
432  * \f[
433  * Q(\mathbf{u}) = \frac{1}{\mu(R)} \int_R \mathbf{u} \cdot d\mathbf{s}
434  * \f]
435  *
436  * through the boundary region \f$ R \f$.
437  */
439  const Array<OneD, Array<OneD, NekDouble>> &inarray)
440 {
441  NekDouble flowrate = 0.0;
442 
443  if (m_flowrateBnd && m_flowrateBndID >= 0)
444  {
445  // If we're an actual boundary, calculate the vector flux through
446  // the boundary.
448 
449  if (!m_homd1DFlowinPlane)
450  {
451  // General case
452  for (int i = 0; i < m_spacedim; ++i)
453  {
454  m_fields[i]->ExtractPhysToBnd(m_flowrateBndID, inarray[i],
455  boundary[i]);
456  }
457  flowrate = m_flowrateBnd->VectorFlux(boundary);
458  }
459  else if (m_planeID == 0)
460  {
461  // Homogeneous with forcing in plane. Calculate flux only on
462  // the meanmode - calculateFlux necessary for hybrid
463  // parallelisation.
464  for (int i = 0; i < m_spacedim; ++i)
465  {
466  m_fields[i]->GetPlane(m_planeID)->ExtractPhysToBnd(
467  m_flowrateBndID, inarray[i], boundary[i]);
468  }
469 
470  // the flowrate is calculated on the mean mode so it needs to be
471  // multiplied by LZ to be consistent with the general case.
472  flowrate = m_flowrateBnd->VectorFlux(boundary) *
473  m_session->GetParameter("LZ");
474  }
475  }
476  else if (m_flowrateBnd && !m_homd1DFlowinPlane)
477  {
478  // 3DH1D case with no Flowrate boundary defined: compute flux
479  // through the zero-th (mean) plane.
480  flowrate = m_flowrateBnd->Integral(inarray[2]);
481  }
482 
483  // Communication to obtain the total flowrate
485  {
486  m_comm->GetColumnComm()->AllReduce(flowrate, LibUtilities::ReduceSum);
487  }
488  else
489  {
490  m_comm->AllReduce(flowrate, LibUtilities::ReduceSum);
491  }
492  return flowrate / m_flowrateArea;
493 }
494 
496 {
497  if (m_flowrateSteps > 0)
498  {
499  if (m_comm->GetRank() == 0 && (step + 1) % m_flowrateSteps == 0)
500  {
501  m_flowrateStream << setw(8) << step << setw(16) << m_time
502  << setw(16) << m_alpha << endl;
503  }
504  }
505 
507 }
508 
509 /**
510  * Destructor
511  */
513 {
514 }
515 
516 /**
517  *
518  */
520 {
522  SolverUtils::AddSummaryItem(s, "Splitting Scheme",
523  "Velocity correction (strong press. form)");
524 
525  if (m_extrapolation->GetSubStepName().size())
526  {
527  SolverUtils::AddSummaryItem(s, "Substepping",
528  m_extrapolation->GetSubStepName());
529  }
530 
531  string dealias = m_homogen_dealiasing ? "Homogeneous1D" : "";
533  {
534  dealias += (dealias == "" ? "" : " + ") + string("spectral/hp");
535  }
536  if (dealias != "")
537  {
538  SolverUtils::AddSummaryItem(s, "Dealiasing", dealias);
539  }
540 
541  string smoothing = m_useSpecVanVisc ? "spectral/hp" : "";
542  if (smoothing != "")
543  {
545  {
547  s, "Smoothing-SpecHP",
548  "SVV (" + smoothing + " Exp Kernel(cut-off = " +
549  boost::lexical_cast<string>(m_sVVCutoffRatio) +
550  ", diff coeff = " +
551  boost::lexical_cast<string>(m_sVVDiffCoeff) + "))");
552  }
553  else
554  {
555  if (m_IsSVVPowerKernel)
556  {
558  s, "Smoothing-SpecHP",
559  "SVV (" + smoothing + " Power Kernel (Power ratio =" +
560  boost::lexical_cast<string>(m_sVVCutoffRatio) +
561  ", diff coeff = " +
562  boost::lexical_cast<string>(m_sVVDiffCoeff) +
563  "*Uh/p))");
564  }
565  else
566  {
568  s, "Smoothing-SpecHP",
569  "SVV (" + smoothing + " DG Kernel (diff coeff = " +
570  boost::lexical_cast<string>(m_sVVDiffCoeff) +
571  "*Uh/p))");
572  }
573  }
574  }
575 
577  {
579  s, "Smoothing-Homo1D",
580  "SVV (Homogeneous1D - Exp Kernel(cut-off = " +
581  boost::lexical_cast<string>(m_sVVCutoffRatioHomo1D) +
582  ", diff coeff = " +
583  boost::lexical_cast<string>(m_sVVDiffCoeffHomo1D) + "))");
584  }
585 
587  {
589  s, "GJP Stab. Impl. ",
590  m_session->GetSolverInfo("GJPStabilisation"));
591  SolverUtils::AddSummaryItem(s, "GJP Stab. JumpScale", m_GJPJumpScale);
592 
593  if (boost::iequals(m_session->GetSolverInfo("GJPStabilisation"),
594  "Explicit"))
595  {
597  s, "GJP Normal Velocity",
598  m_session->GetSolverInfo("GJPNormalVelocity"));
599  }
600  }
601 }
602 
603 /**
604  *
605  */
607 {
609 
610  for (int i = 0; i < m_nConvectiveFields; ++i)
611  {
613  }
614 
615  m_flowrateAiidt = 0.0;
616 
618 
619  // Set up Field Meta Data for output files
620  m_fieldMetaDataMap["Kinvis"] = boost::lexical_cast<std::string>(m_kinvis);
621  m_fieldMetaDataMap["TimeStep"] =
622  boost::lexical_cast<std::string>(m_timestep);
623 
624  // set boundary conditions here so that any normal component
625  // correction are imposed before they are imposed on initial
626  // field below
628 
629  // Ensure the initial conditions have correct BCs
630  for (int i = 0; i < m_fields.size(); ++i)
631  {
632  m_fields[i]->ImposeDirichletConditions(m_fields[i]->UpdateCoeffs());
633  m_fields[i]->LocalToGlobal();
634  m_fields[i]->GlobalToLocal();
635  m_fields[i]->BwdTrans(m_fields[i]->GetCoeffs(),
636  m_fields[i]->UpdatePhys());
637  }
638 }
639 
640 /**
641  *
642  */
644 {
645  int nfields = m_fields.size() - 1;
646  for (int k = 0; k < nfields; ++k)
647  {
648  // Backward Transformation in physical space for time evolution
649  m_fields[k]->BwdTrans(m_fields[k]->GetCoeffs(),
650  m_fields[k]->UpdatePhys());
651  }
652 }
653 
654 /**
655  *
656  */
658 {
659 
660  int nfields = m_fields.size() - 1;
661  for (int k = 0; k < nfields; ++k)
662  {
663  // Forward Transformation in physical space for time evolution
664  m_fields[k]->FwdTransLocalElmt(m_fields[k]->GetPhys(),
665  m_fields[k]->UpdateCoeffs());
666  }
667 }
668 
669 /**
670  *
671  */
673 {
674  int vVar = m_session->GetVariables().size();
675  Array<OneD, bool> vChecks(vVar, false);
676  vChecks[vVar - 1] = true;
677  return vChecks;
678 }
679 
680 /**
681  *
682  */
684 {
685  return m_session->GetVariables().size() - 1;
686 }
687 
688 /**
689  * Explicit part of the method - Advection, Forcing + HOPBCs
690  */
692  const Array<OneD, const Array<OneD, NekDouble>> &inarray,
693  Array<OneD, Array<OneD, NekDouble>> &outarray, const NekDouble time)
694 {
695  LibUtilities::Timer timer;
696  timer.Start();
697  EvaluateAdvectionTerms(inarray, outarray, time);
698  timer.Stop();
699  timer.AccumulateRegion("Advection Terms");
700 
701  // Smooth advection
702  if (m_SmoothAdvection)
703  {
704  for (int i = 0; i < m_nConvectiveFields; ++i)
705  {
706  m_pressure->SmoothField(outarray[i]);
707  }
708  }
709 
710  // Add forcing terms
711  for (auto &x : m_forcing)
712  {
713  x->Apply(m_fields, inarray, outarray, time);
714  }
715 
716  // Calculate High-Order pressure boundary conditions
717  timer.Start();
718  m_extrapolation->EvaluatePressureBCs(inarray, outarray, m_kinvis);
719  timer.Stop();
720  timer.AccumulateRegion("Pressure BCs");
721 }
722 
723 /**
724  * Implicit part of the method - Poisson + nConv*Helmholtz
725  */
727  const Array<OneD, const Array<OneD, NekDouble>> &inarray,
728  Array<OneD, Array<OneD, NekDouble>> &outarray, const NekDouble time,
729  const NekDouble aii_Dt)
730 {
731  // Set up flowrate if we're starting for the first time or the value of
732  // aii_Dt has changed.
733  if (m_flowrate > 0.0 && (aii_Dt != m_flowrateAiidt))
734  {
735  SetupFlowrate(aii_Dt);
736  }
737 
738  int physTot = m_fields[0]->GetTotPoints();
739 
740  // Substep the pressure boundary condition if using substepping
741  m_extrapolation->SubStepSetPressureBCs(inarray, aii_Dt, m_kinvis);
742 
743  // Set up forcing term for pressure Poisson equation
744  LibUtilities::Timer timer;
745  timer.Start();
746  SetUpPressureForcing(inarray, m_F, aii_Dt);
747  timer.Stop();
748  timer.AccumulateRegion("Pressure Forcing");
749 
750  // Solve Pressure System
751  timer.Start();
752  SolvePressure(m_F[0]);
753  timer.Stop();
754  timer.AccumulateRegion("Pressure Solve");
755 
756  // Set up forcing term for Helmholtz problems
757  timer.Start();
758  SetUpViscousForcing(inarray, m_F, aii_Dt);
759  timer.Stop();
760  timer.AccumulateRegion("Viscous Forcing");
761 
762  // Solve velocity system
763  timer.Start();
764  SolveViscous(m_F, inarray, outarray, aii_Dt);
765  timer.Stop();
766  timer.AccumulateRegion("Viscous Solve");
767 
768  // Apply flowrate correction
769  if (m_flowrate > 0.0 && m_greenFlux != numeric_limits<NekDouble>::max())
770  {
771  NekDouble currentFlux = MeasureFlowrate(outarray);
772  m_alpha = (m_flowrate - currentFlux) / m_greenFlux;
773 
774  for (int i = 0; i < m_spacedim; ++i)
775  {
776  Vmath::Svtvp(physTot, m_alpha, m_flowrateStokes[i], 1, outarray[i],
777  1, outarray[i], 1);
778  // Enusre coeff space is updated for next time step
779  m_fields[i]->FwdTransLocalElmt(outarray[i],
780  m_fields[i]->UpdateCoeffs());
781  // Impsoe symmetry of flow on coeff space (good to enfore
782  // periodicity).
783  m_fields[i]->LocalToGlobal();
784  m_fields[i]->GlobalToLocal();
785  }
786  }
787 }
788 
789 /**
790  * Forcing term for Poisson solver solver
791  */
793  const Array<OneD, const Array<OneD, NekDouble>> &fields,
795 {
796  int i;
797  int physTot = m_fields[0]->GetTotPoints();
798  int nvel = m_velocity.size();
799 
800  m_fields[0]->PhysDeriv(eX, fields[0], Forcing[0]);
801 
802  for (i = 1; i < nvel; ++i)
803  {
804  // Use Forcing[1] as storage since it is not needed for the pressure
805  m_fields[i]->PhysDeriv(DirCartesianMap[i], fields[i], Forcing[1]);
806  Vmath::Vadd(physTot, Forcing[1], 1, Forcing[0], 1, Forcing[0], 1);
807  }
808 
809  Vmath::Smul(physTot, 1.0 / aii_Dt, Forcing[0], 1, Forcing[0], 1);
810 }
811 
812 /**
813  * Forcing term for Helmholtz solver
814  */
816  const Array<OneD, const Array<OneD, NekDouble>> &inarray,
818 {
819  NekDouble aii_dtinv = 1.0 / aii_Dt;
820  int phystot = m_fields[0]->GetTotPoints();
821 
822  // Grad p
823  m_pressure->BwdTrans(m_pressure->GetCoeffs(), m_pressure->UpdatePhys());
824 
825  int nvel = m_velocity.size();
826  if (nvel == 2)
827  {
828  m_pressure->PhysDeriv(m_pressure->GetPhys(), Forcing[m_velocity[0]],
829  Forcing[m_velocity[1]]);
830  }
831  else
832  {
833  m_pressure->PhysDeriv(m_pressure->GetPhys(), Forcing[m_velocity[0]],
835  }
836 
837  // zero convective fields.
838  for (int i = nvel; i < m_nConvectiveFields; ++i)
839  {
840  Vmath::Zero(phystot, Forcing[i], 1);
841  }
842 
843  // Subtract inarray/(aii_dt) and divide by kinvis. Kinvis will
844  // need to be updated for the convected fields.
845  for (int i = 0; i < m_nConvectiveFields; ++i)
846  {
847  Blas::Daxpy(phystot, -aii_dtinv, inarray[i], 1, Forcing[i], 1);
848  Blas::Dscal(phystot, 1.0 / m_diffCoeff[i], &(Forcing[i])[0], 1);
849  }
850 }
851 
852 /**
853  * Solve pressure system
854  */
857 {
859  // Setup coefficient for equation
860  factors[StdRegions::eFactorLambda] = 0.0;
861 
862  // Solver Pressure Poisson Equation
863  m_pressure->HelmSolve(Forcing, m_pressure->UpdateCoeffs(), factors);
864 
865  // Add presure to outflow bc if using convective like BCs
866  m_extrapolation->AddPressureToOutflowBCs(m_kinvis);
867 }
868 
869 /**
870  * Solve velocity system
871  */
873  const Array<OneD, const Array<OneD, NekDouble>> &Forcing,
874  const Array<OneD, const Array<OneD, NekDouble>> &inarray,
875  Array<OneD, Array<OneD, NekDouble>> &outarray, const NekDouble aii_Dt)
876 {
880 
881  AppendSVVFactors(factors, varFactorsMap);
882 
883  // Calculate Normal velocity at Trace for GJP explicit stabiliation
884  if (m_useGJPNormalVel)
885  {
887  std::dynamic_pointer_cast<MultiRegions::ContField>(m_fields[0]);
888 
890  cfield->GetGJPForcing();
891 
892  int nTracePts = GJPData->GetNumTracePts();
893  Array<OneD, NekDouble> unorm(nTracePts, 1.0);
894  Array<OneD, NekDouble> Fwd(nTracePts), Bwd(nTracePts);
895  Array<OneD, Array<OneD, NekDouble>> traceNormals =
896  GJPData->GetTraceNormals();
897 
898  m_fields[0]->GetFwdBwdTracePhys(inarray[0], Fwd, Bwd, true, true);
899  Vmath::Vmul(nTracePts, Fwd, 1, traceNormals[0], 1, unorm, 1);
900 
901  // Evaluate u.n on trace
902  for (int f = 1; f < m_fields[0]->GetCoordim(0); ++f)
903  {
904  m_fields[0]->GetFwdBwdTracePhys(inarray[f], Fwd, Bwd, true, true);
905  Vmath::Vvtvp(nTracePts, Fwd, 1, traceNormals[f], 1, unorm, 1, unorm,
906  1);
907  }
908  Vmath::Vabs(nTracePts, unorm, 1, unorm, 1);
909  varCoeffMap[StdRegions::eVarCoeffGJPNormVel] = unorm;
910  }
911 
912  // Solve Helmholtz system and put in Physical space
913  for (int i = 0; i < m_nConvectiveFields; ++i)
914  {
915  // test by adding GJP implicit
917  {
919  }
920 
921  // Setup coefficients for equation
922  factors[StdRegions::eFactorLambda] = 1.0 / aii_Dt / m_diffCoeff[i];
923  m_fields[i]->HelmSolve(Forcing[i], m_fields[i]->UpdateCoeffs(), factors,
924  varCoeffMap, varFactorsMap);
925  m_fields[i]->BwdTrans(m_fields[i]->GetCoeffs(), outarray[i]);
926  }
927 }
928 
930 {
931 
932  m_session->MatchSolverInfo("SpectralVanishingViscosity", "PowerKernel",
933  m_useSpecVanVisc, false);
934 
935  if (m_useSpecVanVisc)
936  {
937  m_useHomo1DSpecVanVisc = true;
938  }
939  else
940  {
941  m_session->MatchSolverInfo("SpectralVanishingViscositySpectralHP",
942  "PowerKernel", m_useSpecVanVisc, false);
943  }
944 
945  if (m_useSpecVanVisc)
946  {
947  m_IsSVVPowerKernel = true;
948  }
949  else
950  {
951  m_session->MatchSolverInfo("SpectralVanishingViscosity", "DGKernel",
952  m_useSpecVanVisc, false);
953  if (m_useSpecVanVisc)
954  {
955  m_useHomo1DSpecVanVisc = true;
956  }
957  else
958  {
959  m_session->MatchSolverInfo("SpectralVanishingViscositySpectralHP",
960  "DGKernel", m_useSpecVanVisc, false);
961  }
962 
963  if (m_useSpecVanVisc)
964  {
965  m_IsSVVPowerKernel = false;
966  }
967  }
968 
969  // set up varcoeff kernel if PowerKernel or DG is specified
970  if (m_useSpecVanVisc)
971  {
972  Array<OneD, Array<OneD, NekDouble>> SVVVelFields =
974  if (m_session->DefinesFunction("SVVVelocityMagnitude"))
975  {
976  if (m_comm->GetRank() == 0)
977  {
978  cout << "Seting up SVV velocity from "
979  "SVVVelocityMagnitude section in session file"
980  << endl;
981  }
982  int nvel = m_velocity.size();
983  int phystot = m_fields[0]->GetTotPoints();
984  SVVVelFields = Array<OneD, Array<OneD, NekDouble>>(nvel);
985  vector<string> vars;
986  for (int i = 0; i < nvel; ++i)
987  {
988  SVVVelFields[i] = Array<OneD, NekDouble>(phystot);
989  vars.push_back(m_session->GetVariable(m_velocity[i]));
990  }
991 
992  // Load up files into m_fields;
993  GetFunction("SVVVelocityMagnitude")->Evaluate(vars, SVVVelFields);
994  }
995 
997  SVVVarDiffCoeff(1.0, m_svvVarDiffCoeff, SVVVelFields);
998  m_session->LoadParameter("SVVDiffCoeff", m_sVVDiffCoeff, 1.0);
999  }
1000  else
1001  {
1003  m_session->LoadParameter("SVVDiffCoeff", m_sVVDiffCoeff, 0.1);
1004  }
1005 
1006  // Load parameters for Spectral Vanishing Viscosity
1007  if (m_useSpecVanVisc == false)
1008  {
1009  m_session->MatchSolverInfo("SpectralVanishingViscosity", "True",
1010  m_useSpecVanVisc, false);
1011  if (m_useSpecVanVisc == false)
1012  {
1013  m_session->MatchSolverInfo("SpectralVanishingViscosity",
1014  "ExpKernel", m_useSpecVanVisc, false);
1015  }
1017 
1018  if (m_useSpecVanVisc == false)
1019  {
1020  m_session->MatchSolverInfo("SpectralVanishingViscositySpectralHP",
1021  "True", m_useSpecVanVisc, false);
1022  if (m_useSpecVanVisc == false)
1023  {
1024  m_session->MatchSolverInfo(
1025  "SpectralVanishingViscositySpectralHP", "ExpKernel",
1026  m_useSpecVanVisc, false);
1027  }
1028  }
1029  }
1030 
1031  // Case of only Homo1D kernel
1032  if (m_session->DefinesSolverInfo("SpectralVanishingViscosityHomo1D"))
1033  {
1034  m_session->MatchSolverInfo("SpectralVanishingViscosityHomo1D", "True",
1035  m_useHomo1DSpecVanVisc, false);
1036  if (m_useHomo1DSpecVanVisc == false)
1037  {
1038  m_session->MatchSolverInfo("SpectralVanishingViscosityHomo1D",
1039  "ExpKernel", m_useHomo1DSpecVanVisc,
1040  false);
1041  }
1042  }
1043 
1044  m_session->LoadParameter("SVVCutoffRatio", m_sVVCutoffRatio, 0.75);
1045  m_session->LoadParameter("SVVCutoffRatioHomo1D", m_sVVCutoffRatioHomo1D,
1047  m_session->LoadParameter("SVVDiffCoeffHomo1D", m_sVVDiffCoeffHomo1D,
1048  m_sVVDiffCoeff);
1049 
1051  {
1052  ASSERTL0(
1053  m_nConvectiveFields > 2,
1054  "Expect to have three velocity fields with homogenous expansion");
1055 
1057  {
1059  planes = m_fields[0]->GetZIDs();
1060 
1061  int num_planes = planes.size();
1062  Array<OneD, NekDouble> SVV(num_planes, 0.0);
1063  NekDouble fac;
1064  int kmodes = m_fields[0]->GetHomogeneousBasis()->GetNumModes();
1065  int pstart;
1066 
1067  pstart = m_sVVCutoffRatioHomo1D * kmodes;
1068 
1069  for (int n = 0; n < num_planes; ++n)
1070  {
1071  if (planes[n] > pstart)
1072  {
1073  fac = (NekDouble)((planes[n] - kmodes) *
1074  (planes[n] - kmodes)) /
1075  ((NekDouble)((planes[n] - pstart) *
1076  (planes[n] - pstart)));
1077  SVV[n] = m_sVVDiffCoeffHomo1D * exp(-fac) / m_kinvis;
1078  }
1079  }
1080 
1081  for (int i = 0; i < m_velocity.size(); ++i)
1082  {
1083  m_fields[m_velocity[i]]->SetHomo1DSpecVanVisc(SVV);
1084  }
1085  }
1086  }
1087 }
1088 
1090  const NekDouble velmag, Array<OneD, NekDouble> &diffcoeff,
1091  const Array<OneD, Array<OneD, NekDouble>> &vel)
1092 {
1093  int phystot = m_fields[0]->GetTotPoints();
1094  int nel = m_fields[0]->GetNumElmts();
1095  int nvel, cnt;
1096 
1098 
1099  Vmath::Fill(nel, velmag, diffcoeff, 1);
1100 
1101  if (vel != NullNekDoubleArrayOfArray)
1102  {
1103  Array<OneD, NekDouble> Velmag(phystot);
1104  nvel = vel.size();
1105  // calculate magnitude of v
1106  Vmath::Vmul(phystot, vel[0], 1, vel[0], 1, Velmag, 1);
1107  for (int n = 1; n < nvel; ++n)
1108  {
1109  Vmath::Vvtvp(phystot, vel[n], 1, vel[n], 1, Velmag, 1, Velmag, 1);
1110  }
1111  Vmath::Vsqrt(phystot, Velmag, 1, Velmag, 1);
1112 
1113  cnt = 0;
1115  // calculate mean value of vel mag.
1116  for (int i = 0; i < nel; ++i)
1117  {
1118  int nq = m_fields[0]->GetExp(i)->GetTotPoints();
1119  tmp = Velmag + cnt;
1120  diffcoeff[i] = m_fields[0]->GetExp(i)->Integral(tmp);
1121  Vmath::Fill(nq, 1.0, tmp, 1);
1122  NekDouble area = m_fields[0]->GetExp(i)->Integral(tmp);
1123  diffcoeff[i] = diffcoeff[i] / area;
1124  cnt += nq;
1125  }
1126  }
1127  else
1128  {
1129  nvel = m_expdim;
1130  }
1131 
1132  for (int e = 0; e < nel; e++)
1133  {
1134  LocalRegions::ExpansionSharedPtr exp = m_fields[0]->GetExp(e);
1135  NekDouble h = 0;
1136 
1137  // Find maximum length of edge.
1138  for (int i = 0; i < exp->GetGeom()->GetNumEdges(); ++i)
1139  {
1140  h = max(h, exp->GetGeom()->GetEdge(i)->GetVertex(0)->dist(
1141  *(exp->GetGeom()->GetEdge(i)->GetVertex(1))));
1142  }
1143 
1144  int p = 0;
1145  for (int i = 0; i < m_expdim; ++i)
1146  {
1147  p = max(p, exp->GetBasisNumModes(i) - 1);
1148  }
1149 
1150  diffcoeff[e] *= h / p;
1151  }
1152 }
1153 
1155  StdRegions::ConstFactorMap &factors,
1156  MultiRegions::VarFactorsMap &varFactorsMap)
1157 {
1158 
1159  if (m_useSpecVanVisc)
1160  {
1164  {
1165  if (m_IsSVVPowerKernel)
1166  {
1169  }
1170  else
1171  {
1172  varFactorsMap[StdRegions::eFactorSVVDGKerDiffCoeff] =
1174  }
1175  }
1176  }
1177 }
1178 } // namespace Nektar
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:215
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
Definition: ErrorUtil.hpp:249
This class is the base class for Navier Stokes problems.
void SetBoundaryConditions(NekDouble time)
time dependent boundary conditions updating
MultiRegions::ExpListSharedPtr m_pressure
Pointer to field holding pressure field.
NekDouble m_kinvis
Kinematic viscosity.
bool m_SmoothAdvection
bool to identify if advection term smoothing is requested
virtual void v_InitObject(bool DeclareField=true) override
Init object for UnsteadySystem class.
ExtrapolateSharedPtr m_extrapolation
Array< OneD, int > m_velocity
int which identifies which components of m_fields contains the velocity (u,v,w);
void EvaluateAdvectionTerms(const Array< OneD, const Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray, const NekDouble time)
EquationType m_equationType
equation type;
int m_nConvectiveFields
Number of fields to be convected;.
std::vector< SolverUtils::ForcingSharedPtr > m_forcing
Forcing terms.
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:198
tBaseSharedPtr CreateInstance(tKey idKey, tParam... args)
Create an instance of the class referred to by idKey.
Definition: NekFactory.hpp:144
void DefineOdeRhs(FuncPointerT func, ObjectPointerT obj)
void DefineImplicitSolve(FuncPointerT func, ObjectPointerT obj)
void AccumulateRegion(std::string, int iolevel=0)
Accumulate elapsed time for a region.
Definition: Timer.cpp:73
SolverUtils::AdvectionSharedPtr m_advObject
Advection term.
virtual SOLVER_UTILS_EXPORT bool v_PostIntegrate(int step)
int m_spacedim
Spatial dimension (>= expansion dim).
int m_expdim
Expansion dimension.
LibUtilities::CommSharedPtr m_comm
Communicator.
NekDouble m_timestep
Time step size.
NekDouble m_time
Current time of simulation.
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
bool m_specHP_dealiasing
Flag to determine if dealisising is usde for the Spectral/hp element discretisation.
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
enum HomogeneousType m_HomogeneousType
bool m_homogen_dealiasing
Flag to determine if dealiasing is used for homogeneous simulations.
SOLVER_UTILS_EXPORT int GetNcoeffs()
LibUtilities::FieldMetaDataMap m_fieldMetaDataMap
Map to identify relevant solver info to dump in output fields.
SOLVER_UTILS_EXPORT SessionFunctionSharedPtr GetFunction(std::string name, const MultiRegions::ExpListSharedPtr &field=MultiRegions::NullExpListSharedPtr, bool cache=false)
Get a SessionFunction by name.
SOLVER_UTILS_EXPORT int GetTotPoints()
Defines a forcing term to be explicitly applied.
Definition: Forcing.h:73
Base class for unsteady solvers.
LibUtilities::TimeIntegrationSchemeOperators m_ode
The time integration scheme operators to use.
virtual SOLVER_UTILS_EXPORT void v_DoInitialise()
Sets up initial conditions.
LibUtilities::TimeIntegrationSchemeSharedPtr m_intScheme
Wrapper to the time integration scheme.
bool m_explicitDiffusion
Indicates if explicit or implicit treatment of diffusion is used.
virtual SOLVER_UTILS_EXPORT void v_GenerateSummary(SummaryList &s)
Print a summary of time stepping parameters.
NekDouble m_greenFlux
Flux of the Stokes function solution.
virtual std::string v_GetExtrapolateStr(void)
MultiRegions::ExpListSharedPtr m_flowrateBnd
Flowrate reference surface.
Array< OneD, NekDouble > m_svvVarDiffCoeff
Array of coefficient if power kernel is used in SVV.
void SetUpPressureForcing(const Array< OneD, const Array< OneD, NekDouble >> &fields, Array< OneD, Array< OneD, NekDouble >> &Forcing, const NekDouble aii_Dt)
void EvaluateAdvection_SetPressureBCs(const Array< OneD, const Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray, const NekDouble time)
NekDouble m_sVVDiffCoeffHomo1D
Diffusion coefficient of SVV modes in homogeneous 1D Direction.
NekDouble m_flowrate
Desired volumetric flowrate.
bool m_IsSVVPowerKernel
Identifier for Power Kernel otherwise DG kernel.
NekDouble MeasureFlowrate(const Array< OneD, Array< OneD, NekDouble >> &inarray)
Measure the volumetric flow rate through the volumetric flow rate reference surface.
virtual void v_GenerateSummary(SolverUtils::SummaryList &s)
Print a summary of time stepping parameters.
virtual void v_DoInitialise(void)
Sets up initial conditions.
NekDouble m_sVVCutoffRatio
cutt off ratio from which to start decayhing modes
NekDouble m_sVVDiffCoeff
Diffusion coefficient of SVV modes.
void SolveUnsteadyStokesSystem(const Array< OneD, const Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray, const NekDouble time, const NekDouble a_iixDt)
Array< OneD, Array< OneD, NekDouble > > m_flowrateStokes
Stokes solution used to impose flowrate.
NekDouble m_alpha
Current flowrate correction.
bool m_useGJPStabilisation
bool to identify if GJP semi-implicit is active.
virtual Array< OneD, bool > v_GetSystemSingularChecks()
virtual void v_TransPhysToCoeff(void)
Virtual function for transformation to coefficient space.
virtual void v_SolvePressure(const Array< OneD, NekDouble > &Forcing)
void SolveViscous(const Array< OneD, const Array< OneD, NekDouble >> &Forcing, const Array< OneD, const Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray, const NekDouble aii_Dt)
void SetUpViscousForcing(const Array< OneD, const Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &Forcing, const NekDouble aii_Dt)
Array< OneD, Array< OneD, NekDouble > > m_F
virtual void v_EvaluateAdvection_SetPressureBCs(const Array< OneD, const Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray, const NekDouble time)
virtual void v_TransCoeffToPhys(void)
Virtual function for transformation to physical space.
virtual void v_InitObject(bool DeclareField=true)
Init object for UnsteadySystem class.
void SetupFlowrate(NekDouble aii_dt)
Set up the Stokes solution used to impose constant flowrate through a boundary.
bool m_useSpecVanVisc
bool to identify if spectral vanishing viscosity is active.
virtual void v_SetUpViscousForcing(const Array< OneD, const Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &Forcing, const NekDouble aii_Dt)
virtual void v_SetUpPressureForcing(const Array< OneD, const Array< OneD, NekDouble >> &fields, Array< OneD, Array< OneD, NekDouble >> &Forcing, const NekDouble aii_Dt)
void SVVVarDiffCoeff(const NekDouble velmag, Array< OneD, NekDouble > &diffcoeff, const Array< OneD, Array< OneD, NekDouble >> &vel=NullNekDoubleArrayOfArray)
Array< OneD, NekDouble > m_diffCoeff
Diffusion coefficients (will be kinvis for velocities)
int m_flowrateSteps
Interval at which to record flowrate data.
NekDouble m_flowrateArea
Area of the boundary through which we are measuring the flowrate.
std::ofstream m_flowrateStream
Output stream to record flowrate.
void AppendSVVFactors(StdRegions::ConstFactorMap &factors, MultiRegions::VarFactorsMap &varFactorsMap)
NekDouble m_flowrateAiidt
Value of aii_dt used to compute Stokes flowrate solution.
virtual void v_SolveViscous(const Array< OneD, const Array< OneD, NekDouble >> &Forcing, const Array< OneD, const Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray, const NekDouble aii_Dt)
void SolvePressure(const Array< OneD, NekDouble > &Forcing)
int m_planeID
Plane ID for cases with homogeneous expansion.
bool m_useHomo1DSpecVanVisc
bool to identify if spectral vanishing viscosity is active.
bool m_useGJPNormalVel
bool to identify if GJP normal Velocity should be applied in explicit approach
int m_flowrateBndID
Boundary ID of the flowrate reference surface.
virtual std::string v_GetSubSteppingExtrapolateStr(const std::string &instr)
static void Dscal(const int &n, const double &alpha, double *x, const int &incx)
BLAS level 1: x = alpha x.
Definition: Blas.hpp:168
static void Daxpy(const int &n, const double &alpha, const double *x, const int &incx, const double *y, const int &incy)
BLAS level 1: y = alpha x plus y.
Definition: Blas.hpp:154
std::shared_ptr< SessionReader > SessionReaderSharedPtr
std::shared_ptr< Equation > EquationSharedPtr
Definition: Equation.h:130
std::shared_ptr< Expansion > ExpansionSharedPtr
Definition: Expansion.h:68
MultiRegions::Direction const DirCartesianMap[]
Definition: ExpList.h:89
static VarFactorsMap NullVarFactorsMap
std::map< StdRegions::ConstFactorType, Array< OneD, NekDouble > > VarFactorsMap
std::shared_ptr< GJPStabilisation > GJPStabilisationSharedPtr
std::shared_ptr< ContField > ContFieldSharedPtr
Definition: ContField.h:277
std::vector< std::pair< std::string, std::string > > SummaryList
Definition: Misc.h:48
EquationSystemFactory & GetEquationSystemFactory()
void AddSummaryItem(SummaryList &l, const std::string &name, const std::string &value)
Adds a summary item to the summary info list.
Definition: Misc.cpp:49
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:172
std::map< StdRegions::VarCoeffType, Array< OneD, NekDouble > > VarCoeffMap
Definition: StdRegions.hpp:240
std::map< ConstFactorType, NekDouble > ConstFactorMap
Definition: StdRegions.hpp:282
static VarCoeffMap NullVarCoeffMap
Definition: StdRegions.hpp:241
The above copyright notice and this permission notice shall be included.
Definition: CoupledSolver.h:1
@ eUnsteadyStokes
@ eUnsteadyNavierStokes
static Array< OneD, Array< OneD, NekDouble > > NullNekDoubleArrayOfArray
ExtrapolateFactory & GetExtrapolateFactory()
Definition: Extrapolate.cpp:48
static Array< OneD, NekDouble > NullNekDouble1DArray
double NekDouble
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
Definition: Vmath.cpp:534
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 Vabs(int n, const T *x, const int incx, T *y, const int incy)
vabs: y = |x|
Definition: Vmath.cpp:553
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
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 Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:492
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition: Vmath.cpp:45