38 #include <boost/algorithm/string.hpp>
61 m_SmoothAdvection(false),
71 int numfields =
m_fields.num_elements();
72 std::string velids[] = {
"u",
"v",
"w"};
79 for(j = 0; j < numfields; ++j)
82 if(boost::iequals(velids[i], var))
88 ASSERTL0(j != numfields,
"Failed to find field: " + var);
125 ASSERTL0(
false,
"Unknown or undefined equation type");
131 std::string vConvectiveType;
136 vConvectiveType =
"NoAdvection";
140 vConvectiveType =
"Convective";
143 vConvectiveType =
"Linearised";
153 vConvectiveType =
m_session->GetTag(
"AdvectiveType");
170 for (i = 0; i <
m_fields.num_elements(); ++i)
178 BndConds =
m_fields[i]->GetBndConditions();
179 BndExp =
m_fields[i]->GetBndCondExpansions();
180 for(
int n = 0; n < BndConds.num_elements(); ++n)
182 if(boost::iequals(BndConds[n]->GetUserDefined(),
"Radiation"))
185 "Radiation boundary condition must be of type Robin <R>");
192 radpts += BndExp[n]->GetTotPoints();
194 if(boost::iequals(BndConds[n]->GetUserDefined(),
"ZeroNormalComponent"))
197 "Zero Normal Component boundary condition option must be of type Dirichlet <D>");
211 for(
int n = 0; n < BndConds.num_elements(); ++n)
213 if(boost::iequals(BndConds[n]->GetUserDefined(),
"Radiation"))
216 int npoints = BndExp[n]->GetNpoints();
222 BndExp[n]->GetCoords(x0,x1,x2);
225 boost::static_pointer_cast<
227 >(BndConds[n])->m_robinPrimitiveCoeff;
230 tmpArray = m_fieldsRadiationFactor[i]+ radpts);
257 ASSERTL1(flux.num_elements() ==
m_velocity.num_elements(),
"Dimension of flux array and velocity array do not match");
259 for(
int j = 0; j < flux.num_elements(); ++j)
290 for(i = 0; i < nDimensions; ++i)
297 for(i = 0; i < numflux.num_elements(); ++i)
300 m_fields[i]->GetFwdBwdTracePhys(physfield[i], Fwd, Bwd);
303 m_fields[i]->GetTrace()->Upwind(Vn, Fwd, Bwd, numflux[i]);
307 Vmath::Vmul(nTracePts, numflux[i], 1, Vn, 1, numflux[i], 1);
319 int nqtot =
m_fields[0]->GetTotPoints();
324 for(i = 0; i < VelDim; ++i)
338 if(wk.num_elements())
340 ASSERTL0(wk.num_elements() >= nqtot*VelDim,
341 "Workspace is not sufficient");
350 velocity, inarray, outarray,
m_time);
360 int nvariables =
m_fields.num_elements();
362 for (i = 0; i < nvariables; ++i)
364 for(n = 0; n <
m_fields[i]->GetBndConditions().num_elements(); ++n)
366 if(
m_fields[i]->GetBndConditions()[n]->IsTimeDependent())
369 m_fields[i]->EvaluateBoundaryConditions(time, varName);
391 BndConds =
m_fields[fieldid]->GetBndConditions();
392 BndExp =
m_fields[fieldid]->GetBndCondExpansions();
398 int elmtid,nq,offset, boundary;
402 for(cnt = n = 0; n < BndConds.num_elements(); ++n)
404 std::string type = BndConds[n]->GetUserDefined();
406 if((BndConds[n]->GetBoundaryConditionType() ==
SpatialDomains::eRobin)&&(boost::iequals(type,
"Radiation")))
408 for(i = 0; i < BndExp[n]->GetExpSize(); ++i,cnt++)
411 elmt =
m_fields[fieldid]->GetExp(elmtid);
412 offset =
m_fields[fieldid]->GetPhys_Offset(elmtid);
414 U =
m_fields[fieldid]->UpdatePhys() + offset;
415 Bc = BndExp[n]->GetExp(i);
420 nq = Bc->GetTotPoints();
422 elmt->GetTracePhysVals(boundary,Bc,U,ubc);
427 Bvals = BndExp[n]->UpdateCoeffs()+BndExp[n]->GetCoeff_Offset(i);
429 Bc->IProductWRTBase(ubc,Bvals);
431 cnt1 += BndExp[n]->GetTotPoints();
435 cnt += BndExp[n]->GetExpSize();
445 static bool Setup =
false;
470 int elmtid,nq, boundary;
477 for(cnt = n = 0; n < BndConds[0].num_elements(); ++n)
479 if((BndConds[0][n]->GetBoundaryConditionType() ==
SpatialDomains::eDirichlet)&& (boost::iequals(BndConds[0][n]->GetUserDefined(),
"ZeroNormalComponent")))
481 for(i = 0; i < BndExp[0][n]->GetExpSize(); ++i,cnt++)
487 normals = elmt->GetSurfaceNormal(boundary);
489 nq = BndExp[0][n]->GetExp(i)->GetTotPoints();
494 Bphys = BndExp[k][n]->UpdatePhys()+
495 BndExp[k][n]->GetPhys_Offset(i);
496 Bc = BndExp[k][n]->GetExp(i);
506 Bphys = BndExp[k][n]->UpdatePhys()+
507 BndExp[k][n]->GetPhys_Offset(i);
508 Bcoeffs = BndExp[k][n]->UpdateCoeffs()+
509 BndExp[k][n]->GetCoeff_Offset(i);
510 Bc = BndExp[k][n]->GetExp(i);
513 Bc->FwdTrans_BndConstrained(Bphys,Bcoeffs);
519 cnt += BndExp[0][n]->GetExpSize();
542 bool returnval =
false;
547 int ncoeffs =
m_fields[0]->GetNcoeffs();
549 for(
int i = 0; i <
m_fields.num_elements(); ++i)
570 int n_element =
m_fields[0]->GetExpSize();
585 for(
int i = 0; i < 2; ++i)
594 for(
int i = 0; i < n_vel; ++i)
602 for(
int el = 0; el < n_element; ++el)
604 cfl[el] =
m_timestep*(stdVelocity[el] * cLambda *
605 (ExpOrder[el]-1) * (ExpOrder[el]-1));
616 int n_element =
m_fields[0]->GetExpSize();
623 CFL = CFL_loc = cfl[elmtid];
627 elmtid =
m_fields[0]->GetExp(elmtid)->GetGeom()->GetGlobalID();
638 elmtid = elmtid%
m_fields[0]->GetPlane(0)->GetExpSize();
665 if(
m_comm->GetRank() == 0)
667 cout <<
"CFL (zero plane): "<< cfl <<
" (in elmt "
668 << elmtid <<
")" << endl;
676 cout <<
"Reached Steady State to tolerance "
EquationType m_equationType
equation type;
bool m_singleMode
Flag to determine if single homogeneous mode is used.
void SetZeroNormalVelocity()
Set Normal Velocity Component to Zero.
#define ASSERTL0(condition, msg)
void SetBoundaryConditions(NekDouble time)
time dependent boundary conditions updating
Array< OneD, Array< OneD, int > > m_fieldsBCToTraceID
Mapping from BCs to Elmt Edge IDs.
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.
virtual bool v_PreIntegrate(int step)
NekDouble m_time
Current time of simulation.
void SetRadiationBoundaryForcing(int fieldid)
Set Radiation forcing term.
NekDouble m_kinvis
Kinematic viscosity.
bool CalcSteadyState(void)
evaluate steady state
NekDouble m_timestep
Time step size.
bool m_halfMode
Flag to determine if half homogeneous mode is used.
SOLVER_UTILS_EXPORT const Array< OneD, int > GetNumExpModesPerExp()
Array< OneD, int > m_velocity
int which identifies which components of m_fields contains the velocity (u,v,w);
ExtrapolateSharedPtr m_extrapolation
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
SOLVER_UTILS_EXPORT typedef boost::shared_ptr< Forcing > ForcingSharedPtr
A shared pointer to an EquationSystem object.
static SOLVER_UTILS_EXPORT std::vector< ForcingSharedPtr > Load(const LibUtilities::SessionReaderSharedPtr &pSession, const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const unsigned int &pNumForcingFields=0)
virtual ~IncNavierStokes()
virtual bool v_PostIntegrate(int step)
int m_nConvectiveFields
Number of fields to be convected;.
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
int Imax(int n, const T *x, const int incx)
Return the index of the maximum element in x.
LibUtilities::CommSharedPtr m_comm
Communicator.
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Array holding trace normals for DG simulations in the forwards direction.
Array< OneD, Array< OneD, NekDouble > > m_fieldsRadiationFactor
RHS Factor for Radiation Condition.
int m_cflsteps
dump cfl estimate
std::vector< SolverUtils::ForcingSharedPtr > m_forcing
Forcing terms.
Base class for unsteady solvers.
NekDouble Evaluate() const
int m_steadyStateSteps
Check for steady state at step interval.
SpatialDomains::BoundaryConditionsSharedPtr m_boundaryConditions
Pointer to boundary conditions object.
int m_spacedim
Spatial dimension (>= expansion dim).
LibUtilities::FieldMetaDataMap m_fieldMetaDataMap
Map to identify relevant solver info to dump in output fields.
AdvectionFactory & GetAdvectionFactory()
Gets the factory for initialising advection objects.
void Neg(int n, T *x, const int incx)
Negate x = -x.
virtual void v_NumericalFlux(Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &numflux)
T Dot(int n, const T *w, const T *x)
vvtvp (vector times vector times vector): z = w*x*y
SOLVER_UTILS_EXPORT int GetPhys_Offset(int n)
SOLVER_UTILS_EXPORT int GetNpoints()
void EvaluateAdvectionTerms(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, Array< OneD, NekDouble > &wk=NullNekDouble1DArray)
void AddForcing(const SolverUtils::ForcingSharedPtr &pForce)
Array< OneD, Array< OneD, int > > m_fieldsBCToElmtID
Mapping from BCs to Elmt IDs.
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
SOLVER_UTILS_EXPORT int GetTraceNpoints()
const std::string kEquationTypeStr[]
int m_infosteps
Number of time steps between outputting status information.
NekDouble m_steadyStateTol
Tolerance to which steady state should be evaluated at.
virtual SOLVER_UTILS_EXPORT void v_InitObject()
Init object for UnsteadySystem class.
A base class for PDEs which include an advection component.
boost::shared_ptr< StdExpansion > StdExpansionSharedPtr
virtual void v_GetFluxVector(const int i, Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &flux)
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
LibUtilities::TimeIntegrationSolutionSharedPtr m_intSoln
virtual void v_InitObject()
Init object for UnsteadySystem class.
virtual int v_GetForceDimension()=0
enum HomogeneousType m_HomogeneousType
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.
NekDouble GetCFLEstimate(int &elmtid)
Array< OneD, NekDouble > GetElmtCFLVals(void)