39#include <boost/algorithm/string.hpp>
72 std::string vExtrapolation =
"Mapping";
81 size_t physTot =
m_fields[0]->GetTotPoints();
92 "Integration method not suitable: "
93 "Options include IMEXGear or IMEXOrder{1,2,3,4}");
105 m_session->MatchSolverInfo(
"MappingImplicitPressure",
"True",
107 m_session->MatchSolverInfo(
"MappingImplicitViscous",
"True",
109 m_session->MatchSolverInfo(
"MappingNeglectViscous",
"True",
150 size_t physTot =
m_fields[0]->GetTotPoints();
185 x->Apply(
m_fields, inarray, outarray, time);
221 size_t physTot =
m_fields[0]->GetTotPoints();
230 for (
size_t i = 0; i < nvel; ++i)
234 m_fields[i]->HomogeneousBwdTrans(physTot, fields[i], wk);
243 m_fields[i]->HomogeneousFwdTrans(physTot, wk, wk);
255 bool wavespace =
m_fields[0]->GetWaveSpace();
262 for (
size_t i = 0; i < tmp.size(); i++)
269 physTot,
m_fields[i]->GetPhys(), velocity[i]);
278 m_mapping->DotGradJacobian(velocity, wk);
282 for (
size_t i = 0; i < nvel; ++i)
286 Vmath::Vmul(physTot, Jac, 1, tmp[i], 1, tmp[i], 1);
290 for (
size_t i = 0; i < nvel; ++i)
300 m_mapping->DotGradJacobian(tmp, velocity[0]);
303 Vmath::Vadd(physTot, velocity[0], 1, wk, 1, wk, 1);
327 m_fields[0]->SetWaveSpace(wavespace);
340 size_t physTot =
m_fields[0]->GetTotPoints();
359 for (
size_t i = 0; i < nvel; i++)
366 for (
size_t i = 0; i < nvel; i++)
384 for (
size_t i = 0; i < nvel; i++)
392 for (
size_t i = 0; i < nvel; i++)
402 for (
size_t i = 0; i < nvel; ++i)
420 size_t physTot =
m_fields[0]->GetTotPoints();
422 bool converged =
false;
428 m_session->LoadParameter(
"MappingMaxIter", maxIter, 5000);
438 for (
size_t i = 0; i < nvel; ++i)
461 "VCSMapping exceeded maximum number of iterations.");
471 for (
size_t i = 0; i < nvel; ++i)
474 previous_iter, gradP[i]);
477 m_pressure->HomogeneousBwdTrans(physTot, gradP[i], wk1[i]);
489 Vmath::Vmul(physTot, F_corrected, 1, Jac, 1, F_corrected, 1);
496 m_pressure->HomogeneousFwdTrans(physTot, F_corrected,
504 Vmath::Vsub(physTot, F_corrected, 1, wk1[0], 1, F_corrected, 1);
510 Vmath::Vadd(physTot, wk1[0], 1, F_corrected, 1, F_corrected, 1);
541 std::cout <<
" Pressure system (mapping) converged in " << s
542 <<
" iterations with error = " << error << std::endl;
562 size_t physTot =
m_fields[0]->GetTotPoints();
564 bool converged =
false;
569 m_session->LoadParameter(
"MappingMaxIter", maxIter, 5000);
580 for (
size_t i = 0; i < nvel; ++i)
599 for (
size_t i = 0; i < nvel; ++i)
611 "VCSMapping exceeded maximum number of iterations.");
622 for (
size_t i = 0; i < nvel; ++i)
624 m_fields[0]->HomogeneousBwdTrans(physTot, previous_iter[i],
630 for (
size_t i = 0; i < nvel; ++i)
637 m_mapping->VelocityLaplacian(wk, F_corrected,
642 for (
size_t i = 0; i < nvel; ++i)
644 m_fields[0]->HomogeneousFwdTrans(physTot, F_corrected[i],
650 for (
size_t i = 0; i < nvel; ++i)
652 Vmath::Vcopy(physTot, F_corrected[i], 1, F_corrected[i], 1);
657 for (
size_t i = 0; i < nvel; ++i)
661 1, F_corrected[i], 1);
671 m_fields[i]->HelmSolve(F_corrected[i],
678 error =
m_fields[i]->L2(outarray[i], previous_iter[i]);
680 if (forcing_L2[i] != 0)
694 if (error > max_error)
700 Vmath::Vcopy(physTot, outarray[i], 1, previous_iter[i], 1);
705 std::cout <<
" Velocity system (mapping) converged in " << s
706 <<
" iterations with error = " << max_error << std::endl;
718 size_t physTot =
m_fields[0]->GetTotPoints();
736 m_fields[0]->HomogeneousBwdTrans(physTot, vel[i], velPhys[i]);
801 size_t physTot =
m_fields[0]->GetTotPoints();
807 m_mapping->ApplyChristoffelContravar(velPhys, wk);
810 for (
size_t i = 0; i < nvel; i++)
813 for (
size_t j = 0; j < nvel; j++)
815 Vmath::Vvtvp(physTot, wk[i * nvel + j], 1, velPhys[j], 1,
816 outarray[i], 1, outarray[i], 1);
827 size_t physTot =
m_fields[0]->GetTotPoints();
833 for (
size_t i = 0; i < nvel; i++)
840 m_mapping->ContravarFromCartesian(tmp, coordVel);
843 m_mapping->ApplyChristoffelContravar(velPhys, wk);
844 for (
size_t i = 0; i < nvel; i++)
848 m_fields[0]->PhysDeriv(velPhys[i], tmp[0], tmp[1]);
849 for (
size_t j = 0; j < nvel; j++)
857 m_fields[0]->HomogeneousBwdTrans(physTot, tmp[2], tmp[2]);
861 Vmath::Vadd(physTot, wk[i * nvel + j], 1, tmp[j], 1,
862 wk[i * nvel + j], 1);
864 Vmath::Vvtvp(physTot, coordVel[j], 1, wk[i * nvel + j], 1,
865 outarray[i], 1, outarray[i], 1);
870 bool wavespace =
m_fields[0]->GetWaveSpace();
874 m_mapping->ApplyChristoffelContravar(coordVel, wk);
875 for (
size_t i = 0; i < nvel; i++)
879 m_fields[0]->PhysDeriv(coordVel[i], tmp[0], tmp[1]);
883 m_fields[0]->PhysDeriv(coordVel[i], tmp[0], tmp[1], tmp[2]);
886 for (
size_t j = 0; j < nvel; j++)
888 Vmath::Vadd(physTot, wk[i * nvel + j], 1, tmp[j], 1,
889 wk[i * nvel + j], 1);
892 Vmath::Vvtvp(physTot, velPhys[j], 1, wk[i * nvel + j], 1,
893 outarray[i], 1, outarray[i], 1);
898 m_fields[0]->SetWaveSpace(wavespace);
904 size_t physTot =
m_fields[0]->GetTotPoints();
916 for (
size_t i = 0; i < nvel; ++i)
921 for (
size_t i = 0; i < nvel; ++i)
932 m_mapping->VelocityLaplacian(velPhys, outarray, 1.0);
#define ASSERTL0(condition, msg)
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mode...
static GLOBAL_MAPPING_EXPORT MappingSharedPtr Load(const LibUtilities::SessionReaderSharedPtr &pSession, const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields)
Return a pointer to the mapping, creating it on first call.
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
ExtrapolateSharedPtr m_extrapolation
Array< OneD, int > m_velocity
int which identifies which components of m_fields contains the velocity (u,v,w);
int m_intSteps
Number of time integration steps AND Order of extrapolation for pressure boundary conditions.
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.
SolverUtils::AdvectionSharedPtr m_advObject
Advection term.
NekDouble m_timestep
Time step size.
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
SOLVER_UTILS_EXPORT int GetTotPoints()
LibUtilities::FieldMetaDataMap m_fieldMetaDataMap
Map to identify relevant solver info to dump in output fields.
Defines a forcing term to be explicitly applied.
Base class for unsteady solvers.
SOLVER_UTILS_EXPORT void v_DoInitialise(bool dumpInitialConditions=true) override
Sets up initial conditions.
LibUtilities::TimeIntegrationSchemeSharedPtr m_intScheme
Wrapper to the time integration scheme.
void v_SolvePressure(const Array< OneD, NekDouble > &Forcing) override
void v_SetUpPressureForcing(const Array< OneD, const Array< OneD, NekDouble > > &fields, Array< OneD, Array< OneD, NekDouble > > &Forcing, const NekDouble aii_Dt) override
void MappingViscousCorrection(const Array< OneD, const Array< OneD, NekDouble > > &velPhys, Array< OneD, Array< OneD, NekDouble > > &outarray)
void v_DoInitialise(bool dumpInitialConditions=true) override
Sets up initial conditions.
void v_SetUpViscousForcing(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &Forcing, const NekDouble aii_Dt) override
static std::string className
Name of class.
NekDouble m_pressureTolerance
NekDouble m_viscousRelaxation
NekDouble m_pressureRelaxation
static SolverUtils::EquationSystemSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
Creates an instance of this class.
NekDouble m_viscousTolerance
void v_EvaluateAdvection_SetPressureBCs(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time) override
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
GlobalMapping::MappingSharedPtr m_mapping
void MappingPressureCorrection(Array< OneD, Array< OneD, NekDouble > > &outarray)
void ApplyIncNSMappingForcing(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
Array< OneD, Array< OneD, NekDouble > > m_presForcingCorrection
VCSMapping(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
static std::string solverTypeLookupId
void v_InitObject(bool DeclareField=true) override
Initialisation object for EquationSystem.
void MappingAccelerationCorrection(const Array< OneD, const Array< OneD, NekDouble > > &vel, const Array< OneD, const Array< OneD, NekDouble > > &velPhys, Array< OneD, Array< OneD, NekDouble > > &outarray)
Array< OneD, Array< OneD, NekDouble > > m_gradP
void MappingAdvectionCorrection(const Array< OneD, const Array< OneD, NekDouble > > &velPhys, Array< OneD, Array< OneD, NekDouble > > &outarray)
NekDouble m_sVVCutoffRatio
cutt off ratio from which to start decayhing modes
NekDouble m_sVVDiffCoeff
Diffusion coefficient of SVV modes.
virtual void v_SolvePressure(const Array< OneD, NekDouble > &Forcing)
Array< OneD, Array< OneD, NekDouble > > m_F
virtual void v_SetUpPressureForcing(const Array< OneD, const Array< OneD, NekDouble > > &fields, Array< OneD, Array< OneD, NekDouble > > &Forcing, const NekDouble aii_Dt)
bool m_useSpecVanVisc
bool to identify if spectral vanishing viscosity is active.
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 v_InitObject(bool DeclareField=true) override
Initialisation object for EquationSystem.
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
MultiRegions::Direction const DirCartesianMap[]
EquationSystemFactory & GetEquationSystemFactory()
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
std::map< ConstFactorType, NekDouble > ConstFactorMap
StdRegions::ConstFactorMap factors
ExtrapolateFactory & GetExtrapolateFactory()
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 Neg(int n, T *x, const int incx)
Negate x = -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 Vdiv(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 Zero(int n, T *x, const int incx)
Zero vector.
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
void Vsub(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Subtract vector z = x-y.