Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
DriverAdaptive.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File DriverAdaptive.cpp
4 //
5 // For more information, please see: http://www.nektar.info
6 //
7 // The MIT License
8 //
9 // Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10 // Department of Aeronautics, Imperial College London (UK), and Scientific
11 // Computing and Imaging Institute, University of Utah (USA).
12 //
13 // License for the specific language governing rights and limitations under
14 // Permission is hereby granted, free of charge, to any person obtaining a
15 // copy of this software and associated documentation files (the "Software"),
16 // to deal in the Software without restriction, including without limitation
17 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
18 // and/or sell copies of the Software, and to permit persons to whom the
19 // Software is furnished to do so, subject to the following conditions:
20 //
21 // The above copyright notice and this permission notice shall be included
22 // in all copies or substantial portions of the Software.
23 //
24 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
25 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
26 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
27 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
28 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
29 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
30 // DEALINGS IN THE SOFTWARE.
31 //
32 // Description: Driver class for adaptive solver
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
36 #include <iomanip>
37 
39 #include <StdRegions/StdQuadExp.h>
40 #include <StdRegions/StdTriExp.h>
41 #include <GlobalMapping/Mapping.h>
42 
43 using namespace std;
44 
45 namespace Nektar
46 {
47 namespace SolverUtils
48 {
49 
50 string DriverAdaptive::className = GetDriverFactory().RegisterCreatorFunction(
51  "Adaptive", DriverAdaptive::create);
52 string DriverAdaptive::driverLookupId =
53  LibUtilities::SessionReader::RegisterEnumValue("Driver", "Adaptive", 0);
54 
55 /**
56  *
57  */
58 DriverAdaptive::DriverAdaptive(
60  : Driver(pSession)
61 {
62 }
63 
64 /**
65  *
66  */
68 {
69 }
70 
71 /**
72  *
73  */
74 void DriverAdaptive::v_InitObject(ostream &out)
75 {
77 }
78 
79 void DriverAdaptive::v_Execute(ostream &out)
80 {
81  time_t starttime, endtime;
82  NekDouble CPUtime;
83 
84  m_equ[0]->PrintSummary(out);
85 
86  // First run using original order
87  time(&starttime);
88  m_equ[0]->DoInitialise();
89 
90  // Obtain initial time in case a restart was used
91  NekDouble startTime = m_equ[0]->GetFinalTime();
92  m_equ[0]->DoSolve();
93 
94  // Load session parameters and solver info
95  bool isHomogeneous1D;
96  int nRuns, minP, maxP, sensorVar;
97  NekDouble lowerTol, upperTol;
98  m_session->LoadParameter ("NumRuns", nRuns, 1);
99  m_session->LoadParameter ("AdaptiveMinModes", minP, 4);
100  m_session->LoadParameter ("AdaptiveMaxModes", maxP, 12);
101  m_session->LoadParameter ("AdaptiveLowerTolerance", lowerTol, 1e-8);
102  m_session->LoadParameter ("AdaptiveUpperTolerance", upperTol, 1e-6);
103  m_session->LoadParameter ("AdaptiveSensorVariable", sensorVar, 0);
104  m_session->MatchSolverInfo("Homogeneous", "1D", isHomogeneous1D, false);
105 
106  // Get number of elements and planes
107  int nExp, nPlanes;
108  if (isHomogeneous1D)
109  {
110  nExp = m_equ[0]->UpdateFields()[0]->GetPlane(0)->GetExpSize();
111  nPlanes = m_equ[0]->UpdateFields()[0]->GetZIDs().num_elements();
112  }
113  else
114  {
115  nExp = m_equ[0]->UpdateFields()[0]->GetExpSize();
116  nPlanes = 1;
117  }
118 
119  int nFields = m_equ[0]->UpdateFields().num_elements();
120  int numSteps = m_session->GetParameter("NumSteps");
121  NekDouble period = m_session->GetParameter("TimeStep") * numSteps;
122 
123  Array<OneD, NekDouble> coeffs, phys, physReduced, tmpArray;
124 
125  // Get mapping
128  m_equ[0]->UpdateFields());
129 
130  // Adaptive loop
132  for (int i = 1; i < nRuns; i++)
133  {
134  // Get field expansions
136  m_equ[0]->UpdateFields();
137 
138  // Determine the change to be applied in the order
139  map<int, int> deltaP;
140  int offset = 0;
141  for (int n = 0; n < nExp; n++)
142  {
143  offset = fields[sensorVar]->GetPhys_Offset(n);
144  Exp = fields[sensorVar]->GetExp(n);
145 
146  int P = Exp->GetBasis(0)->GetNumModes();
147  int Q = Exp->GetBasis(1)->GetNumModes();
148  int qa = Exp->GetBasis(0)->GetNumPoints();
149  int qb = Exp->GetBasis(1)->GetNumPoints();
150 
151  // Declare orthogonal basis.
152  LibUtilities::PointsKey pa(qa, Exp->GetBasis(0)->GetPointsType());
153  LibUtilities::PointsKey pb(qb, Exp->GetBasis(1)->GetPointsType());
154 
156  switch (Exp->GetGeom()->GetShapeType())
157  {
159  {
161  pa);
163  pb);
164  OrthoExp = MemoryManager<
165  StdRegions::StdQuadExp>::AllocateSharedPtr(Ba, Bb);
166  break;
167  }
169  {
171  pa);
173  pb);
174  OrthoExp =
176  Ba, Bb);
177  break;
178  }
179  default:
180  ASSERTL0(false, "Shape not supported.");
181  break;
182  }
183 
184  int nq = OrthoExp->GetTotPoints();
185 
186  NekDouble error = 0;
187  NekDouble fac = 0;
188  NekDouble tmp = 0;
189 
190  coeffs = Array<OneD, NekDouble>(OrthoExp->GetNcoeffs());
191  physReduced = Array<OneD, NekDouble>(OrthoExp->GetTotPoints());
192  tmpArray = Array<OneD, NekDouble>(OrthoExp->GetTotPoints(), 0.0);
193 
194  // Refinement based only on one variable
195  for (int plane = 0; plane < nPlanes; plane++)
196  {
197  if (isHomogeneous1D)
198  {
199  phys =
200  fields[sensorVar]->GetPlane(plane)->GetPhys() + offset;
201  }
202  else
203  {
204  phys = fields[sensorVar]->GetPhys() + offset;
205  }
206 
207  // Project solution to lower order
208  OrthoExp->FwdTrans(phys, coeffs);
209  OrthoExp->BwdTrans(coeffs, physReduced);
210 
211  // Calculate error =||phys-physReduced||^2 / ||phys||^2
212  Vmath::Vsub(nq, phys, 1, physReduced, 1, tmpArray, 1);
213  Vmath::Vmul(nq, tmpArray, 1, tmpArray, 1, tmpArray, 1);
214  tmp = Exp->Integral(tmpArray);
215 
216  Vmath::Vmul(nq, phys, 1, phys, 1, tmpArray, 1);
217  fac = Exp->Integral(tmpArray);
218 
219  tmp = abs(tmp / fac);
220 
221  if (tmp != tmp)
222  {
223  ASSERTL0(false,
224  "Adaptive procedure encountered NaN value.");
225  }
226 
227  // Get maximum value along planes
228  error = (tmp > error) ? tmp : error;
229  }
230 
231  // Combine planes from different processes
232  m_session->GetComm()->GetColumnComm()->AllReduce(
233  error, LibUtilities::ReduceMax);
234 
235  // Override tolerances if function is defined
236  if (m_session->DefinesFunction("AdaptiveLowerTolerance"))
237  {
238  int nq = Exp->GetTotPoints();
239 
240  // Obtain points from the the element
241  Array<OneD, NekDouble> xc0, xc1, xc2;
242  xc0 = Array<OneD, NekDouble>(nq, 0.0);
243  xc1 = Array<OneD, NekDouble>(nq, 0.0);
244  xc2 = Array<OneD, NekDouble>(nq, 0.0);
245  Exp->GetCoords(xc0, xc1, xc2);
246 
247  // Evaluate function from session file
248  Array<OneD, NekDouble> tolerance(nq, 0.0);
250  m_session->GetFunction("AdaptiveLowerTolerance", 0);
251  ffunc->Evaluate(xc0, xc1, xc2, tolerance);
252  lowerTol = Vmath::Vsum(nq, tolerance, 1) / nq;
253  }
254 
255  if (m_session->DefinesFunction("AdaptiveUpperTolerance"))
256  {
257  int nq = Exp->GetTotPoints();
258 
259  // Obtain points from the the element
260  Array<OneD, NekDouble> xc0, xc1, xc2;
261  xc0 = Array<OneD, NekDouble>(nq, 0.0);
262  xc1 = Array<OneD, NekDouble>(nq, 0.0);
263  xc2 = Array<OneD, NekDouble>(nq, 0.0);
264  Exp->GetCoords(xc0, xc1, xc2);
265 
266  // Evaluate function from session file
267  Array<OneD, NekDouble> tolerance(nq, 0.0);
269  m_session->GetFunction("AdaptiveUpperTolerance", 0);
270  ffunc->Evaluate(xc0, xc1, xc2, tolerance);
271  upperTol = Vmath::Vsum(nq, tolerance, 1) / nq;
272  }
273 
274  // Determine what to do with the polynomial order
275  if ((error > upperTol) && (P < maxP))
276  {
277  deltaP[Exp->GetGeom()->GetGlobalID()] = 1;
278  }
279  else if ((error < lowerTol) && P > minP)
280  {
281  deltaP[Exp->GetGeom()->GetGlobalID()] = -1;
282  }
283  else
284  {
285  deltaP[Exp->GetGeom()->GetGlobalID()] = 0;
286  }
287  }
288 
289  // Write new expansion section to the session reader
290  ReplaceExpansion(fields, deltaP);
291 
292  // Reset GlobalLinSys Manager to avoid using too much memory
293  //
294  // @todo This could be made better by replacing individual matrices
295  // within the linear system.
298  ClearManager(std::string("GlobalLinSys"));
299 
300  int chkNumber = m_equ[0]->GetCheckpointNumber();
301  int chkSteps = m_equ[0]->GetCheckpointSteps();
302 
303  // Initialise driver again
305 
306  // Update mapping (must be before m_equ[0]->DoInitialise();)
307  mapping->ReplaceField(m_equ[0]->UpdateFields());
308 
309  // Set chkSteps to zero to avoid writing initial condition
310  m_equ[0]->SetCheckpointSteps(0);
311 
312  // Initialise equation
313  m_equ[0]->DoInitialise();
314  m_equ[0]->SetInitialStep(i * numSteps);
315  m_equ[0]->SetSteps(i * numSteps + numSteps);
316  m_equ[0]->SetTime(startTime + i * period);
317  m_equ[0]->SetBoundaryConditions(startTime + i * period);
318  m_equ[0]->SetCheckpointNumber(chkNumber);
319  m_equ[0]->SetCheckpointSteps(chkSteps);
320 
321  // Project solution to new expansion
322  for (int n = 0; n < nFields; n++)
323  {
324  m_equ[0]->UpdateFields()[n]->ExtractCoeffsToCoeffs(
325  fields[n], fields[n]->GetCoeffs(),
326  m_equ[0]->UpdateFields()[n]->UpdateCoeffs());
327  m_equ[0]->UpdateFields()[n]->BwdTrans_IterPerExp(
328  m_equ[0]->UpdateFields()[n]->GetCoeffs(),
329  m_equ[0]->UpdateFields()[n]->UpdatePhys());
330  }
331 
332  // Solve equation
333  m_equ[0]->DoSolve();
334  }
335 
336  time(&endtime);
337 
338  m_equ[0]->Output();
339 
340  if (m_comm->GetRank() == 0)
341  {
342  CPUtime = difftime(endtime, starttime);
343  cout << "-------------------------------------------" << endl;
344  cout << "Total Computation Time = " << CPUtime << "s" << endl;
345  cout << "-------------------------------------------" << endl;
346  }
347 
348  // Evaluate and output computation time and solution accuracy.
349  // The specific format of the error output is essential for the
350  // regression tests to work.
351 
352  // Evaluate L2 Error
353  for (int i = 0; i < m_equ[0]->GetNvariables(); ++i)
354  {
355  Array<OneD, NekDouble> exactsoln(m_equ[0]->GetTotPoints(), 0.0);
356 
357  // Evaluate "ExactSolution" function, or zero array
358  m_equ[0]->EvaluateExactSolution(i, exactsoln, m_equ[0]->GetFinalTime());
359 
360  NekDouble vL2Error = m_equ[0]->L2Error (i, exactsoln);
361  NekDouble vLinfError = m_equ[0]->LinfError(i, exactsoln);
362 
363  if (m_comm->GetRank() == 0)
364  {
365  out << "L 2 error (variable " << m_equ[0]->GetVariable(i)
366  << ") : " << vL2Error << endl;
367  out << "L inf error (variable " << m_equ[0]->GetVariable(i)
368  << ") : " << vLinfError << endl;
369  }
370  }
371 }
372 
373 /**
374  * @brief Update EXPANSIONS tag inside XML schema to reflect new polynomial
375  * order distribution.
376  *
377  * @param fields Input fields.
378  * @param deltaP Map of polynomial order expansions
379  */
382  map<int, int> deltaP)
383 {
384  int nExp, nDim;
385 
386  // Get field definitions
387  std::vector<LibUtilities::FieldDefinitionsSharedPtr> fielddefs =
388  fields[0]->GetFieldDefinitions();
389 
390  int expDim = 2;
391 
392  if (fielddefs[0]->m_numHomogeneousDir == 1)
393  {
394  nDim = 3;
395  }
396  else
397  {
398  nDim = 2;
399  }
400 
401  // Add variables to field definition
402  for (int i = 0; i < fielddefs.size(); ++i)
403  {
404  for (int j = 0; j < fields.num_elements(); ++j)
405  {
406  fielddefs[i]->m_fields.push_back(m_session->GetVariable(j));
407  }
408  }
409 
410  // Get tinyxml objects
411  TiXmlElement *exp_tag = m_session->GetElement("NEKTAR/EXPANSIONS");
412 
413  // Clear current expansions
414  exp_tag->Clear();
415 
416  // Write new expansion information
417  for (int f = 0; f < fielddefs.size(); ++f)
418  {
419  nExp = fielddefs[f]->m_elementIDs.size();
420 
421  //---------------------------------------------
422  // Write ELEMENTS
423  TiXmlElement *elemTag = new TiXmlElement("ELEMENTS");
424  exp_tag->LinkEndChild(elemTag);
425 
426  // Write FIELDS
427  std::string fieldsString;
428  {
429  std::stringstream fieldsStringStream;
430  bool first = true;
431  for (std::vector<int>::size_type i = 0;
432  i < fielddefs[f]->m_fields.size(); i++)
433  {
434  if (!first)
435  {
436  fieldsStringStream << ",";
437  }
438  fieldsStringStream << fielddefs[f]->m_fields[i];
439  first = false;
440  }
441  fieldsString = fieldsStringStream.str();
442  }
443  elemTag->SetAttribute("FIELDS", fieldsString);
444 
445  // Write SHAPE
446  std::string shapeString;
447  {
448  std::stringstream shapeStringStream;
449  shapeStringStream
450  << LibUtilities::ShapeTypeMap[fielddefs[f]->m_shapeType];
451  if (fielddefs[f]->m_numHomogeneousDir == 1)
452  {
453  shapeStringStream << "-HomogenousExp1D";
454  }
455  else if (fielddefs[f]->m_numHomogeneousDir == 2)
456  {
457  shapeStringStream << "-HomogenousExp2D";
458  }
459 
460  shapeString = shapeStringStream.str();
461  }
462  elemTag->SetAttribute("SHAPE", shapeString);
463 
464  // Write BASIS
465  std::string basisString;
466  {
467  std::stringstream basisStringStream;
468  bool first = true;
469  for (std::vector<LibUtilities::BasisType>::size_type i = 0;
470  i < fielddefs[f]->m_basis.size(); i++)
471  {
472  if (!first)
473  {
474  basisStringStream << ",";
475  }
476  basisStringStream
477  << LibUtilities::BasisTypeMap[fielddefs[f]->m_basis[i]];
478  first = false;
479  }
480  basisString = basisStringStream.str();
481  }
482  elemTag->SetAttribute("BASIS", basisString);
483 
484  // Write homogeneuous length details
485  if (fielddefs[f]->m_numHomogeneousDir)
486  {
487  std::string homoLenString;
488  {
489  std::stringstream homoLenStringStream;
490  bool first = true;
491  for (int i = 0; i < fielddefs[f]->m_numHomogeneousDir; ++i)
492  {
493  if (!first)
494  {
495  homoLenStringStream << ",";
496  }
497  homoLenStringStream
498  << fielddefs[f]->m_homogeneousLengths[i];
499  first = false;
500  }
501  homoLenString = homoLenStringStream.str();
502  }
503  elemTag->SetAttribute("HOMOGENEOUSLENGTHS", homoLenString);
504  }
505 
506  // Write homogeneuous planes/lines details
507  if (fielddefs[f]->m_numHomogeneousDir)
508  {
509  if (fielddefs[f]->m_homogeneousYIDs.size() > 0)
510  {
511  std::string homoYIDsString;
512  {
513  std::stringstream homoYIDsStringStream;
514  bool first = true;
515  for (int i = 0; i < fielddefs[f]->m_homogeneousYIDs.size();
516  i++)
517  {
518  if (!first)
519  {
520  homoYIDsStringStream << ",";
521  }
522  homoYIDsStringStream
523  << fielddefs[f]->m_homogeneousYIDs[i];
524  first = false;
525  }
526  homoYIDsString = homoYIDsStringStream.str();
527  }
528  elemTag->SetAttribute("HOMOGENEOUSYIDS", homoYIDsString);
529  }
530 
531  if (fielddefs[f]->m_homogeneousZIDs.size() > 0)
532  {
533  std::string homoZIDsString;
534  {
535  std::stringstream homoZIDsStringStream;
536  bool first = true;
537  for (int i = 0; i < fielddefs[f]->m_homogeneousZIDs.size();
538  i++)
539  {
540  if (!first)
541  {
542  homoZIDsStringStream << ",";
543  }
544  homoZIDsStringStream
545  << fielddefs[f]->m_homogeneousZIDs[i];
546  first = false;
547  }
548  homoZIDsString = homoZIDsStringStream.str();
549  }
550  elemTag->SetAttribute("HOMOGENEOUSZIDS", homoZIDsString);
551  }
552  }
553 
554  // Write NUMMODESPERDIR
555  std::string numModesString;
556  {
557  std::stringstream numModesStringStream;
558 
559  numModesStringStream << "MIXORDER:";
560  bool first = true;
561  int eId;
562  Array<OneD, int> order(nDim, 0);
563  for (int n = 0; n < nExp; n++)
564  {
565  eId = fielddefs[f]->m_elementIDs[n];
566 
567  for (int i = 0; i < expDim; i++)
568  {
569  order[i] = deltaP[eId];
570  }
571 
572  for (int i = 0; i < nDim; i++)
573  {
574  if (fielddefs[f]->m_uniOrder)
575  {
576  order[i] += fielddefs[f]->m_numModes[i];
577  }
578  else
579  {
580  order[i] += fielddefs[f]->m_numModes[n * nDim + i];
581  }
582 
583  if (!first)
584  {
585  numModesStringStream << ",";
586  }
587 
588  numModesStringStream << order[i];
589  first = false;
590  }
591  }
592  numModesString = numModesStringStream.str();
593  }
594  elemTag->SetAttribute("NUMMODESPERDIR", numModesString);
595 
596  // Write IDs. Should ideally look at ways of compressing this stream if
597  // just sequential.
598  elemTag->SetAttribute(
599  "ID", ParseUtils::GenerateSeqString(fielddefs[f]->m_elementIDs));
600  }
601 }
602 
603 }
604 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
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.
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
const char *const BasisTypeMap[]
Definition: Foundations.hpp:47
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
static std::string GenerateSeqString(const std::vector< unsigned int > &elmtids)
Definition: ParseUtils.hpp:159
STL namespace.
virtual SOLVER_UTILS_EXPORT ~DriverAdaptive()
Destructor.
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:51
const char *const ShapeTypeMap[]
Definition: ShapeType.hpp:66
Principle Orthogonal Functions .
Definition: BasisType.h:47
virtual SOLVER_UTILS_EXPORT void v_Execute(std::ostream &out=std::cout)
Virtual function for solve implementation.
virtual SOLVER_UTILS_EXPORT void v_InitObject(std::ostream &out=std::cout)
Definition: Driver.cpp:92
Principle Orthogonal Functions .
Definition: BasisType.h:46
Defines a specification for a set of points.
Definition: Points.h:58
double NekDouble
boost::shared_ptr< Expansion > ExpansionSharedPtr
Definition: Expansion.h:68
GLOBAL_MAPPING_EXPORT typedef boost::shared_ptr< Mapping > MappingSharedPtr
A shared pointer to a Mapping object.
Definition: Mapping.h:51
Describe a linear system.
boost::shared_ptr< Equation > EquationSharedPtr
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.
Definition: Vmath.cpp:329
A global linear system.
Definition: GlobalLinSys.h:74
Array< OneD, EquationSystemSharedPtr > m_equ
Equation system to solve.
Definition: Driver.h:94
virtual SOLVER_UTILS_EXPORT void v_InitObject(std::ostream &out=std::cout)
Second-stage initialisation.
DriverFactory & GetDriverFactory()
Definition: Driver.cpp:66
LibUtilities::CommSharedPtr m_comm
Communication object.
Definition: Driver.h:85
Base class for the development of solvers.
Definition: Driver.h:65
T Vsum(int n, const T *x, const int incx)
Subtract return sum(x)
Definition: Vmath.cpp:723
boost::shared_ptr< StdExpansion > StdExpansionSharedPtr
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.
Definition: Mapping.cpp:266
Describes the specification for a Basis.
Definition: Basis.h:50
LibUtilities::SessionReaderSharedPtr m_session
Session reader object.
Definition: Driver.h:88
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.
Definition: Vmath.cpp:169
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215