Nektar++
VelocityCorrectionSchemeImplicit.cpp
Go to the documentation of this file.
1///////////////////////////////////////////////////////////////////////////////
2//
3// File VelocityCorrectionSchemeImplicit.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: Implicit Velocity Correction Scheme for the Incompressible
32// Navier Stokes equations similar to Dong and Shen 2010
33///////////////////////////////////////////////////////////////////////////////
34
39
40#include <boost/algorithm/string.hpp>
41
42using namespace std;
43
44namespace Nektar
45{
46using namespace MultiRegions;
47
50 "VCSImplicit", VCSImplicit::create);
51
53 LibUtilities::SessionReader::RegisterEnumValue("SolverType", "VCSImplicit",
55
56/**
57 * Constructor. No special instructions here.
58 * v_DoInitialise sets the scheme up.
59 *
60 * \param
61 * \param
62 */
65 : UnsteadySystem(pSession, pGraph),
66 VelocityCorrectionScheme(pSession, pGraph)
67{
68}
69
70/**
71 *
72 */
74{
76 SolverUtils::AddSummaryItem(s, "Splitting Scheme",
77 "Implicit velocity correction");
78
79 if (m_extrapolation->GetSubStepName().size())
80 {
81 SolverUtils::AddSummaryItem(s, "Substepping",
82 m_extrapolation->GetSubStepName());
83 }
84
85 string dealias = m_homogen_dealiasing ? "Homogeneous1D" : "";
87 {
88 dealias += (dealias == "" ? "" : " + ") + string("spectral/hp");
89 }
90 if (dealias != "")
91 {
92 SolverUtils::AddSummaryItem(s, "Dealiasing", dealias);
93 }
94
95 string smoothing = m_useSpecVanVisc ? "spectral/hp" : "";
96 if (smoothing != "")
97 {
99 {
101 s, "Smoothing-SpecHP",
102 "SVV (" + smoothing + " Exp Kernel(cut-off = " +
103 boost::lexical_cast<string>(m_sVVCutoffRatio) +
104 ", diff coeff = " +
105 boost::lexical_cast<string>(m_sVVDiffCoeff) + "))");
106 }
107 else
108 {
110 {
112 s, "Smoothing-SpecHP",
113 "SVV (" + smoothing + " Power Kernel (Power ratio =" +
114 boost::lexical_cast<string>(m_sVVCutoffRatio) +
115 ", diff coeff = " +
116 boost::lexical_cast<string>(m_sVVDiffCoeff) +
117 "*Uh/p))");
118 }
119 else
120 {
122 s, "Smoothing-SpecHP",
123 "SVV (" + smoothing + " DG Kernel (diff coeff = " +
124 boost::lexical_cast<string>(m_sVVDiffCoeff) +
125 "*Uh/p))");
126 }
127 }
128 }
129
131 {
133 s, "Smoothing-Homo1D",
134 "SVV (Homogeneous1D - Exp Kernel(cut-off = " +
135 boost::lexical_cast<string>(m_sVVCutoffRatioHomo1D) +
136 ", diff coeff = " +
137 boost::lexical_cast<string>(m_sVVDiffCoeffHomo1D) + "))");
138 }
139
141 {
143 s, "GJP Stab. Impl. ",
144 m_session->GetSolverInfo("GJPStabilisation"));
145 SolverUtils::AddSummaryItem(s, "GJP Stab. JumpScale", m_GJPJumpScale);
146
147 if (boost::iequals(m_session->GetSolverInfo("GJPStabilisation"),
148 "Explicit"))
149 {
151 s, "GJP Normal Velocity",
152 m_session->GetSolverInfo("GJPNormalVelocity"));
153 }
154 }
155}
156
157/**
158 *
159 */
160void VCSImplicit::v_DoInitialise(bool dumpInitialConditions)
161{
162 VelocityCorrectionScheme::v_DoInitialise(dumpInitialConditions);
163
164 // Get integration order for extrapolation
165 m_intOrder = m_intScheme->GetOrder();
166
167 // Initialise Advection velocity
169 for (int i = 0; i < m_velocity.size(); i++)
170 {
172 }
173
174 // Initialise storage for explicit advection terms
176 for (int i = 0; i < m_nConvectiveFields; i++)
177 {
178 m_advection[i] =
180 }
181
182 // Check advection velocity: Extrapolated or Updated scheme
183 // Default to Extrapolated scheme
184 m_session->MatchSolverInfo("AdvectionVelocity", "Extrapolated",
185 m_advectionVelocity, true);
186
187 // Extrapolated velocity (for 2nd order scheme)
188 // Only initiliase for Extrapolated scheme
190 {
191 m_extVel =
193 for (int i = 0; i < m_velocity.size(); i++)
194 {
196 for (int j = 0; j < m_intOrder; j++)
197 {
198 m_extVel[i][j] =
200 }
201 }
202 }
203
204 // Check whether multiple or single GlobalLinSys are
205 // required and when to delete them
206 // Initialise to always reset GlobalLinSys
209
210 // Check skew-symmetric advection operator
211 // Default to Convective operator = false
212 m_session->MatchSolverInfo("AdvectionOperator", "SkewSymmetric",
214}
215
216/**
217 * Computes the forcing term for the pressure solve in
218 * VCSImplicit::v_SolvePressure(). It uses the weak pressure forcing similar to
219 * VCSWeakPressure (coefficient space).
220 *
221 * @param fields Holds the BDF formula of previous solutions \f$\sum_q
222 * \frac{\alpha_q}{\gamma} [u^{n-q}, v^{n-q}, \ldots]\f$.
223 * @param Forcing Array for the forcing term of the pressure solve. May
224 * contain temporary values on input. On output, Forcing = \f$ \int_\Omega
225 * \nabla \phi \cdot \left( \frac{1}{\Delta t} \sum_q \alpha_q \mathbf{u}^{n-q}
226 * + \mathbf{f}^{n+1} - N(\mathbf{u})^n \right) \f$.
227 *
228 */
230 const Array<OneD, const Array<OneD, NekDouble>> &fields,
232{
233 int ncoeffs = m_pressure->GetNcoeffs();
234 int phystot = m_fields[0]->GetTotPoints();
235 int nvel = m_velocity.size();
236
237 // Add 1/dt \sum_q \alpha_q u^{n-q} - N(u^n)
238 for (int i = 0; i < nvel; i++)
239 {
240 Vmath::Svtvp(phystot, 1.0 / aii_Dt, fields[i], 1, m_advection[i], 1,
241 Forcing[i], 1);
242 }
243
244 // Add forcing terms
245 for (auto &x : m_forcing)
246 {
248 }
249
250 // Do Innerproduct with derivative base
251 m_pressure->IProductWRTDerivBase(Forcing, Forcing[0]);
252
253 // Negate to reverse negation in HelmSolve
254 Vmath::Neg(ncoeffs, Forcing[0], 1);
255}
256
257/**
258 * Computes the forcing term and advection velocity for ADR solves.
259 *
260 * @param inarray Holds the BDF formula of previous solutions \f$\sum_q
261 \frac{\alpha_q}{\gamma} [u^{n-q}, v^{n-q}, \ldots]\f$.
262 * @param Forcing Shared forcing array m_F. Used in pressure and velocity
263 solve. On input, holds temporary values from pressure solve. On output: Forcing
264 = \f$ - \frac{1}{\nu} \left( \frac{1}{\Delta t} \sum_q \alpha_q
265 \mathbf{u}^{n-q} - \nabla p^{n+1} + \mathbf{f}^{n+1} \right) \f$.
266
267 * @param aii_Dt \f$\frac{\Delta t}{\gamma}\f$.
268 *
269 * Additionally, the advection velocity \f$ \tilde{\mathbf{u}}^{n+1} \f$ is
270 computed and handed to the velocity solve via the array m_AdvVel. It is
271 computed as \f$ \tilde{\mathbf{u}}^{n+1} = \sum_q \frac{\alpha_q}{\gamma}
272 \mathbf{u}^{n-q} - \frac{\Delta t}{\gamma} \left[ \nabla p^{n+1} -
273 N(\mathbf{u})^n - \nu \nabla \times \omega^n + \mathbf{f}^{n+1} \right] \f$.
274 *
275 */
277 const Array<OneD, const Array<OneD, NekDouble>> &inarray,
279{
280 int phystot = m_fields[0]->GetTotPoints();
281 int nvel = m_velocity.size();
282
283 // Update (PhysSpace) pressure to n+1
284 m_pressure->BwdTrans(m_pressure->GetCoeffs(), m_pressure->UpdatePhys());
285
286 // Compute gradient \nabla p^{n+1}
287 if (nvel == 2)
288 {
289 m_pressure->PhysDeriv(m_pressure->GetPhys(), Forcing[m_velocity[0]],
290 Forcing[m_velocity[1]]);
291 }
292 else
293 {
294 m_pressure->PhysDeriv(m_pressure->GetPhys(), Forcing[m_velocity[0]],
296 }
297
298 // Negate pressure gradient: -1 * \nabla p
299 for (int i = 0; i < nvel; ++i)
300 {
301 Vmath::Neg(phystot, Forcing[i], 1);
302 }
303
304 // Zero convective fields
305 for (int i = nvel; i < m_nConvectiveFields; ++i)
306 {
307 Vmath::Zero(phystot, Forcing[i], 1);
308 }
309
310 // Add forcing terms: += f^{n+1}
311 for (auto &x : m_forcing)
312 {
314 }
315
316 // Build Advection Velocity
317 for (int i = 0; i < nvel; ++i)
318 {
319 // Extrapolated Advection Velocity = u^*n+1
321 {
322 Vmath::Smul(phystot, 1.0 / m_diffCoeff[i],
323 m_extVel[i][m_intOrder - 1], 1, m_AdvVel[i], 1);
324 }
325 // Updated Advection Velocity = u^*n+1 - (...)
326 else
327 {
328 // Get velocity fields [u^n, v^n, ..] for curl of vorticity
330 for (int i = 0; i < nvel; i++)
331 {
332 velocity[i] = Array<OneD, NekDouble>(m_fields[i]->GetPhys());
333 }
334 // Curl of vorticity \nabla \times \nabla \times [u^n, v^n, ..]
335 m_fields[0]->CurlCurl(velocity, velocity);
336
337 // \frac{\Delta t}{\gamma}(-\nabla p + f) + inarray
338 Vmath::Svtvp(phystot, aii_Dt, Forcing[i], 1, inarray[i], 1,
339 m_AdvVel[i], 1);
340 // += \frac{\Delta t}{\gamma} -N(u)
341 Vmath::Svtvp(phystot, aii_Dt, m_advection[i], 1, m_AdvVel[i], 1,
342 m_AdvVel[i], 1);
343 // += -\frac{\Delta t}{\gamma} CurlCurl(u)
344 Vmath::Svtvp(phystot, -aii_Dt * m_diffCoeff[i], velocity[i], 1,
345 m_AdvVel[i], 1, m_AdvVel[i], 1);
346 // *= 1/\nu
347 Vmath::Smul(phystot, 1.0 / m_diffCoeff[i], m_AdvVel[i], 1,
348 m_AdvVel[i], 1);
349 }
350 }
351
352 // Build Forcing term
353 for (int i = 0; i < nvel; ++i)
354 {
355 // \frac{\gamma}{\Delta t} * inarray - \nabla p + f
356 Vmath::Svtvp(phystot, 1.0 / aii_Dt, inarray[i], 1, Forcing[i], 1,
357 Forcing[i], 1);
358
359 // *= -1/kinvis
360 Vmath::Smul(phystot, -1.0 / m_diffCoeff[i], Forcing[i], 1, Forcing[i],
361 1);
362 }
363}
364
365/**
366 * Solve pressure system via a Poisson problem with ContField::HelmSolve().
367 * Uses coefficient space forcing and hence sets the
368 * argument PhysSpaceForcing=false in ContFieldHelmSolve().
369 *
370 * @param Forcing See output description in
371 * VCSImplicit::v_SetUpPressureForcing().
372 *
373 */
375{
377 // Setup coefficient for equation, Lambda = 0 ensures Poisson instead of
378 // Helmholtz problem
380
381 // Solve Pressure Poisson Equation (with Weak Forcing)
382 m_pressure->HelmSolve(Forcing, m_pressure->UpdateCoeffs(), factors,
385 false);
386
387 // Add presure to outflow bc if using convective like BCs
388 m_extrapolation->AddPressureToOutflowBCs(m_kinvis);
389}
390
391/**
392 * Solve velocity system via the
393 * ContField::LinearAdvectionDiffusionReactionSolve().
394 *
395 * @param Forcing Holds the forcing term for each velocity component. See output
396 * description in VCSImplicit::v_SetUpViscousForcing().
397 * @param inarray Unused in this routine. Still, holds \f$\sum_q
398 * \frac{\alpha_q}{\gamma} [u^{n-q}, v^{n-q}, \ldots] \f$.
399 * @param outarray Holds vector of previous solution \f$ [u^n, v^n, \ldots] \f$.
400 * Overwritten upon output with new solution [u^{n+1}, v^{n+1}, \ldots].
401 *
402 * Additionally, the advection velocity \f$ \tilde{\mathbf{u}}^{n+1} \f$ is
403 * handed to the velocity solve via the array m_AdvVel. It is computed in
404 * VCSImplicit::v_SetUpViscousForcing(). It is defined as \f$
405 * \tilde{\mathbf{u}}^{n+1} = \sum_q \frac{\alpha_q}{\gamma} \mathbf{u}^{n-q} -
406 * \frac{\Delta t}{\gamma} \left[ \nabla p^{n+1} - N(\mathbf{u})^n - \nu \nabla
407 * \times \omega^n + \mathbf{f}^{n+1} \right] \f$.
408 *
409 */
412 const Array<OneD, const Array<OneD, NekDouble>> &inarray,
413 Array<OneD, Array<OneD, NekDouble>> &outarray, const NekDouble aii_Dt)
414{
416 StdRegions::VarCoeffMap varcoeffs;
418
419 AppendSVVFactors(factors, varfactors);
420 ComputeGJPNormalVelocity(inarray, varcoeffs);
421
422 // Compute mass varcoeff for implicit Skew-symmetric operator
423 AddImplicitSkewSymAdvection(varcoeffs, aii_Dt);
424
425 // Set advection velocities
429 for (int i = 0; i < m_velocity.size(); i++)
430 {
431 varcoeffs[varcoefftypes[i]] = m_AdvVel[i];
432 }
433
434 // Solve Advection-Diffusion-Reaction system
435 for (int i = 0; i < m_nConvectiveFields; ++i)
436 {
437 // Add diffusion coefficient to GJP matrix operator (Implicit part)
439 {
441 }
442
443 // \lambda = - /frac{\gamma}{\Delta t \nu}
444 factors[StdRegions::eFactorLambda] = -1.0 / aii_Dt / m_diffCoeff[i];
445
446 auto gkey = m_fields[i]->LinearAdvectionDiffusionReactionSolve(
447 Forcing[i], m_fields[i]->UpdateCoeffs(), factors, varcoeffs,
448 varfactors);
449
450 // Nuke GlobalLinSys to avoid memory leak
451 // The condition nukes, if:
452 // - Remove after last velocity solve (w-velocity in 3D,
453 // assumes same matrix for each velocity)
454 // - Catches Null return from 3DH1D solve by checking matrix type
455 // is a LinearADR matrix (variant)
456 if (m_unsetGlobalLinSys[i] &&
457 (gkey.GetMatrixType() ==
459 gkey.GetMatrixType() ==
461 {
462 m_fields[i]->UnsetGlobalLinSys(gkey, true);
463 }
464
465 // Transform solution to PhysSpace
466 m_fields[i]->BwdTrans(m_fields[i]->GetCoeffs(), outarray[i]);
467 }
468}
469
470/**
471 * Explicit part of the method - HOPBCs based on VCSWeakPressure.
472 * For the implicit scheme, we do not extrapolate the advection and rotational
473 * diffusion terms, see also S. Dong and J. Shen (2010).
474 *
475 * @param inarray Holds the vector of previous solutions \f$ [u^n, v^n, \ldots]
476 * \f$.
477 * @param outarray In the implicit scheme, currently an empty array. TODO Don't
478 * waste memory here.
479 * @param time Holds the new time \f$ t^{n+1} = (n+1) \Delta t \f$.
480 *
481 */
483 const Array<OneD, const Array<OneD, NekDouble>> &inarray,
484 Array<OneD, Array<OneD, NekDouble>> &outarray, const NekDouble time)
485{
486 // Zero RHS derivatives, avoids undefined values
487 for (int i = 0; i < m_velocity.size(); ++i)
488 {
489 Vmath::Zero(outarray[i].size(), outarray[i], 1);
490 }
491
492 // Calculate High-Order pressure boundary conditions
494 timer.Start();
495 m_extrapolation->EvaluatePressureBCs(inarray, outarray, m_kinvis);
496 timer.Stop();
497 timer.AccumulateRegion("Pressure BCs");
498
499 // Extrapolate advection velocity with StifflyStableBeta coefficients
501 {
502 for (int i = 0; i < m_velocity.size(); ++i)
503 {
504 Vmath::Vcopy(inarray[i].size(), inarray[i], 1,
505 m_extVel[i][m_intOrder - 1], 1);
506 m_extrapolation->ExtrapolateArray(m_extVel[i]);
507 }
508 }
509
510 // Compute explicit advection terms
511 EvaluateAdvectionTerms(inarray, m_advection, time);
512}
513
515 NekDouble aii_Dt)
516{
518 {
519 return;
520 }
521
522 // Setup temporary storages
523 int nvel = m_velocity.size();
524 int phystot = m_fields[0]->GetTotPoints();
525 Array<OneD, NekDouble> divAdvVel(phystot, 0.0);
526 Array<OneD, NekDouble> tmpstore(phystot);
527 Array<OneD, NekDouble> ones(phystot, 1.0);
528
529 // Compute divergence of advection velocity: 1/\nu \tilde{u}
530 m_fields[0]->PhysDeriv(eX, m_AdvVel[0], divAdvVel);
531 for (int i = 1; i < nvel; ++i)
532 {
533 m_fields[i]->PhysDeriv(DirCartesianMap[i], m_AdvVel[i], tmpstore);
534 Vmath::Vadd(phystot, divAdvVel, 1, tmpstore, 1, divAdvVel, 1);
535 }
536 // Smul: 1/(2*lambda) * divAdvVel
537 NekDouble lambda = 1.0 / aii_Dt / m_diffCoeff[0];
538 Vmath::Smul(phystot, 1.0 / 2.0 / lambda, divAdvVel, 1, divAdvVel, 1);
539
540 // Vadd: 1 + 1/(2*lambda) * divAdvVel
541 Vmath::Vadd(phystot, ones, 1, divAdvVel, 1, divAdvVel, 1);
542
543 // Assign VarcoeffMass
544 varcoeffs[StdRegions::eVarCoeffMass] = divAdvVel;
545}
546
548 Array<OneD, int32_t> &unsetGlobalLinSys)
549{
550 // Loop over all convected Fields (excludes pressure)
551 for (int i = 0; i < m_nConvectiveFields; ++i)
552 {
553 // Always unset after last solve
554 if (i == m_nConvectiveFields - 1)
555 {
556 unsetGlobalLinSys[i] = 1;
557 continue;
558 }
559
560 auto contField =
561 std::dynamic_pointer_cast<MultiRegions::ContField>(m_fields[i]);
562 // Check against all solves after this one, whether same
563 // Global Linear system is required
564 for (int j = i + 1; j < m_nConvectiveFields; ++j)
565 {
566 // If same BCs (ie GlobalLinSys) in j-solves, save the GlobalLinSys
567 // for re-use
568 if (std::dynamic_pointer_cast<MultiRegions::ContField>(m_fields[j])
569 ->SameTypeOfBoundaryConditions(*contField))
570 {
571 unsetGlobalLinSys[i] = 0;
572 }
573 // If no same BCs (ie GlobalLinSys) in j-solves, UnsetGlobalLinSys()
574 // for i-solve
575 else
576 {
577 unsetGlobalLinSys[i] = 1;
578 }
579 }
580 }
581}
582
583} // namespace Nektar
MultiRegions::ExpListSharedPtr m_pressure
Pointer to field holding pressure field.
NekDouble m_kinvis
Kinematic viscosity.
ExtrapolateSharedPtr m_extrapolation
Array< OneD, int > m_velocity
int which identifies which components of m_fields contains the velocity (u,v,w);
int m_nConvectiveFields
Number of fields to be convected;.
std::vector< SolverUtils::ForcingSharedPtr > m_forcing
Forcing terms.
void EvaluateAdvectionTerms(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
static std::string RegisterEnumValue(std::string pEnum, std::string pString, int pEnumValue)
Registers an enumeration value.
void AccumulateRegion(std::string, int iolevel=0)
Accumulate elapsed time for a region.
Definition: Timer.cpp:70
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 GetTotPoints()
Defines a forcing term to be explicitly applied.
Definition: Forcing.h:71
Base class for unsteady solvers.
LibUtilities::TimeIntegrationSchemeSharedPtr m_intScheme
Wrapper to the time integration scheme.
SOLVER_UTILS_EXPORT void v_GenerateSummary(SummaryList &s) override
Print a summary of time stepping parameters.
Array< OneD, int32_t > m_unsetGlobalLinSys
Array checking whether GlobalLinSys needs to be unset.
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_extVel
3D Array for extrapolated Advection Velocities [dir][time-levle][dof]
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) override
void v_EvaluateAdvection_SetPressureBCs(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time) override
void v_SetUpPressureForcing(const Array< OneD, const Array< OneD, NekDouble > > &fields, Array< OneD, Array< OneD, NekDouble > > &Forcing, const NekDouble aii_Dt) override
void v_GenerateSummary(SolverUtils::SummaryList &s) override
Print a summary of time stepping parameters.
Array< OneD, Array< OneD, NekDouble > > m_AdvVel
void v_DoInitialise(bool dumpInitialConditions) override
Sets up initial conditions.
Array< OneD, Array< OneD, NekDouble > > m_advection
2D Array for Advection Velocities [dir][dof]
bool m_advectionVelocity
bool to identify implicit scheme ie advection velocity
void AddImplicitSkewSymAdvection(StdRegions::VarCoeffMap varcoeffs, NekDouble aii_Dt)
bool m_implicitSkewSymAdvection
bool to identify advection operator
VCSImplicit(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
void CheckUnsetGlobalLinSys(Array< OneD, int32_t > &unsetGlobalLinSys)
int m_intOrder
integer for advection velocity
static std::string className
Name of class.
void v_SetUpViscousForcing(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &Forcing, const NekDouble aii_Dt) override
static SolverUtils::EquationSystemSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
Creates an instance of this class.
void v_SolvePressure(const Array< OneD, NekDouble > &Forcing) override
Array< OneD, NekDouble > m_svvVarDiffCoeff
Array of coefficient if power kernel is used in SVV.
NekDouble m_sVVDiffCoeffHomo1D
Diffusion coefficient of SVV modes in homogeneous 1D Direction.
bool m_IsSVVPowerKernel
Identifier for Power Kernel otherwise DG kernel.
NekDouble m_sVVCutoffRatio
cutt off ratio from which to start decayhing modes
NekDouble m_sVVDiffCoeff
Diffusion coefficient of SVV modes.
bool m_useGJPStabilisation
bool to identify if GJP semi-implicit is active.
bool m_useSpecVanVisc
bool to identify if spectral vanishing viscosity is active.
Array< OneD, NekDouble > m_diffCoeff
Diffusion coefficients (will be kinvis for velocities)
void AppendSVVFactors(StdRegions::ConstFactorMap &factors, MultiRegions::VarFactorsMap &varFactorsMap)
void v_DoInitialise(bool dumpInitialConditions=true) override
Sets up initial conditions.
bool m_useHomo1DSpecVanVisc
bool to identify if spectral vanishing viscosity is active.
void ComputeGJPNormalVelocity(const Array< OneD, const Array< OneD, NekDouble > > &inarray, StdRegions::VarCoeffMap &varcoeffs)
std::shared_ptr< SessionReader > SessionReaderSharedPtr
MultiRegions::Direction const DirCartesianMap[]
Definition: ExpList.h:87
static VarFactorsMap NullVarFactorsMap
std::map< StdRegions::ConstFactorType, Array< OneD, NekDouble > > VarFactorsMap
std::vector< std::pair< std::string, std::string > > SummaryList
Definition: Misc.h:46
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:47
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:174
@ eLinearAdvectionDiffusionReactionGJP
Definition: StdRegions.hpp:105
std::map< ConstFactorType, NekDouble > ConstFactorMap
Definition: StdRegions.hpp:430
static VarCoeffMap NullVarCoeffMap
Definition: StdRegions.hpp:376
std::map< StdRegions::VarCoeffType, VarCoeffEntry > VarCoeffMap
Definition: StdRegions.hpp:375
StdRegions::ConstFactorMap factors
static Array< OneD, NekDouble > NullNekDouble1DArray
double NekDouble
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.hpp:396
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.hpp:292
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.hpp:180
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.hpp:100
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.hpp:273
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.hpp:825
STL namespace.