44#include <boost/algorithm/string.hpp>
50using namespace MultiRegions;
89 ASSERTL0(
false,
"Need to set up pressure field definition");
96 std::string varName =
m_session->GetVariable(n);
97 if (
m_session->DefinesFunction(
"DiffusionCoefficient", varName))
100 m_session->GetFunction(
"DiffusionCoefficient", varName);
115 m_session->MatchSolverInfo(
"GJPStabilisation",
"False",
129 "Can only specify GJPNormalVelocity with"
130 " GJPStabilisation set to Explicit currently");
158 if (
m_session->DefinesSolverInfo(
"Extrapolation"))
161 m_session->GetSolverInfo(
"Extrapolation"));
230 std::string forces[] = {
"X",
"Y",
"Z"};
237 std::string varName = std::string(
"Force") + forces[i];
238 defined =
m_session->DefinesFunction(
"FlowrateForce", varName);
246 "A 'FlowrateForce' function must defined with components "
247 "[ForceX, ...] to define direction of flowrate forcing");
250 m_session->GetFunction(
"FlowrateForce", varName);
251 flowrateForce[i] = ffunc->Evaluate();
266 flowrateForce[2] = 1.0;
270 for (
size_t i = 0; i < bcs.size(); ++i)
272 if (boost::iequals(bcs[i]->GetUserDefined(),
"Flowrate"))
282 "One boundary region must be marked using the 'Flowrate' "
283 "user-defined type to monitor the volumetric flowrate.");
298 for (
size_t i = 0; i < zIDs.size(); ++i)
307 ASSERTL1(tmpId <= 0,
"Should be either at location 0 or -1 if not "
339 for (
size_t i = 0; i < zIDs.size(); ++i)
362 int nqTot =
m_fields[0]->GetNpoints();
374 m_fields[i]->HomogeneousFwdTrans(nqTot, inTmp[i], inTmp2);
390 for (
size_t n = 0; n < PBndConds.size(); ++n)
392 if (PBndConds[n]->GetBoundaryConditionType() ==
419 std::string filename =
m_session->GetSessionName();
424 <<
"# -------------------------------------------"
508 << setw(16) <<
m_alpha << endl;
529 "Velocity correction (strong press. form)");
540 dealias += (dealias ==
"" ?
"" :
" + ") +
string(
"spectral/hp");
553 s,
"Smoothing-SpecHP",
554 "SVV (" + smoothing +
" Exp Kernel(cut-off = " +
564 s,
"Smoothing-SpecHP",
565 "SVV (" + smoothing +
" Power Kernel (Power ratio =" +
574 s,
"Smoothing-SpecHP",
575 "SVV (" + smoothing +
" DG Kernel (diff coeff = " +
585 s,
"Smoothing-Homo1D",
586 "SVV (Homogeneous1D - Exp Kernel(cut-off = " +
595 s,
"GJP Stab. Impl. ",
596 m_session->GetSolverInfo(
"GJPStabilisation"));
599 if (boost::iequals(
m_session->GetSolverInfo(
"GJPStabilisation"),
603 s,
"GJP Normal Velocity",
604 m_session->GetSolverInfo(
"GJPNormalVelocity"));
636 for (
size_t i = 0; i <
m_fields.size(); ++i)
651 size_t nfields =
m_fields.size() - 1;
652 for (
size_t k = 0; k < nfields; ++k)
666 size_t nfields =
m_fields.size() - 1;
667 for (
size_t k = 0; k < nfields; ++k)
680 int vVar =
m_session->GetVariables().size();
682 vChecks[vVar - 1] =
true;
691 return m_session->GetVariables().size() - 1;
719 x->Apply(
m_fields, inarray, outarray, time);
744 size_t physTot =
m_fields[0]->GetTotPoints();
785 m_fields[i]->FwdTransLocalElmt(outarray[i],
803 size_t physTot =
m_fields[0]->GetTotPoints();
808 for (i = 1; i < nvel; ++i)
826 size_t phystot =
m_fields[0]->GetTotPoints();
902 varCoeffMap, varFactorsMap);
910 m_session->MatchSolverInfo(
"SpectralVanishingViscosity",
"PowerKernel",
919 m_session->MatchSolverInfo(
"SpectralVanishingViscositySpectralHP",
929 m_session->MatchSolverInfo(
"SpectralVanishingViscosity",
"DGKernel",
937 m_session->MatchSolverInfo(
"SpectralVanishingViscositySpectralHP",
952 if (
m_session->DefinesFunction(
"SVVVelocityMagnitude"))
954 if (
m_comm->GetRank() == 0)
956 cout <<
"Seting up SVV velocity from "
957 "SVVVelocityMagnitude section in session file"
961 size_t phystot =
m_fields[0]->GetTotPoints();
964 for (
size_t i = 0; i < nvel; ++i)
971 GetFunction(
"SVVVelocityMagnitude")->Evaluate(vars, SVVVelFields);
987 m_session->MatchSolverInfo(
"SpectralVanishingViscosity",
"True",
991 m_session->MatchSolverInfo(
"SpectralVanishingViscosity",
998 m_session->MatchSolverInfo(
"SpectralVanishingViscositySpectralHP",
1003 "SpectralVanishingViscositySpectralHP",
"ExpKernel",
1010 if (
m_session->DefinesSolverInfo(
"SpectralVanishingViscosityHomo1D"))
1012 m_session->MatchSolverInfo(
"SpectralVanishingViscosityHomo1D",
"True",
1016 m_session->MatchSolverInfo(
"SpectralVanishingViscosityHomo1D",
1032 "Expect to have three velocity fields with homogenous expansion");
1039 size_t num_planes = planes.size();
1042 size_t kmodes =
m_fields[0]->GetHomogeneousBasis()->GetNumModes();
1047 for (
size_t n = 0; n < num_planes; ++n)
1049 if (planes[n] > pstart)
1051 fac = (
NekDouble)((planes[n] - kmodes) *
1052 (planes[n] - kmodes)) /
1054 (planes[n] - pstart)));
1059 for (
size_t i = 0; i <
m_velocity.size(); ++i)
1071 size_t phystot =
m_fields[0]->GetTotPoints();
1072 size_t nel =
m_fields[0]->GetNumElmts();
1084 Vmath::Vmul(phystot, vel[0], 1, vel[0], 1, Velmag, 1);
1085 for (
size_t n = 1; n < nvel; ++n)
1087 Vmath::Vvtvp(phystot, vel[n], 1, vel[n], 1, Velmag, 1, Velmag, 1);
1094 for (
size_t i = 0; i < nel; ++i)
1096 size_t nq =
m_fields[0]->GetExp(i)->GetTotPoints();
1098 diffcoeff[i] =
m_fields[0]->GetExp(i)->Integral(tmp);
1101 diffcoeff[i] = diffcoeff[i] / area;
1110 for (
size_t e = 0; e < nel; e++)
1116 size_t nEdge = exp->GetGeom()->GetNumEdges();
1117 for (
size_t i = 0; i < nEdge; ++i)
1119 h = max(h, exp->GetGeom()->GetEdge(i)->GetVertex(0)->dist(
1120 *(exp->GetGeom()->GetEdge(i)->GetVertex(1))));
1126 p = max(
p, exp->GetBasisNumModes(i) - 1);
1129 diffcoeff[e] *= h /
p;
1171 std::dynamic_pointer_cast<MultiRegions::ContField>(
m_fields[0]);
1174 cfield->GetGJPForcing();
1176 int nTracePts = GJPData->GetNumTracePts();
1180 GJPData->GetTraceNormals();
1182 m_fields[0]->GetFwdBwdTracePhys(inarray[0], Fwd, Bwd,
true,
true);
1183 Vmath::Vmul(nTracePts, Fwd, 1, traceNormals[0], 1, unorm, 1);
1186 for (
int f = 1; f <
m_fields[0]->GetCoordim(0); ++f)
1188 m_fields[0]->GetFwdBwdTracePhys(inarray[f], Fwd, Bwd,
true,
true);
1189 Vmath::Vvtvp(nTracePts, Fwd, 1, traceNormals[f], 1, unorm, 1, unorm,
#define ASSERTL0(condition, msg)
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
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
void v_InitObject(bool DeclareField=true) override
Initialisation object for EquationSystem.
ExtrapolateSharedPtr m_extrapolation
Array< OneD, int > m_velocity
int which identifies which components of m_fields contains the velocity (u,v,w);
EquationType m_equationType
equation type;
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.
tBaseSharedPtr CreateInstance(tKey idKey, tParam... args)
Create an instance of the class referred to by idKey.
static std::string RegisterEnumValue(std::string pEnum, std::string pString, int pEnumValue)
Registers an enumeration value.
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.
SolverUtils::AdvectionSharedPtr m_advObject
Advection term.
SOLVER_UTILS_EXPORT bool v_PostIntegrate(int step) override
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.
Base class for unsteady solvers.
LibUtilities::TimeIntegrationSchemeOperators m_ode
The time integration scheme operators to use.
SOLVER_UTILS_EXPORT void v_DoInitialise(bool dumpInitialConditions=true) override
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.
SOLVER_UTILS_EXPORT void v_GenerateSummary(SummaryList &s) override
Print a summary of time stepping parameters.
std::vector< int > m_intVariables
NekDouble m_sVVCutoffRatioHomo1D
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 v_TransCoeffToPhys(void) override
Virtual function for transformation to physical space.
NekDouble MeasureFlowrate(const Array< OneD, Array< OneD, NekDouble > > &inarray)
Measure the volumetric flow rate through the volumetric flow rate reference surface.
NekDouble m_sVVDiffCoeffHomo1D
Diffusion coefficient of SVV modes in homogeneous 1D Direction.
void v_TransPhysToCoeff(void) override
Virtual function for transformation to coefficient space.
static std::string className
Name of class.
NekDouble m_flowrate
Desired volumetric flowrate.
~VelocityCorrectionScheme() override
bool m_IsSVVPowerKernel
Identifier for Power Kernel otherwise DG kernel.
bool v_PostIntegrate(int step) override
static SolverUtils::EquationSystemSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
Creates an instance of this class.
NekDouble m_sVVCutoffRatio
cutt off ratio from which to start decayhing modes
NekDouble m_sVVDiffCoeff
Diffusion coefficient of SVV modes.
void v_GenerateSummary(SolverUtils::SummaryList &s) override
Print a summary of time stepping parameters.
void EvaluateAdvection_SetPressureBCs(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
Array< OneD, Array< OneD, NekDouble > > m_flowrateStokes
Stokes solution used to impose flowrate.
void SolveUnsteadyStokesSystem(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time, const NekDouble a_iixDt)
NekDouble m_alpha
Current flowrate correction.
bool m_useGJPStabilisation
bool to identify if GJP semi-implicit is active.
virtual void v_SolvePressure(const Array< OneD, NekDouble > &Forcing)
VelocityCorrectionScheme(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
Constructor.
Array< OneD, Array< OneD, NekDouble > > m_F
void SetUpViscousForcing(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &Forcing, const NekDouble aii_Dt)
virtual void v_SolveUnsteadyStokesSystem(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time, const NekDouble a_iixDt)
void SetupFlowrate(NekDouble aii_dt)
Set up the Stokes solution used to impose constant flowrate through a boundary.
void SetUpPressureForcing(const Array< OneD, const Array< OneD, NekDouble > > &fields, Array< OneD, Array< OneD, NekDouble > > &Forcing, const NekDouble aii_Dt)
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)
virtual void v_SetUpPressureForcing(const Array< OneD, const Array< OneD, NekDouble > > &fields, Array< OneD, Array< OneD, NekDouble > > &Forcing, const NekDouble aii_Dt)
static std::string solverTypeLookupId
bool m_useSpecVanVisc
bool to identify if spectral vanishing viscosity is active.
void SetUpExtrapolation(void)
Array< OneD, NekDouble > m_diffCoeff
Diffusion coefficients (will be kinvis for velocities)
int m_flowrateSteps
Interval at which to record flowrate data.
void SVVVarDiffCoeff(const NekDouble velmag, Array< OneD, NekDouble > &diffcoeff, const Array< OneD, Array< OneD, NekDouble > > &vel=NullNekDoubleArrayOfArray)
int v_GetForceDimension() override
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)
NekDouble m_flowrateArea
Area of the boundary through which we are measuring the flowrate.
std::ofstream m_flowrateStream
Output stream to record flowrate.
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_EvaluateAdvection_SetPressureBCs(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
void v_InitObject(bool DeclareField=true) override
Initialisation object for EquationSystem.
void AppendSVVFactors(StdRegions::ConstFactorMap &factors, MultiRegions::VarFactorsMap &varFactorsMap)
NekDouble m_flowrateAiidt
Value of aii_dt used to compute Stokes flowrate solution.
void SolvePressure(const Array< OneD, NekDouble > &Forcing)
void v_DoInitialise(bool dumpInitialConditions=true) override
Sets up initial conditions.
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
Array< OneD, bool > v_GetSystemSingularChecks() override
int m_flowrateBndID
Boundary ID of the flowrate reference surface.
virtual std::string v_GetSubSteppingExtrapolateStr(const std::string &instr)
void ComputeGJPNormalVelocity(const Array< OneD, const Array< OneD, NekDouble > > &inarray, StdRegions::VarCoeffMap &varcoeffs)
static void Dscal(const int &n, const double &alpha, double *x, const int &incx)
BLAS level 1: x = alpha x.
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.
std::shared_ptr< SessionReader > SessionReaderSharedPtr
std::shared_ptr< Equation > EquationSharedPtr
std::shared_ptr< Expansion > ExpansionSharedPtr
MultiRegions::Direction const DirCartesianMap[]
static VarFactorsMap NullVarFactorsMap
std::map< StdRegions::ConstFactorType, Array< OneD, NekDouble > > VarFactorsMap
std::shared_ptr< GJPStabilisation > GJPStabilisationSharedPtr
std::shared_ptr< ContField > ContFieldSharedPtr
std::vector< std::pair< std::string, std::string > > SummaryList
EquationSystemFactory & GetEquationSystemFactory()
void AddSummaryItem(SummaryList &l, const std::string &name, const std::string &value)
Adds a summary item to the summary info list.
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
@ eFactorSVVDGKerDiffCoeff
@ eFactorSVVPowerKerDiffCoeff
std::map< ConstFactorType, NekDouble > ConstFactorMap
static VarCoeffMap NullVarCoeffMap
std::map< StdRegions::VarCoeffType, VarCoeffEntry > VarCoeffMap
StdRegions::ConstFactorMap factors
@ eVelocityCorrectionScheme
static Array< OneD, Array< OneD, NekDouble > > NullNekDoubleArrayOfArray
ExtrapolateFactory & GetExtrapolateFactory()
static Array< OneD, NekDouble > NullNekDouble1DArray
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
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.
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.
void Vabs(int n, const T *x, const int incx, T *y, const int incy)
vabs: y = |x|
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
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.
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
void Zero(int n, T *x, const int incx)
Zero vector.
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.