40 #include <boost/algorithm/string.hpp>
46 string VCSMapping::className =
57 VCSMapping::VCSMapping(
71 "Could not create mapping in VCSMapping.");
73 std::string vExtrapolation =
"Mapping";
86 int physTot =
m_fields[0]->GetTotPoints();
88 int intMethod =
m_intScheme->GetIntegrationMethod();
116 m_session->MatchSolverInfo(
"MappingImplicitPressure",
"True",
118 m_session->MatchSolverInfo(
"MappingImplicitViscous",
"True",
120 m_session->MatchSolverInfo(
"MappingNeglectViscous",
"True",
129 m_session->LoadParameter(
"MappingPressureTolerance",
131 m_session->LoadParameter(
"MappingViscousTolerance",
133 m_session->LoadParameter(
"MappingPressureRelaxation",
135 m_session->LoadParameter(
"MappingViscousRelaxation",
153 boost::lexical_cast<std::string>(
m_kinvis);
172 int physTot =
m_fields[0]->GetTotPoints();
207 std::vector<SolverUtils::ForcingSharedPtr>::const_iterator x;
210 (*x)->Apply(
m_fields, inarray, outarray, time);
248 int physTot =
m_fields[0]->GetTotPoints();
257 for(
int i = 0; i < nvel; ++i)
261 m_fields[i]->HomogeneousBwdTrans(fields[i],wk);
270 m_fields[i]->HomogeneousFwdTrans(wk,wk);
282 bool wavespace =
m_fields[0]->GetWaveSpace();
289 for(
int i = 0; i < tmp.num_elements(); i++)
305 m_mapping->DotGradJacobian(velocity, wk);
309 for(
int i = 0; i < nvel; ++i)
313 Vmath::Vmul(physTot, Jac, 1, tmp[i], 1, tmp[i], 1);
317 for(
int i = 0; i < nvel; ++i)
327 m_mapping->DotGradJacobian(tmp, velocity[0]);
330 Vmath::Vadd(physTot, velocity[0], 1, wk, 1, wk, 1);
354 m_fields[0]->SetWaveSpace(wavespace);
368 int physTot =
m_fields[0]->GetTotPoints();
387 for (
int i=0; i<nvel; i++)
394 for (
int i=0; i<nvel; i++)
412 for (
int i=0; i<nvel; i++)
420 for (
int i=0; i<nvel; i++)
429 for(
int i = 0; i < nvel; ++i)
431 Blas::Daxpy(physTot,-aii_dtinv,inarray[i],1,
Forcing[i],1);
448 int physTot =
m_fields[0]->GetTotPoints();
450 bool converged =
false;
456 m_session->LoadParameter(
"MappingMaxIter",maxIter,5000);
466 for(
int i = 0; i < nvel; ++i)
490 "VCSMapping exceeded maximum number of iterations.");
501 for(
int i = 0; i < nvel; ++i)
504 previous_iter, gradP[i]);
507 m_pressure->HomogeneousBwdTrans(gradP[i], wk1[i]);
528 m_pressure->HomogeneousFwdTrans(F_corrected, F_corrected);
543 Vmath::Vadd(physTot, wk1[0], 1, F_corrected, 1, F_corrected, 1);
557 if ( forcing_L2 != 0)
574 std::cout <<
" Pressure system (mapping) converged in " << s <<
575 " iterations with error = " << error << std::endl;
594 int physTot =
m_fields[0]->GetTotPoints();
596 bool converged =
false;
601 m_session->LoadParameter(
"MappingMaxIter",maxIter,5000);
612 for(
int i = 0; i < nvel; ++i)
631 for(
int i = 0; i < nvel; ++i)
644 "VCSMapping exceeded maximum number of iterations.");
655 for (
int i = 0; i < nvel; ++i)
657 m_fields[0]->HomogeneousBwdTrans(previous_iter[i],
663 for (
int i = 0; i < nvel; ++i)
670 m_mapping->VelocityLaplacian(wk, F_corrected,
675 for (
int i = 0; i < nvel; ++i)
677 m_fields[0]->HomogeneousFwdTrans(F_corrected[i],
683 for (
int i = 0; i < nvel; ++i)
691 for (
int i = 0; i < nvel; ++i)
706 m_fields[i]->HelmSolve(F_corrected[i],
714 error =
m_fields[i]->L2(outarray[i], previous_iter[i]);
716 if ( forcing_L2[i] != 0)
730 if (error > max_error)
736 Vmath::Vcopy(physTot, outarray[i], 1, previous_iter[i], 1);
741 std::cout <<
" Velocity system (mapping) converged in " << s <<
742 " iterations with error = " << max_error << std::endl;
754 int physTot =
m_fields[0]->GetTotPoints();
772 m_fields[0]->HomogeneousBwdTrans(vel[i],velPhys[i]);
793 Vmath::Vadd(physTot, tmp[i], 1, Forcing[i], 1, Forcing[i], 1);
803 Vmath::Vadd(physTot, tmp[i], 1, Forcing[i], 1, Forcing[i], 1);
813 Vmath::Vadd(physTot, tmp[i], 1, Forcing[i], 1, Forcing[i], 1);
822 m_fields[0]->HomogeneousFwdTrans(Forcing[i],Forcing[i]);
829 Vmath::Vadd(physTot, outarray[i], 1, Forcing[i], 1, outarray[i], 1);
837 int physTot =
m_fields[0]->GetTotPoints();
843 m_mapping->ApplyChristoffelContravar(velPhys, wk);
846 for (
int i = 0; i< nvel; i++)
849 for (
int j = 0; j< nvel; j++)
852 outarray[i],1,outarray[i],1);
863 int physTot =
m_fields[0]->GetTotPoints();
869 for (
int i = 0; i< nvel; i++)
876 m_mapping->ContravarFromCartesian(tmp, coordVel);
879 m_mapping->ApplyChristoffelContravar(velPhys, wk);
880 for (
int i=0; i< nvel; i++)
884 m_fields[0]->PhysDeriv(velPhys[i], tmp[0], tmp[1]);
885 for (
int j=0; j< nvel; j++)
893 m_fields[0]->HomogeneousBwdTrans(tmp[2],tmp[2]);
901 outarray[i], 1, outarray[i], 1);
906 bool wavespace =
m_fields[0]->GetWaveSpace();
910 m_mapping->ApplyChristoffelContravar(coordVel, wk);
911 for (
int i=0; i< nvel; i++)
915 m_fields[0]->PhysDeriv(coordVel[i], tmp[0], tmp[1]);
919 m_fields[0]->PhysDeriv(coordVel[i], tmp[0], tmp[1], tmp[2]);
922 for (
int j=0; j< nvel; j++)
929 outarray[i], 1, outarray[i], 1);
934 m_fields[0]->SetWaveSpace(wavespace);
940 int physTot =
m_fields[0]->GetTotPoints();
952 for(
int i = 0; i < nvel; ++i)
957 for(
int i = 0; i < nvel; ++i)
969 m_mapping->VelocityLaplacian(velPhys, outarray, 1.0);
#define ASSERTL0(condition, msg)
virtual void v_InitObject()
Init object for UnsteadySystem class.
NekDouble m_pressureTolerance
virtual void v_DoInitialise(void)
Sets up initial conditions.
void ApplyIncNSMappingForcing(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
void MappingAccelerationCorrection(const Array< OneD, const Array< OneD, NekDouble > > &vel, const Array< OneD, const Array< OneD, NekDouble > > &velPhys, Array< OneD, Array< OneD, NekDouble > > &outarray)
SolverUtils::AdvectionSharedPtr m_advObject
Advection term.
tBaseSharedPtr CreateInstance(tKey idKey BOOST_PP_COMMA_IF(MAX_PARAM) BOOST_PP_ENUM_BINARY_PARAMS(MAX_PARAM, tParam, x))
Create an instance of the class referred to by idKey.
NekDouble m_kinvis
Kinematic viscosity.
ExtrapolateFactory & GetExtrapolateFactory()
NekDouble m_timestep
Time step size.
virtual void v_SolveViscous(const Array< OneD, const Array< OneD, NekDouble > > &Forcing, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble aii_Dt)
Array< OneD, int > m_velocity
int which identifies which components of m_fields contains the velocity (u,v,w);
NekDouble m_sVVDiffCoeff
Diffusion coefficient of SVV modes.
void MappingAdvectionCorrection(const Array< OneD, const Array< OneD, NekDouble > > &velPhys, Array< OneD, Array< OneD, NekDouble > > &outarray)
ExtrapolateSharedPtr m_extrapolation
GlobalMapping::MappingSharedPtr m_mapping
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
virtual void v_SetUpViscousForcing(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &Forcing, const NekDouble aii_Dt)
void EvaluateAdvectionTerms(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
std::map< ConstFactorType, NekDouble > ConstFactorMap
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.
int m_nConvectiveFields
Number of fields to be convected;.
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
NekDouble m_viscousRelaxation
Array< OneD, Array< OneD, NekDouble > > m_gradP
void MappingViscousCorrection(const Array< OneD, const Array< OneD, NekDouble > > &velPhys, Array< OneD, Array< OneD, NekDouble > > &outarray)
IMEX 2nd order scheme using Backward Different Formula & Extrapolation.
virtual void v_SolvePressure(const Array< OneD, NekDouble > &Forcing)
virtual void v_SetUpPressureForcing(const Array< OneD, const Array< OneD, NekDouble > > &fields, Array< OneD, Array< OneD, NekDouble > > &Forcing, const NekDouble aii_Dt)
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
std::vector< SolverUtils::ForcingSharedPtr > m_forcing
Forcing terms.
Base class for unsteady solvers.
bool m_SmoothAdvection
bool to identify if advection term smoothing is requested
LibUtilities::FieldMetaDataMap m_fieldMetaDataMap
Map to identify relevant solver info to dump in output fields.
void Neg(int n, T *x, const int incx)
Negate x = -x.
virtual void v_SolvePressure(const Array< OneD, NekDouble > &Forcing)
Array< OneD, Array< OneD, NekDouble > > m_presForcingCorrection
EquationSystemFactory & GetEquationSystemFactory()
Array< OneD, Array< OneD, NekDouble > > m_F
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.
MultiRegions::Direction const DirCartesianMap[]
virtual void v_SolveViscous(const Array< OneD, const Array< OneD, NekDouble > > &Forcing, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble aii_Dt)
void MappingPressureCorrection(Array< OneD, Array< OneD, NekDouble > > &outarray)
NekDouble m_sVVCutoffRatio
cutt off ratio from which to start decayhing modes
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
bool m_useSpecVanVisc
bool to identify if spectral vanishing viscosity is active.
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
IMEX 3rd order scheme using Backward Different Formula & Extrapolation.
LibUtilities::TimeIntegrationWrapperSharedPtr m_intScheme
Wrapper to the time integration scheme.
NekDouble m_viscousTolerance
virtual void v_InitObject()
Init object for UnsteadySystem class.
MultiRegions::ExpListSharedPtr m_pressure
Pointer to field holding pressure field.
IMEX 1st order scheme using Euler Backwards/Euler Forwards.
Defines a forcing term to be explicitly applied.
void Zero(int n, T *x, const int incx)
Zero vector.
virtual void v_SetUpPressureForcing(const Array< OneD, const Array< OneD, NekDouble > > &fields, Array< OneD, Array< OneD, NekDouble > > &Forcing, const NekDouble aii_Dt)
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
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.
virtual SOLVER_UTILS_EXPORT void v_DoInitialise()
Sets up initial conditions.
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 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.
virtual void v_EvaluateAdvection_SetPressureBCs(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
NekDouble m_pressureRelaxation
static FlagList NullFlagList
An empty flag list.
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.