59 std::vector<std::string> vel;
74 ASSERTL0(
false,
"Implicit unsteady Advection not set up.");
88 int nvariables = inarray.num_elements();
95 int ncoeffs = inarray[0].num_elements();
99 for(j = 1; j < nvariables; ++j)
101 WeakAdv[j] = WeakAdv[j-1] + ncoeffs;
106 for(j = 0; j < nvariables; ++j)
108 m_fields[j]->MultiplyByElmtInvMass(WeakAdv[j], WeakAdv[j]);
109 m_fields[j]->BwdTrans(WeakAdv[j],outarray[j]);
118 for(j = 0; j < nvariables; ++j)
137 int nvariables = inarray.num_elements();
147 for(j = 0; j < nvariables; ++j)
158 for(j = 0; j < nvariables; ++j)
160 m_fields[j]->FwdTrans(inarray[j],coeffs);
161 m_fields[j]->BwdTrans_IterPerExp(coeffs,outarray[j]);
166 ASSERTL0(
false,
"Unknown projection scheme");
178 "Dimension of flux array and velocity array do not match");
180 for (
int j = 0; j < flux.num_elements(); ++j)
202 for (i = 0; i < nvel; ++i)
207 Fwd, 1, Vn, 1, Vn, 1);
210 for (i = 0; i < numflux.num_elements(); ++i)
212 m_fields[i]->GetFwdBwdTracePhys(physfield[i], Fwd, Bwd);
214 m_fields[i]->GetTrace()->Upwind(Vn, Fwd, Bwd, numflux[i]);
216 Vmath::Vmul(nTraceNumPoints, numflux[i], 1, Vn, 1, numflux[i], 1);
235 int n_element =
m_fields[0]->GetExpSize();
245 for (
int el = 0; el < n_element; ++el)
247 int npoints =
m_fields[0]->GetExp(el)->GetTotPoints();
250 if(boost::dynamic_pointer_cast<LocalRegions::TriExp>(
m_fields[0]->GetExp(el)))
254 tstep[el] = CFL[el]/(stdVelocity[el]);
256 else if(boost::dynamic_pointer_cast<LocalRegions::QuadExp>(
m_fields[0]->GetExp(el)))
260 tstep[el] = CFL[el]/(stdVelocity[el]);
284 int n_elements =
m_fields[0]->GetExpSize();
301 if (TimeStability == 1.0)
305 else if (TimeStability == 2.0)
309 else if (TimeStability == 2.784)
315 ASSERTL0(
false,
"Dominant eigenvalues database not present for this time-stepping scheme")
320 TimeStep = (TimeStability/SpatialStability)*CFL;
336 int n_element =
m_fields[0]->GetExpSize();
337 int nvel = inarray.num_elements();
345 for (
int i = 0; i < nvel; ++i)
352 for (
int el = 0; el < n_element; ++el)
359 el2D->GetGeom2D()->GetMetricInfo()->GetJac();
361 el2D->GetGeom2D()->GetMetricInfo()->GetDerivFactors();
363 if (el2D->GetGeom2D()->GetMetricInfo()->GetGtype()
366 for (
int i = 0; i < n_points; i++)
368 stdVelocity[0][i] = gmat[0][i]*inarray[0][i]
369 + gmat[2][i]*inarray[1][i];
371 stdVelocity[1][i] = gmat[1][i]*inarray[0][i]
372 + gmat[3][i]*inarray[1][i];
377 for (
int i = 0; i < n_points; i++)
379 stdVelocity[0][i] = gmat[0][0]*inarray[0][i]
380 + gmat[2][0]*inarray[1][i];
382 stdVelocity[1][i] = gmat[1][0]*inarray[0][i]
383 + gmat[3][0]*inarray[1][i];
388 for (
int i = 0; i < n_points; i++)
390 pntVelocity = sqrt(stdVelocity[0][i]*stdVelocity[0][i]
391 + stdVelocity[1][i]*stdVelocity[1][i]);
393 if (pntVelocity>stdV[el])
395 stdV[el] = pntVelocity;
402 for (
int el = 0; el < n_element; ++el)
410 el3D->GetGeom3D()->GetMetricInfo()->GetJac();
412 el3D->GetGeom3D()->GetMetricInfo()->GetDerivFactors();
414 if (el3D->GetGeom3D()->GetMetricInfo()->GetGtype()
417 for (
int i = 0; i < n_points; i++)
419 stdVelocity[0][i] = gmat[0][i]*inarray[0][i]
420 + gmat[3][i]*inarray[1][i]
421 + gmat[6][i]*inarray[2][i];
423 stdVelocity[1][i] = gmat[1][i]*inarray[0][i]
424 + gmat[4][i]*inarray[1][i]
425 + gmat[7][i]*inarray[2][i];
427 stdVelocity[2][i] = gmat[2][i]*inarray[0][i]
428 + gmat[5][i]*inarray[1][i]
429 + gmat[8][i]*inarray[2][i];
435 el3D->GetGeom3D()->GetMetricInfo()->GetJac();
437 el3D->GetGeom3D()->GetMetricInfo()->GetDerivFactors();
439 for (
int i = 0; i < n_points; i++)
441 stdVelocity[0][i] = gmat[0][0]*inarray[0][i]
442 + gmat[3][0]*inarray[1][i]
443 + gmat[6][0]*inarray[2][i];
445 stdVelocity[1][i] = gmat[1][0]*inarray[0][i]
446 + gmat[4][0]*inarray[1][i]
447 + gmat[7][0]*inarray[2][i];
449 stdVelocity[2][i] = gmat[2][0]*inarray[0][i]
450 + gmat[5][0]*inarray[1][i]
451 + gmat[8][0]*inarray[2][i];
455 for (
int i = 0; i < n_points; i++)
457 pntVelocity = sqrt(stdVelocity[0][i]*stdVelocity[0][i]
458 + stdVelocity[1][i]*stdVelocity[1][i]
459 + stdVelocity[2][i]*stdVelocity[2][i]);
461 if (pntVelocity > stdV[el])
463 stdV[el] = pntVelocity;
f SOLVER_UTILS_EXPORT void AdvectionNonConservativeForm(const Array< OneD, Array< OneD, NekDouble > > &V, const Array< OneD, const NekDouble > &u, Array< OneD, NekDouble > &outarray, Array< OneD, NekDouble > &wk=NullNekDouble1DArray)
Compute the non-conservative advection.
#define ASSERTL0(condition, msg)
static NekDouble EigenvaluesAnaMeshesRK2[6][14]
virtual void v_InitObject()
Init object for UnsteadySystem class.
LibUtilities::TimeIntegrationSchemeOperators m_ode
The time integration scheme operators to use.
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.
int m_expdim
Expansion dimension.
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.
Array< OneD, NekDouble > GetStdVelocity(const Array< OneD, Array< OneD, NekDouble > > inarray)
virtual void v_GenerateSummary(SummaryList &s)
Print a summary of time stepping parameters.
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
SOLVER_UTILS_EXPORT int GetTotPoints()
static NekDouble EigenvaluesAnaMeshesRK4[6][14]
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Array holding trace normals for DG simulations in the forwards direction.
int GetTotPoints() const
This function returns the total number of quadrature points used in the element.
void DefineProjection(FuncPointerT func, ObjectPointerT obj)
virtual SOLVER_UTILS_EXPORT void v_GenerateSummary(SummaryList &s)
Print a summary of time stepping parameters.
static std::string className
bool m_explicitAdvection
Indicates if explicit or implicit treatment of advection is used.
void DoOdeRhs(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
void DefineOdeRhs(FuncPointerT func, ObjectPointerT obj)
Base class for unsteady solvers.
boost::shared_ptr< Expansion3D > Expansion3DSharedPtr
virtual NekDouble v_GetTimeStep(const Array< OneD, int > ExpOrder, const Array< OneD, NekDouble > CFL, NekDouble timeCFL)
int m_spacedim
Spatial dimension (>= expansion dim).
virtual void v_GetFluxVector(const int i, Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &flux)
void Neg(int n, T *x, const int incx)
Negate x = -x.
virtual SOLVER_UTILS_EXPORT void v_InitObject()
Init object for UnsteadySystem class.
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 DoOdeProjection(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
SOLVER_UTILS_EXPORT void SetBoundaryConditions(NekDouble time)
Evaluates the boundary conditions at the given time.
SOLVER_UTILS_EXPORT int GetNpoints()
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
SOLVER_UTILS_EXPORT int GetTraceNpoints()
SOLVER_UTILS_EXPORT int GetNcoeffs()
Array< OneD, Array< OneD, NekDouble > > m_velocity
static NekDouble EigenvaluesAnaMeshesAB2[6][14]
CFLtester(const LibUtilities::SessionReaderSharedPtr &pSession)
SOLVER_UTILS_EXPORT void WeakDGAdvection(const Array< OneD, Array< OneD, NekDouble > > &InField, Array< OneD, Array< OneD, NekDouble > > &OutField, bool NumericalFluxIncludesNormal=true, bool InFieldIsInPhysSpace=false, int nvariables=0)
Calculate the weak discontinuous Galerkin advection.
static EquationSystemSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession)
#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.
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_NumericalFlux(Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &numflux)
boost::shared_ptr< Expansion2D > Expansion2DSharedPtr
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.