38 #include <boost/algorithm/string.hpp>
73 m_SmoothAdvection(false),
83 int numfields =
m_fields.num_elements();
84 std::string velids[] = {
"u",
"v",
"w"};
91 for(j = 0; j < numfields; ++j)
94 if(boost::iequals(velids[i], var))
100 ASSERTL0(j != numfields,
"Failed to find field: " + var);
137 ASSERTL0(
false,
"Unknown or undefined equation type");
143 std::string vConvectiveType;
148 vConvectiveType =
"NoAdvection";
152 vConvectiveType =
"Convective";
155 vConvectiveType =
"Linearised";
165 vConvectiveType =
m_session->GetTag(
"AdvectiveType");
182 for (i = 0; i <
m_fields.num_elements(); ++i)
190 BndConds =
m_fields[i]->GetBndConditions();
191 BndExp =
m_fields[i]->GetBndCondExpansions();
192 for(
int n = 0; n < BndConds.num_elements(); ++n)
194 if(boost::iequals(BndConds[n]->GetUserDefined(),
"Radiation"))
197 "Radiation boundary condition must be of type Robin <R>");
204 radpts += BndExp[n]->GetTotPoints();
206 if(boost::iequals(BndConds[n]->GetUserDefined(),
"ZeroNormalComponent"))
209 "Zero Normal Component boundary condition option must be of type Dirichlet <D>");
223 for(
int n = 0; n < BndConds.num_elements(); ++n)
225 if(boost::iequals(BndConds[n]->GetUserDefined(),
"Radiation"))
228 int npoints = BndExp[n]->GetNpoints();
234 BndExp[n]->GetCoords(x0,x1,x2);
237 boost::static_pointer_cast<
239 >(BndConds[n])->m_robinPrimitiveCoeff;
242 tmpArray = m_fieldsRadiationFactor[i]+ radpts);
250 for (
int i = 0; i <
m_fields.num_elements(); ++i)
252 for(
int n = 0; n <
m_fields[i]->GetBndConditions().num_elements(); ++n)
254 if(boost::istarts_with(
m_fields[i]->GetBndConditions()[n]->GetUserDefined(),
"Womersley"))
283 m_fields[i]->GetBndConditions()[n]->GetUserDefined());
311 ASSERTL1(flux.num_elements() ==
m_velocity.num_elements(),
"Dimension of flux array and velocity array do not match");
313 for(
int j = 0; j < flux.num_elements(); ++j)
344 for(i = 0; i < nDimensions; ++i)
351 for(i = 0; i < numflux.num_elements(); ++i)
354 m_fields[i]->GetFwdBwdTracePhys(physfield[i], Fwd, Bwd);
357 m_fields[i]->GetTrace()->Upwind(Vn, Fwd, Bwd, numflux[i]);
361 Vmath::Vmul(nTracePts, numflux[i], 1, Vn, 1, numflux[i], 1);
375 for(i = 0; i < VelDim; ++i)
381 velocity, inarray, outarray,
m_time);
391 int nvariables =
m_fields.num_elements();
393 for (i = 0; i < nvariables; ++i)
395 for(n = 0; n <
m_fields[i]->GetBndConditions().num_elements(); ++n)
397 if(
m_fields[i]->GetBndConditions()[n]->IsTimeDependent())
400 m_fields[i]->EvaluateBoundaryConditions(time, varName);
402 else if(boost::istarts_with(
m_fields[i]->GetBndConditions()[n]->GetUserDefined(),
"Womersley"))
425 BndConds =
m_fields[fieldid]->GetBndConditions();
426 BndExp =
m_fields[fieldid]->GetBndCondExpansions();
432 int elmtid,nq,offset, boundary;
436 for(cnt = n = 0; n < BndConds.num_elements(); ++n)
438 std::string type = BndConds[n]->GetUserDefined();
440 if((BndConds[n]->GetBoundaryConditionType() ==
SpatialDomains::eRobin)&&(boost::iequals(type,
"Radiation")))
442 for(i = 0; i < BndExp[n]->GetExpSize(); ++i,cnt++)
445 elmt =
m_fields[fieldid]->GetExp(elmtid);
446 offset =
m_fields[fieldid]->GetPhys_Offset(elmtid);
448 U =
m_fields[fieldid]->UpdatePhys() + offset;
449 Bc = BndExp[n]->GetExp(i);
454 nq = Bc->GetTotPoints();
456 elmt->GetTracePhysVals(boundary,Bc,U,ubc);
461 Bvals = BndExp[n]->UpdateCoeffs()+BndExp[n]->GetCoeff_Offset(i);
463 Bc->IProductWRTBase(ubc,Bvals);
465 cnt1 += BndExp[n]->GetTotPoints();
469 cnt += BndExp[n]->GetExpSize();
479 static bool Setup =
false;
504 int elmtid,nq, boundary;
511 for(cnt = n = 0; n < BndConds[0].num_elements(); ++n)
513 if((BndConds[0][n]->GetBoundaryConditionType() ==
SpatialDomains::eDirichlet)&& (boost::iequals(BndConds[0][n]->GetUserDefined(),
"ZeroNormalComponent")))
515 for(i = 0; i < BndExp[0][n]->GetExpSize(); ++i,cnt++)
521 normals = elmt->GetSurfaceNormal(boundary);
523 nq = BndExp[0][n]->GetExp(i)->GetTotPoints();
528 Bphys = BndExp[k][n]->UpdatePhys()+
529 BndExp[k][n]->GetPhys_Offset(i);
530 Bc = BndExp[k][n]->GetExp(i);
540 Bphys = BndExp[k][n]->UpdatePhys()+
541 BndExp[k][n]->GetPhys_Offset(i);
542 Bcoeffs = BndExp[k][n]->UpdateCoeffs()+
543 BndExp[k][n]->GetCoeff_Offset(i);
544 Bc = BndExp[k][n]->GetExp(i);
547 Bc->FwdTrans_BndConstrained(Bphys,Bcoeffs);
553 cnt += BndExp[0][n]->GetExpSize();
567 std::complex<NekDouble> za, zar, zJ0, zJ0r, zq, zvel, zJ0rJ0;
570 int M = WomParam->m_wom_vel_r.size();
582 std::complex<NekDouble> z1 (1.0,0.0);
583 std::complex<NekDouble> zi (0.0,1.0);
584 std::complex<NekDouble> comp_conj (-1.0,1.0);
588 BndExp =
m_fields[fldid]->GetBndCondExpansions();
593 int elmtid,offset, boundary,nfq;
598 for(i = 0; i < BndExp[bndid]->GetExpSize(); ++i,cnt++)
602 elmt =
m_fields[fldid]->GetExp(elmtid);
603 offset =
m_fields[fldid]->GetPhys_Offset(elmtid);
606 bc = BndExp[bndid]->GetExp(i);
609 nfq=bc->GetTotPoints();
610 w =
m_fields[fldid]->UpdatePhys() + offset;
616 bc->GetCoords(x,y,z);
619 elmt->GetTracePhysVals(boundary,bc,w,wbc);
626 r = sqrt((x[j]-x0[0])*(x[j]-x0[0]) +
627 (y[j]-x0[1])*(y[j]-x0[1]) +
628 (z[j]-x0[2])*(z[j]-x0[2]))/R;
630 wbc[j] = WomParam->m_wom_vel_r[0]*(1. - r*r);
634 kt = 2.0 * M_PI * k *
m_time / T;
635 za = alpha * sqrt((
NekDouble)k/2.0) * comp_conj;
640 zq = std::exp(zi * kt) * std::complex<NekDouble>(
641 WomParam->m_wom_vel_r[k],
642 WomParam->m_wom_vel_i[k]);
643 zvel = zq * (z1 - zJ0rJ0);
644 wbc[j] = wbc[j] + zvel.real();
651 Bvals = BndExp[bndid]->UpdateCoeffs()+
652 BndExp[bndid]->GetCoeff_Offset(i);
654 bc->FwdTrans(wbc,Bvals);
661 std::string::size_type indxBeg = womStr.find_first_of(
':') + 1;
662 string filename = womStr.substr(indxBeg,string::npos);
664 std::complex<NekDouble> coef;
667 TiXmlDocument doc(filename);
669 bool loadOkay = doc.LoadFile();
670 ASSERTL0(loadOkay,(std::string(
"Failed to load file: ") +
673 TiXmlHandle docHandle(&doc);
677 TiXmlElement *nektar = doc.FirstChildElement(
"NEKTAR");
678 ASSERTL0(nektar,
"Unable to find NEKTAR tag in file.");
680 TiXmlElement *wombc = nektar->FirstChildElement(
"WOMERSLEYBC");
681 ASSERTL0(wombc,
"Unable to find WOMERSLEYBC tag in file.");
684 TiXmlElement *womparam = wombc->FirstChildElement(
"WOMPARAMS");
685 ASSERTL0(womparam,
"Unable to find WOMPARAMS tag in file.");
688 TiXmlElement *params = womparam->FirstChildElement(
"W");
689 map<std::string,std::string> Wparams;
696 propstr = params->Attribute(
"PROPERTY");
698 ASSERTL0(!propstr.empty(),
"Failed to read PROPERTY value Womersley BC Parameter");
702 valstr = params->Attribute(
"VALUE");
704 ASSERTL0(!valstr.empty(),
"Failed to read VALUE value Womersley BC Parameter");
706 std::transform(propstr.begin(),propstr.end(),propstr.begin(),
708 Wparams[propstr] = valstr;
710 params = params->NextSiblingElement(
"W");
715 ASSERTL0(Wparams.count(
"RADIUS") == 1,
716 "Failed to find Radius parameter in Womersley boundary conditions");
717 std::vector<NekDouble>
rad;
719 Wparams[
"RADIUS"].c_str(),rad);
722 ASSERTL0(Wparams.count(
"PERIOD") == 1,
723 "Failed to find period parameter in Womersley boundary conditions");
724 std::vector<NekDouble> period;
726 Wparams[
"PERIOD"].c_str(),period);
730 ASSERTL0(Wparams.count(
"AXISNORMAL") == 1,
731 "Failed to find axisnormal parameter in Womersley boundary conditions");
732 std::vector<NekDouble> anorm;
734 Wparams[
"AXISNORMAL"].c_str(),anorm);
740 ASSERTL0(Wparams.count(
"AXISPOINT") == 1,
741 "Failed to find axispoint parameter in Womersley boundary conditions");
742 std::vector<NekDouble> apt;
744 Wparams[
"AXISPOINT"].c_str(),apt);
752 TiXmlElement *coeff = wombc->FirstChildElement(
"FOURIERCOEFFS");
755 TiXmlElement *fval = coeff->FirstChildElement(
"F");
758 int nextFourierCoeff = -1;
764 TiXmlAttribute *fvalAttr = fval->FirstAttribute();
765 std::string attrName(fvalAttr->Name());
767 ASSERTL0(attrName ==
"ID", (std::string(
"Unknown attribute name: ") + attrName).c_str());
769 err = fvalAttr->QueryIntValue(&indx);
770 ASSERTL0(err == TIXML_SUCCESS,
"Unable to read attribute ID.");
772 std::string coeffStr = fval->FirstChild()->ToText()->ValueStr();
773 vector<NekDouble> coeffvals;
776 ASSERTL0(parseGood,(std::string(
"Problem reading value of fourier coefficient, ID=") + boost::lexical_cast<string>(indx)).c_str());
777 ASSERTL1(coeffvals.size() == 2,(std::string(
"Have not read two entries of Fourier coefficicent from ID="+ boost::lexical_cast<string>(indx)).c_str()));
781 fval = fval->NextSiblingElement(
"F");
785 std::ifstream file(filename);
788 ASSERTL1(file.is_open(),(std::string(
"Missing file ") + filename).c_str());
790 while(std::getline(file,line))
792 std::stringstream
stream(line);
820 bool returnval =
false;
825 int ncoeffs =
m_fields[0]->GetNcoeffs();
827 for(
int i = 0; i <
m_fields.num_elements(); ++i)
848 int n_element =
m_fields[0]->GetExpSize();
863 for(
int i = 0; i < 2; ++i)
872 for(
int i = 0; i < n_vel; ++i)
880 for(
int el = 0; el < n_element; ++el)
882 cfl[el] =
m_timestep*(stdVelocity[el] * cLambda *
883 (ExpOrder[el]-1) * (ExpOrder[el]-1));
894 int n_element =
m_fields[0]->GetExpSize();
901 CFL = CFL_loc = cfl[elmtid];
905 elmtid =
m_fields[0]->GetExp(elmtid)->GetGeom()->GetGlobalID();
916 elmtid = elmtid%
m_fields[0]->GetPlane(0)->GetExpSize();
944 if(
m_comm->GetRank() == 0)
946 cout <<
"CFL (zero plane): "<< cfl <<
" (in elmt "
947 << elmtid <<
")" << endl;
955 cout <<
"Reached Steady State to tolerance "
EquationType m_equationType
equation type;
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)
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
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.
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.
void EvaluateAdvectionTerms(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
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.
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
void SetUpWomersley(const int bndid, std::string womstr)
Set Up Womersley details.
int m_cflsteps
dump cfl estimate
std::map< int, WomersleyParamsSharedPtr > m_womersleyParams
Womersley parameters if required.
std::vector< SolverUtils::ForcingSharedPtr > m_forcing
Forcing terms.
Base class for unsteady solvers.
boost::shared_ptr< WomersleyParams > WomersleyParamsSharedPtr
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)
static NekDouble rad(NekDouble x, NekDouble y)
SOLVER_UTILS_EXPORT int GetNpoints()
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()
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...
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.
void SetWomersleyBoundary(const int fldid, const int bndid)
Set Womersley Profile if specified.
virtual SOLVER_UTILS_EXPORT void v_InitObject()
Init object for UnsteadySystem class.
static bool GenerateUnOrderedVector(const char *const str, std::vector< NekDouble > &vec)
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)