83 m_equ[0]->PrintSummary(out);
87 m_equ[0]->DoInitialise();
95 int nRuns, minP, maxP, sensorVar;
97 m_session->LoadParameter(
"NumRuns", nRuns, 1);
98 m_session->LoadParameter(
"AdaptiveMinModes", minP, 4);
99 m_session->LoadParameter(
"AdaptiveMaxModes", maxP, 12);
100 m_session->LoadParameter(
"AdaptiveLowerTolerance", lowerTol, 1e-8);
101 m_session->LoadParameter(
"AdaptiveUpperTolerance", upperTol, 1e-6);
102 m_session->LoadParameter(
"AdaptiveSensorVariable", sensorVar, 0);
103 m_session->MatchSolverInfo(
"Homogeneous",
"1D", isHomogeneous1D,
false);
109 nExp =
m_equ[0]->UpdateFields()[0]->GetPlane(0)->GetExpSize();
110 nPlanes =
m_equ[0]->UpdateFields()[0]->GetZIDs().size();
114 nExp =
m_equ[0]->UpdateFields()[0]->GetExpSize();
117 int expdim =
m_equ[0]->UpdateFields()[0]->GetGraph()->GetMeshDimension();
119 int nFields =
m_equ[0]->UpdateFields().size();
120 int numSteps =
m_session->GetParameter(
"NumSteps");
134 for (
int i = 1; i < nRuns; i++)
138 m_equ[0]->UpdateFields();
141 std::map<int, int> deltaP;
143 for (
int n = 0; n < nExp; n++)
145 offset = fields[sensorVar]->GetPhys_Offset(n);
146 Exp = fields[sensorVar]->GetExp(n);
148 for (
int k = 0; k < expdim; ++k)
150 P[k] = Exp->GetBasis(k)->GetNumModes();
151 numPoints[k] = Exp->GetBasis(k)->GetNumPoints();
153 numPoints[k], Exp->GetBasis(k)->GetPointsType());
158 switch (Exp->GetGeom()->GetShapeType())
233 ASSERTL0(
false,
"Shape not supported.");
237 int nq = OrthoExp->GetTotPoints();
248 for (
int plane = 0; plane < nPlanes; plane++)
253 fields[sensorVar]->GetPlane(plane)->GetPhys() + offset;
257 phys = fields[sensorVar]->GetPhys() + offset;
261 OrthoExp->FwdTrans(phys, coeffs);
262 OrthoExp->BwdTrans(coeffs, physReduced);
265 Vmath::Vsub(nq, phys, 1, physReduced, 1, tmpArray, 1);
266 Vmath::Vmul(nq, tmpArray, 1, tmpArray, 1, tmpArray, 1);
267 tmp = Exp->Integral(tmpArray);
270 fac = Exp->Integral(tmpArray);
272 tmp =
abs(tmp / fac);
277 "Adaptive procedure encountered NaN value.");
281 error = (tmp > error) ? tmp : error;
285 m_session->GetComm()->GetColumnComm()->AllReduce(
289 if (
m_session->DefinesFunction(
"AdaptiveLowerTolerance"))
291 int nq = Exp->GetTotPoints();
298 Exp->GetCoords(xc0, xc1, xc2);
303 m_session->GetFunction(
"AdaptiveLowerTolerance", 0);
304 ffunc->Evaluate(xc0, xc1, xc2, tolerance);
308 if (
m_session->DefinesFunction(
"AdaptiveUpperTolerance"))
310 int nq = Exp->GetTotPoints();
317 Exp->GetCoords(xc0, xc1, xc2);
322 m_session->GetFunction(
"AdaptiveUpperTolerance", 0);
323 ffunc->Evaluate(xc0, xc1, xc2, tolerance);
328 if ((error > upperTol) && (
P[0] < maxP))
330 deltaP[Exp->GetGeom()->GetGlobalID()] = 1;
332 else if ((error < lowerTol) &&
P[0] > minP)
334 deltaP[Exp->GetGeom()->GetGlobalID()] = -1;
338 deltaP[Exp->GetGeom()->GetGlobalID()] = 0;
352 PoolCreated(std::string(
"GlobalLinSys")))
356 ClearManager(std::string(
"GlobalLinSys"));
359 int chkNumber =
m_equ[0]->GetCheckpointNumber();
360 int chkSteps =
m_equ[0]->GetCheckpointSteps();
366 mapping->ReplaceField(
m_equ[0]->UpdateFields());
369 m_equ[0]->SetCheckpointSteps(0);
372 m_equ[0]->DoInitialise();
373 m_equ[0]->SetInitialStep(i * numSteps);
374 m_equ[0]->SetSteps(i * numSteps + numSteps);
375 m_equ[0]->SetTime(startTime + i * period);
376 m_equ[0]->SetBoundaryConditions(startTime + i * period);
377 m_equ[0]->SetCheckpointNumber(chkNumber);
378 m_equ[0]->SetCheckpointSteps(chkSteps);
381 for (
int n = 0; n < nFields; n++)
383 m_equ[0]->UpdateFields()[n]->ExtractCoeffsToCoeffs(
384 fields[n], fields[n]->GetCoeffs(),
385 m_equ[0]->UpdateFields()[n]->UpdateCoeffs());
386 m_equ[0]->UpdateFields()[n]->BwdTrans(
387 m_equ[0]->UpdateFields()[n]->GetCoeffs(),
388 m_equ[0]->UpdateFields()[n]->UpdatePhys());
399 if (
m_comm->GetRank() == 0)
401 CPUtime = timer.
Elapsed().count();
402 std::cout <<
"-------------------------------------------" << std::endl;
403 std::cout <<
"Total Computation Time = " << CPUtime <<
"s" << std::endl;
404 std::cout <<
"-------------------------------------------" << std::endl;
412 for (
int i = 0; i <
m_equ[0]->GetNvariables(); ++i)
417 m_equ[0]->EvaluateExactSolution(i, exactsoln,
m_equ[0]->GetTime());
422 if (
m_comm->GetRank() == 0)
424 out <<
"L 2 error (variable " <<
m_equ[0]->GetVariable(i)
425 <<
") : " << vL2Error << std::endl;
426 out <<
"L inf error (variable " <<
m_equ[0]->GetVariable(i)
427 <<
") : " << vLinfError << std::endl;
441 std::map<int, int> deltaP)
444 int expdim =
m_equ[0]->UpdateFields()[0]->GetGraph()->GetMeshDimension();
447 std::vector<LibUtilities::FieldDefinitionsSharedPtr> fielddefs =
448 fields[0]->GetFieldDefinitions();
450 if (fielddefs[0]->m_numHomogeneousDir == 1)
460 for (
int i = 0; i < fielddefs.size(); ++i)
462 for (
int j = 0; j < fields.size(); ++j)
464 fielddefs[i]->m_fields.push_back(
m_session->GetVariable(j));
469 TiXmlElement *exp_tag =
m_session->GetElement(
"NEKTAR/EXPANSIONS");
475 for (
int f = 0; f < fielddefs.size(); ++f)
477 nExp = fielddefs[f]->m_elementIDs.size();
481 TiXmlElement *elemTag =
new TiXmlElement(
"ELEMENTS");
482 exp_tag->LinkEndChild(elemTag);
485 std::string fieldsString;
487 std::stringstream fieldsStringStream;
489 for (std::vector<int>::size_type i = 0;
490 i < fielddefs[f]->m_fields.size(); i++)
494 fieldsStringStream <<
",";
496 fieldsStringStream << fielddefs[f]->m_fields[i];
499 fieldsString = fieldsStringStream.str();
501 elemTag->SetAttribute(
"FIELDS", fieldsString);
504 std::string shapeString;
506 std::stringstream shapeStringStream;
509 if (fielddefs[f]->m_numHomogeneousDir == 1)
511 shapeStringStream <<
"-HomogenousExp1D";
513 else if (fielddefs[f]->m_numHomogeneousDir == 2)
515 shapeStringStream <<
"-HomogenousExp2D";
518 shapeString = shapeStringStream.str();
520 elemTag->SetAttribute(
"SHAPE", shapeString);
523 std::string basisString;
525 std::stringstream basisStringStream;
527 for (std::vector<LibUtilities::BasisType>::size_type i = 0;
528 i < fielddefs[f]->m_basis.size(); i++)
532 basisStringStream <<
",";
538 basisString = basisStringStream.str();
540 elemTag->SetAttribute(
"BASIS", basisString);
543 if (fielddefs[f]->m_numHomogeneousDir)
545 std::string homoLenString;
547 std::stringstream homoLenStringStream;
549 for (
int i = 0; i < fielddefs[f]->m_numHomogeneousDir; ++i)
553 homoLenStringStream <<
",";
556 << fielddefs[f]->m_homogeneousLengths[i];
559 homoLenString = homoLenStringStream.str();
561 elemTag->SetAttribute(
"HOMOGENEOUSLENGTHS", homoLenString);
565 if (fielddefs[f]->m_numHomogeneousDir)
567 if (fielddefs[f]->m_homogeneousYIDs.size() > 0)
569 std::string homoYIDsString;
571 std::stringstream homoYIDsStringStream;
573 for (
int i = 0; i < fielddefs[f]->m_homogeneousYIDs.size();
578 homoYIDsStringStream <<
",";
581 << fielddefs[f]->m_homogeneousYIDs[i];
584 homoYIDsString = homoYIDsStringStream.str();
586 elemTag->SetAttribute(
"HOMOGENEOUSYIDS", homoYIDsString);
589 if (fielddefs[f]->m_homogeneousZIDs.size() > 0)
591 std::string homoZIDsString;
593 std::stringstream homoZIDsStringStream;
595 for (
int i = 0; i < fielddefs[f]->m_homogeneousZIDs.size();
600 homoZIDsStringStream <<
",";
603 << fielddefs[f]->m_homogeneousZIDs[i];
606 homoZIDsString = homoZIDsStringStream.str();
608 elemTag->SetAttribute(
"HOMOGENEOUSZIDS", homoZIDsString);
613 std::string numModesString;
615 std::stringstream numModesStringStream;
617 numModesStringStream <<
"MIXORDER:";
621 for (
int n = 0; n < nExp; n++)
623 eId = fielddefs[f]->m_elementIDs[n];
625 for (
int i = 0; i < expdim; i++)
627 order[i] = deltaP[eId];
630 for (
int i = 0; i < nDim; i++)
632 if (fielddefs[f]->m_uniOrder)
634 order[i] += fielddefs[f]->m_numModes[i];
638 order[i] += fielddefs[f]->m_numModes[n * nDim + i];
643 numModesStringStream <<
",";
646 numModesStringStream << order[i];
650 numModesString = numModesStringStream.str();
652 elemTag->SetAttribute(
"NUMMODESPERDIR", numModesString);
656 elemTag->SetAttribute(
#define ASSERTL0(condition, msg)
static GLOBAL_MAPPING_EXPORT MappingSharedPtr Load(const LibUtilities::SessionReaderSharedPtr &pSession, const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields)
Return a pointer to the mapping, creating it on first call.
Describes the specification for a Basis.
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
Defines a specification for a set of points.
static std::string RegisterEnumValue(std::string pEnum, std::string pString, int pEnumValue)
Registers an enumeration value.
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
Describe a linear system.
static std::string GenerateSeqString(const std::vector< T > &v)
Generate a compressed comma-separated string representation of a vector of unsigned integers.
static std::string driverLookupId
static DriverSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
Creates an instance of this class.
SOLVER_UTILS_EXPORT void ReplaceExpansion(Array< OneD, MultiRegions::ExpListSharedPtr > &fields, std::map< int, int > deltaP)
Update EXPANSIONS tag inside XML schema to reflect new polynomial order distribution.
SOLVER_UTILS_EXPORT DriverAdaptive(const LibUtilities::SessionReaderSharedPtr pSession, const SpatialDomains::MeshGraphSharedPtr pGraph)
Constructor.
SOLVER_UTILS_EXPORT void v_InitObject(std::ostream &out=std::cout) override
Virtual function for initialisation implementation.
SOLVER_UTILS_EXPORT void v_Execute(std::ostream &out=std::cout) override
Virtual function for solve implementation.
static std::string className
Name of the class.
Base class for the development of solvers.
LibUtilities::SessionReaderSharedPtr m_session
Session reader object.
virtual SOLVER_UTILS_EXPORT void v_InitObject(std::ostream &out=std::cout)
Virtual function for initialisation implementation.
LibUtilities::CommSharedPtr m_comm
Communication object.
SpatialDomains::MeshGraphSharedPtr m_graph
MeshGraph object.
Array< OneD, EquationSystemSharedPtr > m_equ
Equation system to solve.
Class representing a prismatic element in reference space.
GLOBAL_MAPPING_EXPORT typedef std::shared_ptr< Mapping > MappingSharedPtr
A shared pointer to a Mapping object.
const char *const ShapeTypeMap[SIZE_ShapeType]
const char *const BasisTypeMap[]
std::shared_ptr< SessionReader > SessionReaderSharedPtr
std::shared_ptr< Equation > EquationSharedPtr
@ P
Monomial polynomials .
@ eOrtho_A
Principle Orthogonal Functions .
@ eOrtho_C
Principle Orthogonal Functions .
@ eOrtho_B
Principle Orthogonal Functions .
std::shared_ptr< Expansion > ExpansionSharedPtr
DriverFactory & GetDriverFactory()
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
std::shared_ptr< StdExpansion > StdExpansionSharedPtr
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.
T Vsum(int n, const T *x, const int incx)
Subtract return sum(x)
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 > abs(scalarT< T > in)