Virtual function for solve implementation.
87{
90
91 m_equ[0]->PrintSummary(out);
92
93
95 m_equ[0]->DoInitialise();
96
97
100
101
102 bool isHomogeneous1D;
103 int nRuns, minP, maxP, sensorVar;
105 m_session->LoadParameter(
"NumRuns", nRuns, 1);
106 m_session->LoadParameter(
"AdaptiveMinModes", minP, 4);
107 m_session->LoadParameter(
"AdaptiveMaxModes", maxP, 12);
108 m_session->LoadParameter(
"AdaptiveLowerTolerance", lowerTol, 1e-8);
109 m_session->LoadParameter(
"AdaptiveUpperTolerance", upperTol, 1e-6);
110 m_session->LoadParameter(
"AdaptiveSensorVariable", sensorVar, 0);
111 m_session->MatchSolverInfo(
"Homogeneous",
"1D", isHomogeneous1D,
false);
112
113
114 int nExp, nPlanes;
115 if (isHomogeneous1D)
116 {
117 nExp =
m_equ[0]->UpdateFields()[0]->GetPlane(0)->GetExpSize();
118 nPlanes =
m_equ[0]->UpdateFields()[0]->GetZIDs().size();
119 }
120 else
121 {
122 nExp =
m_equ[0]->UpdateFields()[0]->GetExpSize();
123 nPlanes = 1;
124 }
125 int expdim =
m_equ[0]->UpdateFields()[0]->GetGraph()->GetMeshDimension();
126
127 int nFields =
m_equ[0]->UpdateFields().size();
128 int numSteps =
m_session->GetParameter(
"NumSteps");
130
131 Array<OneD, NekDouble> coeffs, phys, physReduced, tmpArray;
132
133
136
137
138 Array<OneD, int>
P(expdim);
139 Array<OneD, int> numPoints(expdim);
140 Array<OneD, LibUtilities::PointsKey> ptsKey(expdim);
142 for (int i = 1; i < nRuns; i++)
143 {
144
145 Array<OneD, MultiRegions::ExpListSharedPtr> fields =
146 m_equ[0]->UpdateFields();
147
148
149 map<int, int> deltaP;
150 int offset = 0;
151 for (int n = 0; n < nExp; n++)
152 {
153 offset = fields[sensorVar]->GetPhys_Offset(n);
154 Exp = fields[sensorVar]->GetExp(n);
155
156 for (int k = 0; k < expdim; ++k)
157 {
158 P[k] = Exp->GetBasis(k)->GetNumModes();
159 numPoints[k] = Exp->GetBasis(k)->GetNumPoints();
160 ptsKey[k] = LibUtilities::PointsKey(
161 numPoints[k], Exp->GetBasis(k)->GetPointsType());
162 }
163
164
166 switch (Exp->GetGeom()->GetShapeType())
167 {
169 {
171 ptsKey[0]);
173 ptsKey[1]);
174 OrthoExp = MemoryManager<
175 StdRegions::StdQuadExp>::AllocateSharedPtr(Ba, Bb);
176 break;
177 }
179 {
181 ptsKey[0]);
183 ptsKey[1]);
184 OrthoExp =
186 Ba, Bb);
187 break;
188 }
190 {
192 ptsKey[0]);
194 ptsKey[1]);
196 ptsKey[2]);
197 OrthoExp =
199 Ba, Bb, Bc);
200 break;
201 }
203 {
205 ptsKey[0]);
207 ptsKey[1]);
209 ptsKey[2]);
210 OrthoExp =
212 Ba, Bb, Bc);
213 break;
214 }
216 {
218 ptsKey[0]);
220 ptsKey[1]);
222 ptsKey[2]);
223 OrthoExp = MemoryManager<
224 StdRegions::StdPrismExp>::AllocateSharedPtr(Ba, Bb, Bc);
225 break;
226 }
228 {
230 ptsKey[0]);
232 ptsKey[1]);
234 ptsKey[2]);
235 OrthoExp =
237 Ba, Bb, Bc);
238 break;
239 }
240 default:
241 ASSERTL0(
false,
"Shape not supported.");
242 break;
243 }
244
245 int nq = OrthoExp->GetTotPoints();
246
250
251 coeffs = Array<OneD, NekDouble>(OrthoExp->GetNcoeffs());
252 physReduced = Array<OneD, NekDouble>(OrthoExp->GetTotPoints());
253 tmpArray = Array<OneD, NekDouble>(OrthoExp->GetTotPoints(), 0.0);
254
255
256 for (int plane = 0; plane < nPlanes; plane++)
257 {
258 if (isHomogeneous1D)
259 {
260 phys =
261 fields[sensorVar]->GetPlane(plane)->GetPhys() + offset;
262 }
263 else
264 {
265 phys = fields[sensorVar]->GetPhys() + offset;
266 }
267
268
269 OrthoExp->FwdTrans(phys, coeffs);
270 OrthoExp->BwdTrans(coeffs, physReduced);
271
272
273 Vmath::Vsub(nq, phys, 1, physReduced, 1, tmpArray, 1);
274 Vmath::Vmul(nq, tmpArray, 1, tmpArray, 1, tmpArray, 1);
275 tmp = Exp->Integral(tmpArray);
276
278 fac = Exp->Integral(tmpArray);
279
280 tmp =
abs(tmp / fac);
281
282 if (tmp != tmp)
283 {
285 "Adaptive procedure encountered NaN value.");
286 }
287
288
289 error = (tmp > error) ? tmp : error;
290 }
291
292
293 m_session->GetComm()->GetColumnComm()->AllReduce(
295
296
297 if (
m_session->DefinesFunction(
"AdaptiveLowerTolerance"))
298 {
299 int nq = Exp->GetTotPoints();
300
301
302 Array<OneD, NekDouble> xc0, xc1, xc2;
303 xc0 = Array<OneD, NekDouble>(nq, 0.0);
304 xc1 = Array<OneD, NekDouble>(nq, 0.0);
305 xc2 = Array<OneD, NekDouble>(nq, 0.0);
306 Exp->GetCoords(xc0, xc1, xc2);
307
308
309 Array<OneD, NekDouble> tolerance(nq, 0.0);
311 m_session->GetFunction(
"AdaptiveLowerTolerance", 0);
312 ffunc->Evaluate(xc0, xc1, xc2, tolerance);
314 }
315
316 if (
m_session->DefinesFunction(
"AdaptiveUpperTolerance"))
317 {
318 int nq = Exp->GetTotPoints();
319
320
321 Array<OneD, NekDouble> xc0, xc1, xc2;
322 xc0 = Array<OneD, NekDouble>(nq, 0.0);
323 xc1 = Array<OneD, NekDouble>(nq, 0.0);
324 xc2 = Array<OneD, NekDouble>(nq, 0.0);
325 Exp->GetCoords(xc0, xc1, xc2);
326
327
328 Array<OneD, NekDouble> tolerance(nq, 0.0);
330 m_session->GetFunction(
"AdaptiveUpperTolerance", 0);
331 ffunc->Evaluate(xc0, xc1, xc2, tolerance);
333 }
334
335
336 if ((error > upperTol) && (
P[0] < maxP))
337 {
338 deltaP[Exp->GetGeom()->GetGlobalID()] = 1;
339 }
340 else if ((error < lowerTol) &&
P[0] > minP)
341 {
342 deltaP[Exp->GetGeom()->GetGlobalID()] = -1;
343 }
344 else
345 {
346 deltaP[Exp->GetGeom()->GetGlobalID()] = 0;
347 }
348 }
349
350
353
354
355
356
357
358 if (LibUtilities::NekManager<MultiRegions::GlobalLinSysKey,
359 MultiRegions::GlobalLinSys>::
360 PoolCreated(std::string("GlobalLinSys")))
361 {
362 LibUtilities::NekManager<MultiRegions::GlobalLinSysKey,
363 MultiRegions::GlobalLinSys>::
364 ClearManager(std::string("GlobalLinSys"));
365 }
366
367 int chkNumber =
m_equ[0]->GetCheckpointNumber();
368 int chkSteps =
m_equ[0]->GetCheckpointSteps();
369
370
372
373
374 mapping->ReplaceField(
m_equ[0]->UpdateFields());
375
376
377 m_equ[0]->SetCheckpointSteps(0);
378
379
380 m_equ[0]->DoInitialise();
381 m_equ[0]->SetInitialStep(i * numSteps);
382 m_equ[0]->SetSteps(i * numSteps + numSteps);
383 m_equ[0]->SetTime(startTime + i * period);
384 m_equ[0]->SetBoundaryConditions(startTime + i * period);
385 m_equ[0]->SetCheckpointNumber(chkNumber);
386 m_equ[0]->SetCheckpointSteps(chkSteps);
387
388
389 for (int n = 0; n < nFields; n++)
390 {
391 m_equ[0]->UpdateFields()[n]->ExtractCoeffsToCoeffs(
392 fields[n], fields[n]->GetCoeffs(),
393 m_equ[0]->UpdateFields()[n]->UpdateCoeffs());
394 m_equ[0]->UpdateFields()[n]->BwdTrans(
395 m_equ[0]->UpdateFields()[n]->GetCoeffs(),
396 m_equ[0]->UpdateFields()[n]->UpdatePhys());
397 }
398
399
401 }
402
404
406
407 if (
m_comm->GetRank() == 0)
408 {
409 CPUtime = timer.
Elapsed().count();
410 cout << "-------------------------------------------" << endl;
411 cout << "Total Computation Time = " << CPUtime << "s" << endl;
412 cout << "-------------------------------------------" << endl;
413 }
414
415
416
417
418
419
420 for (
int i = 0; i <
m_equ[0]->GetNvariables(); ++i)
421 {
422 Array<OneD, NekDouble> exactsoln(
m_equ[0]->GetTotPoints(), 0.0);
423
424
425 m_equ[0]->EvaluateExactSolution(i, exactsoln,
m_equ[0]->GetTime());
426
429
430 if (
m_comm->GetRank() == 0)
431 {
432 out <<
"L 2 error (variable " <<
m_equ[0]->GetVariable(i)
433 << ") : " << vL2Error << endl;
434 out <<
"L inf error (variable " <<
m_equ[0]->GetVariable(i)
435 << ") : " << vLinfError << endl;
436 }
437 }
438}
#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.
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.
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.
GLOBAL_MAPPING_EXPORT typedef std::shared_ptr< Mapping > MappingSharedPtr
A shared pointer to a Mapping object.
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
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)