43 #include <boost/algorithm/string.hpp>
49 using namespace MultiRegions;
51 string VelocityCorrectionScheme::className =
53 "VelocityCorrectionScheme", VelocityCorrectionScheme::create);
61 VelocityCorrectionScheme::VelocityCorrectionScheme(
84 ASSERTL0(
false,
"Need to set up pressure field definition");
91 std::string varName =
m_session->GetVariable(n);
92 if (
m_session->DefinesFunction(
"DiffusionCoefficient", varName))
95 m_session->GetFunction(
"DiffusionCoefficient", varName);
110 m_session->MatchSolverInfo(
"GJPStabilisation",
"False",
124 "Can only specify GJPNormalVelocity with"
125 " GJPStabilisation set to Explicit currently");
153 if (
m_session->DefinesSolverInfo(
"Extrapolation"))
156 m_session->GetSolverInfo(
"Extrapolation"));
224 std::string forces[] = {
"X",
"Y",
"Z"};
231 std::string varName = std::string(
"Force") + forces[i];
232 defined =
m_session->DefinesFunction(
"FlowrateForce", varName);
240 "A 'FlowrateForce' function must defined with components "
241 "[ForceX, ...] to define direction of flowrate forcing");
244 m_session->GetFunction(
"FlowrateForce", varName);
245 flowrateForce[i] = ffunc->Evaluate();
260 flowrateForce[2] = 1.0;
264 for (
int i = 0; i < bcs.size(); ++i)
266 if (boost::iequals(bcs[i]->GetUserDefined(),
"Flowrate"))
276 "One boundary region must be marked using the 'Flowrate' "
277 "user-defined type to monitor the volumetric flowrate.");
292 for (
int i = 0; i < zIDs.size(); ++i)
301 ASSERTL1(tmpId <= 0,
"Should be either at location 0 or -1 if not "
333 for (
int i = 0; i < zIDs.size(); ++i)
356 int nqTot =
m_fields[0]->GetNpoints();
368 m_fields[i]->HomogeneousFwdTrans(inTmp[i], inTmp2);
384 for (
int n = 0; n < PBndConds.size(); ++n)
386 if (PBndConds[n]->GetBoundaryConditionType() ==
413 std::string filename =
m_session->GetSessionName();
418 <<
"# -------------------------------------------"
502 << setw(16) <<
m_alpha << endl;
523 "Velocity correction (strong press. form)");
534 dealias += (dealias ==
"" ?
"" :
" + ") +
string(
"spectral/hp");
547 s,
"Smoothing-SpecHP",
548 "SVV (" + smoothing +
" Exp Kernel(cut-off = " +
558 s,
"Smoothing-SpecHP",
559 "SVV (" + smoothing +
" Power Kernel (Power ratio =" +
568 s,
"Smoothing-SpecHP",
569 "SVV (" + smoothing +
" DG Kernel (diff coeff = " +
579 s,
"Smoothing-Homo1D",
580 "SVV (Homogeneous1D - Exp Kernel(cut-off = " +
589 s,
"GJP Stab. Impl. ",
590 m_session->GetSolverInfo(
"GJPStabilisation"));
593 if (boost::iequals(
m_session->GetSolverInfo(
"GJPStabilisation"),
597 s,
"GJP Normal Velocity",
598 m_session->GetSolverInfo(
"GJPNormalVelocity"));
630 for (
int i = 0; i <
m_fields.size(); ++i)
646 for (
int k = 0; k < nfields; ++k)
661 for (
int k = 0; k < nfields; ++k)
674 int vVar =
m_session->GetVariables().size();
676 vChecks[vVar - 1] =
true;
685 return m_session->GetVariables().size() - 1;
713 x->Apply(
m_fields, inarray, outarray, time);
738 int physTot =
m_fields[0]->GetTotPoints();
779 m_fields[i]->FwdTransLocalElmt(outarray[i],
797 int physTot =
m_fields[0]->GetTotPoints();
802 for (i = 1; i < nvel; ++i)
820 int phystot =
m_fields[0]->GetTotPoints();
887 std::dynamic_pointer_cast<MultiRegions::ContField>(
m_fields[0]);
890 cfield->GetGJPForcing();
892 int nTracePts = GJPData->GetNumTracePts();
896 GJPData->GetTraceNormals();
898 m_fields[0]->GetFwdBwdTracePhys(inarray[0], Fwd, Bwd,
true,
true);
899 Vmath::Vmul(nTracePts, Fwd, 1, traceNormals[0], 1, unorm, 1);
902 for (
int f = 1; f <
m_fields[0]->GetCoordim(0); ++f)
904 m_fields[0]->GetFwdBwdTracePhys(inarray[f], Fwd, Bwd,
true,
true);
905 Vmath::Vvtvp(nTracePts, Fwd, 1, traceNormals[f], 1, unorm, 1, unorm,
924 varCoeffMap, varFactorsMap);
932 m_session->MatchSolverInfo(
"SpectralVanishingViscosity",
"PowerKernel",
941 m_session->MatchSolverInfo(
"SpectralVanishingViscositySpectralHP",
951 m_session->MatchSolverInfo(
"SpectralVanishingViscosity",
"DGKernel",
959 m_session->MatchSolverInfo(
"SpectralVanishingViscositySpectralHP",
974 if (
m_session->DefinesFunction(
"SVVVelocityMagnitude"))
976 if (
m_comm->GetRank() == 0)
978 cout <<
"Seting up SVV velocity from "
979 "SVVVelocityMagnitude section in session file"
983 int phystot =
m_fields[0]->GetTotPoints();
986 for (
int i = 0; i < nvel; ++i)
993 GetFunction(
"SVVVelocityMagnitude")->Evaluate(vars, SVVVelFields);
1009 m_session->MatchSolverInfo(
"SpectralVanishingViscosity",
"True",
1013 m_session->MatchSolverInfo(
"SpectralVanishingViscosity",
1020 m_session->MatchSolverInfo(
"SpectralVanishingViscositySpectralHP",
1025 "SpectralVanishingViscositySpectralHP",
"ExpKernel",
1032 if (
m_session->DefinesSolverInfo(
"SpectralVanishingViscosityHomo1D"))
1034 m_session->MatchSolverInfo(
"SpectralVanishingViscosityHomo1D",
"True",
1038 m_session->MatchSolverInfo(
"SpectralVanishingViscosityHomo1D",
1054 "Expect to have three velocity fields with homogenous expansion");
1061 int num_planes = planes.size();
1064 int kmodes =
m_fields[0]->GetHomogeneousBasis()->GetNumModes();
1069 for (
int n = 0; n < num_planes; ++n)
1071 if (planes[n] > pstart)
1073 fac = (
NekDouble)((planes[n] - kmodes) *
1074 (planes[n] - kmodes)) /
1076 (planes[n] - pstart)));
1093 int phystot =
m_fields[0]->GetTotPoints();
1094 int nel =
m_fields[0]->GetNumElmts();
1106 Vmath::Vmul(phystot, vel[0], 1, vel[0], 1, Velmag, 1);
1107 for (
int n = 1; n < nvel; ++n)
1109 Vmath::Vvtvp(phystot, vel[n], 1, vel[n], 1, Velmag, 1, Velmag, 1);
1116 for (
int i = 0; i < nel; ++i)
1118 int nq =
m_fields[0]->GetExp(i)->GetTotPoints();
1120 diffcoeff[i] =
m_fields[0]->GetExp(i)->Integral(tmp);
1123 diffcoeff[i] = diffcoeff[i] / area;
1132 for (
int e = 0; e < nel; e++)
1138 for (
int i = 0; i < exp->GetGeom()->GetNumEdges(); ++i)
1140 h = max(h, exp->GetGeom()->GetEdge(i)->GetVertex(0)->dist(
1141 *(exp->GetGeom()->GetEdge(i)->GetVertex(1))));
1147 p = max(
p, exp->GetBasisNumModes(i) - 1);
1150 diffcoeff[e] *= h /
p;
#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
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.
tBaseSharedPtr CreateInstance(tKey idKey, tParam... args)
Create an instance of the class referred to by idKey.
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.
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.
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.
std::vector< int > m_intVariables
NekDouble m_sVVCutoffRatioHomo1D
NekDouble m_greenFlux
Flux of the Stokes function solution.
virtual std::string v_GetExtrapolateStr(void)
virtual ~VelocityCorrectionScheme()
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)
virtual bool v_PostIntegrate(int step)
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.
void SetUpExtrapolation(void)
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.
virtual int v_GetForceDimension()
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.
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< StdRegions::VarCoeffType, Array< OneD, NekDouble > > VarCoeffMap
std::map< ConstFactorType, NekDouble > ConstFactorMap
static VarCoeffMap NullVarCoeffMap
The above copyright notice and this permission notice shall be included.
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.