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);
170 int physTot =
m_fields[0]->GetTotPoints();
217 std::vector<SolverUtils::ForcingSharedPtr>::const_iterator x;
220 (*x)->Apply(
m_fields, inarray, outarray, time);
246 int physTot =
m_fields[0]->GetTotPoints();
255 for(
int i = 0; i < nvel; ++i)
259 m_fields[i]->HomogeneousBwdTrans(fields[i],wk);
268 m_fields[i]->HomogeneousFwdTrans(wk,wk);
280 bool wavespace =
m_fields[0]->GetWaveSpace();
287 for(
int i = 0; i < tmp.num_elements(); i++)
303 m_mapping->DotGradJacobian(velocity, wk);
307 for(
int i = 0; i < nvel; ++i)
311 Vmath::Vmul(physTot, Jac, 1, tmp[i], 1, tmp[i], 1);
315 for(
int i = 0; i < nvel; ++i)
325 m_mapping->DotGradJacobian(tmp, velocity[0]);
328 Vmath::Vadd(physTot, velocity[0], 1, wk, 1, wk, 1);
340 m_fields[0]->HomogeneousFwdTrans(wk,wk);
345 m_fields[0]->SetWaveSpace(wavespace);
359 int physTot =
m_fields[0]->GetTotPoints();
378 for (
int i=0; i<nvel; i++)
385 for (
int i=0; i<nvel; i++)
403 for (
int i=0; i<nvel; i++)
411 for (
int i=0; i<nvel; i++)
420 for(
int i = 0; i < nvel; ++i)
422 Blas::Daxpy(physTot,-aii_dtinv,inarray[i],1,
Forcing[i],1);
439 int physTot =
m_fields[0]->GetTotPoints();
441 bool converged =
false;
447 m_session->LoadParameter(
"MappingMaxIter",maxIter,5000);
457 for(
int i = 0; i < nvel; ++i)
481 "VCSMapping exceeded maximum number of iterations.");
492 for(
int i = 0; i < nvel; ++i)
495 previous_iter, gradP[i]);
498 m_pressure->HomogeneousBwdTrans(gradP[i], wk1[i]);
519 m_pressure->HomogeneousFwdTrans(F_corrected, F_corrected);
534 Vmath::Vadd(physTot, wk1[0], 1, F_corrected, 1, F_corrected, 1);
548 if ( forcing_L2 != 0)
565 std::cout <<
" Pressure system (mapping) converged in " << s <<
566 " iterations with error = " << error << std::endl;
585 int physTot =
m_fields[0]->GetTotPoints();
587 bool converged =
false;
592 m_session->LoadParameter(
"MappingMaxIter",maxIter,5000);
603 for(
int i = 0; i < nvel; ++i)
622 for(
int i = 0; i < nvel; ++i)
635 "VCSMapping exceeded maximum number of iterations.");
646 for (
int i = 0; i < nvel; ++i)
648 m_fields[0]->HomogeneousBwdTrans(previous_iter[i],
654 for (
int i = 0; i < nvel; ++i)
661 m_mapping->VelocityLaplacian(wk, F_corrected,
666 for (
int i = 0; i < nvel; ++i)
668 m_fields[0]->HomogeneousFwdTrans(F_corrected[i],
674 for (
int i = 0; i < nvel; ++i)
682 for (
int i = 0; i < nvel; ++i)
697 m_fields[i]->HelmSolve(F_corrected[i],
705 error =
m_fields[i]->L2(outarray[i], previous_iter[i]);
707 if ( forcing_L2[i] != 0)
721 if (error > max_error)
727 Vmath::Vcopy(physTot, outarray[i], 1, previous_iter[i], 1);
732 std::cout <<
" Velocity system (mapping) converged in " << s <<
733 " iterations with error = " << max_error << std::endl;
744 int physTot =
m_fields[0]->GetTotPoints();
762 m_fields[0]->HomogeneousBwdTrans(vel[i],velPhys[i]);
783 Vmath::Vadd(physTot, tmp[i], 1, Forcing[i], 1, Forcing[i], 1);
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);
812 m_fields[0]->HomogeneousFwdTrans(Forcing[i],Forcing[i]);
819 Vmath::Vadd(physTot, outarray[i], 1, Forcing[i], 1, outarray[i], 1);
827 int physTot =
m_fields[0]->GetTotPoints();
833 m_mapping->ApplyChristoffelContravar(velPhys, wk);
836 for (
int i = 0; i< nvel; i++)
839 for (
int j = 0; j< nvel; j++)
842 outarray[i],1,outarray[i],1);
853 int physTot =
m_fields[0]->GetTotPoints();
859 for (
int i = 0; i< nvel; i++)
866 m_mapping->ContravarFromCartesian(tmp, coordVel);
869 m_mapping->ApplyChristoffelContravar(velPhys, wk);
870 for (
int i=0; i< nvel; i++)
874 m_fields[0]->PhysDeriv(velPhys[i], tmp[0], tmp[1]);
875 for (
int j=0; j< nvel; j++)
883 m_fields[0]->HomogeneousBwdTrans(tmp[2],tmp[2]);
891 outarray[i], 1, outarray[i], 1);
896 bool wavespace =
m_fields[0]->GetWaveSpace();
900 m_mapping->ApplyChristoffelContravar(coordVel, wk);
901 for (
int i=0; i< nvel; i++)
905 m_fields[0]->PhysDeriv(coordVel[i], tmp[0], tmp[1]);
909 m_fields[0]->PhysDeriv(coordVel[i], tmp[0], tmp[1], tmp[2]);
912 for (
int j=0; j< nvel; j++)
919 outarray[i], 1, outarray[i], 1);
924 m_fields[0]->SetWaveSpace(wavespace);
930 int physTot =
m_fields[0]->GetTotPoints();
942 for(
int i = 0; i < nvel; ++i)
947 for(
int i = 0; i < nvel; ++i)
959 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(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)
void MappingAccelerationCorrection(const Array< OneD, Array< OneD, NekDouble > > &vel, const Array< OneD, Array< OneD, NekDouble > > &velPhys, Array< OneD, Array< OneD, NekDouble > > &outarray)
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.
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)
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
void MappingViscousCorrection(const Array< OneD, Array< OneD, NekDouble > > &velPhys, Array< OneD, Array< OneD, NekDouble > > &outarray)
Array< OneD, Array< OneD, NekDouble > > m_gradP
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
void MappingAdvectionCorrection(const Array< OneD, Array< OneD, NekDouble > > &velPhys, Array< OneD, Array< OneD, NekDouble > > &outarray)
EquationSystemFactory & GetEquationSystemFactory()
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)
void EvaluateAdvectionTerms(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, Array< OneD, NekDouble > &wk=NullNekDouble1DArray)
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.