44 string UnsteadyAdvectionDiffusion::className
46 "UnsteadyAdvectionDiffusion",
47 UnsteadyAdvectionDiffusion::create);
49 UnsteadyAdvectionDiffusion::UnsteadyAdvectionDiffusion(
69 m_session->MatchSolverInfo(
"Extrapolation",
"SubStepping",
74 std::vector<std::string> vel;
103 m_session->LoadSolverInfo(
"AdvectionType", advName,
"WeakDG");
105 CreateInstance(advName, advName);
108 m_session->LoadSolverInfo(
"UpwindType", riemName,
"Upwind");
110 CreateInstance(riemName);
119 std::string diffName;
120 m_session->LoadSolverInfo(
"DiffusionType", diffName,
"LDG");
122 CreateInstance(diffName, diffName);
137 m_session->LoadSolverInfo(
"AdvectionType", advName,
140 CreateInstance(advName, advName);
144 if(advName.compare(
"WeakDG") == 0)
147 m_session->LoadSolverInfo(
"UpwindType", riemName,
"Upwind");
149 CreateInstance(riemName);
160 ASSERTL0(
false,
"Explicit Galerkin diffusion not set up.");
167 ASSERTL0(
false,
"Unsupported projection type.");
181 "Projection must be set to Mixed_CG_Discontinuous for "
219 for (i = 0; i < velfield.num_elements(); ++i)
221 m_fields[0]->ExtractTracePhys(velfield[i], tmp);
247 int nVariables = inarray.num_elements();
254 inarray, outarray, time);
257 for (
int i = 0; i < nVariables; ++i)
265 for (
int i = 0; i < nVariables; ++i)
272 for (
int i = 0; i < nVariables; ++i)
296 int nvariables = inarray.num_elements();
305 for(i = 0; i < nvariables; ++i)
316 for(i = 0; i < nvariables; ++i)
318 m_fields[i]->FwdTrans(inarray[i], coeffs);
319 m_fields[i]->BwdTrans_IterPerExp(coeffs, outarray[i]);
325 ASSERTL0(
false,
"Unknown projection scheme");
345 int nvariables = inarray.num_elements();
362 for (
int n = 0; n < nvariables; ++n)
371 for (
int i = 0; i < nvariables; ++i)
375 inarray[i], 1, F[i], 1);
381 for (
int i = 0; i < nvariables; ++i)
402 "Dimension of flux array and velocity array do not match");
404 const int nq =
m_fields[0]->GetNpoints();
406 for (
int i = 0; i < flux.num_elements(); ++i)
408 for (
int j = 0; j < flux[0].num_elements(); ++j)
432 for (
int k = 0; k < flux.num_elements(); ++k)
475 static int ncalls = 1;
493 cout <<
"Sub-integrating using "<< nsubsteps
498 for (
int m = 0; m < nint; ++m)
501 fields = integrationSoln->UpdateSolutionVector()[m];
510 for(n = 0; n < nsubsteps; ++n)
517 integrationSoln->SetSolVector(m,fields);
527 int n_element =
m_fields[0]->GetExpSize();
539 for(
int el = 0; el < n_element; ++el)
542 (stdVelocity[el] * cLambda *
543 (ExpOrder[el]-1) * (ExpOrder[el]-1));
573 ASSERTL0(0,
"Integration method not suitable: Options include BackwardEuler or BDFImplicitOrder1");
576 m_intSteps = IntegrationScheme->GetIntegrationSteps();
592 int nVariables = inarray.num_elements();
595 int ncoeffs =
m_fields[0]->GetNcoeffs();
600 for(i = 1; i < nVariables; ++i)
602 WeakAdv[i] = WeakAdv[i-1] + ncoeffs;
609 inarray, outarray, time);
611 for(i = 0; i < nVariables; ++i)
613 m_fields[i]->IProductWRTBase(outarray[i],WeakAdv[i]);
622 for(i = 0; i < nVariables; ++i)
628 m_fields[i]->MultiplyByElmtInvMass(WeakAdv[i], WeakAdv[i]);
631 m_fields[i]->BwdTrans(WeakAdv[i], outarray[i]);
643 ASSERTL1(inarray.num_elements() == outarray.num_elements(),
"Inarray and outarray of different sizes ");
645 for(
int i = 0; i < inarray.num_elements(); ++i)
647 Vmath::Vcopy(inarray[i].num_elements(),inarray[i],1,outarray[i],1);
656 ASSERTL1(physfield.num_elements() == Outarray.num_elements(),
657 "Physfield and outarray are of different dimensions");
662 int nTracePts =
m_fields[0]->GetTrace()->GetNpoints();
676 for(i = 0; i < physfield.num_elements(); ++i)
680 m_fields[i]->GetFwdBwdTracePhys(physfield[i], Fwd, Bwd);
683 m_fields[0]->GetTrace()->Upwind(Vn, Fwd, Bwd, numflux);
686 Vmath::Vsub(nTracePts, numflux, 1, Fwd, 1, Fwd, 1);
687 Vmath::Vsub(nTracePts, numflux, 1, Bwd, 1, Bwd, 1);
694 m_fields[0]->AddFwdBwdTraceIntegral(Fwd,Bwd,Outarray[i]);
703 int n_points_0 =
m_fields[0]->GetExp(0)->GetTotPoints();
704 int n_element =
m_fields[0]->GetExpSize();
705 int nvel = inarray.num_elements();
708 ASSERTL0(nvel >= 2,
"Method not implemented for 1D");
717 for (
int i = 0; i < nvel; ++i)
725 for (
int el = 0; el < n_element; ++el)
727 int n_points =
m_fields[0]->GetExp(el)->GetTotPoints();
728 ptsKeys =
m_fields[0]->GetExp(el)->GetPointsKeys();
731 if(n_points != n_points_0)
733 for (
int j = 0; j < nvel; ++j)
737 n_points_0 = n_points;
741 m_fields[0]->GetExp(el)->GetGeom()->GetMetricInfo()->GetDerivFactors(ptsKeys);
743 if (
m_fields[0]->GetExp(el)->GetGeom()->GetMetricInfo()->GetGtype()
746 for (
int i = 0; i < n_points; i++)
748 stdVelocity[0][i] = gmat[0][i]*inarray[0][i+cnt]
749 + gmat[2][i]*inarray[1][i+cnt];
751 stdVelocity[1][i] = gmat[1][i]*inarray[0][i+cnt]
752 + gmat[3][i]*inarray[1][i+cnt];
757 for (
int i = 0; i < n_points; i++)
759 stdVelocity[0][i] = gmat[0][0]*inarray[0][i+cnt]
760 + gmat[2][0]*inarray[1][i+cnt];
762 stdVelocity[1][i] = gmat[1][0]*inarray[0][i+cnt]
763 + gmat[3][0]*inarray[1][i+cnt];
770 for (
int i = 0; i < n_points; i++)
772 pntVelocity = stdVelocity[0][i]*stdVelocity[0][i]
773 + stdVelocity[1][i]*stdVelocity[1][i];
775 if (pntVelocity>maxV[el])
777 maxV[el] = pntVelocity;
780 maxV[el] = sqrt(maxV[el]);
786 for (
int el = 0; el < n_element; ++el)
789 int n_points =
m_fields[0]->GetExp(el)->GetTotPoints();
790 ptsKeys =
m_fields[0]->GetExp(el)->GetPointsKeys();
793 if(n_points != n_points_0)
795 for (
int j = 0; j < nvel; ++j)
799 n_points_0 = n_points;
803 m_fields[0]->GetExp(el)->GetGeom()->GetMetricInfo()->GetDerivFactors(ptsKeys);
805 if (
m_fields[0]->GetExp(el)->GetGeom()->GetMetricInfo()->GetGtype()
808 for (
int i = 0; i < n_points; i++)
810 stdVelocity[0][i] = gmat[0][i]*inarray[0][i+cnt]
811 + gmat[3][i]*inarray[1][i+cnt]
812 + gmat[6][i]*inarray[2][i+cnt];
814 stdVelocity[1][i] = gmat[1][i]*inarray[0][i+cnt]
815 + gmat[4][i]*inarray[1][i+cnt]
816 + gmat[7][i]*inarray[2][i+cnt];
818 stdVelocity[2][i] = gmat[2][i]*inarray[0][i+cnt]
819 + gmat[5][i]*inarray[1][i+cnt]
820 + gmat[8][i]*inarray[2][i+cnt];
825 for (
int i = 0; i < n_points; i++)
827 stdVelocity[0][i] = gmat[0][0]*inarray[0][i+cnt]
828 + gmat[3][0]*inarray[1][i+cnt]
829 + gmat[6][0]*inarray[2][i+cnt];
831 stdVelocity[1][i] = gmat[1][0]*inarray[0][i+cnt]
832 + gmat[4][0]*inarray[1][i+cnt]
833 + gmat[7][0]*inarray[2][i+cnt];
835 stdVelocity[2][i] = gmat[2][0]*inarray[0][i+cnt]
836 + gmat[5][0]*inarray[1][i+cnt]
837 + gmat[8][0]*inarray[2][i+cnt];
843 for (
int i = 0; i < n_points; i++)
845 pntVelocity = stdVelocity[0][i]*stdVelocity[0][i]
846 + stdVelocity[1][i]*stdVelocity[1][i]
847 + stdVelocity[2][i]*stdVelocity[2][i];
849 if (pntVelocity > maxV[el])
851 maxV[el] = pntVelocity;
855 maxV[el] = sqrt(maxV[el]);
void SubStepProjection(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
#define ASSERTL0(condition, msg)
void GetFluxVectorAdv(const Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &flux)
Evaluate the flux at each solution point for the advection part.
Array< OneD, NekDouble > GetMaxStdVelocity(const Array< OneD, Array< OneD, NekDouble > > inarray)
std::vector< PointsKey > PointsKeyVector
BDF multi-step scheme of order 1 (implicit)
bool m_homoInitialFwd
Flag to determine if simulation should start in homogeneous forward transformed state.
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.
SolverUtils::AdvectionSharedPtr m_advObject
Advection term.
NekDouble m_time
Current time of simulation.
bool m_explicitDiffusion
Indicates if explicit or implicit treatment of diffusion is used.
void DefineImplicitSolve(FuncPointerT func, ObjectPointerT obj)
LibUtilities::TimeIntegrationWrapperSharedPtr m_subStepIntegrationScheme
LibUtilities::TimeIntegrationSchemeOperators m_ode
The time integration scheme operators to use.
NekDouble m_timestep
Time step size.
std::vector< std::pair< std::string, std::string > > SummaryList
T Vmin(int n, const T *x, const int incx)
Return the minimum element in x - called vmin to avoid conflict with min.
Array< OneD, Array< OneD, NekDouble > > m_velocity
boost::shared_ptr< TimeIntegrationWrapper > TimeIntegrationWrapperSharedPtr
DiffusionFactory & GetDiffusionFactory()
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 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
enum MultiRegions::ProjectionType m_projectionType
Type of projection; e.g continuous or discontinuous.
std::map< ConstFactorType, NekDouble > ConstFactorMap
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
void SubStepAdvection(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
Array< OneD, NekDouble > & GetNormalVel(const Array< OneD, const Array< OneD, NekDouble > > &velfield)
Get the normal velocity based on input velfield.
Array< OneD, NekDouble > & GetNormalVelocity()
Get the normal velocity based on m_velocity.
SolverUtils::RiemannSolverSharedPtr m_riemannSolver
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Array holding trace normals for DG simulations in the forwards direction.
void DoOdeProjection(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
Perform the projection.
void DefineProjection(FuncPointerT func, ObjectPointerT obj)
virtual SOLVER_UTILS_EXPORT void v_GenerateSummary(SummaryList &s)
Print a summary of time stepping parameters.
boost::shared_ptr< Comm > CommSharedPtr
Pointer to a Communicator object.
virtual void DoOdeRhs(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
Compute the RHS.
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
NekDouble m_sVVCutoffRatio
virtual bool v_PreIntegrate(int step)
PreIntegration step for substepping.
void DefineOdeRhs(FuncPointerT func, ObjectPointerT obj)
Base class for unsteady solvers.
Array< OneD, NekDouble > m_traceVn
virtual void v_GenerateSummary(SolverUtils::SummaryList &s)
Print Summary.
RiemannSolverFactory & GetRiemannSolverFactory()
int m_spacedim
Spatial dimension (>= expansion dim).
void SetUpSubSteppingTimeIntegration(int intMethod, const LibUtilities::TimeIntegrationWrapperSharedPtr &IntegrationScheme)
virtual ~UnsteadyAdvectionDiffusion()
Destructor.
SolverUtils::DiffusionSharedPtr m_diffusion
AdvectionFactory & GetAdvectionFactory()
Gets the factory for initialising advection objects.
void Neg(int n, T *x, const int incx)
Negate x = -x.
TimeIntegrationWrapperFactory & GetTimeIntegrationWrapperFactory()
SOLVER_UTILS_EXPORT void EvaluateFunction(Array< OneD, Array< OneD, NekDouble > > &pArray, std::string pFunctionName, const NekDouble pTime=0.0, const int domain=0)
Evaluates a function as specified in the session file.
EquationSystemFactory & GetEquationSystemFactory()
void SubStepAdvance(const LibUtilities::TimeIntegrationSolutionSharedPtr &integrationSoln, int nstep, NekDouble time)
SOLVER_UTILS_EXPORT void SetBoundaryConditions(NekDouble time)
Evaluates the boundary conditions at the given time.
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.
BDF multi-step scheme of order 2 (implicit)
virtual void DoImplicitSolve(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, NekDouble time, NekDouble lambda)
Solve implicitly the diffusion term.
SOLVER_UTILS_EXPORT int GetNpoints()
LibUtilities::TimeIntegrationSchemeOperators m_subStepIntegrationOps
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
SOLVER_UTILS_EXPORT int GetExpSize()
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
SOLVER_UTILS_EXPORT int GetTraceNpoints()
void AddAdvectionPenaltyFlux(const Array< OneD, const Array< OneD, NekDouble > > &velfield, const Array< OneD, const Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &Outarray)
SOLVER_UTILS_EXPORT int GetNcoeffs()
boost::shared_ptr< TimeIntegrationSolution > TimeIntegrationSolutionSharedPtr
LibUtilities::TimeIntegrationWrapperSharedPtr m_intScheme
Wrapper to the time integration scheme.
virtual SOLVER_UTILS_EXPORT void v_InitObject()
Init object for UnsteadySystem class.
void GetFluxVectorDiff(const int i, const int j, const Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &derivatives, Array< OneD, Array< OneD, NekDouble > > &flux)
Evaluate the flux at each solution point for the diffusion part.
void Zero(int n, T *x, const int incx)
Zero vector.
A base class for PDEs which include an advection component.
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Geometry is curved or has non-constant factors.
NekDouble m_cflSafetyFactor
virtual void v_InitObject()
Initialise the object.
LibUtilities::TimeIntegrationSolutionSharedPtr m_intSoln
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.
static FlagList NullFlagList
An empty flag list.
NekDouble GetSubstepTimeStep()
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.