39 #include <boost/algorithm/string.hpp>
45 string VCSMapping::className =
47 "VCSMapping", VCSMapping::create);
69 std::string vExtrapolation =
"Mapping";
77 int physTot =
m_fields[0]->GetTotPoints();
88 "Integration method not suitable: "
89 "Options include IMEXGear or IMEXOrder{1,2,3,4}");
101 m_session->MatchSolverInfo(
"MappingImplicitPressure",
"True",
103 m_session->MatchSolverInfo(
"MappingImplicitViscous",
"True",
105 m_session->MatchSolverInfo(
"MappingNeglectViscous",
"True",
153 int physTot =
m_fields[0]->GetTotPoints();
188 x->Apply(
m_fields, inarray, outarray, time);
224 int physTot =
m_fields[0]->GetTotPoints();
233 for (
int i = 0; i < nvel; ++i)
237 m_fields[i]->HomogeneousBwdTrans(fields[i], wk);
246 m_fields[i]->HomogeneousFwdTrans(wk, wk);
258 bool wavespace =
m_fields[0]->GetWaveSpace();
265 for (
int i = 0; i < tmp.size(); i++)
281 m_mapping->DotGradJacobian(velocity, wk);
285 for (
int i = 0; i < nvel; ++i)
289 Vmath::Vmul(physTot, Jac, 1, tmp[i], 1, tmp[i], 1);
293 for (
int i = 0; i < nvel; ++i)
303 m_mapping->DotGradJacobian(tmp, velocity[0]);
306 Vmath::Vadd(physTot, velocity[0], 1, wk, 1, wk, 1);
330 m_fields[0]->SetWaveSpace(wavespace);
343 int physTot =
m_fields[0]->GetTotPoints();
362 for (
int i = 0; i < nvel; i++)
369 for (
int i = 0; i < nvel; i++)
387 for (
int i = 0; i < nvel; i++)
395 for (
int i = 0; i < nvel; i++)
404 for (
int i = 0; i < nvel; ++i)
422 int physTot =
m_fields[0]->GetTotPoints();
424 bool converged =
false;
430 m_session->LoadParameter(
"MappingMaxIter", maxIter, 5000);
440 for (
int i = 0; i < nvel; ++i)
463 "VCSMapping exceeded maximum number of iterations.");
473 for (
int i = 0; i < nvel; ++i)
476 previous_iter, gradP[i]);
479 m_pressure->HomogeneousBwdTrans(gradP[i], wk1[i]);
491 Vmath::Vmul(physTot, F_corrected, 1, Jac, 1, F_corrected, 1);
498 m_pressure->HomogeneousFwdTrans(F_corrected, F_corrected);
505 Vmath::Vsub(physTot, F_corrected, 1, wk1[0], 1, F_corrected, 1);
511 Vmath::Vadd(physTot, wk1[0], 1, F_corrected, 1, F_corrected, 1);
542 std::cout <<
" Pressure system (mapping) converged in " << s
543 <<
" iterations with error = " << error << std::endl;
556 boost::ignore_unused(inarray);
565 int physTot =
m_fields[0]->GetTotPoints();
567 bool converged =
false;
572 m_session->LoadParameter(
"MappingMaxIter", maxIter, 5000);
583 for (
int i = 0; i < nvel; ++i)
602 for (
int i = 0; i < nvel; ++i)
614 "VCSMapping exceeded maximum number of iterations.");
625 for (
int i = 0; i < nvel; ++i)
627 m_fields[0]->HomogeneousBwdTrans(previous_iter[i], wk[i]);
632 for (
int i = 0; i < nvel; ++i)
639 m_mapping->VelocityLaplacian(wk, F_corrected,
644 for (
int i = 0; i < nvel; ++i)
646 m_fields[0]->HomogeneousFwdTrans(F_corrected[i],
652 for (
int i = 0; i < nvel; ++i)
654 Vmath::Vcopy(physTot, F_corrected[i], 1, F_corrected[i], 1);
659 for (
int i = 0; i < nvel; ++i)
663 1, F_corrected[i], 1);
673 m_fields[i]->HelmSolve(F_corrected[i],
674 m_fields[i]->UpdateCoeffs(), factors);
680 error =
m_fields[i]->L2(outarray[i], previous_iter[i]);
682 if (forcing_L2[i] != 0)
696 if (error > max_error)
702 Vmath::Vcopy(physTot, outarray[i], 1, previous_iter[i], 1);
707 std::cout <<
" Velocity system (mapping) converged in " << s
708 <<
" iterations with error = " << max_error << std::endl;
720 int physTot =
m_fields[0]->GetTotPoints();
738 m_fields[0]->HomogeneousBwdTrans(vel[i], velPhys[i]);
803 int physTot =
m_fields[0]->GetTotPoints();
809 m_mapping->ApplyChristoffelContravar(velPhys, wk);
812 for (
int i = 0; i < nvel; i++)
815 for (
int j = 0; j < nvel; j++)
817 Vmath::Vvtvp(physTot, wk[i * nvel + j], 1, velPhys[j], 1,
818 outarray[i], 1, outarray[i], 1);
829 int physTot =
m_fields[0]->GetTotPoints();
835 for (
int i = 0; i < nvel; i++)
842 m_mapping->ContravarFromCartesian(tmp, coordVel);
845 m_mapping->ApplyChristoffelContravar(velPhys, wk);
846 for (
int i = 0; i < nvel; i++)
850 m_fields[0]->PhysDeriv(velPhys[i], tmp[0], tmp[1]);
851 for (
int j = 0; j < nvel; j++)
859 m_fields[0]->HomogeneousBwdTrans(tmp[2], tmp[2]);
863 Vmath::Vadd(physTot, wk[i * nvel + j], 1, tmp[j], 1,
864 wk[i * nvel + j], 1);
866 Vmath::Vvtvp(physTot, coordVel[j], 1, wk[i * nvel + j], 1,
867 outarray[i], 1, outarray[i], 1);
872 bool wavespace =
m_fields[0]->GetWaveSpace();
876 m_mapping->ApplyChristoffelContravar(coordVel, wk);
877 for (
int i = 0; i < nvel; i++)
881 m_fields[0]->PhysDeriv(coordVel[i], tmp[0], tmp[1]);
885 m_fields[0]->PhysDeriv(coordVel[i], tmp[0], tmp[1], tmp[2]);
888 for (
int j = 0; j < nvel; j++)
890 Vmath::Vadd(physTot, wk[i * nvel + j], 1, tmp[j], 1,
891 wk[i * nvel + j], 1);
894 Vmath::Vvtvp(physTot, velPhys[j], 1, wk[i * nvel + j], 1,
895 outarray[i], 1, outarray[i], 1);
900 m_fields[0]->SetWaveSpace(wavespace);
906 int physTot =
m_fields[0]->GetTotPoints();
918 for (
int i = 0; i < nvel; ++i)
923 for (
int i = 0; i < nvel; ++i)
934 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);
void EvaluateAdvectionTerms(const Array< OneD, const Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray, const NekDouble time)
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.
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.
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.
LibUtilities::FieldMetaDataMap m_fieldMetaDataMap
Map to identify relevant solver info to dump in output fields.
SOLVER_UTILS_EXPORT int GetTotPoints()
Defines a forcing term to be explicitly applied.
Base class for unsteady solvers.
virtual SOLVER_UTILS_EXPORT void v_DoInitialise()
Sets up initial conditions.
LibUtilities::TimeIntegrationSchemeSharedPtr m_intScheme
Wrapper to the time integration scheme.
virtual void v_SetUpViscousForcing(const Array< OneD, const Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &Forcing, const NekDouble aii_Dt)
void MappingAccelerationCorrection(const Array< OneD, const Array< OneD, NekDouble >> &vel, const Array< OneD, const Array< OneD, NekDouble >> &velPhys, Array< OneD, Array< OneD, NekDouble >> &outarray)
void MappingViscousCorrection(const Array< OneD, const Array< OneD, NekDouble >> &velPhys, Array< OneD, Array< OneD, NekDouble >> &outarray)
virtual void v_SolvePressure(const Array< OneD, NekDouble > &Forcing)
void ApplyIncNSMappingForcing(const Array< OneD, const Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray)
NekDouble m_pressureTolerance
NekDouble m_viscousRelaxation
NekDouble m_pressureRelaxation
NekDouble m_viscousTolerance
void MappingPressureCorrection(Array< OneD, Array< OneD, NekDouble >> &outarray)
GlobalMapping::MappingSharedPtr m_mapping
void MappingAdvectionCorrection(const Array< OneD, const Array< OneD, NekDouble >> &velPhys, Array< OneD, Array< OneD, NekDouble >> &outarray)
virtual void v_InitObject(bool DeclareField=true)
Init object for UnsteadySystem class.
Array< OneD, Array< OneD, NekDouble > > m_presForcingCorrection
virtual void v_SetUpPressureForcing(const Array< OneD, const Array< OneD, NekDouble >> &fields, Array< OneD, Array< OneD, NekDouble >> &Forcing, const NekDouble aii_Dt)
virtual void v_DoInitialise(void)
Sets up initial conditions.
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_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)
Array< OneD, Array< OneD, NekDouble > > m_gradP
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_InitObject(bool DeclareField=true)
Init object for UnsteadySystem class.
bool m_useSpecVanVisc
bool to identify if spectral vanishing viscosity is active.
virtual void v_SetUpPressureForcing(const Array< OneD, const Array< OneD, NekDouble >> &fields, Array< OneD, Array< OneD, NekDouble >> &Forcing, const NekDouble aii_Dt)
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)
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
The above copyright notice and this permission notice shall be included.
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.