36#include <boost/algorithm/string.hpp>
85 m_SmoothAdvection(false)
95 std::string velids[] = {
"u",
"v",
"w"};
102 for (j = 0; j < numfields; ++j)
105 if (boost::iequals(velids[i], var))
111 ASSERTL0(j != numfields,
"Failed to find field: " + var);
131 std::string vConvectiveType;
136 vConvectiveType =
"NoAdvection";
140 vConvectiveType =
"Convective";
143 vConvectiveType =
"Linearised";
150 if (
m_session->DefinesTag(
"AdvectiveType") &&
154 vConvectiveType =
m_session->GetTag(
"AdvectiveType");
159 vConvectiveType, vConvectiveType);
171 "X",
"Y",
"Z",
"Theta_x",
"Theta_y",
"Theta_z",
172 "U",
"V",
"W",
"Omega_x",
"Omega_y",
"Omega_z",
173 "A_x",
"A_y",
"A_z",
"DOmega_x",
"DOmega_y",
"DOmega_z",
189 for (
size_t i = 0; i <
m_fields.size(); ++i)
197 BndConds =
m_fields[i]->GetBndConditions();
198 BndExp =
m_fields[i]->GetBndCondExpansions();
199 for (
size_t n = 0; n < BndConds.size(); ++n)
201 if (boost::iequals(BndConds[n]->GetUserDefined(),
"Radiation"))
204 BndConds[n]->GetBoundaryConditionType() ==
206 "Radiation boundary condition must be of type Robin <R>");
214 radpts += BndExp[n]->GetTotPoints();
216 if (boost::iequals(BndConds[n]->GetUserDefined(),
217 "ZeroNormalComponent"))
219 ASSERTL0(BndConds[n]->GetBoundaryConditionType() ==
221 "Zero Normal Component boundary condition option must "
222 "be of type Dirichlet <D>");
237 for (
size_t n = 0; n < BndConds.size(); ++n)
239 if (boost::iequals(BndConds[n]->GetUserDefined(),
"Radiation"))
242 int npoints = BndExp[n]->GetNpoints();
248 BndExp[n]->GetCoords(x0, x1, x2);
251 std::static_pointer_cast<
253 ->m_robinPrimitiveCoeff;
265 for (
size_t i = 0; i <
m_fields.size(); ++i)
267 for (
size_t n = 0; n <
m_fields[i]->GetBndConditions().size(); ++n)
269 if (boost::istarts_with(
270 m_fields[i]->GetBndConditions()[n]->GetUserDefined(),
281 i, n,
m_fields[i]->GetBndConditions()[n]->GetUserDefined());
305 size_t npoints =
m_fields[0]->GetNpoints();
306 for (
size_t i = 0; i < VelDim; ++i)
313 x->PreApply(
m_fields, velocity, velocity, time);
327 size_t nvariables =
m_fields.size();
329 for (i = 0; i < nvariables; ++i)
331 for (n = 0; n <
m_fields[i]->GetBndConditions().size(); ++n)
333 if (
m_fields[i]->GetBndConditions()[n]->IsTimeDependent())
336 m_fields[i]->EvaluateBoundaryConditions(time, varName);
338 else if (boost::istarts_with(
339 m_fields[i]->GetBndConditions()[n]->GetUserDefined(),
365 BndConds =
m_fields[fieldid]->GetBndConditions();
366 BndExp =
m_fields[fieldid]->GetBndCondExpansions();
372 size_t elmtid, nq, offset, boundary;
376 for (cnt = n = 0; n < BndConds.size(); ++n)
378 std::string type = BndConds[n]->GetUserDefined();
380 if ((BndConds[n]->GetBoundaryConditionType() ==
382 (boost::iequals(type,
"Radiation")))
384 size_t nExp = BndExp[n]->GetExpSize();
385 for (i = 0; i < nExp; ++i, cnt++)
388 elmt =
m_fields[fieldid]->GetExp(elmtid);
389 offset =
m_fields[fieldid]->GetPhys_Offset(elmtid);
391 U =
m_fields[fieldid]->UpdatePhys() + offset;
392 Bc = BndExp[n]->GetExp(i);
397 nq = Bc->GetTotPoints();
399 elmt->GetTracePhysVals(boundary, Bc, U, ubc);
404 1, &ubc[0], 1, &ubc[0], 1);
407 BndExp[n]->UpdateCoeffs() + BndExp[n]->GetCoeff_Offset(i);
409 Bc->IProductWRTBase(ubc, Bvals);
411 cnt1 += BndExp[n]->GetTotPoints();
415 cnt += BndExp[n]->GetExpSize();
424 static bool Setup =
false;
447 size_t elmtid, nq, boundary;
454 for (cnt = n = 0; n < BndConds[0].size(); ++n)
456 if ((BndConds[0][n]->GetBoundaryConditionType() ==
458 (boost::iequals(BndConds[0][n]->GetUserDefined(),
459 "ZeroNormalComponent")))
461 size_t nExp = BndExp[0][n]->GetExpSize();
462 for (i = 0; i < nExp; ++i, cnt++)
468 normals = elmt->GetTraceNormal(boundary);
470 nq = BndExp[0][n]->GetExp(i)->GetTotPoints();
475 Bphys = BndExp[k][n]->UpdatePhys() +
476 BndExp[k][n]->GetPhys_Offset(i);
477 Bc = BndExp[k][n]->GetExp(i);
487 Bphys = BndExp[k][n]->UpdatePhys() +
488 BndExp[k][n]->GetPhys_Offset(i);
489 Bcoeffs = BndExp[k][n]->UpdateCoeffs() +
490 BndExp[k][n]->GetCoeff_Offset(i);
491 Bc = BndExp[k][n]->GetExp(i);
492 Vmath::Vvtvp(nq, normvel, 1, normals[k], 1, Bphys, 1, Bphys,
494 Bc->FwdTransBndConstrained(Bphys, Bcoeffs);
500 cnt += BndExp[0][n]->GetExpSize();
511 "Womersley parameters for this boundary have not been set up");
517 size_t M_coeffs = WomParam->m_wom_vel.size();
520 NekDouble axis_normal = WomParam->m_axisnormal[fldid];
530 BndCondExp =
m_fields[fldid]->GetBndCondExpansions()[bndid];
536 size_t exp_npts = BndCondExp->GetExpSize();
542 for (k = 1; k < M_coeffs; k++)
545 zt[k] = std::exp(zi * omega_c * k_c * m_time_c);
549 for (i = 0; i < exp_npts; ++i, cnt++)
552 bc = BndCondExp->GetExp(i);
553 nfq = bc->GetTotPoints();
557 for (j = 0; j < nfq; j++)
559 wbc[j] = WomParam->m_poiseuille[i][j];
560 for (k = 1; k < M_coeffs; k++)
562 zvel = WomParam->m_zvel[i][j][k] * zt[k];
563 wbc[j] = wbc[j] + zvel.real();
570 Bvals = BndCondExp->UpdateCoeffs() + BndCondExp->GetCoeff_Offset(i);
573 bc->FwdTrans(wbc, Bvals);
580 std::string::size_type indxBeg = womStr.find_first_of(
':') + 1;
581 std::string filename = womStr.substr(indxBeg, std::string::npos);
583 TiXmlDocument doc(filename);
585 bool loadOkay = doc.LoadFile();
587 (std::string(
"Failed to load file: ") + filename).c_str());
589 TiXmlHandle docHandle(&doc);
593 TiXmlElement *nektar = doc.FirstChildElement(
"NEKTAR");
594 ASSERTL0(nektar,
"Unable to find NEKTAR tag in file.");
596 TiXmlElement *wombc = nektar->FirstChildElement(
"WOMERSLEYBC");
597 ASSERTL0(wombc,
"Unable to find WOMERSLEYBC tag in file.");
600 TiXmlElement *womparam = wombc->FirstChildElement(
"WOMPARAMS");
601 ASSERTL0(womparam,
"Unable to find WOMPARAMS tag in file.");
604 TiXmlElement *params = womparam->FirstChildElement(
"W");
605 std::map<std::string, std::string> Wparams;
612 propstr = params->Attribute(
"PROPERTY");
615 "Failed to read PROPERTY value Womersley BC Parameter");
618 valstr = params->Attribute(
"VALUE");
621 "Failed to read VALUE value Womersley BC Parameter");
623 std::transform(propstr.begin(), propstr.end(), propstr.begin(),
625 Wparams[propstr] = valstr;
627 params = params->NextSiblingElement(
"W");
634 Wparams.count(
"RADIUS") == 1,
635 "Failed to find Radius parameter in Womersley boundary conditions");
636 std::vector<NekDouble>
rad;
641 Wparams.count(
"PERIOD") == 1,
642 "Failed to find period parameter in Womersley boundary conditions");
643 std::vector<NekDouble> period;
648 Wparams.count(
"AXISNORMAL") == 1,
649 "Failed to find axisnormal parameter in Womersley boundary conditions");
650 std::vector<NekDouble> anorm;
657 Wparams.count(
"AXISPOINT") == 1,
658 "Failed to find axispoint parameter in Womersley boundary conditions");
659 std::vector<NekDouble> apt;
668 TiXmlElement *coeff = wombc->FirstChildElement(
"FOURIERCOEFFS");
671 TiXmlElement *fval = coeff->FirstChildElement(
"F");
677 TiXmlAttribute *fvalAttr = fval->FirstAttribute();
678 std::string attrName(fvalAttr->Name());
681 (std::string(
"Unknown attribute name: ") + attrName).c_str());
683 err = fvalAttr->QueryIntValue(&indx);
684 ASSERTL0(err == TIXML_SUCCESS,
"Unable to read attribute ID.");
686 std::string coeffStr = fval->FirstChild()->ToText()->ValueStr();
687 std::vector<NekDouble> coeffvals;
692 (std::string(
"Problem reading value of fourier coefficient, ID=") +
693 boost::lexical_cast<std::string>(indx))
696 coeffvals.size() == 2,
698 "Have not read two entries of Fourier coefficicent from ID=" +
699 boost::lexical_cast<std::string>(indx))
705 fval = fval->NextSiblingElement(
"F");
724 BndCondExp =
m_fields[fldid]->GetBndCondExpansions()[bndid];
731 size_t exp_npts = BndCondExp->GetExpSize();
745 for (k = 1; k < M_coeffs; k++)
748 lamda_n[k] = i_pow_3q2 * alpha_c *
sqrt(k_c[k]);
753 for (i = 0; i < exp_npts; ++i, cnt++)
756 bc = BndCondExp->GetExp(i);
757 nfq = bc->GetTotPoints();
762 bc->GetCoords(x, y,
z);
770 for (j = 0; j < nfq; j++)
773 (y[j] - x0[1]) * (y[j] - x0[1]) +
774 (
z[j] - x0[2]) * (
z[j] - x0[2])) /
781 (1. - rqR.real() * rqR.real());
787 for (k = 1; k < M_coeffs; k++)
792 (z1 - (zJ0r / zJ0[k]));
810 [[maybe_unused]]
const NekDouble SpeedSoundFactor)
813 size_t nelmt =
m_fields[0]->GetExpSize();
822 for (
size_t i = 0; i < 2; ++i)
831 for (
size_t i = 0; i < nvel; ++i)
849 if (physfield.size())
851 pressure = physfield[physfield.size() - 1];
862 int nPts = physfield[0].size();
875 velocity[i] = physfield[i];
891 "Arrays have different dimensions, cannot set moving frame "
894 Vmath::Vcopy(vFrameVels.size(), vFrameVels, 1, temp, 1);
904 "Arrays have different dimensions, cannot get moving frame "
926 vFrameDisp.size() == 6,
927 "Arrays have different size, cannot set moving frame displacement");
929 Vmath::Vcopy(vFrameDisp.size(), vFrameDisp, 1, temp, 1);
943 vFrameDisp.size() == 6,
944 "Arrays have different size, cannot get moving frame displacement");
959 "Arrays have different size, cannot set moving frame pivot");
961 Vmath::Vcopy(vFramePivot.size(), vFramePivot, 1, temp, 1);
963 Vmath::Vcopy(vFramePivot.size(), vFramePivot, 1, temp, 1);
987 std::vector<std::string> vForceList;
988 bool hasForce{
false};
990 if (!
m_session->DefinesElement(
"Nektar/Forcing"))
995 TiXmlElement *vForcing =
m_session->GetElement(
"Nektar/Forcing");
998 TiXmlElement *vForce = vForcing->FirstChildElement(
"FORCE");
1001 std::string vType = vForce->Attribute(
"TYPE");
1003 vForceList.push_back(vType);
1004 vForce = vForce->NextSiblingElement(
"FORCE");
1008 for (
auto &f : vForceList)
1010 if (boost::iequals(f, sForce))
#define ASSERTL0(condition, msg)
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
Array< OneD, Array< OneD, int > > m_fieldsBCToTraceID
Mapping from BCs to Elmt Edge IDs.
Array< OneD, NekDouble > v_GetMaxStdVelocity(const NekDouble SpeedSoundFactor) override
std::map< int, std::map< int, WomersleyParamsSharedPtr > > m_womersleyParams
Womersley parameters if required.
void v_SetAeroForce(Array< OneD, NekDouble > forces) override
bool v_GetMovingFrameVelocities(Array< OneD, NekDouble > &vFrameVels, const int step) override
virtual int v_GetForceDimension()=0
void SetWomersleyBoundary(const int fldid, const int bndid)
Set Womersley Profile if specified.
void v_GetAeroForce(Array< OneD, NekDouble > forces) override
void SetZeroNormalVelocity()
Set Normal Velocity Component to Zero.
void v_SetMovingFrameDisp(const Array< OneD, NekDouble > &vFrameDisp, const int step) override
MultiRegions::ExpListSharedPtr v_GetPressure() override
Array< OneD, NekDouble > m_aeroForces
void SetBoundaryConditions(NekDouble time)
time dependent boundary conditions updating
bool v_GetMovingFrameDisp(Array< OneD, NekDouble > &vFrameDisp, const int step) override
NekDouble m_kinvis
Kinematic viscosity.
void v_GetDensity(const Array< OneD, const Array< OneD, NekDouble > > &physfield, Array< OneD, NekDouble > &density) override
void SetRadiationBoundaryForcing(int fieldid)
Set Radiation forcing term.
Array< OneD, Array< OneD, NekDouble > > m_fieldsRadiationFactor
RHS Factor for Radiation Condition.
bool DefinedForcing(const std::string &sForce)
bool v_PreIntegrate(int step) override
void v_SetMovingFrameVelocities(const Array< OneD, NekDouble > &vFrameVels, const int step) override
void SetUpWomersley(const int fldid, const int bndid, std::string womstr)
Set Up Womersley details.
void v_InitObject(bool DeclareField=true) override
Initialisation object for EquationSystem.
ExtrapolateSharedPtr m_extrapolation
Array< OneD, int > m_velocity
int which identifies which components of m_fields contains the velocity (u,v,w);
static std::string eqTypeLookupIds[]
void v_SetMovingFramePivot(const Array< OneD, NekDouble > &vFramePivot) override
Array< OneD, Array< OneD, int > > m_fieldsBCToElmtID
Mapping from BCs to Elmt IDs.
void v_GetVelocity(const Array< OneD, const Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &velocity) override
EquationType m_equationType
equation type;
int m_nConvectiveFields
Number of fields to be convected;.
void AddForcing(const SolverUtils::ForcingSharedPtr &pForce)
std::vector< SolverUtils::ForcingSharedPtr > m_forcing
Forcing terms.
void EvaluateAdvectionTerms(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
IncNavierStokes(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
Constructor.
NekDouble Evaluate() const
tBaseSharedPtr CreateInstance(tKey idKey, tParam... args)
Create an instance of the class referred to by idKey.
static std::string RegisterEnumValue(std::string pEnum, std::string pString, int pEnumValue)
Registers an enumeration value.
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
static bool GenerateVector(const std::string &str, std::vector< T > &out)
Takes a comma-separated string and converts it to entries in a vector.
A base class for PDEs which include an advection component.
SolverUtils::AdvectionSharedPtr m_advObject
Advection term.
SOLVER_UTILS_EXPORT void v_InitObject(bool DeclareField=true) override
Initialisation object for EquationSystem.
int m_spacedim
Spatial dimension (>= expansion dim).
NekDouble m_timestep
Time step size.
NekDouble m_time
Current time of simulation.
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
std::vector< std::string > m_strFrameData
variable name in m_movingFrameData
SOLVER_UTILS_EXPORT int GetPhys_Offset(int n)
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
enum HomogeneousType m_HomogeneousType
LibUtilities::FieldMetaDataMap m_fieldMetaDataMap
Map to identify relevant solver info to dump in output fields.
SpatialDomains::BoundaryConditionsSharedPtr m_boundaryConditions
Pointer to boundary conditions object.
Array< OneD, NekDouble > m_movingFrameData
Moving reference frame status in the inertial frame X, Y, Z, Theta_x, Theta_y, Theta_z,...
static SOLVER_UTILS_EXPORT std::vector< ForcingSharedPtr > Load(const LibUtilities::SessionReaderSharedPtr &pSession, const std::weak_ptr< EquationSystem > &pEquation, const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const unsigned int &pNumForcingFields=0)
Base class for unsteady solvers.
static NekDouble rad(NekDouble x, NekDouble y)
std::shared_ptr< SessionReader > SessionReaderSharedPtr
std::shared_ptr< Expansion > ExpansionSharedPtr
std::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
AdvectionFactory & GetAdvectionFactory()
Gets the factory for initialising advection objects.
SOLVER_UTILS_EXPORT typedef std::shared_ptr< Forcing > ForcingSharedPtr
A shared pointer to an EquationSystem object.
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
std::shared_ptr< StdExpansion > StdExpansionSharedPtr
std::vector< double > z(NPUPPER)
std::complex< double > NekComplexDouble
std::shared_ptr< WomersleyParams > WomersleyParamsSharedPtr
const std::string kEquationTypeStr[]
std::complex< Nektar::NekDouble > ImagBesselComp(int n, std::complex< Nektar::NekDouble > y)
Calcualte the bessel function of the first kind with complex double input y. Taken from Numerical Rec...
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.
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
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
scalarT< T > sqrt(scalarT< T > in)