40using namespace LibUtilities;
45 "Unsteady Advection-Diffusion equation");
64 m_session->MatchSolverInfo(
"SpectralVanishingViscosity",
"True",
73 m_session->MatchSolverInfo(
"Extrapolation",
"SubStepping",
88 m_session->LoadSolverInfo(
"DiffusionType", diffName,
"LDG");
97 "SubSteppingScheme is not set up for DG projection");
106 m_session->LoadSolverInfo(
"AdvectionType", advName,
108 if (advName.compare(
"WeakDG") == 0)
116 std::string riemName;
117 m_session->LoadSolverInfo(
"UpwindType", riemName,
"Upwind");
130 ASSERTL0(
false,
"Explicit Galerkin diffusion not set up.");
137 ASSERTL0(
false,
"Unsupported projection type.");
149 "Projection must be set to Mixed_CG_Discontinuous for "
168 int nVariables = inarray.size();
178 for (
int i = 0; i < nVariables; ++i)
185 for (
int i = 0; i < nVariables; ++i)
187 Vmath::Vadd(nSolutionPts, &outarrayDiff[i][0], 1, &outarray[i][0],
188 1, &outarray[i][0], 1);
205 int nvariables = inarray.size();
229 for (
int n = 0; n < nvariables; ++n)
238 for (
int i = 0; i < nvariables; ++i)
248 for (
int i = 0; i < nvariables; ++i)
266 unsigned int nDim = qfield.size();
267 unsigned int nConvectiveFields = qfield[0].size();
268 unsigned int nPts = qfield[0][0].size();
269 for (
unsigned int j = 0; j < nDim; ++j)
271 for (
unsigned int i = 0; i < nConvectiveFields; ++i)
284 std::stringstream ss;
315 static int ncalls = 1;
332 std::cout <<
"Sub-integrating using " << nsubsteps
339 for (
int m = 0; m < nint; ++m)
342 fields = solutionVector[m];
350 for (
int n = 0; n < nsubsteps; ++n)
365 int n_element =
m_fields[0]->GetExpSize();
377 for (
int el = 0; el < n_element; ++el)
381 (ExpOrder[el] - 1) * (ExpOrder[el] - 1));
395 unsigned int order = IntegrationScheme->GetOrder();
399 if ((IntegrationScheme->GetName() ==
"Euler" &&
400 IntegrationScheme->GetVariant() ==
"Backward") ||
401 (IntegrationScheme->GetName() ==
"BDFImplicit" &&
402 (order == 1 || order == 2)))
407 "RungeKutta",
"SSP", order, std::vector<NekDouble>());
412 "Integration method not suitable: "
413 "Options include BackwardEuler or BDFImplicitOrder1");
416 m_intSteps = IntegrationScheme->GetNumIntegrationPhases();
433 int nVariables = inarray.size();
436 int ncoeffs =
m_fields[0]->GetNcoeffs();
441 for (i = 1; i < nVariables; ++i)
443 WeakAdv[i] = WeakAdv[i - 1] + ncoeffs;
452 for (i = 0; i < nVariables; ++i)
454 m_fields[i]->IProductWRTBase(outarray[i], WeakAdv[i]);
462 for (i = 0; i < nVariables; ++i)
468 m_fields[i]->MultiplyByElmtInvMass(WeakAdv[i], WeakAdv[i]);
471 m_fields[i]->BwdTrans(WeakAdv[i], outarray[i]);
483 ASSERTL1(inarray.size() == outarray.size(),
484 "Inarray and outarray of different sizes ");
486 for (
int i = 0; i < inarray.size(); ++i)
488 Vmath::Vcopy(inarray[i].size(), inarray[i], 1, outarray[i], 1);
497 ASSERTL1(physfield.size() == Outarray.size(),
498 "Physfield and outarray are of different dimensions");
503 int nTracePts =
m_fields[0]->GetTrace()->GetNpoints();
517 for (i = 0; i < physfield.size(); ++i)
521 m_fields[i]->GetFwdBwdTracePhys(physfield[i], Fwd, Bwd);
524 m_fields[0]->GetTrace()->Upwind(Vn, Fwd, Bwd, numflux);
527 Vmath::Vsub(nTracePts, numflux, 1, Fwd, 1, Fwd, 1);
528 Vmath::Vsub(nTracePts, numflux, 1, Bwd, 1, Bwd, 1);
535 m_fields[0]->AddFwdBwdTraceIntegral(Fwd, Bwd, Outarray[i]);
543 int n_points_0 =
m_fields[0]->GetExp(0)->GetTotPoints();
544 int n_element =
m_fields[0]->GetExpSize();
545 int nvel = inarray.size();
548 ASSERTL0(nvel >= 2,
"Method not implemented for 1D");
557 for (
int i = 0; i < nvel; ++i)
565 for (
int el = 0; el < n_element; ++el)
567 int n_points =
m_fields[0]->GetExp(el)->GetTotPoints();
568 ptsKeys =
m_fields[0]->GetExp(el)->GetPointsKeys();
571 if (n_points != n_points_0)
573 for (
int j = 0; j < nvel; ++j)
577 n_points_0 = n_points;
584 ->GetDerivFactors(ptsKeys);
592 for (
int i = 0; i < n_points; i++)
594 stdVelocity[0][i] = gmat[0][i] * inarray[0][i + cnt] +
595 gmat[2][i] * inarray[1][i + cnt];
597 stdVelocity[1][i] = gmat[1][i] * inarray[0][i + cnt] +
598 gmat[3][i] * inarray[1][i + cnt];
603 for (
int i = 0; i < n_points; i++)
605 stdVelocity[0][i] = gmat[0][0] * inarray[0][i + cnt] +
606 gmat[2][0] * inarray[1][i + cnt];
608 stdVelocity[1][i] = gmat[1][0] * inarray[0][i + cnt] +
609 gmat[3][0] * inarray[1][i + cnt];
615 for (
int i = 0; i < n_points; i++)
617 pntVelocity = stdVelocity[0][i] * stdVelocity[0][i] +
618 stdVelocity[1][i] * stdVelocity[1][i];
620 if (pntVelocity > maxV[el])
622 maxV[el] = pntVelocity;
625 maxV[el] =
sqrt(maxV[el]);
631 for (
int el = 0; el < n_element; ++el)
634 int n_points =
m_fields[0]->GetExp(el)->GetTotPoints();
635 ptsKeys =
m_fields[0]->GetExp(el)->GetPointsKeys();
638 if (n_points != n_points_0)
640 for (
int j = 0; j < nvel; ++j)
644 n_points_0 = n_points;
651 ->GetDerivFactors(ptsKeys);
659 for (
int i = 0; i < n_points; i++)
661 stdVelocity[0][i] = gmat[0][i] * inarray[0][i + cnt] +
662 gmat[3][i] * inarray[1][i + cnt] +
663 gmat[6][i] * inarray[2][i + cnt];
665 stdVelocity[1][i] = gmat[1][i] * inarray[0][i + cnt] +
666 gmat[4][i] * inarray[1][i + cnt] +
667 gmat[7][i] * inarray[2][i + cnt];
669 stdVelocity[2][i] = gmat[2][i] * inarray[0][i + cnt] +
670 gmat[5][i] * inarray[1][i + cnt] +
671 gmat[8][i] * inarray[2][i + cnt];
676 for (
int i = 0; i < n_points; i++)
678 stdVelocity[0][i] = gmat[0][0] * inarray[0][i + cnt] +
679 gmat[3][0] * inarray[1][i + cnt] +
680 gmat[6][0] * inarray[2][i + cnt];
682 stdVelocity[1][i] = gmat[1][0] * inarray[0][i + cnt] +
683 gmat[4][0] * inarray[1][i + cnt] +
684 gmat[7][0] * inarray[2][i + cnt];
686 stdVelocity[2][i] = gmat[2][0] * inarray[0][i + cnt] +
687 gmat[5][0] * inarray[1][i + cnt] +
688 gmat[8][0] * inarray[2][i + cnt];
694 for (
int i = 0; i < n_points; i++)
696 pntVelocity = stdVelocity[0][i] * stdVelocity[0][i] +
697 stdVelocity[1][i] * stdVelocity[1][i] +
698 stdVelocity[2][i] * stdVelocity[2][i];
700 if (pntVelocity > maxV[el])
702 maxV[el] = pntVelocity;
706 maxV[el] =
sqrt(maxV[el]);
#define ASSERTL0(condition, msg)
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mode...
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
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.
void DefineProjection(FuncPointerT func, ObjectPointerT obj)
void DefineOdeRhs(FuncPointerT func, ObjectPointerT obj)
void DefineImplicitSolve(FuncPointerT func, ObjectPointerT obj)
SolverUtils::AdvectionSharedPtr m_advObject
Advection term.
SOLVER_UTILS_EXPORT int GetTraceNpoints()
NekDouble m_timestep
Time step size.
int m_infosteps
Number of time steps between outputting status information.
NekDouble m_time
Current time of simulation.
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 GetNpoints()
enum MultiRegions::ProjectionType m_projectionType
Type of projection; e.g continuous or discontinuous.
SOLVER_UTILS_EXPORT void SetBoundaryConditions(NekDouble time)
Evaluates the boundary conditions at the given time.
Base class for unsteady solvers.
LibUtilities::TimeIntegrationSchemeOperators m_ode
The time integration scheme operators to use.
NekDouble m_cflSafetyFactor
CFL safety factor (comprise between 0 to 1).
LibUtilities::TimeIntegrationSchemeSharedPtr m_intScheme
Wrapper to the time integration scheme.
bool m_explicitDiffusion
Indicates if explicit or implicit treatment of diffusion is used.
UnsteadyAdvectionDiffusion(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
Session reader.
SolverUtils::DiffusionSharedPtr m_diffusion
bool v_PreIntegrate(int step) override
PreIntegration step for substepping.
void v_InitObject(bool DeclareFields=true) override
Initialise the object.
void DoOdeRhs(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
Compute the RHS.
void SetUpSubSteppingTimeIntegration(const LibUtilities::TimeIntegrationSchemeSharedPtr &IntegrationScheme)
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.
NekDouble m_sVVCutoffRatio
void AddAdvectionPenaltyFlux(const Array< OneD, const Array< OneD, NekDouble > > &velfield, const Array< OneD, const Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &Outarray)
static SolverUtils::EquationSystemSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
Creates an instance of this class.
void SubStepAdvection(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
void SubStepAdvance(int nstep, NekDouble time)
static std::string className
Name of class.
Array< OneD, NekDouble > GetMaxStdVelocity(const Array< OneD, Array< OneD, NekDouble > > inarray)
void GetFluxVectorDiff(const Array< OneD, Array< OneD, NekDouble > > &inarray, const Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &qfield, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &viscousTensor)
Evaluate the flux at each solution point for the diffusion part.
void v_GenerateSummary(SolverUtils::SummaryList &s) override
Print Summary.
void SubStepProjection(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
LibUtilities::TimeIntegrationSchemeSharedPtr m_subStepIntegrationScheme
NekDouble GetSubstepTimeStep()
LibUtilities::TimeIntegrationSchemeOperators m_subStepIntegrationOps
Array< OneD, NekDouble > m_traceVn
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.
void v_GenerateSummary(SolverUtils::SummaryList &s) override
Print Summary.
Array< OneD, Array< OneD, NekDouble > > m_velocity
Advection velocity.
void DoOdeRhs(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
Compute the RHS.
void v_InitObject(bool DeclareFields=true) override
Initialise the object.
bool m_useGJPStabilisation
SolverUtils::RiemannSolverSharedPtr m_riemannSolver
TimeIntegrationSchemeFactory & GetTimeIntegrationSchemeFactory()
std::shared_ptr< SessionReader > SessionReaderSharedPtr
std::vector< PointsKey > PointsKeyVector
std::shared_ptr< TimeIntegrationScheme > TimeIntegrationSchemeSharedPtr
std::shared_ptr< Comm > CommSharedPtr
Pointer to a Communicator object.
@ eMixed_CG_Discontinuous
std::vector< std::pair< std::string, std::string > > SummaryList
DiffusionFactory & GetDiffusionFactory()
EquationSystemFactory & GetEquationSystemFactory()
void AddSummaryItem(SummaryList &l, const std::string &name, const std::string &value)
Adds a summary item to the summary info list.
RiemannSolverFactory & GetRiemannSolverFactory()
@ eDeformed
Geometry is curved or has non-constant factors.
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
std::map< ConstFactorType, NekDouble > ConstFactorMap
StdRegions::ConstFactorMap factors
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.
T Vmin(int n, const T *x, const int incx)
Return the minimum element in x - called vmin to avoid conflict with min.
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 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.
scalarT< T > sqrt(scalarT< T > in)