Nektar++
EquationSystem.cpp
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////////
2 //
3 // File EquationSystem.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 // Permission is hereby granted, free of charge, to any person obtaining a
14 // copy of this software and associated documentation files (the "Software"),
15 // to deal in the Software without restriction, including without limitation
16 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
17 // and/or sell copies of the Software, and to permit persons to whom the
18 // Software is furnished to do so, subject to the following conditions:
19 //
20 // The above copyright notice and this permission notice shall be included
21 // in all copies or substantial portions of the Software.
22 //
23 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29 // DEALINGS IN THE SOFTWARE.
30 //
31 // Description: Main wrapper class for Advection Diffusion Reaction Solver
32 //
33 ///////////////////////////////////////////////////////////////////////////////
34 
35 #include <boost/core/ignore_unused.hpp>
36 
39 
41 #include <LocalRegions/MatrixKey.h>
42 #include <MultiRegions/ContField.h>
45 
46 #include <MultiRegions/ExpList.h>
49 
52 
53 #include <GlobalMapping/Mapping.h>
54 
55 #include <boost/format.hpp>
56 
57 #include <iostream>
58 #include <string>
59 
60 using namespace std;
61 
62 namespace Nektar
63 {
64 namespace SolverUtils
65 {
66 
67 std::string EquationSystem::equationSystemTypeLookupIds[2] = {
68  LibUtilities::SessionReader::RegisterEnumValue("DEALIASING", "True", 0),
69  LibUtilities::SessionReader::RegisterEnumValue("DEALIASING", "False", 1)};
70 
71 /**
72  * @class EquationSystem
73  *
74  * This class is a base class for all solver implementations. It
75  * provides the underlying generic functionality and interface for
76  * solving equations.
77  *
78  * To solve a steady-state equation, create a derived class from this
79  * class and reimplement the virtual functions to provide custom
80  * implementation for the problem.
81  *
82  * To solve unsteady problems, derive from the UnsteadySystem class
83  * instead which provides general time integration.
84  */
86 {
87  static EquationSystemFactory instance;
88  return instance;
89 }
90 
91 /**
92  * This constructor is protected as the objects of this class are never
93  * instantiated directly.
94  * @param pSession The session reader holding problem parameters.
95  */
96 EquationSystem::EquationSystem(
99  : m_comm(pSession->GetComm()), m_session(pSession), m_graph(pGraph),
100  m_lambda(0), m_fieldMetaDataMap(LibUtilities::NullFieldMetaDataMap)
101 {
102  // set up session names in fieldMetaDataMap
103  const vector<std::string> filenames = m_session->GetFilenames();
104 
105  for (int i = 0; i < filenames.size(); ++i)
106  {
107  string sessionname = "SessionName";
108  sessionname += boost::lexical_cast<std::string>(i);
109  m_fieldMetaDataMap[sessionname] = filenames[i];
110  m_fieldMetaDataMap["ChkFileNum"] = boost::lexical_cast<std::string>(0);
111  }
112 }
113 
114 /**
115  * @brief Initialisation object for EquationSystem.
116  */
117 void EquationSystem::v_InitObject(bool DeclareFields)
118 {
119  // Save the basename of input file name for output details
120  m_sessionName = m_session->GetSessionName();
121 
122  // Instantiate a field reader/writer
124 
125  // Also read and store the boundary conditions
128  m_session, m_graph);
129 
130  // Set space dimension for use in class
131  m_spacedim = m_graph->GetSpaceDimension();
132 
133  // Setting parameteres for homogenous problems
134  m_HomoDirec = 0;
135  m_useFFT = false;
136  m_homogen_dealiasing = false;
137  m_singleMode = false;
138  m_halfMode = false;
139  m_multipleModes = false;
141 
142  m_verbose = m_session->DefinesCmdLineArgument("verbose");
143  m_root = false;
144  if (0 == m_comm->GetRank())
145  {
146  m_root = true;
147  }
148 
149  if (m_session->DefinesSolverInfo("HOMOGENEOUS"))
150  {
151  std::string HomoStr = m_session->GetSolverInfo("HOMOGENEOUS");
152  m_spacedim = 3;
153 
154  if ((HomoStr == "HOMOGENEOUS1D") || (HomoStr == "Homogeneous1D") ||
155  (HomoStr == "1D") || (HomoStr == "Homo1D"))
156  {
158  m_session->LoadParameter("LZ", m_LhomZ);
159  m_HomoDirec = 1;
160 
161  if (m_session->DefinesSolverInfo("ModeType"))
162  {
163  m_session->MatchSolverInfo("ModeType", "SingleMode",
164  m_singleMode, false);
165  m_session->MatchSolverInfo("ModeType", "HalfMode", m_halfMode,
166  false);
167  m_session->MatchSolverInfo("ModeType", "MultipleModes",
168  m_multipleModes, false);
169  }
170 
171  // Stability Analysis flags
172  if (m_session->DefinesSolverInfo("ModeType"))
173  {
174  if (m_singleMode)
175  {
176  m_npointsZ = 2;
177  }
178  else if (m_halfMode)
179  {
180  m_npointsZ = 1;
181  }
182  else if (m_multipleModes)
183  {
184  m_npointsZ = m_session->GetParameter("HomModesZ");
185  }
186  else
187  {
188  ASSERTL0(false, "SolverInfo ModeType not valid");
189  }
190  }
191  else
192  {
193  m_npointsZ = m_session->GetParameter("HomModesZ");
194  }
195  }
196 
197  if ((HomoStr == "HOMOGENEOUS2D") || (HomoStr == "Homogeneous2D") ||
198  (HomoStr == "2D") || (HomoStr == "Homo2D"))
199  {
201  m_session->LoadParameter("HomModesY", m_npointsY);
202  m_session->LoadParameter("LY", m_LhomY);
203  m_session->LoadParameter("HomModesZ", m_npointsZ);
204  m_session->LoadParameter("LZ", m_LhomZ);
205  m_HomoDirec = 2;
206  }
207 
208  if ((HomoStr == "HOMOGENEOUS3D") || (HomoStr == "Homogeneous3D") ||
209  (HomoStr == "3D") || (HomoStr == "Homo3D"))
210  {
212  m_session->LoadParameter("HomModesY", m_npointsY);
213  m_session->LoadParameter("LY", m_LhomY);
214  m_session->LoadParameter("HomModesZ", m_npointsZ);
215  m_session->LoadParameter("LZ", m_LhomZ);
216  m_HomoDirec = 2;
217  }
218 
219  m_session->MatchSolverInfo("USEFFT", "FFTW", m_useFFT, false);
220 
221  m_session->MatchSolverInfo("DEALIASING", "True", m_homogen_dealiasing,
222  false);
223  }
224  else
225  {
226  // set to default value so can use to identify 2d or 3D
227  // (homogeneous) expansions
228  m_npointsZ = 1;
229  }
230 
231  m_session->MatchSolverInfo("SPECTRALHPDEALIASING", "True",
232  m_specHP_dealiasing, false);
233  if (m_specHP_dealiasing == false)
234  {
235  m_session->MatchSolverInfo("SPECTRALHPDEALIASING", "On",
236  m_specHP_dealiasing, false);
237  }
238 
239  // Options to determine type of projection from file or directly
240  // from constructor
241  if (m_session->DefinesSolverInfo("PROJECTION"))
242  {
243  std::string ProjectStr = m_session->GetSolverInfo("PROJECTION");
244 
245  if ((ProjectStr == "Continuous") || (ProjectStr == "Galerkin") ||
246  (ProjectStr == "CONTINUOUS") || (ProjectStr == "GALERKIN"))
247  {
249  }
250  else if ((ProjectStr == "MixedCGDG") ||
251  (ProjectStr == "Mixed_CG_Discontinuous"))
252  {
254  }
255  else if (ProjectStr == "DisContinuous")
256  {
258  }
259  else
260  {
261  ASSERTL0(false, "PROJECTION value not recognised");
262  }
263  }
264  else
265  {
266  cerr << "Projection type not specified in SOLVERINFO,"
267  "defaulting to continuous Galerkin"
268  << endl;
270  }
271 
272  // Enforce singularity check for some problems
274 
275  int i;
276  int nvariables = m_session->GetVariables().size();
277  bool DeclareCoeffPhysArrays = true;
278 
280  m_spacedim = m_graph->GetSpaceDimension() + m_HomoDirec;
281  m_expdim = m_graph->GetMeshDimension();
282 
283  if (DeclareFields) // declare field if required
284  {
285  /// Continuous field
288  {
289  switch (m_expdim)
290  {
291  case 1:
292  {
295  {
296  const LibUtilities::PointsKey PkeyY(
298  const LibUtilities::BasisKey BkeyY(
300  const LibUtilities::PointsKey PkeyZ(
302  const LibUtilities::BasisKey BkeyZ(
304 
305  for (i = 0; i < m_fields.size(); i++)
306  {
307  m_fields[i] = MemoryManager<
308  MultiRegions ::ContField3DHomogeneous2D>::
309  AllocateSharedPtr(m_session, BkeyY, BkeyZ,
312  m_session->GetVariable(i));
313  }
314  }
315  else
316  {
317  for (i = 0; i < m_fields.size(); i++)
318  {
319  m_fields[i] =
323  m_session->GetVariable(i));
324  }
325  }
326  break;
327  }
328  case 2:
329  {
331  {
332  // Fourier single mode stability analysis
333  if (m_singleMode)
334  {
335  const LibUtilities::PointsKey PkeyZ(
336  m_npointsZ,
338 
339  const LibUtilities::BasisKey BkeyZ(
341  PkeyZ);
342 
343  for (i = 0; i < m_fields.size(); i++)
344  {
345  m_fields[i] = MemoryManager<
346  MultiRegions ::ContField3DHomogeneous1D>::
347  AllocateSharedPtr(
348  m_session, BkeyZ, m_LhomZ, m_useFFT,
350  m_session->GetVariable(i),
352  }
353  }
354  // Half mode stability analysis
355  else if (m_halfMode)
356  {
357  const LibUtilities::PointsKey PkeyZ(
358  m_npointsZ,
360 
361  const LibUtilities::BasisKey BkeyZR(
363  PkeyZ);
364 
365  const LibUtilities::BasisKey BkeyZI(
367  PkeyZ);
368 
369  for (i = 0; i < m_fields.size(); i++)
370  {
371  if (m_session->GetVariable(i).compare("w") == 0)
372  {
373  m_fields[i] = MemoryManager<
374  MultiRegions ::
375  ContField3DHomogeneous1D>::
376  AllocateSharedPtr(
377  m_session, BkeyZI, m_LhomZ,
379  m_graph, m_session->GetVariable(i),
381  }
382  else
383  {
384  m_fields[i] = MemoryManager<
385  MultiRegions ::
386  ContField3DHomogeneous1D>::
387  AllocateSharedPtr(
388  m_session, BkeyZR, m_LhomZ,
390  m_graph, m_session->GetVariable(i),
392  }
393  }
394  }
395  // Normal homogeneous 1D
396  else
397  {
398  const LibUtilities::PointsKey PkeyZ(
400  const LibUtilities::BasisKey BkeyZ(
402 
403  for (i = 0; i < m_fields.size(); i++)
404  {
405  m_fields[i] = MemoryManager<
406  MultiRegions ::ContField3DHomogeneous1D>::
407  AllocateSharedPtr(
408  m_session, BkeyZ, m_LhomZ, m_useFFT,
410  m_session->GetVariable(i),
412  }
413  }
414  }
415  else
416  {
417  i = 0;
421  m_session->GetVariable(i),
422  DeclareCoeffPhysArrays,
424  m_fields[0] = firstfield;
425  for (i = 1; i < m_fields.size(); i++)
426  {
427  if (m_graph->SameExpansionInfo(
428  m_session->GetVariable(0),
429  m_session->GetVariable(i)))
430  {
431  m_fields[i] =
434  *firstfield, m_graph,
435  m_session->GetVariable(i),
436  DeclareCoeffPhysArrays,
438  }
439  else
440  {
441  m_fields[i] =
445  m_session->GetVariable(i),
446  DeclareCoeffPhysArrays,
448  }
449  }
450 
451  if (m_projectionType ==
453  {
454  /// Setting up the normals
457 
458  for (i = 0; i < m_spacedim; ++i)
459  {
460  m_traceNormals[i] =
462  }
463 
464  m_fields[0]->GetTrace()->GetNormals(m_traceNormals);
465  }
466  }
467 
468  break;
469  }
470  case 3:
471  {
472  i = 0;
476  m_session->GetVariable(i),
477  DeclareCoeffPhysArrays,
479 
480  m_fields[0] = firstfield;
481  for (i = 1; i < m_fields.size(); i++)
482  {
483  if (m_graph->SameExpansionInfo(
484  m_session->GetVariable(0),
485  m_session->GetVariable(i)))
486  {
487  m_fields[i] =
490  *firstfield, m_graph,
491  m_session->GetVariable(i),
492  DeclareCoeffPhysArrays,
494  }
495  else
496  {
497  m_fields[i] =
501  m_session->GetVariable(i),
502  DeclareCoeffPhysArrays,
504  }
505  }
506 
507  if (m_projectionType ==
509  {
510  /// Setting up the normals
513  for (i = 0; i < m_spacedim; ++i)
514  {
515  m_traceNormals[i] =
517  }
518 
519  m_fields[0]->GetTrace()->GetNormals(m_traceNormals);
520  // Call the trace on all fields to ensure DG setup.
521  for (i = 1; i < m_fields.size(); ++i)
522  {
523  m_fields[i]->GetTrace();
524  }
525  }
526  break;
527  }
528  default:
529  ASSERTL0(false, "Expansion dimension not recognised");
530  break;
531  }
532  }
533  // Discontinuous field
534  else
535  {
536  switch (m_expdim)
537  {
538  case 1:
539  {
542  {
543  const LibUtilities::PointsKey PkeyY(
545  const LibUtilities::BasisKey BkeyY(
547  const LibUtilities::PointsKey PkeyZ(
549  const LibUtilities::BasisKey BkeyZ(
551 
552  for (i = 0; i < m_fields.size(); i++)
553  {
554  m_fields[i] = MemoryManager<
555  MultiRegions ::DisContField3DHomogeneous2D>::
556  AllocateSharedPtr(m_session, BkeyY, BkeyZ,
559  m_session->GetVariable(i));
560  }
561  }
562  else
563  {
564  for (i = 0; i < m_fields.size(); i++)
565  {
566  m_fields[i] =
570  m_session->GetVariable(i));
571  }
572  }
573 
574  break;
575  }
576  case 2:
577  {
579  {
580  const LibUtilities::PointsKey PkeyZ(
582  const LibUtilities::BasisKey BkeyZ(
584 
585  for (i = 0; i < m_fields.size(); i++)
586  {
587  m_fields[i] = MemoryManager<
588  MultiRegions ::DisContField3DHomogeneous1D>::
589  AllocateSharedPtr(m_session, BkeyZ, m_LhomZ,
590  m_useFFT,
592  m_session->GetVariable(i));
593  }
594  }
595  else
596  {
597  for (i = 0; i < m_fields.size(); i++)
598  {
599  m_fields[i] =
603  m_session->GetVariable(i));
604  }
605  }
606 
607  break;
608  }
609  case 3:
610  {
612  {
613  ASSERTL0(
614  false,
615  "3D fully periodic problems not implemented yet");
616  }
617  else
618  {
619  for (i = 0; i < m_fields.size(); i++)
620  {
621  m_fields[i] =
625  m_session->GetVariable(i));
626  }
627  }
628  break;
629  }
630  default:
631  ASSERTL0(false, "Expansion dimension not recognised");
632  break;
633  }
634 
635  // Setting up the normals
637 
638  for (i = 0; i < m_spacedim; ++i)
639  {
640  m_traceNormals[i] =
642  }
643 
644  m_fields[0]->GetTrace()->GetNormals(m_traceNormals);
645  }
646  // Zero all physical fields initially
647  ZeroPhysFields();
648  }
649 
650  // Set Default Parameter
651  m_session->LoadParameter("Time", m_time, 0.0);
652  m_session->LoadParameter("TimeStep", m_timestep, 0.0);
653  m_session->LoadParameter("NumSteps", m_steps, 0);
654  m_session->LoadParameter("IO_CheckSteps", m_checksteps, 0);
655  m_session->LoadParameter("IO_CheckTime", m_checktime, 0.0);
656  m_session->LoadParameter("FinTime", m_fintime, 0);
657  m_session->LoadParameter("NumQuadPointsError", m_NumQuadPointsError, 0);
658 
659  // Check uniqueness of checkpoint output
660  ASSERTL0((m_checktime == 0.0 && m_checksteps == 0) ||
661  (m_checktime > 0.0 && m_checksteps == 0) ||
662  (m_checktime == 0.0 && m_checksteps > 0),
663  "Only one of IO_CheckTime and IO_CheckSteps "
664  "should be set!");
665  m_session->LoadParameter("TimeIncrementFactor", m_TimeIncrementFactor, 1.0);
666 
667  m_nchk = 0;
668 }
669 
670 /**
671  * @brief Destructor for class EquationSystem.
672  */
674 {
676  LocalRegions::MatrixKey::opLess>::ClearManager();
678  LocalRegions::MatrixKey::opLess>::ClearManager();
679 }
680 
682  std::string name, const MultiRegions::ExpListSharedPtr &field, bool cache)
683 {
684  MultiRegions::ExpListSharedPtr vField = field;
685  if (!field)
686  {
687  vField = m_fields[0];
688  }
689 
690  if (cache)
691  {
692  if ((m_sessionFunctions.find(name) == m_sessionFunctions.end()) ||
693  (m_sessionFunctions[name]->GetSession() != m_session) ||
694  (m_sessionFunctions[name]->GetExpansion() != vField))
695  {
698  m_session, vField, name, cache);
699  }
700 
701  return m_sessionFunctions[name];
702  }
703  else
704  {
706  new SessionFunction(m_session, vField, name, cache));
707  }
708 }
709 
710 /**
711  * If boundary conditions are time-dependent, they will be evaluated at
712  * the time specified.
713  * @param time The time at which to evaluate the BCs
714  */
716 {
717  std::string varName;
718  int nvariables = m_fields.size();
719  for (int i = 0; i < nvariables; ++i)
720  {
721  varName = m_session->GetVariable(i);
722  m_fields[i]->EvaluateBoundaryConditions(time, varName);
723  }
724 }
725 
726 /**
727  * Compute the error in the L2-norm.
728  * @param field The field to compare.
729  * @param exactsoln The exact solution to compare with.
730  * @param Normalised Normalise L2-error.
731  * @returns Error in the L2-norm.
732  */
734  const Array<OneD, NekDouble> &exactsoln,
735  bool Normalised)
736 {
737  NekDouble L2error = -1.0;
738 
739  if (m_NumQuadPointsError == 0)
740  {
741  if (m_fields[field]->GetPhysState() == false)
742  {
743  m_fields[field]->BwdTrans(m_fields[field]->GetCoeffs(),
744  m_fields[field]->UpdatePhys());
745  }
746 
747  if (exactsoln.size())
748  {
749  L2error =
750  m_fields[field]->L2(m_fields[field]->GetPhys(), exactsoln);
751  }
752  else if (m_session->DefinesFunction("ExactSolution"))
753  {
754  Array<OneD, NekDouble> exactsoln(m_fields[field]->GetNpoints());
755 
756  GetFunction("ExactSolution")
757  ->Evaluate(m_session->GetVariable(field), exactsoln, m_time);
758 
759  L2error =
760  m_fields[field]->L2(m_fields[field]->GetPhys(), exactsoln);
761  }
762  else
763  {
764  L2error = m_fields[field]->L2(m_fields[field]->GetPhys());
765  }
766 
767  if (Normalised == true)
768  {
769  Array<OneD, NekDouble> one(m_fields[field]->GetNpoints(), 1.0);
770 
771  NekDouble Vol = m_fields[field]->Integral(one);
772  L2error = sqrt(L2error * L2error / Vol);
773  }
774  }
775  else
776  {
777  Array<OneD, NekDouble> L2INF(2);
778  L2INF = ErrorExtraPoints(field);
779  L2error = L2INF[0];
780  }
781  return L2error;
782 }
783 
784 /**
785  * Compute the error in the L_inf-norm
786  * @param field The field to compare.
787  * @param exactsoln The exact solution to compare with.
788  * @returns Error in the L_inft-norm.
789  */
791  const Array<OneD, NekDouble> &exactsoln)
792 {
793  NekDouble Linferror = -1.0;
794 
795  if (m_NumQuadPointsError == 0)
796  {
797  if (m_fields[field]->GetPhysState() == false)
798  {
799  m_fields[field]->BwdTrans(m_fields[field]->GetCoeffs(),
800  m_fields[field]->UpdatePhys());
801  }
802 
803  if (exactsoln.size())
804  {
805  Linferror =
806  m_fields[field]->Linf(m_fields[field]->GetPhys(), exactsoln);
807  }
808  else if (m_session->DefinesFunction("ExactSolution"))
809  {
810  Array<OneD, NekDouble> exactsoln(m_fields[field]->GetNpoints());
811 
812  GetFunction("ExactSolution")
813  ->Evaluate(m_session->GetVariable(field), exactsoln, m_time);
814 
815  Linferror =
816  m_fields[field]->Linf(m_fields[field]->GetPhys(), exactsoln);
817  }
818  else
819  {
820  Linferror = m_fields[field]->Linf(m_fields[field]->GetPhys());
821  }
822  }
823  else
824  {
825  Array<OneD, NekDouble> L2INF(2);
826  L2INF = ErrorExtraPoints(field);
827  Linferror = L2INF[1];
828  }
829 
830  return Linferror;
831 }
832 
833 /**
834  * Compute the error in the L2-norm, L-inf for a larger number of
835  * quadrature points.
836  * @param field The field to compare.
837  * @returns Error in the L2-norm and L-inf norm.
838  */
840 {
841  int NumModes = GetNumExpModes();
842  Array<OneD, NekDouble> L2INF(2);
843 
847  LibUtilities::eGaussRadauMAlpha1Beta0);
852  const LibUtilities::BasisKey BkeyT1(LibUtilities::eModified_A, NumModes,
853  PkeyT1);
854  const LibUtilities::BasisKey BkeyT2(LibUtilities::eModified_B, NumModes,
855  PkeyT2);
856  const LibUtilities::BasisKey BkeyQ1(LibUtilities::eModified_A, NumModes,
857  PkeyQ1);
858  const LibUtilities::BasisKey BkeyQ2(LibUtilities::eModified_A, NumModes,
859  PkeyQ2);
860 
861  LibUtilities::BasisKeyVector Tkeys, Qkeys;
862 
863  // make a copy of the ExpansionInfoMap
864  SpatialDomains::ExpansionInfoMap NewExpInfo = m_graph->GetExpansionInfo();
867  NewExpInfo);
868 
869  // reset new graph with new keys
870  Tkeys.push_back(BkeyT1);
871  Tkeys.push_back(BkeyT2);
872  m_graph->ResetExpansionInfoToBasisKey(ExpInfo, LibUtilities::eTriangle,
873  Tkeys);
874  Qkeys.push_back(BkeyQ1);
875  Qkeys.push_back(BkeyQ2);
876  m_graph->ResetExpansionInfoToBasisKey(ExpInfo, LibUtilities::eQuadrilateral,
877  Qkeys);
878 
881  NewExpInfo);
882 
883  int ErrorCoordim = ErrorExp->GetCoordim(0);
884  int ErrorNq = ErrorExp->GetTotPoints();
885 
886  Array<OneD, NekDouble> ErrorXc0(ErrorNq, 0.0);
887  Array<OneD, NekDouble> ErrorXc1(ErrorNq, 0.0);
888  Array<OneD, NekDouble> ErrorXc2(ErrorNq, 0.0);
889 
890  switch (ErrorCoordim)
891  {
892  case 1:
893  ErrorExp->GetCoords(ErrorXc0);
894  break;
895  case 2:
896  ErrorExp->GetCoords(ErrorXc0, ErrorXc1);
897  break;
898  case 3:
899  ErrorExp->GetCoords(ErrorXc0, ErrorXc1, ErrorXc2);
900  break;
901  }
903  m_session->GetFunction("ExactSolution", field);
904 
905  // Evaluate the exact solution
906  Array<OneD, NekDouble> ErrorSol(ErrorNq);
907 
908  exSol->Evaluate(ErrorXc0, ErrorXc1, ErrorXc2, m_time, ErrorSol);
909 
910  // Calcualte spectral/hp approximation on the quadrature points
911  // of this new expansion basis
912  ErrorExp->BwdTrans(m_fields[field]->GetCoeffs(), ErrorExp->UpdatePhys());
913 
914  L2INF[0] = ErrorExp->L2(ErrorExp->GetPhys(), ErrorSol);
915  L2INF[1] = ErrorExp->Linf(ErrorExp->GetPhys(), ErrorSol);
916 
917  return L2INF;
918 }
919 
920 /**
921  * Set the physical fields based on a restart file, or a function
922  * describing the initial condition given in the session.
923  * @param initialtime Time at which to evaluate the function.
924  * @param dumpInitialConditions Write the initial condition to file?
925  */
927  bool dumpInitialConditions,
928  const int domain)
929 {
930  boost::ignore_unused(initialtime);
931 
932  if (m_session->GetComm()->GetRank() == 0)
933  {
934  cout << "Initial Conditions:" << endl;
935  }
936 
937  if (m_session->DefinesFunction("InitialConditions"))
938  {
939  GetFunction("InitialConditions")
940  ->Evaluate(m_session->GetVariables(), m_fields, m_time, domain);
941  // Enforce C0 Continutiy of initial condiiton
944  {
945  for (int i = 0; i < m_fields.size(); ++i)
946  {
947  m_fields[i]->LocalToGlobal();
948  m_fields[i]->GlobalToLocal();
949  m_fields[i]->BwdTrans(m_fields[i]->GetCoeffs(),
950  m_fields[i]->UpdatePhys());
951  }
952  }
953 
954  if (m_session->GetComm()->GetRank() == 0)
955  {
956 
957  for (int i = 0; i < m_fields.size(); ++i)
958  {
959  std::string varName = m_session->GetVariable(i);
960  cout << " - Field " << varName << ": "
961  << GetFunction("InitialConditions")
962  ->Describe(varName, domain)
963  << endl;
964  }
965  }
966  }
967  else
968  {
969  int nq = m_fields[0]->GetNpoints();
970  for (int i = 0; i < m_fields.size(); i++)
971  {
972  Vmath::Zero(nq, m_fields[i]->UpdatePhys(), 1);
973  m_fields[i]->SetPhysState(true);
974  Vmath::Zero(m_fields[i]->GetNcoeffs(), m_fields[i]->UpdateCoeffs(),
975  1);
976  if (m_session->GetComm()->GetRank() == 0)
977  {
978  cout << " - Field " << m_session->GetVariable(i)
979  << ": 0 (default)" << endl;
980  }
981  }
982  }
983 
984  if (dumpInitialConditions && m_checksteps && m_nchk == 0)
985  {
987  }
988  ++m_nchk;
989 }
990 
992  Array<OneD, NekDouble> &outfield,
993  const NekDouble time)
994 {
995  ASSERTL0(outfield.size() == m_fields[field]->GetNpoints(),
996  "ExactSolution array size mismatch.");
997  Vmath::Zero(outfield.size(), outfield, 1);
998  if (m_session->DefinesFunction("ExactSolution"))
999  {
1000  GetFunction("ExactSolution")
1001  ->Evaluate(m_session->GetVariable(field), outfield, time);
1002  }
1003 }
1004 
1005 /**
1006  * By default, nothing needs initialising at the EquationSystem level.
1007  */
1009 {
1010 }
1011 
1012 /**
1013  *
1014  */
1016 {
1017 }
1018 
1019 /**
1020  * Virtual function to define if operator in DoSolve is
1021  * negated with regard to the strong form. This is currently
1022  * only used in Arnoldi solves. Default is false.
1023  */
1025 {
1026  return false;
1027 }
1028 
1029 /**
1030  *
1031  */
1033 {
1034 }
1035 
1036 /**
1037  *
1038  */
1040 {
1041 }
1042 
1043 /// Virtual function for generating summary information.
1045 {
1046  SessionSummary(l);
1047 }
1048 
1049 /**
1050  * Write the field data to file. The file is named according to the
1051  * session name with the extension .fld appended.
1052  */
1054 {
1055  WriteFld(m_sessionName + ".fld");
1056 }
1057 
1058 /**
1059  * Zero the physical fields.
1060  */
1062 {
1063  for (int i = 0; i < m_fields.size(); i++)
1064  {
1065  Vmath::Zero(m_fields[i]->GetNpoints(), m_fields[i]->UpdatePhys(), 1);
1066  }
1067 }
1068 
1069 /**
1070  * FwdTrans the m_fields members
1071  */
1073 {
1074  for (int i = 0; i < m_fields.size(); i++)
1075  {
1076  m_fields[i]->FwdTrans(m_fields[i]->GetPhys(),
1077  m_fields[i]->UpdateCoeffs());
1078  m_fields[i]->SetPhysState(false);
1079  }
1080 }
1081 
1082 /**
1083  * Write the n-th checkpoint file.
1084  * @param n The index of the checkpoint file.
1085  */
1087 {
1088  std::string outname =
1089  m_sessionName + "_" + boost::lexical_cast<std::string>(n);
1090  WriteFld(outname + ".chk");
1091 }
1092 
1093 /**
1094  * Write the n-th checkpoint file.
1095  * @param n The index of the checkpoint file.
1096  */
1098  const int n, MultiRegions::ExpListSharedPtr &field,
1099  std::vector<Array<OneD, NekDouble>> &fieldcoeffs,
1100  std::vector<std::string> &variables)
1101 {
1102  std::string outname =
1103  m_sessionName + "_" + boost::lexical_cast<std::string>(n);
1104  WriteFld(outname, field, fieldcoeffs, variables);
1105 }
1106 
1107 /**
1108  * Write the n-th base flow into a .chk file
1109  * @param n The index of the base flow file.
1110  */
1112 {
1113  std::string outname =
1114  m_sessionName + "_BaseFlow_" + boost::lexical_cast<std::string>(n);
1115 
1116  WriteFld(outname + ".chk");
1117 }
1118 
1119 /**
1120  * Writes the field data to a file with the given filename.
1121  * @param outname Filename to write to.
1122  */
1123 void EquationSystem::WriteFld(const std::string &outname)
1124 {
1125  std::vector<Array<OneD, NekDouble>> fieldcoeffs(m_fields.size());
1126  std::vector<std::string> variables(m_fields.size());
1127 
1128  for (int i = 0; i < m_fields.size(); ++i)
1129  {
1130  if (m_fields[i]->GetNcoeffs() == m_fields[0]->GetNcoeffs())
1131  {
1132  fieldcoeffs[i] = m_fields[i]->UpdateCoeffs();
1133  }
1134  else
1135  {
1136  fieldcoeffs[i] = Array<OneD, NekDouble>(m_fields[0]->GetNcoeffs());
1137  m_fields[0]->ExtractCoeffsToCoeffs(
1138  m_fields[i], m_fields[i]->GetCoeffs(), fieldcoeffs[i]);
1139  }
1140  variables[i] = m_boundaryConditions->GetVariable(i);
1141  }
1142 
1143  ExtraFldOutput(fieldcoeffs, variables);
1144 
1145  WriteFld(outname, m_fields[0], fieldcoeffs, variables);
1146 }
1147 
1148 /**
1149  * Writes the field data to a file with the given filename.
1150  * @param outname Filename to write to.
1151  * @param field ExpList on which data is based.
1152  * @param fieldcoeffs An array of array of expansion coefficients.
1153  * @param variables An array of variable names.
1154  */
1155 void EquationSystem::WriteFld(const std::string &outname,
1157  std::vector<Array<OneD, NekDouble>> &fieldcoeffs,
1158  std::vector<std::string> &variables)
1159 {
1160  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef =
1161  field->GetFieldDefinitions();
1162  std::vector<std::vector<NekDouble>> FieldData(FieldDef.size());
1163 
1164  // Copy Data into FieldData and set variable
1165  for (int j = 0; j < fieldcoeffs.size(); ++j)
1166  {
1167  for (int i = 0; i < FieldDef.size(); ++i)
1168  {
1169  // Could do a search here to find correct variable
1170  FieldDef[i]->m_fields.push_back(variables[j]);
1171  field->AppendFieldData(FieldDef[i], FieldData[i], fieldcoeffs[j]);
1172  }
1173  }
1174 
1175  // Update time in field info if required
1176  if (m_fieldMetaDataMap.find("Time") != m_fieldMetaDataMap.end())
1177  {
1178  m_fieldMetaDataMap["Time"] = boost::lexical_cast<std::string>(m_time);
1179  }
1180 
1181  // Update step in field info if required
1182  if (m_fieldMetaDataMap.find("ChkFileNum") != m_fieldMetaDataMap.end())
1183  {
1184  m_fieldMetaDataMap["ChkFileNum"] =
1185  boost::lexical_cast<std::string>(m_nchk);
1186  }
1187 
1188  // If necessary, add mapping information to metadata
1189  // and output mapping coordinates
1191  fields[0] = field;
1195  mapping->Output(fieldMetaDataMap, outname);
1196 
1197  // If necessary, add informaton for moving frame reference to metadata
1198  if (m_fieldMetaDataMap.find("Theta_x") != m_fieldMetaDataMap.end())
1199  {
1200  // if one theta exists, add all three thetas
1201  std::vector<std::string> vSuffix = {"_x", "_y", "_z"};
1202  for (int i = 0; i < 3; ++i)
1203  {
1204  std::string sTheta = "Theta" + vSuffix[i];
1205  m_fieldMetaDataMap[sTheta] =
1206  boost::lexical_cast<std::string>(m_movingFrameTheta[i]);
1207  }
1208  }
1209 
1210 #ifdef NEKTAR_DISABLE_BACKUPS
1211  bool backup = false;
1212 #else
1213  bool backup = true;
1214 #endif
1215 
1216  m_fld->Write(outname, FieldDef, FieldData, fieldMetaDataMap, backup);
1217 }
1218 
1219 /**
1220  * Import field from infile and load into \a m_fields. This routine will
1221  * also perform a \a BwdTrans to ensure data is in both the physical and
1222  * coefficient storage.
1223  * @param infile Filename to read.
1224  */
1226  const std::string &infile,
1228 {
1229  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef;
1230  std::vector<std::vector<NekDouble>> FieldData;
1231  LibUtilities::FieldIOSharedPtr field_fld =
1233  field_fld->Import(infile, FieldDef, FieldData);
1234 
1235  // Copy FieldData into m_fields
1236  for (int j = 0; j < pFields.size(); ++j)
1237  {
1238  Vmath::Zero(pFields[j]->GetNcoeffs(), pFields[j]->UpdateCoeffs(), 1);
1239 
1240  for (int i = 0; i < FieldDef.size(); ++i)
1241  {
1242  ASSERTL1(FieldDef[i]->m_fields[j] == m_session->GetVariable(j),
1243  std::string("Order of ") + infile +
1244  std::string(" data and that defined in "
1245  "m_boundaryconditions differs"));
1246 
1247  pFields[j]->ExtractDataToCoeffs(FieldDef[i], FieldData[i],
1248  FieldDef[i]->m_fields[j],
1249  pFields[j]->UpdateCoeffs());
1250  }
1251  pFields[j]->BwdTrans(pFields[j]->GetCoeffs(), pFields[j]->UpdatePhys());
1252  }
1253 }
1254 
1255 /**
1256  * Import field from infile and load into \a m_fields. This routine will
1257  * also perform a \a BwdTrans to ensure data is in both the physical and
1258  * coefficient storage.
1259  * @param infile Filename to read.
1260  * If optionan \a ndomains is specified it assumes we loop over nodmains
1261  * for each nvariables.
1262  */
1264  const std::string &infile,
1265  Array<OneD, MultiRegions::ExpListSharedPtr> &pFields, const int ndomains)
1266 {
1267  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef;
1268  std::vector<std::vector<NekDouble>> FieldData;
1269 
1270  LibUtilities::Import(infile, FieldDef, FieldData);
1271 
1272  int nvariables = GetNvariables();
1273 
1274  ASSERTL0(
1275  ndomains * nvariables == pFields.size(),
1276  "Number of fields does not match the number of variables and domains");
1277 
1278  // Copy FieldData into m_fields
1279  for (int j = 0; j < ndomains; ++j)
1280  {
1281  for (int i = 0; i < nvariables; ++i)
1282  {
1283  Vmath::Zero(pFields[j * nvariables + i]->GetNcoeffs(),
1284  pFields[j * nvariables + i]->UpdateCoeffs(), 1);
1285 
1286  for (int n = 0; n < FieldDef.size(); ++n)
1287  {
1288  ASSERTL1(FieldDef[n]->m_fields[i] == m_session->GetVariable(i),
1289  std::string("Order of ") + infile +
1290  std::string(" data and that defined in "
1291  "m_boundaryconditions differs"));
1292 
1293  pFields[j * nvariables + i]->ExtractDataToCoeffs(
1294  FieldDef[n], FieldData[n], FieldDef[n]->m_fields[i],
1295  pFields[j * nvariables + i]->UpdateCoeffs());
1296  }
1297  pFields[j * nvariables + i]->BwdTrans(
1298  pFields[j * nvariables + i]->GetCoeffs(),
1299  pFields[j * nvariables + i]->UpdatePhys());
1300  }
1301  }
1302 }
1303 
1304 /**
1305  * Import field from infile and load into \a pField. This routine will
1306  * also perform a \a BwdTrans to ensure data is in both the physical and
1307  * coefficient storage.
1308  */
1309 void EquationSystem::ImportFld(const std::string &infile,
1311  std::string &pFieldName)
1312 {
1313  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef;
1314  std::vector<std::vector<NekDouble>> FieldData;
1315 
1316  LibUtilities::FieldIOSharedPtr field_fld =
1318  field_fld->Import(infile, FieldDef, FieldData);
1319  int idx = -1;
1320 
1321  Vmath::Zero(pField->GetNcoeffs(), pField->UpdateCoeffs(), 1);
1322 
1323  for (int i = 0; i < FieldDef.size(); ++i)
1324  {
1325  // find the index of the required field in the file.
1326  for (int j = 0; j < FieldData.size(); ++j)
1327  {
1328  if (FieldDef[i]->m_fields[j] == pFieldName)
1329  {
1330  idx = j;
1331  }
1332  }
1333  ASSERTL1(idx >= 0, "Field " + pFieldName + " not found.");
1334 
1335  pField->ExtractDataToCoeffs(FieldDef[i], FieldData[i],
1336  FieldDef[i]->m_fields[idx],
1337  pField->UpdateCoeffs());
1338  }
1339  pField->BwdTrans(pField->GetCoeffs(), pField->UpdatePhys());
1340 }
1341 
1342 /**
1343  * Import field from infile and load into the array \a coeffs.
1344  *
1345  * @param infile Filename to read.
1346  * @param fieldStr an array of string identifying fields to be imported
1347  * @param coeffs array of array of coefficients to store imported data
1348  */
1349 void EquationSystem::ImportFld(const std::string &infile,
1350  std::vector<std::string> &fieldStr,
1351  Array<OneD, Array<OneD, NekDouble>> &coeffs)
1352 {
1353 
1354  ASSERTL0(fieldStr.size() <= coeffs.size(),
1355  "length of fieldstr should be the same as pFields");
1356 
1357  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef;
1358  std::vector<std::vector<NekDouble>> FieldData;
1359 
1360  LibUtilities::FieldIOSharedPtr field_fld =
1362  field_fld->Import(infile, FieldDef, FieldData);
1363 
1364  // Copy FieldData into m_fields
1365  for (int j = 0; j < fieldStr.size(); ++j)
1366  {
1367  Vmath::Zero(coeffs[j].size(), coeffs[j], 1);
1368  for (int i = 0; i < FieldDef.size(); ++i)
1369  {
1370  m_fields[0]->ExtractDataToCoeffs(FieldDef[i], FieldData[i],
1371  fieldStr[j], coeffs[j]);
1372  }
1373  }
1374 }
1375 
1376 /**
1377  * Write out a summary of the session data.
1378  * @param out Output stream to write data to.
1379  */
1381 {
1382  AddSummaryItem(s, "EquationType", m_session->GetSolverInfo("EQTYPE"));
1383  AddSummaryItem(s, "Session Name", m_sessionName);
1384  AddSummaryItem(s, "Spatial Dim.", m_spacedim);
1385  AddSummaryItem(s, "Max SEM Exp. Order",
1386  m_fields[0]->EvalBasisNumModesMax());
1387 
1388  if (m_session->GetComm()->GetSize() > 1)
1389  {
1390  AddSummaryItem(s, "Num. Processes", m_session->GetComm()->GetSize());
1391  }
1392 
1394  {
1395  AddSummaryItem(s, "Quasi-3D", "Homogeneous in z-direction");
1396  AddSummaryItem(s, "Expansion Dim.", m_expdim + 1);
1397  AddSummaryItem(s, "Num. Hom. Modes (z)", m_npointsZ);
1398  AddSummaryItem(s, "Hom. length (LZ)", m_LhomZ);
1399  AddSummaryItem(s, "FFT Type", m_useFFT ? "FFTW" : "MVM");
1400  if (m_halfMode)
1401  {
1402  AddSummaryItem(s, "ModeType", "Half Mode");
1403  }
1404  else if (m_singleMode)
1405  {
1406  AddSummaryItem(s, "ModeType", "Single Mode");
1407  }
1408  else if (m_multipleModes)
1409  {
1410  AddSummaryItem(s, "ModeType", "Multiple Modes");
1411  }
1412  }
1413  else if (m_HomogeneousType == eHomogeneous2D)
1414  {
1415  AddSummaryItem(s, "Quasi-3D", "Homogeneous in yz-plane");
1416  AddSummaryItem(s, "Expansion Dim.", m_expdim + 2);
1417  AddSummaryItem(s, "Num. Hom. Modes (y)", m_npointsY);
1418  AddSummaryItem(s, "Num. Hom. Modes (z)", m_npointsZ);
1419  AddSummaryItem(s, "Hom. length (LY)", m_LhomY);
1420  AddSummaryItem(s, "Hom. length (LZ)", m_LhomZ);
1421  AddSummaryItem(s, "FFT Type", m_useFFT ? "FFTW" : "MVM");
1422  }
1423  else
1424  {
1425  AddSummaryItem(s, "Expansion Dim.", m_expdim);
1426  }
1427 
1428  if (m_session->DefinesSolverInfo("UpwindType"))
1429  {
1430  AddSummaryItem(s, "Riemann Solver",
1431  m_session->GetSolverInfo("UpwindType"));
1432  }
1433 
1434  if (m_session->DefinesSolverInfo("AdvectionType"))
1435  {
1436  std::string AdvectionType;
1437  AdvectionType = m_session->GetSolverInfo("AdvectionType");
1439  s, "Advection Type",
1440  GetAdvectionFactory().GetClassDescription(AdvectionType));
1441  }
1442 
1444  {
1445  AddSummaryItem(s, "Projection Type", "Continuous Galerkin");
1446  }
1448  {
1449  AddSummaryItem(s, "Projection Type", "Discontinuous Galerkin");
1450  }
1452  {
1453  AddSummaryItem(s, "Projection Type",
1454  "Mixed Continuous Galerkin and Discontinuous");
1455  }
1456 
1457  if (m_session->DefinesSolverInfo("DiffusionType"))
1458  {
1459  std::string DiffusionType;
1460  DiffusionType = m_session->GetSolverInfo("DiffusionType");
1462  s, "Diffusion Type",
1463  GetDiffusionFactory().GetClassDescription(DiffusionType));
1464  }
1465 }
1466 
1468 {
1469  return Array<OneD, bool>(m_session->GetVariables().size(), false);
1470 }
1471 
1473 {
1474  ASSERTL0(false, "This function is not valid for the Base class");
1476  return null;
1477 }
1478 
1480  std::vector<Array<OneD, NekDouble>> &fieldcoeffs,
1481  std::vector<std::string> &variables)
1482 {
1483  boost::ignore_unused(fieldcoeffs, variables);
1484 }
1485 
1486 } // namespace SolverUtils
1487 } // namespace Nektar
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:215
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
Definition: ErrorUtil.hpp:249
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:269
Describes the specification for a Basis.
Definition: Basis.h:50
static std::shared_ptr< FieldIO > CreateForFile(const LibUtilities::SessionReaderSharedPtr session, const std::string &filename)
Construct a FieldIO object for a given input filename.
Definition: FieldIO.cpp:227
static std::shared_ptr< FieldIO > CreateDefault(const LibUtilities::SessionReaderSharedPtr session)
Returns an object for the default FieldIO method.
Definition: FieldIO.cpp:198
Provides a generic Factory class.
Definition: NekFactory.hpp:105
Defines a specification for a set of points.
Definition: Points.h:59
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.
virtual SOLVER_UTILS_EXPORT void v_DoInitialise()
Virtual function for initialisation implementation.
int m_spacedim
Spatial dimension (>= expansion dim).
SOLVER_UTILS_EXPORT int GetTraceNpoints()
bool m_useFFT
Flag to determine if FFT is used for homogeneous transform.
SOLVER_UTILS_EXPORT int GetNvariables()
int m_expdim
Expansion dimension.
LibUtilities::FieldIOSharedPtr m_fld
Field input/output.
SOLVER_UTILS_EXPORT void SessionSummary(SummaryList &vSummary)
Write out a session summary.
virtual SOLVER_UTILS_EXPORT void v_SetInitialConditions(NekDouble initialtime=0.0, bool dumpInitialConditions=true, const int domain=0)
virtual SOLVER_UTILS_EXPORT ~EquationSystem()
Destructor.
SpatialDomains::MeshGraphSharedPtr m_graph
Pointer to graph defining mesh.
LibUtilities::CommSharedPtr m_comm
Communicator.
NekDouble m_timestep
Time step size.
SOLVER_UTILS_EXPORT void ImportFld(const std::string &infile, Array< OneD, MultiRegions::ExpListSharedPtr > &pFields)
Input field data from the given file.
NekDouble m_time
Current time of simulation.
bool m_multipleModes
Flag to determine if use multiple homogenenous modes are used.
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
virtual SOLVER_UTILS_EXPORT void v_TransPhysToCoeff()
Virtual function for transformation to coefficient space.
NekDouble m_fintime
Finish time of the simulation.
virtual SOLVER_UTILS_EXPORT NekDouble v_L2Error(unsigned int field, const Array< OneD, NekDouble > &exactsoln=NullNekDouble1DArray, bool Normalised=false)
Virtual function for the L_2 error computation between fields and a given exact solution.
virtual SOLVER_UTILS_EXPORT MultiRegions::ExpListSharedPtr v_GetPressure(void)
int m_npointsZ
number of points in Z direction (if homogeneous)
SOLVER_UTILS_EXPORT void Checkpoint_Output(const int n)
Write checkpoint file of m_fields.
std::map< std::string, SolverUtils::SessionFunctionSharedPtr > m_sessionFunctions
Map of known SessionFunctions.
NekDouble m_checktime
Time between checkpoints.
virtual SOLVER_UTILS_EXPORT void v_InitObject(bool DeclareFeld=true)
Initialisation object for EquationSystem.
virtual SOLVER_UTILS_EXPORT NekDouble v_LinfError(unsigned int field, const Array< OneD, NekDouble > &exactsoln=NullNekDouble1DArray)
Virtual function for the L_inf error computation between fields and a given exact solution.
virtual SOLVER_UTILS_EXPORT void v_EvaluateExactSolution(unsigned int field, Array< OneD, NekDouble > &outfield, const NekDouble time)
SOLVER_UTILS_EXPORT void WriteFld(const std::string &outname)
Write field data to the given filename.
SOLVER_UTILS_EXPORT void ImportFldToMultiDomains(const std::string &infile, Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const int ndomains)
Input field data from the given file to multiple domains.
NekDouble m_LhomY
physical length in Y direction (if homogeneous)
int m_npointsY
number of points in Y direction (if homogeneous)
std::string m_sessionName
Name of the session.
SOLVER_UTILS_EXPORT void Checkpoint_BaseFlow(const int n)
Write base flow file of m_fields.
bool m_specHP_dealiasing
Flag to determine if dealisising is usde for the Spectral/hp element discretisation.
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
bool m_singleMode
Flag to determine if single homogeneous mode is used.
virtual SOLVER_UTILS_EXPORT bool v_NegatedOp()
Virtual function to identify if operator is negated in DoSolve.
SOLVER_UTILS_EXPORT void ZeroPhysFields()
int m_HomoDirec
number of homogenous directions
virtual SOLVER_UTILS_EXPORT void v_ExtraFldOutput(std::vector< Array< OneD, NekDouble >> &fieldcoeffs, std::vector< std::string > &variables)
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Array holding trace normals for DG simulations in the forwards direction.
enum HomogeneousType m_HomogeneousType
SOLVER_UTILS_EXPORT Array< OneD, NekDouble > ErrorExtraPoints(unsigned int field)
Compute error (L2 and L_inf) over an larger set of quadrature points return [L2 Linf].
SOLVER_UTILS_EXPORT int GetNpoints()
Array< OneD, bool > m_checkIfSystemSingular
Flag to indicate if the fields should be checked for singularity.
bool m_homogen_dealiasing
Flag to determine if dealiasing is used for homogeneous simulations.
SOLVER_UTILS_EXPORT int GetNcoeffs()
int m_nchk
Number of checkpoints written so far.
enum MultiRegions::ProjectionType m_projectionType
Type of projection; e.g continuous or discontinuous.
virtual SOLVER_UTILS_EXPORT void v_DoSolve()
Virtual function for solve implementation.
virtual SOLVER_UTILS_EXPORT Array< OneD, bool > v_GetSystemSingularChecks()
int m_steps
Number of steps to take.
int m_NumQuadPointsError
Number of Quadrature points used to work out the error.
Array< OneD, NekDouble > m_movingFrameTheta
Moving frame of reference angles with respect to the.
LibUtilities::FieldMetaDataMap m_fieldMetaDataMap
Map to identify relevant solver info to dump in output fields.
SOLVER_UTILS_EXPORT int GetNumExpModes()
SpatialDomains::BoundaryConditionsSharedPtr m_boundaryConditions
Pointer to boundary conditions object.
SOLVER_UTILS_EXPORT void SetBoundaryConditions(NekDouble time)
Evaluates the boundary conditions at the given time.
NekDouble m_LhomZ
physical length in Z direction (if homogeneous)
virtual SOLVER_UTILS_EXPORT void v_GenerateSummary(SummaryList &l)
Virtual function for generating summary information.
SOLVER_UTILS_EXPORT void ExtraFldOutput(std::vector< Array< OneD, NekDouble >> &fieldcoeffs, std::vector< std::string > &variables)
bool m_halfMode
Flag to determine if half homogeneous mode is used.
SOLVER_UTILS_EXPORT void FwdTransFields()
int m_checksteps
Number of steps between checkpoints.
SOLVER_UTILS_EXPORT SessionFunctionSharedPtr GetFunction(std::string name, const MultiRegions::ExpListSharedPtr &field=MultiRegions::NullExpListSharedPtr, bool cache=false)
Get a SessionFunction by name.
virtual SOLVER_UTILS_EXPORT void v_Output(void)
virtual SOLVER_UTILS_EXPORT void v_TransCoeffToPhys()
Virtual function for transformation to physical space.
GLOBAL_MAPPING_EXPORT typedef std::shared_ptr< Mapping > MappingSharedPtr
A shared pointer to a Mapping object.
Definition: Mapping.h:50
std::vector< BasisKey > BasisKeyVector
Name for a vector of BasisKeys.
std::shared_ptr< FieldIO > FieldIOSharedPtr
Definition: FieldIO.h:301
void Import(const std::string &infilename, std::vector< FieldDefinitionsSharedPtr > &fielddefs, std::vector< std::vector< NekDouble >> &fielddata, FieldMetaDataMap &fieldinfomap, const Array< OneD, int > &ElementIDs)
This function allows for data to be imported from an FLD file when a session and/or communicator is n...
Definition: FieldIO.cpp:291
std::map< std::string, std::string > FieldMetaDataMap
Definition: FieldIO.h:52
std::shared_ptr< SessionReader > SessionReaderSharedPtr
std::shared_ptr< Equation > EquationSharedPtr
Definition: Equation.h:130
static FieldMetaDataMap NullFieldMetaDataMap
Definition: FieldIO.h:53
@ eFourierEvenlySpaced
1D Evenly-spaced points using Fourier Fit
Definition: PointsType.h:76
@ eGaussLobattoLegendre
1D Gauss-Lobatto-Legendre quadrature points
Definition: PointsType.h:53
@ eFourierSingleModeSpaced
1D Non Evenly-spaced points for Single Mode analysis
Definition: PointsType.h:77
@ eModified_B
Principle Modified Functions .
Definition: BasisType.h:51
@ eFourierSingleMode
Fourier ModifiedExpansion with just the first mode .
Definition: BasisType.h:66
@ eModified_A
Principle Modified Functions .
Definition: BasisType.h:50
@ eFourierHalfModeIm
Fourier Modified expansions with just the imaginary part of the first mode .
Definition: BasisType.h:70
@ eFourierHalfModeRe
Fourier Modified expansions with just the real part of the first mode .
Definition: BasisType.h:68
@ eFourier
Fourier Expansion .
Definition: BasisType.h:57
std::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
std::shared_ptr< ContField > ContFieldSharedPtr
Definition: ContField.h:277
AdvectionFactory & GetAdvectionFactory()
Gets the factory for initialising advection objects.
Definition: Advection.cpp:47
std::vector< std::pair< std::string, std::string > > SummaryList
Definition: Misc.h:48
DiffusionFactory & GetDiffusionFactory()
Definition: Diffusion.cpp:41
EquationSystemFactory & GetEquationSystemFactory()
void AddSummaryItem(SummaryList &l, const std::string &name, const std::string &value)
Adds a summary item to the summary info list.
Definition: Misc.cpp:49
std::shared_ptr< SessionFunction > SessionFunctionSharedPtr
std::shared_ptr< ExpansionInfoMap > ExpansionInfoMapShPtr
Definition: MeshGraph.h:145
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:172
std::map< int, ExpansionInfoShPtr > ExpansionInfoMap
Definition: MeshGraph.h:143
The above copyright notice and this permission notice shall be included.
Definition: CoupledSolver.h:1
NekMatrix< NekMatrix< NekMatrix< NekDouble, StandardMatrixTag >, ScaledMatrixTag >, BlockMatrixTag > DNekScalBlkMat
Definition: NekTypeDefs.hpp:68
NekMatrix< NekMatrix< NekDouble, StandardMatrixTag >, ScaledMatrixTag > DNekScalMat
double NekDouble
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:492
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:291
Used to lookup the create function in NekManager.
Definition: MatrixKey.h:71