Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
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 // 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: Main wrapper class for Advection Diffusion Reaction Solver
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
36 
37 #include <SolverUtils/EquationSystem.h>
38 
39 #include <LocalRegions/MatrixKey.h>
40 #include <LibUtilities/BasicUtils/Equation.h>
41 #include <MultiRegions/ContField1D.h>
42 #include <MultiRegions/ContField2D.h>
43 #include <MultiRegions/ContField3D.h>
44 #include <MultiRegions/ContField3DHomogeneous1D.h>
45 #include <MultiRegions/ContField3DHomogeneous2D.h>
46 
47 #include <MultiRegions/ExpList2D.h> // for ExpList2D, etc
48 #include <MultiRegions/ExpList3D.h> // for ExpList3D
49 #include <MultiRegions/ExpList3DHomogeneous1D.h>
50 #include <MultiRegions/ExpList3DHomogeneous2D.h>
51 
52 #include <SolverUtils/AdvectionSystem.h>
53 #include <SolverUtils/Diffusion/Diffusion.h>
54 
55 #include <GlobalMapping/Mapping.h>
56 
57 #include <boost/format.hpp>
58 # include <boost/function.hpp>
59 
60 #include <iostream>
61 #include <string>
62 
63 using namespace std;
64 
65 namespace Nektar
66 {
67  namespace SolverUtils
68  {
69  /**
70  * @class EquationSystem
71  *
72  * This class is a base class for all solver implementations. It
73  * provides the underlying generic functionality and interface for
74  * solving equations.
75  *
76  * To solve a steady-state equation, create a derived class from this
77  * class and reimplement the virtual functions to provide custom
78  * implementation for the problem.
79  *
80  * To solve unsteady problems, derive from the UnsteadySystem class
81  * instead which provides general time integration.
82  */
83  EquationSystemFactory& GetEquationSystemFactory()
84  {
85  typedef Loki::SingletonHolder<EquationSystemFactory,
86  Loki::CreateUsingNew,
87  Loki::NoDestroy,
88  Loki::ClassLevelLockable> Type;
89  return Type::Instance();
90  }
91 
92  /**
93  * This constructor is protected as the objects of this class are never
94  * instantiated directly.
95  * @param pSession The session reader holding problem parameters.
96  */
97  EquationSystem::EquationSystem(
98  const LibUtilities::SessionReaderSharedPtr& pSession)
99  : m_comm (pSession->GetComm()),
100  m_session (pSession),
101  m_lambda (0),
102  m_fieldMetaDataMap(LibUtilities::NullFieldMetaDataMap)
103  {
104  // set up session names in fieldMetaDataMap
105  const vector<std::string> filenames = m_session->GetFilenames();
106 
107  for(int i = 0; i < filenames.size(); ++i)
108  {
109  string sessionname = "SessionName";
110  sessionname += boost::lexical_cast<std::string>(i);
111  m_fieldMetaDataMap[sessionname] = filenames[i];
112  m_fieldMetaDataMap["ChkFileNum"] = boost::lexical_cast<std::string>(0);
113  }
114 
115  }
116 
117  /**
118  * @brief Initialisation object for EquationSystem.
119  */
120  void EquationSystem::v_InitObject()
121  {
122  // Save the basename of input file name for output details
123  m_sessionName = m_session->GetSessionName();
124 
125  // Instantiate a field reader/writer
126  m_fld = LibUtilities::FieldIO::CreateDefault(m_session);
127 
128  // Read the geometry and the expansion information
129  m_graph = SpatialDomains::MeshGraph::Read(m_session);
130 
131  // Also read and store the boundary conditions
132  m_boundaryConditions =
133  MemoryManager<SpatialDomains::BoundaryConditions>::
134  AllocateSharedPtr(m_session, m_graph);
135 
136  // Set space dimension for use in class
137  m_spacedim = m_graph->GetSpaceDimension();
138 
139  // Setting parameteres for homogenous problems
140  m_HomoDirec = 0;
141  m_useFFT = false;
142  m_homogen_dealiasing = false;
143  m_singleMode = false;
144  m_halfMode = false;
145  m_multipleModes = false;
146  m_HomogeneousType = eNotHomogeneous;
147 
148  if (m_session->DefinesSolverInfo("HOMOGENEOUS"))
149  {
150  std::string HomoStr = m_session->GetSolverInfo("HOMOGENEOUS");
151  m_spacedim = 3;
152 
153  if ((HomoStr == "HOMOGENEOUS1D") || (HomoStr == "Homogeneous1D")
154  || (HomoStr == "1D") || (HomoStr == "Homo1D"))
155  {
156  m_HomogeneousType = eHomogeneous1D;
157  m_session->LoadParameter("LZ", m_LhomZ);
158  m_HomoDirec = 1;
159 
160  if(m_session->DefinesSolverInfo("ModeType"))
161  {
162  m_session->MatchSolverInfo("ModeType", "SingleMode",
163  m_singleMode, false);
164  m_session->MatchSolverInfo("ModeType", "HalfMode",
165  m_halfMode, false);
166  m_session->MatchSolverInfo("ModeType", "MultipleModes",
167  m_multipleModes, false);
168  }
169 
170  // Stability Analysis flags
171  if (m_session->DefinesSolverInfo("ModeType"))
172  {
173  if(m_singleMode)
174  {
175  m_npointsZ = 2;
176  }
177  else if(m_halfMode)
178  {
179  m_npointsZ = 1;
180  }
181  else if(m_multipleModes)
182  {
183  m_npointsZ = m_session->GetParameter("HomModesZ");
184  }
185  else
186  {
187  ASSERTL0(false, "SolverInfo ModeType not valid");
188  }
189  }
190  else
191  {
192  m_npointsZ = m_session->GetParameter("HomModesZ");
193  }
194  }
195 
196  if ((HomoStr == "HOMOGENEOUS2D") || (HomoStr == "Homogeneous2D")
197  || (HomoStr == "2D") || (HomoStr == "Homo2D"))
198  {
199  m_HomogeneousType = eHomogeneous2D;
200  m_session->LoadParameter("HomModesY", m_npointsY);
201  m_session->LoadParameter("LY", m_LhomY);
202  m_session->LoadParameter("HomModesZ", m_npointsZ);
203  m_session->LoadParameter("LZ", m_LhomZ);
204  m_HomoDirec = 2;
205  }
206 
207  if ((HomoStr == "HOMOGENEOUS3D") || (HomoStr == "Homogeneous3D")
208  || (HomoStr == "3D") || (HomoStr == "Homo3D"))
209  {
210  m_HomogeneousType = eHomogeneous3D;
211  m_session->LoadParameter("HomModesY", m_npointsY);
212  m_session->LoadParameter("LY", m_LhomY);
213  m_session->LoadParameter("HomModesZ", m_npointsZ);
214  m_session->LoadParameter("LZ", m_LhomZ);
215  m_HomoDirec = 2;
216  }
217 
218  m_session->MatchSolverInfo("USEFFT", "FFTW", m_useFFT, false);
219 
220  m_session->MatchSolverInfo("DEALIASING", "True",
221  m_homogen_dealiasing, false);
222  if(m_homogen_dealiasing == false)
223  {
224  m_session->MatchSolverInfo("DEALIASING", "On",
225  m_homogen_dealiasing, false);
226  }
227  }
228  else
229  {
230  // set to default value so can use to identify 2d or 3D
231  // (homogeneous) expansions
232  m_npointsZ = 1;
233  }
234 
235  m_session->MatchSolverInfo("SPECTRALHPDEALIASING", "True",
236  m_specHP_dealiasing, false);
237  if (m_specHP_dealiasing == false)
238  {
239  m_session->MatchSolverInfo("SPECTRALHPDEALIASING", "On",
240  m_specHP_dealiasing, false);
241  }
242 
243  // Options to determine type of projection from file or directly
244  // from constructor
245  if (m_session->DefinesSolverInfo("PROJECTION"))
246  {
247  std::string ProjectStr = m_session->GetSolverInfo("PROJECTION");
248 
249  if ((ProjectStr == "Continuous") || (ProjectStr == "Galerkin") ||
250  (ProjectStr == "CONTINUOUS") || (ProjectStr == "GALERKIN"))
251  {
252  m_projectionType = MultiRegions::eGalerkin;
253  }
254  else if ((ProjectStr == "MixedCGDG") ||
255  (ProjectStr == "Mixed_CG_Discontinuous"))
256  {
257  m_projectionType = MultiRegions::eMixed_CG_Discontinuous;
258  }
259  else if(ProjectStr == "DisContinuous")
260  {
261  m_projectionType = MultiRegions::eDiscontinuous;
262  }
263  else
264  {
265  ASSERTL0(false,"PROJECTION value not recognised");
266  }
267  }
268  else
269  {
270  cerr << "Projection type not specified in SOLVERINFO,"
271  "defaulting to continuous Galerkin" << endl;
272  m_projectionType = MultiRegions::eGalerkin;
273  }
274 
275  // Enforce singularity check for some problems
276  m_checkIfSystemSingular = v_GetSystemSingularChecks();
277 
278  int i;
279  int nvariables = m_session->GetVariables().size();
280  bool DeclareCoeffPhysArrays = true;
281 
282 
283  m_fields = Array<OneD, MultiRegions::ExpListSharedPtr>(nvariables);
284  m_spacedim = m_graph->GetSpaceDimension()+m_HomoDirec;
285  m_expdim = m_graph->GetMeshDimension();
286 
287  /// Continuous field
288  if (m_projectionType == MultiRegions::eGalerkin ||
289  m_projectionType == MultiRegions::eMixed_CG_Discontinuous)
290  {
291  switch(m_expdim)
292  {
293  case 1:
294  {
295  if (m_HomogeneousType == eHomogeneous2D
296  || m_HomogeneousType == eHomogeneous3D)
297  {
298  const LibUtilities::PointsKey PkeyY
299  (m_npointsY, LibUtilities::eFourierEvenlySpaced);
300  const LibUtilities::BasisKey BkeyY
301  (LibUtilities::eFourier, m_npointsY, PkeyY);
302  const LibUtilities::PointsKey PkeyZ
303  (m_npointsZ, LibUtilities::eFourierEvenlySpaced);
304  const LibUtilities::BasisKey
305  BkeyZ(LibUtilities::eFourier, m_npointsZ, PkeyZ);
306 
307  for (i = 0; i < m_fields.num_elements(); i++)
308  {
309  m_fields[i] = MemoryManager<MultiRegions
310  ::ContField3DHomogeneous2D>
311  ::AllocateSharedPtr(
312  m_session, BkeyY, BkeyZ, m_LhomY,
313  m_LhomZ, m_useFFT,
314  m_homogen_dealiasing, m_graph,
315  m_session->GetVariable(i));
316  }
317  }
318  else
319  {
320  for (i = 0; i < m_fields.num_elements(); i++)
321  {
322  m_fields[i] = MemoryManager
323  <MultiRegions::ContField1D>::
324  AllocateSharedPtr(
325  m_session, m_graph,
326  m_session->GetVariable(i));
327  }
328  }
329  break;
330  }
331  case 2:
332  {
333  if (m_HomogeneousType == eHomogeneous1D)
334  {
335  // Fourier single mode stability analysis
336  if (m_singleMode)
337  {
338  const LibUtilities::PointsKey PkeyZ(
339  m_npointsZ,
340  LibUtilities::eFourierSingleModeSpaced);
341 
342  const LibUtilities::BasisKey BkeyZ(
343  LibUtilities::eFourierSingleMode,
344  m_npointsZ,
345  PkeyZ);
346 
347  for(i = 0; i < m_fields.num_elements(); i++)
348  {
349  m_fields[i] = MemoryManager<MultiRegions
350  ::ContField3DHomogeneous1D>
351  ::AllocateSharedPtr(
352  m_session, BkeyZ, m_LhomZ,
353  m_useFFT, m_homogen_dealiasing,
354  m_graph,
355  m_session->GetVariable(i),
356  m_checkIfSystemSingular[i]);
357  }
358  }
359  // Half mode stability analysis
360  else if(m_halfMode)
361  {
362  const LibUtilities::PointsKey PkeyZ(
363  m_npointsZ,
364  LibUtilities::eFourierSingleModeSpaced);
365 
366  const LibUtilities::BasisKey BkeyZR(
367  LibUtilities::eFourierHalfModeRe,
368  m_npointsZ, PkeyZ);
369 
370  const LibUtilities::BasisKey BkeyZI(
371  LibUtilities::eFourierHalfModeIm,
372  m_npointsZ, PkeyZ);
373 
374 
375  for (i = 0; i < m_fields.num_elements(); i++)
376  {
377  if(m_session->GetVariable(i).compare("w")
378  == 0)
379  {
380  m_fields[i] = MemoryManager<MultiRegions
381  ::ContField3DHomogeneous1D>
382  ::AllocateSharedPtr(
383  m_session, BkeyZI, m_LhomZ,
384  m_useFFT,
385  m_homogen_dealiasing,
386  m_graph,
387  m_session->GetVariable(i),
388  m_checkIfSystemSingular[i]);
389  }
390  else
391  {
392  m_fields[i] = MemoryManager<MultiRegions
393  ::ContField3DHomogeneous1D>
394  ::AllocateSharedPtr(
395  m_session, BkeyZR, m_LhomZ,
396  m_useFFT, m_homogen_dealiasing,
397  m_graph,
398  m_session->GetVariable(i),
399  m_checkIfSystemSingular[i]);
400  }
401 
402 
403  }
404  }
405  // Normal homogeneous 1D
406  else
407  {
408  const LibUtilities::PointsKey PkeyZ(
409  m_npointsZ,
410  LibUtilities::eFourierEvenlySpaced);
411  const LibUtilities::BasisKey BkeyZ(
412  LibUtilities::eFourier, m_npointsZ, PkeyZ);
413 
414  for (i = 0; i < m_fields.num_elements(); i++)
415  {
416  m_fields[i] = MemoryManager<MultiRegions
417  ::ContField3DHomogeneous1D>
418  ::AllocateSharedPtr(
419  m_session, BkeyZ, m_LhomZ,
420  m_useFFT, m_homogen_dealiasing,
421  m_graph,
422  m_session->GetVariable(i),
423  m_checkIfSystemSingular[i]);
424  }
425  }
426  }
427  else
428  {
429  i = 0;
430  MultiRegions::ContField2DSharedPtr firstfield;
431  firstfield = MemoryManager<MultiRegions::
432  ContField2D>::AllocateSharedPtr(
433  m_session, m_graph,
434  m_session->GetVariable(i),
435  DeclareCoeffPhysArrays,
436  m_checkIfSystemSingular[0]);
437  m_fields[0] = firstfield;
438  for (i = 1; i < m_fields.num_elements(); i++)
439  {
440  if (m_graph->
441  SameExpansions(m_session->GetVariable(0),
442  m_session->GetVariable(i)))
443  {
444  m_fields[i] = MemoryManager<MultiRegions::
445  ContField2D>::AllocateSharedPtr(
446  *firstfield, m_graph,
447  m_session->GetVariable(i),
448  DeclareCoeffPhysArrays,
449  m_checkIfSystemSingular[i]);
450  }
451  else
452  {
453  m_fields[i] = MemoryManager<MultiRegions
454  ::ContField2D>::AllocateSharedPtr(
455  m_session, m_graph,
456  m_session->GetVariable(i),
457  DeclareCoeffPhysArrays,
458  m_checkIfSystemSingular[i]);
459  }
460  }
461 
462  if (m_projectionType ==
463  MultiRegions::eMixed_CG_Discontinuous)
464  {
465  /// Setting up the normals
466  m_traceNormals =
467  Array<OneD, Array<OneD, NekDouble> >
468  (m_spacedim);
469 
470  for (i = 0; i < m_spacedim; ++i)
471  {
472  m_traceNormals[i] = Array<OneD, NekDouble>
473  (GetTraceNpoints());
474  }
475 
476  m_fields[0]->GetTrace()->
477  GetNormals(m_traceNormals);
478  }
479 
480  }
481 
482  break;
483  }
484  case 3:
485  {
486  i = 0;
487  MultiRegions::ContField3DSharedPtr firstfield =
488  MemoryManager<MultiRegions::ContField3D>
489  ::AllocateSharedPtr(m_session, m_graph,
490  m_session->GetVariable(i),
491  m_checkIfSystemSingular[i]);
492 
493  m_fields[0] = firstfield;
494  for (i = 1; i < m_fields.num_elements(); i++)
495  {
496  if(m_graph->SameExpansions(
497  m_session->GetVariable(0),
498  m_session->GetVariable(i)))
499  {
500  m_fields[i] = MemoryManager<MultiRegions
501  ::ContField3D>::AllocateSharedPtr(
502  *firstfield, m_graph,
503  m_session->GetVariable(i),
504  m_checkIfSystemSingular[i]);
505  }
506  else
507  {
508  m_fields[i] = MemoryManager<MultiRegions
509  ::ContField3D>::AllocateSharedPtr(
510  m_session, m_graph,
511  m_session->GetVariable(i),
512  m_checkIfSystemSingular[i]);
513  }
514  }
515 
516  if (m_projectionType ==
517  MultiRegions::eMixed_CG_Discontinuous)
518  {
519  /// Setting up the normals
520  m_traceNormals =
521  Array<OneD, Array<OneD, NekDouble> >
522  (m_spacedim);
523  for(i = 0; i < m_spacedim; ++i)
524  {
525  m_traceNormals[i] =
526  Array<OneD, NekDouble> (GetTraceNpoints());
527  }
528 
529  m_fields[0]->GetTrace()->GetNormals(m_traceNormals);
530  // Call the trace on all fields to ensure DG setup.
531  for(i = 1; i < m_fields.num_elements(); ++i)
532  {
533  m_fields[i]->GetTrace();
534  }
535  }
536  break;
537  }
538  default:
539  ASSERTL0(false,"Expansion dimension not recognised");
540  break;
541  }
542  }
543  // Discontinuous field
544  else
545  {
546  switch(m_expdim)
547  {
548  case 1:
549  {
550  if (m_HomogeneousType == eHomogeneous2D
551  || m_HomogeneousType == eHomogeneous3D)
552  {
553  const LibUtilities::PointsKey PkeyY(
554  m_npointsY, LibUtilities::eFourierEvenlySpaced);
555  const LibUtilities::BasisKey BkeyY(
556  LibUtilities::eFourier, m_npointsY, PkeyY);
557  const LibUtilities::PointsKey PkeyZ(
558  m_npointsZ, LibUtilities::eFourierEvenlySpaced);
559  const LibUtilities::BasisKey BkeyZ(
560  LibUtilities::eFourier, m_npointsZ, PkeyZ);
561 
562  for (i = 0; i < m_fields.num_elements(); i++)
563  {
564  m_fields[i] = MemoryManager<MultiRegions
565  ::DisContField3DHomogeneous2D>
566  ::AllocateSharedPtr(
567  m_session, BkeyY, BkeyZ, m_LhomY,
568  m_LhomZ, m_useFFT,
569  m_homogen_dealiasing, m_graph,
570  m_session->GetVariable(i));
571  }
572  }
573  else
574  {
575  for (i = 0; i < m_fields.num_elements(); i++)
576  {
577  m_fields[i] = MemoryManager<MultiRegions::
578  DisContField1D>::AllocateSharedPtr(
579  m_session, m_graph,
580  m_session->GetVariable(i));
581  }
582  }
583 
584  break;
585  }
586  case 2:
587  {
588  if(m_HomogeneousType == eHomogeneous1D)
589  {
590  const LibUtilities::PointsKey PkeyZ(
591  m_npointsZ,LibUtilities::eFourierEvenlySpaced);
592  const LibUtilities::BasisKey BkeyZ(
593  LibUtilities::eFourier, m_npointsZ,PkeyZ);
594 
595  for (i = 0; i < m_fields.num_elements(); i++)
596  {
597  m_fields[i] = MemoryManager<MultiRegions
598  ::DisContField3DHomogeneous1D>
599  ::AllocateSharedPtr(
600  m_session, BkeyZ, m_LhomZ, m_useFFT,
601  m_homogen_dealiasing, m_graph,
602  m_session->GetVariable(i));
603  }
604  }
605  else
606  {
607  for (i = 0; i < m_fields.num_elements(); i++)
608  {
609  m_fields[i] = MemoryManager<MultiRegions::
610  DisContField2D>::AllocateSharedPtr(
611  m_session, m_graph,
612  m_session->GetVariable(i));
613  }
614  }
615 
616  break;
617  }
618  case 3:
619  {
620  if (m_HomogeneousType == eHomogeneous3D)
621  {
622  ASSERTL0(false,
623  "3D fully periodic problems not implemented yet");
624  }
625  else
626  {
627  for (i = 0; i < m_fields.num_elements(); i++)
628  {
629  m_fields[i] = MemoryManager<MultiRegions::
630  DisContField3D>::AllocateSharedPtr(
631  m_session, m_graph,
632  m_session->GetVariable(i));
633  }
634  }
635  break;
636  }
637  default:
638  ASSERTL0(false, "Expansion dimension not recognised");
639  break;
640  }
641 
642  // Setting up the normals
643  m_traceNormals =
644  Array<OneD, Array<OneD, NekDouble> >(m_spacedim);
645 
646  for (i = 0; i < m_spacedim; ++i)
647  {
648  m_traceNormals[i] =
649  Array<OneD, NekDouble> (GetTraceNpoints(), 0.0);
650  }
651 
652  m_fields[0]->GetTrace()->GetNormals(m_traceNormals);
653  }
654 
655  // Set Default Parameter
656  m_session->LoadParameter("Time", m_time, 0.0);
657  m_session->LoadParameter("TimeStep", m_timestep, 0.01);
658  m_session->LoadParameter("NumSteps", m_steps, 0);
659  m_session->LoadParameter("IO_CheckSteps", m_checksteps, 0);
660  m_session->LoadParameter("IO_CheckTime", m_checktime, 0.0);
661  m_session->LoadParameter("FinTime", m_fintime, 0);
662  m_session->LoadParameter("NumQuadPointsError",
663  m_NumQuadPointsError, 0);
664 
665  m_nchk = 0;
666 
667  // Zero all physical fields initially
668  ZeroPhysFields();
669  }
670 
671  /**
672  * @brief Destructor for class EquationSystem.
673  */
674  EquationSystem::~EquationSystem()
675  {
676  LibUtilities::NekManager<LocalRegions::MatrixKey,
677  DNekScalMat, LocalRegions::MatrixKey::opLess>::ClearManager();
678  LibUtilities::NekManager<LocalRegions::MatrixKey,
679  DNekScalBlkMat, LocalRegions::MatrixKey::opLess>::ClearManager();
680  }
681 
682  /**
683  * Evaluates a physical function at each quadrature point in the domain.
684  *
685  * @param pArray The array into which to write the values.
686  * @param pEqn The equation to evaluate.
687  */
688  void EquationSystem::EvaluateFunction(
689  Array<OneD, Array<OneD, NekDouble> >& pArray,
690  std::string pFunctionName,
691  const NekDouble pTime,
692  const int domain)
693  {
694  ASSERTL0(m_session->DefinesFunction(pFunctionName),
695  "Function '" + pFunctionName + "' does not exist.");
696 
697  std::vector<std::string> vFieldNames = m_session->GetVariables();
698 
699  for(int i = 0 ; i < vFieldNames.size(); i++)
700  {
701  EvaluateFunction(vFieldNames[i], pArray[i], pFunctionName,
702  pTime, domain);
703  }
704  }
705 
706  /**
707  * Populates a forcing function for each of the dependent variables
708  * using the expression provided by the BoundaryConditions object.
709  * @param force Array of fields to assign forcing.
710  */
711  void EquationSystem::EvaluateFunction(
712  std::vector<std::string> pFieldNames,
713  Array<OneD, Array<OneD, NekDouble> > &pFields,
714  const std::string &pFunctionName,
715  const NekDouble &pTime,
716  const int domain)
717  {
718  ASSERTL1(pFieldNames.size() == pFields.num_elements(),
719  "Function '" + pFunctionName +
720  "' variable list size mismatch with array storage.");
721  ASSERTL0(m_session->DefinesFunction(pFunctionName),
722  "Function '" + pFunctionName + "' does not exist.");
723 
724  for (int i = 0; i < pFieldNames.size(); i++)
725  {
726  EvaluateFunction(pFieldNames[i], pFields[i], pFunctionName, pTime,
727  domain);
728  }
729  }
730 
731  /**
732  * Populates a function for each of the dependent variables using
733  * the expression or filenames provided by the SessionReader object.
734  * @param force Array of fields to assign forcing.
735  */
736  void EquationSystem::EvaluateFunction(
737  std::vector<std::string> pFieldNames,
738  Array<OneD, MultiRegions::ExpListSharedPtr> &pFields,
739  const std::string &pFunctionName,
740  const NekDouble &pTime,
741  const int domain)
742  {
743  ASSERTL0(m_session->DefinesFunction(pFunctionName),
744  "Function '" + pFunctionName + "' does not exist.");
745  ASSERTL0(pFieldNames.size() == pFields.num_elements(),
746  "Field list / name list size mismatch.");
747 
748  for (int i = 0; i < pFieldNames.size(); i++)
749  {
750  EvaluateFunction(pFieldNames[i], pFields[i]->UpdatePhys(),
751  pFunctionName, pTime, domain);
752  pFields[i]->FwdTrans_IterPerExp(pFields[i]->GetPhys(),
753  pFields[i]->UpdateCoeffs());
754  }
755  }
756 
757  void EquationSystem::EvaluateFunction(std::string pFieldName,
758  Array<OneD, NekDouble> &pArray,
759  const std::string &pFunctionName,
760  const NekDouble &pTime,
761  const int domain)
762  {
763  ASSERTL0(m_session->DefinesFunction(pFunctionName),
764  "Function '" + pFunctionName + "' does not exist.");
765 
766  LibUtilities::FunctionType vType;
767  vType = m_session->GetFunctionType(pFunctionName, pFieldName, domain);
768  if (vType == LibUtilities::eFunctionTypeExpression)
769  {
770  EvaluateFunctionExp(pFieldName, pArray, pFunctionName, pTime, domain);
771  }
772  else if (vType == LibUtilities::eFunctionTypeFile ||
773  vType == LibUtilities::eFunctionTypeTransientFile)
774  {
775  std::string filename =
776  m_session->GetFunctionFilename(pFunctionName, pFieldName, domain);
777 
778  if (boost::filesystem::path(filename).extension() != ".pts")
779  {
780  EvaluateFunctionFld(pFieldName, pArray, pFunctionName, pTime,
781  domain);
782  }
783  else
784  {
785  EvaluateFunctionPts(pFieldName, pArray, pFunctionName, pTime,
786  domain);
787  }
788  }
789  }
790 
791  void EquationSystem::EvaluateFunctionExp(string pFieldName,
792  Array<OneD, NekDouble> &pArray,
793  const string &pFunctionName,
794  const NekDouble &pTime,
795  const int domain)
796  {
797  ASSERTL0(m_session->DefinesFunction(pFunctionName),
798  "Function '" + pFunctionName + "' does not exist.");
799 
800  unsigned int nq = m_fields[0]->GetNpoints();
801  if (pArray.num_elements() < nq)
802  {
803  pArray = Array<OneD, NekDouble>(nq);
804  }
805 
806  LibUtilities::FunctionType vType;
807  vType = m_session->GetFunctionType(pFunctionName, pFieldName, domain);
808  ASSERTL0(vType == LibUtilities::eFunctionTypeExpression,
809  "vType not eFunctionTypeExpression");
810 
811  Array<OneD, NekDouble> x0(nq);
812  Array<OneD, NekDouble> x1(nq);
813  Array<OneD, NekDouble> x2(nq);
814 
815  // Get the coordinates (assuming all fields have the same
816  // discretisation)
817  m_fields[0]->GetCoords(x0, x1, x2);
818  LibUtilities::EquationSharedPtr ffunc =
819  m_session->GetFunction(pFunctionName, pFieldName, domain);
820 
821  ffunc->Evaluate(x0, x1, x2, pTime, pArray);
822  }
823 
824  void EquationSystem::EvaluateFunctionFld(string pFieldName,
825  Array<OneD, NekDouble> &pArray,
826  const string &pFunctionName,
827  const NekDouble &pTime,
828  const int domain)
829  {
830 
831  ASSERTL0(m_session->DefinesFunction(pFunctionName),
832  "Function '" + pFunctionName + "' does not exist.");
833 
834  unsigned int nq = m_fields[0]->GetNpoints();
835  if (pArray.num_elements() < nq)
836  {
837  pArray = Array<OneD, NekDouble>(nq);
838  }
839 
840  LibUtilities::FunctionType vType;
841  vType = m_session->GetFunctionType(pFunctionName, pFieldName, domain);
842  ASSERTL0(vType == LibUtilities::eFunctionTypeFile ||
843  vType == LibUtilities::eFunctionTypeTransientFile,
844  "vType not eFunctionTypeFile or eFunctionTypeTransientFile");
845 
846  std::string filename =
847  m_session->GetFunctionFilename(pFunctionName, pFieldName, domain);
848  std::string fileVar = m_session->GetFunctionFilenameVariable(
849  pFunctionName, pFieldName, domain);
850 
851  if (fileVar.length() == 0)
852  {
853  fileVar = pFieldName;
854  }
855 
856  // In case of eFunctionTypeTransientFile, generate filename from
857  // format string
858  if (vType == LibUtilities::eFunctionTypeTransientFile)
859  {
860  try
861  {
862 #if (defined _WIN32 && _MSC_VER < 1900)
863  // We need this to make sure boost::format has always
864  // two digits in the exponents of Scientific notation.
865  unsigned int old_exponent_format;
866  old_exponent_format = _set_output_format(_TWO_DIGIT_EXPONENT);
867  filename = boost::str(boost::format(filename) % pTime);
868  _set_output_format(old_exponent_format);
869 #else
870  filename = boost::str(boost::format(filename) % pTime);
871 #endif
872  }
873  catch (...)
874  {
875  ASSERTL0(false, "Invalid Filename in function \"" + pFunctionName +
876  "\", variable \"" + fileVar + "\"")
877  }
878  }
879 
880  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef;
881  std::vector<std::vector<NekDouble> > FieldData;
882 
883  // Define list of global element ids
884  int numexp = m_fields[0]->GetExpSize();
885  Array<OneD, int> ElementGIDs(numexp);
886  for (int i = 0; i < numexp; ++i)
887  {
888  ElementGIDs[i] = m_fields[0]->GetExp(i)->GetGeom()->GetGlobalID();
889  }
890 
891  // check if we already loaded this file. For transient files,
892  // funcFilename != filename so we can make sure we only keep the
893  // latest field per funcFilename.
894  std::string funcFilename =
895  m_session->GetFunctionFilename(pFunctionName, pFieldName, domain);
896  if (m_loadedFldFields.find(funcFilename) != m_loadedFldFields.end())
897  {
898  // found
899  if (m_loadedFldFields[funcFilename].first == filename)
900  {
901  // found
902  FieldDef = m_loadedFldFields[funcFilename].second.fieldDef;
903  FieldData = m_loadedFldFields[funcFilename].second.fieldData;
904  }
905  else
906  {
907  LibUtilities::FieldIOSharedPtr pts_fld =
908  LibUtilities::FieldIO::CreateForFile(m_session, filename);
909  pts_fld->Import(filename, FieldDef, FieldData,
910  LibUtilities::NullFieldMetaDataMap,
911  ElementGIDs);
912  }
913  }
914  else
915  {
916  LibUtilities::FieldIOSharedPtr pts_fld =
917  LibUtilities::FieldIO::CreateForFile(m_session, filename);
918  pts_fld->Import(filename, FieldDef, FieldData,
919  LibUtilities::NullFieldMetaDataMap,
920  ElementGIDs);
921  }
922  // Now we overwrite the field we had in
923  // m_loadedFldFields[funcFilename] before, making sure we only keep
924  // one field per funcFilename in memory
925  m_loadedFldFields[funcFilename].first = filename;
926  m_loadedFldFields[funcFilename].second.fieldDef = FieldDef;
927  m_loadedFldFields[funcFilename].second.fieldData = FieldData;
928 
929  int idx = -1;
930  Array<OneD, NekDouble> vCoeffs(m_fields[0]->GetNcoeffs(), 0.0);
931  // Loop over all the expansions
932  for (int i = 0; i < FieldDef.size(); ++i)
933  {
934  // Find the index of the required field in the
935  // expansion segment
936  for (int j = 0; j < FieldDef[i]->m_fields.size(); ++j)
937  {
938  if (FieldDef[i]->m_fields[j] == fileVar)
939  {
940  idx = j;
941  }
942  }
943 
944  if (idx >= 0)
945  {
946  m_fields[0]->ExtractDataToCoeffs(
947  FieldDef[i], FieldData[i], FieldDef[i]->m_fields[idx], vCoeffs);
948  }
949  else
950  {
951  cout << "Field " + fileVar + " not found." << endl;
952  }
953  }
954 
955  m_fields[0]->BwdTrans_IterPerExp(vCoeffs, pArray);
956  }
957 
958  void EquationSystem::EvaluateFunctionPts(string pFieldName,
959  Array<OneD, NekDouble> &pArray,
960  const string &pFunctionName,
961  const NekDouble &pTime,
962  const int domain)
963  {
964 
965  ASSERTL0(m_session->DefinesFunction(pFunctionName),
966  "Function '" + pFunctionName + "' does not exist.");
967 
968  unsigned int nq = m_fields[0]->GetNpoints();
969  if (pArray.num_elements() < nq)
970  {
971  pArray = Array<OneD, NekDouble>(nq);
972  }
973 
974  LibUtilities::FunctionType vType;
975  vType = m_session->GetFunctionType(pFunctionName, pFieldName, domain);
976  ASSERTL0(vType == LibUtilities::eFunctionTypeFile ||
977  vType == LibUtilities::eFunctionTypeTransientFile,
978  "vType not eFunctionTypeFile or eFunctionTypeTransientFile");
979 
980  std::string filename =
981  m_session->GetFunctionFilename(pFunctionName, pFieldName, domain);
982  std::string fileVar = m_session->GetFunctionFilenameVariable(
983  pFunctionName, pFieldName, domain);
984 
985  if (fileVar.length() == 0)
986  {
987  fileVar = pFieldName;
988  }
989 
990  // In case of eFunctionTypeTransientFile, generate filename from
991  // format string
992  if (vType == LibUtilities::eFunctionTypeTransientFile)
993  {
994  try
995  {
996 #if (defined _WIN32 && _MSC_VER < 1900)
997  // We need this to make sure boost::format has always
998  // two digits in the exponents of Scientific notation.
999  unsigned int old_exponent_format;
1000  old_exponent_format = _set_output_format(_TWO_DIGIT_EXPONENT);
1001  filename = boost::str(boost::format(filename) % pTime);
1002  _set_output_format(old_exponent_format);
1003 #else
1004  filename = boost::str(boost::format(filename) % pTime);
1005 #endif
1006  }
1007  catch (...)
1008  {
1009  ASSERTL0(false, "Invalid Filename in function \"" + pFunctionName +
1010  "\", variable \"" + fileVar + "\"")
1011  }
1012  }
1013 
1014  LibUtilities::PtsFieldSharedPtr outPts;
1015  // check if we already loaded this file. For transient files,
1016  // funcFilename != filename so we can make sure we only keep the
1017  // latest pts field per funcFilename.
1018  std::string funcFilename =
1019  m_session->GetFunctionFilename(pFunctionName, pFieldName, domain);
1020 
1021  if (m_loadedPtsFields.find(funcFilename) != m_loadedPtsFields.end())
1022  {
1023  // found
1024  if (m_loadedPtsFields[funcFilename].first == filename)
1025  {
1026  // found
1027  outPts = m_loadedPtsFields[funcFilename].second;
1028  }
1029  else
1030  {
1031  LoadPts(funcFilename, filename, outPts);
1032  }
1033  }
1034  else
1035  {
1036  LoadPts(funcFilename, filename, outPts);
1037  }
1038  // Now we overwrite the field we had in
1039  // m_loadedPtsFields[funcFilename] before, making sure we only keep
1040  // one field per funcFilename in memory
1041  m_loadedPtsFields[funcFilename].first = filename;
1042  m_loadedPtsFields[funcFilename].second = outPts;
1043 
1044  int fieldInd;
1045  vector<string> fieldNames = outPts->GetFieldNames();
1046  for (fieldInd = 0; fieldInd < fieldNames.size(); ++fieldInd)
1047  {
1048  if (outPts->GetFieldName(fieldInd) == pFieldName)
1049  {
1050  break;
1051  }
1052  }
1053  ASSERTL0(fieldInd != fieldNames.size(), "field not found");
1054 
1055  pArray = outPts->GetPts(fieldInd + outPts->GetDim());
1056  }
1057 
1058  void EquationSystem::LoadPts(std::string funcFilename,
1059  std::string filename,
1060  LibUtilities::PtsFieldSharedPtr &outPts)
1061  {
1062  unsigned int nq = m_fields[0]->GetNpoints();
1063 
1064  LibUtilities::PtsFieldSharedPtr inPts;
1065  LibUtilities::PtsIO ptsIO(m_session->GetComm());
1066  ptsIO.Import(filename, inPts);
1067 
1068  Array<OneD, Array<OneD, NekDouble> > pts(inPts->GetDim() +
1069  inPts->GetNFields());
1070  for (int i = 0; i < inPts->GetDim() + inPts->GetNFields(); ++i)
1071  {
1072  pts[i] = Array<OneD, NekDouble>(nq);
1073  }
1074  if (inPts->GetDim() == 1)
1075  {
1076  m_fields[0]->GetCoords(pts[0]);
1077  }
1078  else if (inPts->GetDim() == 2)
1079  {
1080  m_fields[0]->GetCoords(pts[0], pts[1]);
1081  }
1082  else if (inPts->GetDim() == 3)
1083  {
1084  m_fields[0]->GetCoords(pts[0], pts[1], pts[2]);
1085  }
1086  outPts = MemoryManager<LibUtilities::PtsField>::AllocateSharedPtr(
1087  inPts->GetDim(), inPts->GetFieldNames(), pts);
1088 
1089  // check if we already have an interolator for this funcFilename
1090  if (m_interpolators.find(funcFilename) == m_interpolators.end())
1091  {
1092  m_interpolators[funcFilename] =
1093  FieldUtils::Interpolator(Nektar::FieldUtils::eShepard);
1094  if (m_comm->GetRank() == 0)
1095  {
1096  m_interpolators[funcFilename].SetProgressCallback(
1097  &EquationSystem::PrintProgressbar, this);
1098  }
1099  m_interpolators[funcFilename].CalcWeights(inPts, outPts);
1100  if (m_comm->GetRank() == 0)
1101  {
1102  cout << endl;
1103  if (GetSession()->DefinesCmdLineArgument("verbose"))
1104  {
1105  m_interpolators[funcFilename].PrintStatistics();
1106  }
1107  }
1108  }
1109  m_interpolators[funcFilename].Interpolate(inPts, outPts);
1110  }
1111 
1112 
1113  /**
1114  * @brief Provide a description of a function for a given field name.
1115  *
1116  * @param pFieldName Field name.
1117  * @param pFunctionName Function name.
1118  */
1119  std::string EquationSystem::DescribeFunction(
1120  std::string pFieldName,
1121  const std::string &pFunctionName,
1122  const int domain)
1123  {
1124  ASSERTL0(m_session->DefinesFunction(pFunctionName),
1125  "Function '" + pFunctionName + "' does not exist.");
1126 
1127  std::string retVal;
1128  LibUtilities::FunctionType vType;
1129 
1130  vType = m_session->GetFunctionType(pFunctionName, pFieldName);
1131  if (vType == LibUtilities::eFunctionTypeExpression)
1132  {
1133  LibUtilities::EquationSharedPtr ffunc
1134  = m_session->GetFunction(pFunctionName, pFieldName,domain);
1135  retVal = ffunc->GetExpression();
1136  }
1137  else if (vType == LibUtilities::eFunctionTypeFile)
1138  {
1139  std::string filename
1140  = m_session->GetFunctionFilename(pFunctionName, pFieldName,domain);
1141  retVal = "from file " + filename;
1142  }
1143 
1144  return retVal;
1145  }
1146 
1147  /**
1148  * If boundary conditions are time-dependent, they will be evaluated at
1149  * the time specified.
1150  * @param time The time at which to evaluate the BCs
1151  */
1152  void EquationSystem::SetBoundaryConditions(NekDouble time)
1153  {
1154  std::string varName;
1155  int nvariables = m_fields.num_elements();
1156  for (int i = 0; i < nvariables; ++i)
1157  {
1158  varName = m_session->GetVariable(i);
1159  m_fields[i]->EvaluateBoundaryConditions(time, varName);
1160  }
1161  }
1162 
1163  /**
1164  * Compute the error in the L2-norm.
1165  * @param field The field to compare.
1166  * @param exactsoln The exact solution to compare with.
1167  * @param Normalised Normalise L2-error.
1168  * @returns Error in the L2-norm.
1169  */
1170  NekDouble EquationSystem::v_L2Error(
1171  unsigned int field,
1172  const Array<OneD, NekDouble> &exactsoln,
1173  bool Normalised)
1174  {
1175  NekDouble L2error = -1.0;
1176 
1177  if (m_NumQuadPointsError == 0)
1178  {
1179  if (m_fields[field]->GetPhysState() == false)
1180  {
1181  m_fields[field]->BwdTrans(m_fields[field]->GetCoeffs(),
1182  m_fields[field]->UpdatePhys());
1183  }
1184 
1185  if (exactsoln.num_elements())
1186  {
1187  L2error = m_fields[field]->L2(m_fields[field]->GetPhys(), exactsoln);
1188  }
1189  else if (m_session->DefinesFunction("ExactSolution"))
1190  {
1191  Array<OneD, NekDouble>
1192  exactsoln(m_fields[field]->GetNpoints());
1193 
1194  EvaluateFunction(m_session->GetVariable(field), exactsoln,
1195  "ExactSolution", m_time);
1196 
1197  L2error = m_fields[field]->L2(m_fields[field]->GetPhys(), exactsoln);
1198  }
1199  else
1200  {
1201  L2error = m_fields[field]->L2(m_fields[field]->GetPhys());
1202  }
1203 
1204  if (Normalised == true)
1205  {
1206  Array<OneD, NekDouble> one(m_fields[field]->GetNpoints(),
1207  1.0);
1208 
1209  NekDouble Vol = m_fields[field]->PhysIntegral(one);
1210  m_comm->AllReduce(Vol, LibUtilities::ReduceSum);
1211 
1212  L2error = sqrt(L2error*L2error/Vol);
1213  }
1214  }
1215  else
1216  {
1217  Array<OneD,NekDouble> L2INF(2);
1218  L2INF = ErrorExtraPoints(field);
1219  L2error = L2INF[0];
1220  }
1221  return L2error;
1222  }
1223 
1224  /**
1225  * Compute the error in the L_inf-norm
1226  * @param field The field to compare.
1227  * @param exactsoln The exact solution to compare with.
1228  * @returns Error in the L_inft-norm.
1229  */
1230  NekDouble EquationSystem::v_LinfError (
1231  unsigned int field,
1232  const Array<OneD, NekDouble> &exactsoln)
1233  {
1234  NekDouble Linferror = -1.0;
1235 
1236  if (m_NumQuadPointsError == 0)
1237  {
1238  if (m_fields[field]->GetPhysState() == false)
1239  {
1240  m_fields[field]->BwdTrans(m_fields[field]->GetCoeffs(),
1241  m_fields[field]->UpdatePhys());
1242  }
1243 
1244  if (exactsoln.num_elements())
1245  {
1246  Linferror = m_fields[field]->Linf(m_fields[field]->GetPhys(), exactsoln);
1247  }
1248  else if (m_session->DefinesFunction("ExactSolution"))
1249  {
1250  Array<OneD, NekDouble>
1251  exactsoln(m_fields[field]->GetNpoints());
1252 
1253  EvaluateFunction(m_session->GetVariable(field), exactsoln,
1254  "ExactSolution", m_time);
1255 
1256  Linferror = m_fields[field]->Linf(m_fields[field]->GetPhys(), exactsoln);
1257  }
1258  else
1259  {
1260  Linferror = m_fields[field]->Linf(m_fields[field]->GetPhys());
1261  }
1262  }
1263  else
1264  {
1265  Array<OneD,NekDouble> L2INF(2);
1266  L2INF = ErrorExtraPoints(field);
1267  Linferror = L2INF[1];
1268  }
1269 
1270  return Linferror;
1271  }
1272 
1273  /**
1274  * Compute the error in the L2-norm, L-inf for a larger number of
1275  * quadrature points.
1276  * @param field The field to compare.
1277  * @returns Error in the L2-norm and L-inf norm.
1278  */
1279  Array<OneD,NekDouble> EquationSystem::ErrorExtraPoints(
1280  unsigned int field)
1281  {
1282  int NumModes = GetNumExpModes();
1283  Array<OneD,NekDouble> L2INF(2);
1284 
1285  const LibUtilities::PointsKey PkeyT1(
1286  m_NumQuadPointsError,LibUtilities::eGaussLobattoLegendre);
1287  const LibUtilities::PointsKey PkeyT2(
1288  m_NumQuadPointsError, LibUtilities::eGaussRadauMAlpha1Beta0);
1289  const LibUtilities::PointsKey PkeyQ1(
1290  m_NumQuadPointsError, LibUtilities::eGaussLobattoLegendre);
1291  const LibUtilities::PointsKey PkeyQ2(
1292  m_NumQuadPointsError, LibUtilities::eGaussLobattoLegendre);
1293  const LibUtilities::BasisKey BkeyT1(
1294  LibUtilities::eModified_A,NumModes, PkeyT1);
1295  const LibUtilities::BasisKey BkeyT2(
1296  LibUtilities::eModified_B, NumModes, PkeyT2);
1297  const LibUtilities::BasisKey BkeyQ1(
1298  LibUtilities::eModified_A, NumModes, PkeyQ1);
1299  const LibUtilities::BasisKey BkeyQ2(
1300  LibUtilities::eModified_A, NumModes, PkeyQ2);
1301 
1302  MultiRegions::ExpList2DSharedPtr ErrorExp =
1303  MemoryManager<MultiRegions::ExpList2D>::AllocateSharedPtr(
1304  m_session, BkeyT1, BkeyT2, BkeyQ1, BkeyQ2, m_graph);
1305 
1306  int ErrorCoordim = ErrorExp->GetCoordim(0);
1307  int ErrorNq = ErrorExp->GetTotPoints();
1308 
1309  Array<OneD,NekDouble> ErrorXc0(ErrorNq, 0.0);
1310  Array<OneD,NekDouble> ErrorXc1(ErrorNq, 0.0);
1311  Array<OneD,NekDouble> ErrorXc2(ErrorNq, 0.0);
1312 
1313  switch(ErrorCoordim)
1314  {
1315  case 1:
1316  ErrorExp->GetCoords(ErrorXc0);
1317  break;
1318  case 2:
1319  ErrorExp->GetCoords(ErrorXc0, ErrorXc1);
1320  break;
1321  case 3:
1322  ErrorExp->GetCoords(ErrorXc0, ErrorXc1, ErrorXc2);
1323  break;
1324  }
1325  LibUtilities::EquationSharedPtr exSol =
1326  m_session->GetFunction("ExactSolution", field);
1327 
1328  // Evaluate the exact solution
1329  Array<OneD,NekDouble> ErrorSol(ErrorNq);
1330 
1331  exSol->Evaluate(ErrorXc0,ErrorXc1,ErrorXc2,m_time,ErrorSol);
1332 
1333  // Calcualte spectral/hp approximation on the quadrature points
1334  // of this new expansion basis
1335  ErrorExp->BwdTrans_IterPerExp(m_fields[field]->GetCoeffs(),
1336  ErrorExp->UpdatePhys());
1337 
1338  L2INF[0] = ErrorExp->L2 (ErrorExp->GetPhys(), ErrorSol);
1339  L2INF[1] = ErrorExp->Linf(ErrorExp->GetPhys(), ErrorSol);
1340 
1341  return L2INF;
1342  }
1343 
1344 
1345  /**
1346  * Set the physical fields based on a restart file, or a function
1347  * describing the initial condition given in the session.
1348  * @param initialtime Time at which to evaluate the function.
1349  * @param dumpInitialConditions Write the initial condition to file?
1350  */
1351  void EquationSystem::v_SetInitialConditions(NekDouble initialtime,
1352  bool dumpInitialConditions,
1353  const int domain)
1354  {
1355  if (m_session->GetComm()->GetRank() == 0)
1356  {
1357  cout << "Initial Conditions:" << endl;
1358  }
1359 
1360  if (m_session->DefinesFunction("InitialConditions"))
1361  {
1362  EvaluateFunction(m_session->GetVariables(), m_fields,
1363  "InitialConditions", m_time, domain);
1364 
1365  if (m_session->GetComm()->GetRank() == 0)
1366  {
1367 
1368  for (int i = 0; i < m_fields.num_elements(); ++i)
1369  {
1370  std::string varName = m_session->GetVariable(i);
1371  cout << " - Field " << varName << ": "
1372  << DescribeFunction(varName, "InitialConditions",domain)
1373  << endl;
1374  }
1375  }
1376  }
1377  else
1378  {
1379  int nq = m_fields[0]->GetNpoints();
1380  for (int i = 0; i < m_fields.num_elements(); i++)
1381  {
1382  Vmath::Zero(nq, m_fields[i]->UpdatePhys(), 1);
1383  m_fields[i]->SetPhysState(true);
1384  Vmath::Zero(m_fields[i]->GetNcoeffs(),
1385  m_fields[i]->UpdateCoeffs(), 1);
1386  if (m_session->GetComm()->GetRank() == 0)
1387  {
1388  cout << " - Field " << m_session->GetVariable(i)
1389  << ": 0 (default)" << endl;
1390  }
1391  }
1392 
1393  }
1394 
1395  if (dumpInitialConditions && m_checksteps)
1396  {
1397  Checkpoint_Output(m_nchk);
1398  m_nchk++;
1399  }
1400  }
1401 
1402  void EquationSystem::v_EvaluateExactSolution(
1403  unsigned int field,
1404  Array<OneD, NekDouble> &outfield,
1405  const NekDouble time)
1406  {
1407  ASSERTL0 (outfield.num_elements() == m_fields[field]->GetNpoints(),
1408  "ExactSolution array size mismatch.");
1409  Vmath::Zero(outfield.num_elements(), outfield, 1);
1410  if (m_session->DefinesFunction("ExactSolution"))
1411  {
1412  EvaluateFunction(m_session->GetVariable(field), outfield,
1413  "ExactSolution", time);
1414  }
1415  }
1416 
1417 
1418  /**
1419  * By default, nothing needs initialising at the EquationSystem level.
1420  */
1421  void EquationSystem::v_DoInitialise()
1422  {
1423 
1424  }
1425 
1426 
1427  void EquationSystem::InitialiseBaseFlow(
1428  Array<OneD, Array<OneD, NekDouble> > &base)
1429  {
1430  base = Array<OneD, Array<OneD, NekDouble> >(m_spacedim);
1431  std::vector<std::string> vel;
1432  vel.push_back("Vx");
1433  vel.push_back("Vy");
1434  vel.push_back("Vz");
1435  vel.resize(m_spacedim);
1436  SetUpBaseFields(m_graph);
1437  EvaluateFunction(vel, base, "BaseFlow");
1438  }
1439 
1440  void EquationSystem::SetUpBaseFields(
1441  SpatialDomains::MeshGraphSharedPtr &mesh)
1442  {
1443  int i;
1444 
1445  // The number of variables can be different from the dimension
1446  // of the base flow
1447  int nvariables = m_session->GetVariables().size();
1448  m_base = Array<OneD, MultiRegions::ExpListSharedPtr>(nvariables);
1449  if (m_projectionType == MultiRegions::eGalerkin ||
1450  m_projectionType == MultiRegions::eMixed_CG_Discontinuous)
1451  {
1452  switch (m_expdim)
1453  {
1454  case 1:
1455  {
1456  for(i = 0; i < m_base.num_elements(); i++)
1457  {
1458  m_base[i] = MemoryManager<MultiRegions::ContField1D>
1459  ::AllocateSharedPtr(m_session, mesh,
1460  m_session->GetVariable(0));
1461  }
1462  }
1463  break;
1464  case 2:
1465  {
1466  if (m_HomogeneousType == eHomogeneous1D)
1467  {
1468  if (m_singleMode)
1469  {
1470  const LibUtilities::PointsKey PkeyZ(m_npointsZ,
1471  LibUtilities::eFourierSingleModeSpaced);
1472  const LibUtilities::BasisKey BkeyZ(
1473  LibUtilities::eFourier,
1474  m_npointsZ,PkeyZ);
1475 
1476  for (i = 0 ; i < m_base.num_elements(); i++)
1477  {
1478  m_base[i] = MemoryManager<MultiRegions
1479  ::ContField3DHomogeneous1D>
1480  ::AllocateSharedPtr(
1481  m_session, BkeyZ, m_LhomZ,
1482  m_useFFT, m_homogen_dealiasing,
1483  m_graph,
1484  m_session->GetVariable(i));
1485  m_base[i]->SetWaveSpace(true);
1486  }
1487  }
1488  else if (m_halfMode)
1489  {
1490  //1 plane field (half mode expansion)
1491  const LibUtilities::PointsKey PkeyZ(m_npointsZ,
1492  LibUtilities::eFourierSingleModeSpaced);
1493  const LibUtilities::BasisKey BkeyZ(
1494  LibUtilities::eFourierHalfModeRe,
1495  m_npointsZ,PkeyZ);
1496 
1497  for (i = 0 ; i < m_base.num_elements(); i++)
1498  {
1499  m_base[i] = MemoryManager<MultiRegions
1500  ::ContField3DHomogeneous1D>
1501  ::AllocateSharedPtr(
1502  m_session, BkeyZ, m_LhomZ,
1503  m_useFFT, m_homogen_dealiasing,
1504  m_graph,
1505  m_session->GetVariable(i));
1506  m_base[i]->SetWaveSpace(true);
1507  }
1508  }
1509  else
1510  {
1511  const LibUtilities::PointsKey PkeyZ(m_npointsZ,
1512  LibUtilities::eFourierEvenlySpaced);
1513  const LibUtilities::BasisKey BkeyZ(
1514  LibUtilities::eFourier,m_npointsZ,
1515  PkeyZ);
1516 
1517  for (i = 0 ; i < m_base.num_elements(); i++)
1518  {
1519  m_base[i] = MemoryManager<MultiRegions
1520  ::ContField3DHomogeneous1D>
1521  ::AllocateSharedPtr(
1522  m_session, BkeyZ, m_LhomZ,
1523  m_useFFT, m_homogen_dealiasing,
1524  m_graph,
1525  m_session->GetVariable(i));
1526  m_base[i]->SetWaveSpace(false);
1527  }
1528  }
1529  }
1530  else
1531  {
1532  i = 0;
1533  MultiRegions::ContField2DSharedPtr firstbase =
1534  MemoryManager<MultiRegions::ContField2D>
1535  ::AllocateSharedPtr(m_session,mesh,
1536  m_session->GetVariable(i));
1537  m_base[0]=firstbase;
1538 
1539  for (i = 1 ; i < m_base.num_elements(); i++)
1540  {
1541  m_base[i] = MemoryManager<MultiRegions::
1542  ContField2D>::AllocateSharedPtr(
1543  *firstbase,mesh,
1544  m_session->GetVariable(i));
1545  }
1546  }
1547  }
1548  break;
1549  case 3:
1550  {
1551  MultiRegions::ContField3DSharedPtr firstbase =
1552  MemoryManager<MultiRegions::ContField3D>
1553  ::AllocateSharedPtr(m_session, m_graph,
1554  m_session->GetVariable(0));
1555  m_base[0] = firstbase;
1556  for (i = 1 ; i < m_base.num_elements(); i++)
1557  {
1558  m_base[i] = MemoryManager<MultiRegions::ContField3D>
1559  ::AllocateSharedPtr(*firstbase, m_graph,
1560  m_session->GetVariable(0));
1561  }
1562  }
1563  break;
1564  default:
1565  ASSERTL0(false,"Expansion dimension not recognised");
1566  break;
1567  }
1568  }
1569  else
1570  {
1571  switch(m_expdim)
1572  {
1573  case 1:
1574  {
1575  // need to use zero for variable as may be more base
1576  // flows than variables
1577  for(i = 0 ; i < m_base.num_elements(); i++)
1578  {
1579  m_base[i] = MemoryManager<MultiRegions
1580  ::DisContField1D>
1581  ::AllocateSharedPtr(m_session, m_graph,
1582  m_session->GetVariable(0));
1583  }
1584  break;
1585  }
1586  case 2:
1587  {
1588  for(i = 0 ; i < m_base.num_elements(); i++)
1589  {
1590  m_base[i] = MemoryManager<MultiRegions::
1591  DisContField2D>::AllocateSharedPtr(
1592  m_session,m_graph,
1593  m_session->GetVariable(0));
1594  }
1595  break;
1596  }
1597  case 3:
1598  ASSERTL0(false, "3 D not set up");
1599  default:
1600  ASSERTL0(false, "Expansion dimension not recognised");
1601  break;
1602  }
1603  }
1604  }
1605 
1606  // Import base flow from file and load in m_base
1607  void EquationSystem::ImportFldBase(
1608  std::string pInfile,
1609  SpatialDomains::MeshGraphSharedPtr pGraph)
1610  {
1611  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef;
1612  std::vector<std::vector<NekDouble> > FieldData;
1613 
1614  //Get Homogeneous
1615  LibUtilities::FieldIOSharedPtr base_fld =
1616  LibUtilities::FieldIO::CreateForFile(m_session, pInfile);
1617  base_fld->Import(pInfile,FieldDef,FieldData);
1618 
1619  int nvar = m_session->GetVariables().size();
1620  if (m_session->DefinesSolverInfo("HOMOGENEOUS"))
1621  {
1622  std::string HomoStr = m_session->GetSolverInfo("HOMOGENEOUS");
1623  }
1624  // copy FieldData into m_fields
1625  for (int j = 0; j < nvar; ++j)
1626  {
1627  for(int i = 0; i < FieldDef.size(); ++i)
1628  {
1629  bool flag = FieldDef[i]->m_fields[j] ==
1630  m_session->GetVariable(j);
1631  ASSERTL0(flag, (std::string("Order of ") + pInfile
1632  + std::string(" data and that defined in "
1633  "the session differs")).c_str());
1634 
1635  m_base[j]->ExtractDataToCoeffs(FieldDef[i], FieldData[i],
1636  FieldDef[i]->m_fields[j],
1637  m_base[j]->UpdateCoeffs());
1638  }
1639  }
1640  }
1641 
1642  /**
1643  *
1644  */
1645  void EquationSystem::v_DoSolve()
1646  {
1647 
1648  }
1649 
1650  /**
1651  * Virtual function to define if operator in DoSolve is
1652  * negated with regard to the strong form. This is currently
1653  * only used in Arnoldi solves. Default is false.
1654  */
1655  bool EquationSystem::v_NegatedOp(void)
1656  {
1657  return false;
1658  }
1659 
1660  /**
1661  *
1662  */
1663  void EquationSystem::v_TransCoeffToPhys()
1664  {
1665 
1666  }
1667 
1668  /**
1669  *
1670  */
1671  void EquationSystem::v_TransPhysToCoeff()
1672  {
1673 
1674  }
1675 
1676 
1677  /// Virtual function for generating summary information.
1678  void EquationSystem::v_GenerateSummary(SummaryList& l)
1679  {
1680  SessionSummary(l);
1681  }
1682 
1683 
1684  /**
1685  * Write the field data to file. The file is named according to the session
1686  * name with the extension .fld appended.
1687  */
1688  void EquationSystem::v_Output(void)
1689  {
1690  WriteFld(m_sessionName + ".fld");
1691  }
1692 
1693  /**
1694  * Zero the physical fields.
1695  */
1696  void EquationSystem::ZeroPhysFields(void)
1697  {
1698  for (int i = 0; i < m_fields.num_elements(); i++)
1699  {
1700  Vmath::Zero(m_fields[i]->GetNpoints(),
1701  m_fields[i]->UpdatePhys(),1);
1702  }
1703  }
1704 
1705  /**
1706  * FwdTrans the m_fields members
1707  */
1708  void EquationSystem::FwdTransFields(void)
1709  {
1710  for (int i = 0; i < m_fields.num_elements(); i++)
1711  {
1712  m_fields[i]->FwdTrans(m_fields[i]->GetPhys(),
1713  m_fields[i]->UpdateCoeffs());
1714  m_fields[i]->SetPhysState(false);
1715  }
1716  }
1717 
1718  /**
1719  * Computes the weak Green form of advection terms (without boundary
1720  * integral), i.e. \f$ (\nabla \phi \cdot F) \f$ where for example
1721  * \f$ F=uV \f$.
1722  * @param F Fields.
1723  * @param outarray Storage for result.
1724  *
1725  * \note Assuming all fields are of the same expansion and order so that
1726  * we can use the parameters of m_fields[0].
1727  */
1728  void EquationSystem::WeakAdvectionGreensDivergenceForm(
1729  const Array<OneD, Array<OneD, NekDouble> > &F,
1730  Array<OneD, NekDouble> &outarray)
1731  {
1732  m_fields[0]->IProductWRTDerivBase(F,outarray);
1733  }
1734 
1735  /**
1736  * Calculate Inner product of the divergence advection form
1737  * \f$(\phi, \nabla \cdot F)\f$, where for example \f$ F = uV \f$.
1738  * @param F Fields.
1739  * @param outarray Storage for result.
1740  */
1741  void EquationSystem::WeakAdvectionDivergenceForm(
1742  const Array<OneD, Array<OneD, NekDouble> > &F,
1743  Array<OneD, NekDouble> &outarray)
1744  {
1745  // Use dimension of Velocity vector to dictate dimension of operation
1746  int ndim = F.num_elements();
1747  int nPointsTot = m_fields[0]->GetNpoints();
1748  Array<OneD, NekDouble> tmp(nPointsTot);
1749  Array<OneD, NekDouble> div(nPointsTot, 0.0);
1750 
1751  // Evaluate the divergence
1752  for (int i = 0; i < ndim; ++i)
1753  {
1754  //m_fields[0]->PhysDeriv(i,F[i],tmp);
1755  m_fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[i],F[i],tmp);
1756  Vmath::Vadd(nPointsTot, tmp, 1, div, 1, div, 1);
1757  }
1758 
1759  m_fields[0]->IProductWRTBase(div, outarray);
1760  }
1761 
1762  /**
1763  * Calculate Inner product of the divergence advection form
1764  * \f$ (\phi, V\cdot \nabla u) \f$
1765  * @param V Fields.
1766  * @param u Fields.
1767  * @param outarray Storage for result.
1768  */
1769  void EquationSystem::WeakAdvectionNonConservativeForm(
1770  const Array<OneD, Array<OneD, NekDouble> > &V,
1771  const Array<OneD, const NekDouble> &u,
1772  Array<OneD, NekDouble> &outarray,
1773  bool UseContCoeffs)
1774  {
1775  // use dimension of Velocity vector to dictate dimension of operation
1776  int ndim = V.num_elements();
1777 
1778  int nPointsTot = m_fields[0]->GetNpoints();
1779  Array<OneD, NekDouble> tmp(nPointsTot);
1780  Array<OneD, NekDouble> wk(ndim * nPointsTot, 0.0);
1781 
1782  AdvectionNonConservativeForm(V, u, tmp, wk);
1783 
1784  if (UseContCoeffs)
1785  {
1786  m_fields[0]->IProductWRTBase(tmp, outarray,
1787  MultiRegions::eGlobal);
1788  }
1789  else
1790  {
1791  m_fields[0]->IProductWRTBase_IterPerExp(tmp, outarray);
1792  }
1793  }
1794 
1795  /**
1796  * Calculate the inner product \f$ V\cdot \nabla u \f$
1797  * @param V Fields.
1798  * @param u Fields.
1799  * @param outarray Storage for result.
1800  * @param wk Workspace.
1801  */
1802  void EquationSystem::AdvectionNonConservativeForm(
1803  const Array<OneD, Array<OneD, NekDouble> > &V,
1804  const Array<OneD, const NekDouble> &u,
1805  Array<OneD, NekDouble> &outarray,
1806  Array<OneD, NekDouble> &wk)
1807  {
1808  // Use dimension of Velocity vector to dictate dimension of operation
1809  int ndim = V.num_elements();
1810  //int ndim = m_expdim;
1811 
1812  // ToDo: here we should add a check that V has right dimension
1813 
1814  int nPointsTot = m_fields[0]->GetNpoints();
1815  Array<OneD, NekDouble> grad0,grad1,grad2;
1816 
1817  // Check to see if wk space is defined
1818  if (wk.num_elements())
1819  {
1820  grad0 = wk;
1821  }
1822  else
1823  {
1824  grad0 = Array<OneD, NekDouble> (nPointsTot);
1825  }
1826 
1827  // Evaluate V\cdot Grad(u)
1828  switch(ndim)
1829  {
1830  case 1:
1831  m_fields[0]->PhysDeriv(u,grad0);
1832  Vmath::Vmul(nPointsTot, grad0, 1, V[0], 1, outarray,1);
1833  break;
1834  case 2:
1835  grad1 = Array<OneD, NekDouble> (nPointsTot);
1836  m_fields[0]->PhysDeriv(u, grad0, grad1);
1837  Vmath::Vmul (nPointsTot, grad0, 1, V[0], 1, outarray, 1);
1838  Vmath::Vvtvp(nPointsTot, grad1, 1, V[1], 1,
1839  outarray, 1, outarray, 1);
1840  break;
1841  case 3:
1842  grad1 = Array<OneD, NekDouble> (nPointsTot);
1843  grad2 = Array<OneD, NekDouble> (nPointsTot);
1844  m_fields[0]->PhysDeriv(u,grad0,grad1,grad2);
1845  Vmath::Vmul (nPointsTot, grad0, 1, V[0], 1, outarray, 1);
1846  Vmath::Vvtvp(nPointsTot, grad1, 1, V[1], 1,
1847  outarray, 1, outarray, 1);
1848  Vmath::Vvtvp(nPointsTot, grad2, 1, V[2], 1,
1849  outarray, 1, outarray, 1);
1850  break;
1851  default:
1852  ASSERTL0(false,"dimension unknown");
1853  }
1854  }
1855 
1856  /**
1857  * @brief Calculate weak DG advection in the form \f$ \langle\phi,
1858  * \hat{F}\cdot n\rangle - (\nabla \phi \cdot F) \f$
1859  *
1860  * @param InField Fields.
1861  * @param OutField Storage for result.
1862  * @param NumericalFluxIncludesNormal Default: true.
1863  * @param InFieldIsPhysSpace Default: false.
1864  * @param nvariables Number of fields.
1865  */
1866  void EquationSystem::WeakDGAdvection(
1867  const Array<OneD, Array<OneD, NekDouble> >& InField,
1868  Array<OneD, Array<OneD, NekDouble> >& OutField,
1869  bool NumericalFluxIncludesNormal,
1870  bool InFieldIsInPhysSpace,
1871  int nvariables)
1872  {
1873  int i;
1874  int nVelDim = m_expdim;
1875  int nPointsTot = GetNpoints();
1876  int ncoeffs = GetNcoeffs();
1877  int nTracePointsTot = GetTraceNpoints();
1878 
1879  if (!nvariables)
1880  {
1881  nvariables = m_fields.num_elements();
1882  }
1883 
1884  Array<OneD, Array<OneD, NekDouble> > fluxvector(nVelDim);
1885  Array<OneD, Array<OneD, NekDouble> > physfield (nvariables);
1886 
1887  for(i = 0; i < nVelDim; ++i)
1888  {
1889  fluxvector[i] = Array<OneD, NekDouble>(nPointsTot);
1890  }
1891 
1892  // Get the variables in physical space
1893  // already in physical space
1894  if (InFieldIsInPhysSpace == true)
1895  {
1896  for (i = 0; i < nvariables; ++i)
1897  {
1898  physfield[i] = InField[i];
1899  }
1900  }
1901  // otherwise do a backward transformation
1902  else
1903  {
1904  for(i = 0; i < nvariables; ++i)
1905  {
1906  // Could make this point to m_fields[i]->UpdatePhys();
1907  physfield[i] = Array<OneD, NekDouble>(nPointsTot);
1908  m_fields[i]->BwdTrans(InField[i],physfield[i]);
1909  }
1910  }
1911 
1912  // Get the advection part (without numerical flux)
1913  for (i = 0; i < nvariables; ++i)
1914  {
1915  // Get the ith component of the flux vector in (physical space)
1916  GetFluxVector(i, physfield, fluxvector);
1917 
1918  // Calculate the i^th value of (\grad_i \phi, F)
1919  WeakAdvectionGreensDivergenceForm(fluxvector,OutField[i]);
1920  }
1921 
1922  // Get the numerical flux and add to the modal coeffs
1923  // if the NumericalFluxs function already includes the
1924  // normal in the output
1925  if (NumericalFluxIncludesNormal == true)
1926  {
1927  Array<OneD, Array<OneD, NekDouble> > numflux (nvariables);
1928 
1929  for (i = 0; i < nvariables; ++i)
1930  {
1931  numflux[i] = Array<OneD, NekDouble>(nTracePointsTot);
1932  }
1933 
1934  // Evaluate numerical flux in physical space which may in
1935  // general couple all component of vectors
1936  NumericalFlux(physfield, numflux);
1937 
1938  // Evaulate <\phi, \hat{F}\cdot n> - OutField[i]
1939  for (i = 0; i < nvariables; ++i)
1940  {
1941  Vmath::Neg(ncoeffs,OutField[i],1);
1942  m_fields[i]->AddTraceIntegral(numflux[i],OutField[i]);
1943  m_fields[i]->SetPhysState(false);
1944  }
1945  }
1946  // if the NumericalFlux function does not include the
1947  // normal in the output
1948  else
1949  {
1950  Array<OneD, Array<OneD, NekDouble> > numfluxX (nvariables);
1951  Array<OneD, Array<OneD, NekDouble> > numfluxY (nvariables);
1952 
1953  for (i = 0; i < nvariables; ++i)
1954  {
1955  numfluxX[i] = Array<OneD, NekDouble>(nTracePointsTot);
1956  numfluxY[i] = Array<OneD, NekDouble>(nTracePointsTot);
1957  }
1958 
1959  // Evaluate numerical flux in physical space which may in
1960  // general couple all component of vectors
1961  NumericalFlux(physfield, numfluxX, numfluxY);
1962 
1963  // Evaulate <\phi, \hat{F}\cdot n> - OutField[i]
1964  for(i = 0; i < nvariables; ++i)
1965  {
1966  Vmath::Neg(ncoeffs,OutField[i],1);
1967  m_fields[i]->AddTraceIntegral(numfluxX[i], numfluxY[i],
1968  OutField[i]);
1969  m_fields[i]->SetPhysState(false);
1970  }
1971  }
1972  }
1973 
1974  /**
1975  * Calculate weak DG Diffusion in the LDG form
1976  * \f$ \langle\psi, \hat{u}\cdot n\rangle
1977  * - \langle\nabla\psi \cdot u\rangle
1978  * \langle\phi, \hat{q}\cdot n\rangle - (\nabla \phi \cdot q) \rangle \f$
1979  */
1980  void EquationSystem::WeakDGDiffusion(
1981  const Array<OneD, Array<OneD, NekDouble> >& InField,
1982  Array<OneD, Array<OneD, NekDouble> >& OutField,
1983  bool NumericalFluxIncludesNormal,
1984  bool InFieldIsInPhysSpace)
1985  {
1986  int i, j, k;
1987  int nPointsTot = GetNpoints();
1988  int ncoeffs = GetNcoeffs();
1989  int nTracePointsTot = GetTraceNpoints();
1990  int nvariables = m_fields.num_elements();
1991  int nqvar = 2;
1992 
1993  Array<OneD, NekDouble> qcoeffs (ncoeffs);
1994  Array<OneD, NekDouble> temp (ncoeffs);
1995 
1996  Array<OneD, Array<OneD, NekDouble> > fluxvector (m_spacedim);
1997  Array<OneD, Array<OneD, NekDouble> > ufield (nvariables);
1998 
1999  Array<OneD, Array<OneD, Array<OneD, NekDouble> > > flux (nqvar);
2000  Array<OneD, Array<OneD, Array<OneD, NekDouble> > > qfield (nqvar);
2001 
2002  for (j = 0; j < nqvar; ++j)
2003  {
2004  qfield[j] = Array<OneD, Array<OneD, NekDouble> >(nqvar);
2005  flux[j] = Array<OneD, Array<OneD, NekDouble> >(nqvar);
2006 
2007  for (i = 0; i< nvariables; ++i)
2008  {
2009  ufield[i] = Array<OneD, NekDouble>(nPointsTot, 0.0);
2010  qfield[j][i] = Array<OneD, NekDouble>(nPointsTot, 0.0);
2011  flux[j][i] = Array<OneD, NekDouble>(nTracePointsTot, 0.0);
2012  }
2013  }
2014 
2015  for (k = 0; k < m_spacedim; ++k)
2016  {
2017  fluxvector[k] = Array<OneD, NekDouble>(nPointsTot, 0.0);
2018  }
2019 
2020  // Get the variables in physical space already in physical space
2021  if (InFieldIsInPhysSpace == true)
2022  {
2023  for (i = 0; i < nvariables; ++i)
2024  {
2025  ufield[i] = InField[i];
2026  }
2027  }
2028  // Otherwise do a backward transformation
2029  else
2030  {
2031  for (i = 0; i < nvariables; ++i)
2032  {
2033  // Could make this point to m_fields[i]->UpdatePhys();
2034  ufield[i] = Array<OneD, NekDouble>(nPointsTot);
2035  m_fields[i]->BwdTrans(InField[i],ufield[i]);
2036  }
2037  }
2038 
2039  // ##########################################################
2040  // Compute q_{\eta} and q_{\xi} from su
2041  // Obtain Numerical Fluxes
2042  // ##########################################################
2043  NumFluxforScalar(ufield, flux);
2044 
2045  for (j = 0; j < nqvar; ++j)
2046  {
2047  for (i = 0; i < nvariables; ++i)
2048  {
2049  // Get the ith component of the flux vector in
2050  // (physical space) fluxvector = m_tanbasis * u,
2051  // where m_tanbasis = 2 by m_spacedim by nPointsTot
2052  if (m_tanbasis.num_elements())
2053  {
2054  for (k = 0; k < m_spacedim; ++k)
2055  {
2056  Vmath::Vmul(nPointsTot, m_tanbasis[j][k], 1,
2057  ufield[i], 1, fluxvector[k], 1);
2058  }
2059  }
2060  else
2061  {
2062  GetFluxVector(i, j, ufield, fluxvector);
2063  }
2064 
2065  // Calculate the i^th value of (\grad_i \phi, F)
2066  WeakAdvectionGreensDivergenceForm(fluxvector, qcoeffs);
2067 
2068  Vmath::Neg(ncoeffs,qcoeffs,1);
2069  m_fields[i]->AddTraceIntegral(flux[j][i], qcoeffs);
2070  m_fields[i]->SetPhysState(false);
2071 
2072  // Add weighted mass matrix = M ( \nabla \cdot Tanbasis )
2073 // if(m_gradtan.num_elements())
2074 // {
2075 // MultiRegions::GlobalMatrixKey key(StdRegions::eMass,
2076 // m_gradtan[j]);
2077 // m_fields[i]->MultiRegions::ExpList::GeneralMatrixOp(key,
2078 // InField[i], temp);
2079 // Vmath::Svtvp(ncoeffs, -1.0, temp, 1, qcoeffs, 1,
2080 // qcoeffs, 1);
2081 // }
2082 
2083  //Multiply by the inverse of mass matrix
2084  m_fields[i]->MultiplyByElmtInvMass(qcoeffs, qcoeffs);
2085 
2086  // Back to physical space
2087  m_fields[i]->BwdTrans(qcoeffs, qfield[j][i]);
2088  }
2089  }
2090 
2091 
2092  // ##########################################################
2093  // Compute u from q_{\eta} and q_{\xi}
2094  // ##########################################################
2095 
2096  // Obtain Numerical Fluxes
2097  NumFluxforVector(ufield, qfield, flux[0]);
2098 
2099  for (i = 0; i < nvariables; ++i)
2100  {
2101  // L = L(tan_eta) q_eta + L(tan_xi) q_xi
2102  OutField[i] = Array<OneD, NekDouble>(ncoeffs, 0.0);
2103  temp = Array<OneD, NekDouble>(ncoeffs, 0.0);
2104 
2105  if (m_tanbasis.num_elements())
2106  {
2107  for (j = 0; j < nqvar; ++j)
2108  {
2109  for (k = 0; k < m_spacedim; ++k)
2110  {
2111  Vmath::Vmul(nPointsTot, m_tanbasis[j][k], 1,
2112  qfield[j][i], 1, fluxvector[k], 1);
2113  }
2114 
2115  WeakAdvectionGreensDivergenceForm(fluxvector, temp);
2116  Vmath::Vadd(ncoeffs, temp, 1, OutField[i], 1,
2117  OutField[i], 1);
2118  }
2119  }
2120  else
2121  {
2122  for (k = 0; k < m_spacedim; ++k)
2123  {
2124  Vmath::Vcopy(nPointsTot, qfield[k][i], 1,
2125  fluxvector[k], 1);
2126  }
2127 
2128  WeakAdvectionGreensDivergenceForm(fluxvector, OutField[i]);
2129  }
2130 
2131  // Evaulate <\phi, \hat{F}\cdot n> - OutField[i]
2132  Vmath::Neg(ncoeffs,OutField[i],1);
2133  m_fields[i]->AddTraceIntegral(flux[0][i], OutField[i]);
2134  m_fields[i]->SetPhysState(false);
2135  }
2136  }
2137 
2138  /**
2139  * Write the n-th checkpoint file.
2140  * @param n The index of the checkpoint file.
2141  */
2142  void EquationSystem::Checkpoint_Output(const int n)
2143  {
2144  std::string outname = m_sessionName + "_" +
2145  boost::lexical_cast<std::string>(n);
2146  WriteFld(outname + ".chk");
2147  }
2148 
2149  /**
2150  * Write the n-th checkpoint file.
2151  * @param n The index of the checkpoint file.
2152  */
2153  void EquationSystem::Checkpoint_Output(
2154  const int n,
2155  MultiRegions::ExpListSharedPtr &field,
2156  std::vector<Array<OneD, NekDouble> > &fieldcoeffs,
2157  std::vector<std::string> &variables)
2158  {
2159  std::string outname = m_sessionName + "_" +
2160  boost::lexical_cast<std::string>(n);
2161  WriteFld(outname, field, fieldcoeffs, variables);
2162  }
2163 
2164  /**
2165  * Write the n-th base flow into a .chk file
2166  * @param n The index of the base flow file.
2167  */
2168  void EquationSystem::Checkpoint_BaseFlow(const int n)
2169  {
2170  std::string outname = m_sessionName + "_BaseFlow_" +
2171  boost::lexical_cast<std::string>(n);
2172 
2173  WriteFld(outname + ".chk");
2174  }
2175 
2176  /**
2177  * Writes the field data to a file with the given filename.
2178  * @param outname Filename to write to.
2179  */
2180  void EquationSystem::WriteFld(const std::string &outname)
2181  {
2182  std::vector<Array<OneD, NekDouble> > fieldcoeffs(
2183  m_fields.num_elements());
2184  std::vector<std::string> variables(m_fields.num_elements());
2185 
2186  for (int i = 0; i < m_fields.num_elements(); ++i)
2187  {
2188  if (m_fields[i]->GetNcoeffs() == m_fields[0]->GetNcoeffs())
2189  {
2190  fieldcoeffs[i] = m_fields[i]->UpdateCoeffs();
2191  }
2192  else
2193  {
2194  fieldcoeffs[i] = Array<OneD,NekDouble>(m_fields[0]->
2195  GetNcoeffs());
2196  m_fields[0]->ExtractCoeffsToCoeffs(m_fields[i],
2197  m_fields[i]->GetCoeffs(),
2198  fieldcoeffs[i]);
2199  }
2200  variables[i] = m_boundaryConditions->GetVariable(i);
2201  }
2202 
2203  v_ExtraFldOutput(fieldcoeffs, variables);
2204 
2205  WriteFld(outname, m_fields[0], fieldcoeffs, variables);
2206  }
2207 
2208 
2209 
2210  /**
2211  * Writes the field data to a file with the given filename.
2212  * @param outname Filename to write to.
2213  * @param field ExpList on which data is based.
2214  * @param fieldcoeffs An array of array of expansion coefficients.
2215  * @param variables An array of variable names.
2216  */
2217  void EquationSystem::WriteFld(
2218  const std::string &outname,
2219  MultiRegions::ExpListSharedPtr &field,
2220  std::vector<Array<OneD, NekDouble> > &fieldcoeffs,
2221  std::vector<std::string> &variables)
2222  {
2223  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef
2224  = field->GetFieldDefinitions();
2225  std::vector<std::vector<NekDouble> > FieldData(FieldDef.size());
2226 
2227  // Copy Data into FieldData and set variable
2228  for(int j = 0; j < fieldcoeffs.size(); ++j)
2229  {
2230  for(int i = 0; i < FieldDef.size(); ++i)
2231  {
2232  // Could do a search here to find correct variable
2233  FieldDef[i]->m_fields.push_back(variables[j]);
2234  field->AppendFieldData(FieldDef[i], FieldData[i],
2235  fieldcoeffs[j]);
2236  }
2237  }
2238 
2239  // Update time in field info if required
2240  if(m_fieldMetaDataMap.find("Time") != m_fieldMetaDataMap.end())
2241  {
2242  m_fieldMetaDataMap["Time"] = boost::lexical_cast<std::string>(m_time);
2243  }
2244 
2245  // Update step in field info if required
2246  if(m_fieldMetaDataMap.find("ChkFileNum") != m_fieldMetaDataMap.end())
2247  {
2248  m_fieldMetaDataMap["ChkFileNum"] = boost::lexical_cast<std::string>(m_nchk);
2249  }
2250 
2251  // If necessary, add mapping information to metadata
2252  // and output mapping coordinates
2253  Array<OneD, MultiRegions::ExpListSharedPtr> fields(1);
2254  fields[0] = field;
2255  GlobalMapping::MappingSharedPtr mapping =
2256  GlobalMapping::Mapping::Load(m_session, fields);
2257  LibUtilities::FieldMetaDataMap fieldMetaDataMap(m_fieldMetaDataMap);
2258  mapping->Output( fieldMetaDataMap, outname);
2259 
2260  m_fld->Write(outname, FieldDef, FieldData, fieldMetaDataMap, true);
2261  }
2262 
2263 
2264  /**
2265  * Import field from infile and load into \a m_fields. This routine will
2266  * also perform a \a BwdTrans to ensure data is in both the physical and
2267  * coefficient storage.
2268  * @param infile Filename to read.
2269  */
2270  void EquationSystem::ImportFld(
2271  const std::string &infile,
2272  Array<OneD, MultiRegions::ExpListSharedPtr> &pFields)
2273  {
2274  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef;
2275  std::vector<std::vector<NekDouble> > FieldData;
2276  LibUtilities::FieldIOSharedPtr field_fld =
2277  LibUtilities::FieldIO::CreateForFile(m_session, infile);
2278  field_fld->Import(infile,FieldDef,FieldData);
2279 
2280  // Copy FieldData into m_fields
2281  for(int j = 0; j < pFields.num_elements(); ++j)
2282  {
2283  Vmath::Zero(pFields[j]->GetNcoeffs(),
2284  pFields[j]->UpdateCoeffs(),1);
2285 
2286  for(int i = 0; i < FieldDef.size(); ++i)
2287  {
2288  ASSERTL1(FieldDef[i]->m_fields[j] ==
2289  m_session->GetVariable(j),
2290  std::string("Order of ") + infile
2291  + std::string(" data and that defined in "
2292  "m_boundaryconditions differs"));
2293 
2294  pFields[j]->ExtractDataToCoeffs(FieldDef[i], FieldData[i],
2295  FieldDef[i]->m_fields[j],
2296  pFields[j]->UpdateCoeffs());
2297  }
2298  pFields[j]->BwdTrans(pFields[j]->GetCoeffs(),
2299  pFields[j]->UpdatePhys());
2300  }
2301  }
2302 
2303 
2304 
2305  /**
2306  * Import field from infile and load into \a m_fields. This routine will
2307  * also perform a \a BwdTrans to ensure data is in both the physical and
2308  * coefficient storage.
2309  * @param infile Filename to read.
2310  * If optionan \a ndomains is specified it assumes we loop over nodmains for each nvariables.
2311  */
2312  void EquationSystem::ImportFldToMultiDomains(
2313  const std::string &infile,
2314  Array<OneD, MultiRegions::ExpListSharedPtr> &pFields,
2315  const int ndomains)
2316  {
2317  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef;
2318  std::vector<std::vector<NekDouble> > FieldData;
2319 
2320  LibUtilities::Import(infile,FieldDef,FieldData);
2321 
2322  int nvariables = GetNvariables();
2323 
2324  ASSERTL0(ndomains*nvariables == pFields.num_elements(),"Number of fields does not match the number of variables and domains");
2325 
2326  // Copy FieldData into m_fields
2327  for(int j = 0; j < ndomains; ++j)
2328  {
2329  for(int i = 0; i < nvariables; ++i)
2330  {
2331  Vmath::Zero(pFields[j*nvariables+i]->GetNcoeffs(),pFields[j*nvariables+i]->UpdateCoeffs(),1);
2332 
2333  for(int n = 0; n < FieldDef.size(); ++n)
2334  {
2335  ASSERTL1(FieldDef[n]->m_fields[i] == m_session->GetVariable(i),
2336  std::string("Order of ") + infile
2337  + std::string(" data and that defined in "
2338  "m_boundaryconditions differs"));
2339 
2340  pFields[j*nvariables+i]->ExtractDataToCoeffs(FieldDef[n], FieldData[n],
2341  FieldDef[n]->m_fields[i],
2342  pFields[j*nvariables+i]->UpdateCoeffs());
2343  }
2344  pFields[j*nvariables+i]->BwdTrans(pFields[j*nvariables+i]->GetCoeffs(),
2345  pFields[j*nvariables+i]->UpdatePhys());
2346  }
2347  }
2348  }
2349 
2350  /**
2351  * Import field from infile and load into \a pField. This routine will
2352  * also perform a \a BwdTrans to ensure data is in both the physical and
2353  * coefficient storage.
2354  */
2355  void EquationSystem::ImportFld(
2356  const std::string &infile,
2357  MultiRegions::ExpListSharedPtr &pField,
2358  std::string &pFieldName)
2359  {
2360  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef;
2361  std::vector<std::vector<NekDouble> > FieldData;
2362 
2363  LibUtilities::FieldIOSharedPtr field_fld =
2364  LibUtilities::FieldIO::CreateForFile(m_session, infile);
2365  field_fld->Import(infile,FieldDef,FieldData);
2366  int idx = -1;
2367 
2368  Vmath::Zero(pField->GetNcoeffs(),pField->UpdateCoeffs(),1);
2369 
2370  for(int i = 0; i < FieldDef.size(); ++i)
2371  {
2372  // find the index of the required field in the file.
2373  for(int j = 0; j < FieldData.size(); ++j)
2374  {
2375  if (FieldDef[i]->m_fields[j] == pFieldName)
2376  {
2377  idx = j;
2378  }
2379  }
2380  ASSERTL1(idx >= 0, "Field " + pFieldName + " not found.");
2381 
2382  pField->ExtractDataToCoeffs(FieldDef[i], FieldData[i],
2383  FieldDef[i]->m_fields[idx],
2384  pField->UpdateCoeffs());
2385  }
2386  pField->BwdTrans(pField->GetCoeffs(), pField->UpdatePhys());
2387  }
2388 
2389  /**
2390  * Import field from infile and load into the array \a coeffs.
2391  *
2392  * @param infile Filename to read.
2393  * @param fieldStr an array of string identifying fields to be imported
2394  * @param coeffs and array of array of coefficients to store imported data
2395  */
2396  void EquationSystem::ImportFld(
2397  const std::string &infile,
2398  std::vector< std::string> &fieldStr,
2399  Array<OneD, Array<OneD, NekDouble> > &coeffs)
2400  {
2401 
2402  ASSERTL0(fieldStr.size() <= coeffs.num_elements(),
2403  "length of fieldstr should be the same as pFields");
2404 
2405  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef;
2406  std::vector<std::vector<NekDouble> > FieldData;
2407 
2408  LibUtilities::FieldIOSharedPtr field_fld =
2409  LibUtilities::FieldIO::CreateForFile(m_session, infile);
2410  field_fld->Import(infile,FieldDef,FieldData);
2411 
2412  // Copy FieldData into m_fields
2413  for(int j = 0; j < fieldStr.size(); ++j)
2414  {
2415  Vmath::Zero(coeffs[j].num_elements(),coeffs[j],1);
2416  for(int i = 0; i < FieldDef.size(); ++i)
2417  {
2418  m_fields[0]->ExtractDataToCoeffs(FieldDef[i], FieldData[i],
2419  fieldStr[j], coeffs[j]);
2420  }
2421  }
2422  }
2423 
2424  /**
2425  * Write out a summary of the session data.
2426  * @param out Output stream to write data to.
2427  */
2428  void EquationSystem::SessionSummary(SummaryList& s)
2429  {
2430  AddSummaryItem(s, "EquationType", m_session->GetSolverInfo("EQTYPE"));
2431  AddSummaryItem(s, "Session Name", m_sessionName);
2432  AddSummaryItem(s, "Spatial Dim.", m_spacedim);
2433  AddSummaryItem(s, "Max SEM Exp. Order", m_fields[0]->EvalBasisNumModesMax());
2434 
2435  if (m_session->GetComm()->GetSize() > 1)
2436  {
2437  AddSummaryItem(s, "Num. Processes",
2438  m_session->GetComm()->GetSize());
2439  }
2440 
2441  if(m_HomogeneousType == eHomogeneous1D)
2442  {
2443  AddSummaryItem(s, "Quasi-3D", "Homogeneous in z-direction");
2444  AddSummaryItem(s, "Expansion Dim.", m_expdim + 1);
2445  AddSummaryItem(s, "Num. Hom. Modes (z)", m_npointsZ);
2446  AddSummaryItem(s, "Hom. length (LZ)", m_LhomZ);
2447  AddSummaryItem(s, "FFT Type", m_useFFT ? "FFTW" : "MVM");
2448  if (m_halfMode)
2449  {
2450  AddSummaryItem(s, "ModeType", "Half Mode");
2451  }
2452  else if (m_singleMode)
2453  {
2454  AddSummaryItem(s, "ModeType", "Single Mode");
2455  }
2456  else if (m_multipleModes)
2457  {
2458  AddSummaryItem(s, "ModeType", "Multiple Modes");
2459  }
2460  }
2461  else if(m_HomogeneousType == eHomogeneous2D)
2462  {
2463  AddSummaryItem(s, "Quasi-3D", "Homogeneous in yz-plane");
2464  AddSummaryItem(s, "Expansion Dim.", m_expdim + 2);
2465  AddSummaryItem(s, "Num. Hom. Modes (y)", m_npointsY);
2466  AddSummaryItem(s, "Num. Hom. Modes (z)", m_npointsZ);
2467  AddSummaryItem(s, "Hom. length (LY)", m_LhomY);
2468  AddSummaryItem(s, "Hom. length (LZ)", m_LhomZ);
2469  AddSummaryItem(s, "FFT Type", m_useFFT ? "FFTW" : "MVM");
2470  }
2471  else
2472  {
2473  AddSummaryItem(s, "Expansion Dim.", m_expdim);
2474  }
2475 
2476  if (m_session->DefinesSolverInfo("UpwindType"))
2477  {
2478  AddSummaryItem(s, "Riemann Solver",
2479  m_session->GetSolverInfo("UpwindType"));
2480  }
2481 
2482  if (m_session->DefinesSolverInfo("AdvectionType"))
2483  {
2484  std::string AdvectionType;
2485  AdvectionType = m_session->GetSolverInfo("AdvectionType");
2486  AddSummaryItem(s, "Advection Type", GetAdvectionFactory().
2487  GetClassDescription(AdvectionType));
2488  }
2489 
2490  if (m_projectionType == MultiRegions::eGalerkin)
2491  {
2492  AddSummaryItem(s, "Projection Type", "Continuous Galerkin");
2493  }
2494  else if (m_projectionType == MultiRegions::eDiscontinuous)
2495  {
2496  AddSummaryItem(s, "Projection Type", "Discontinuous Galerkin");
2497  }
2498  else if (m_projectionType == MultiRegions::eMixed_CG_Discontinuous)
2499  {
2500  AddSummaryItem(s, "Projection Type",
2501  "Mixed Continuous Galerkin and Discontinuous");
2502  }
2503 
2504  if (m_session->DefinesSolverInfo("DiffusionType"))
2505  {
2506  std::string DiffusionType;
2507  DiffusionType = m_session->GetSolverInfo("DiffusionType");
2508  AddSummaryItem(s, "Diffusion Type", GetDiffusionFactory().
2509  GetClassDescription(DiffusionType));
2510  }
2511  }
2512 
2513  /**
2514  * Performs a case-insensitive string comparison (from web).
2515  * @param s1 First string to compare.
2516  * @param s2 Second string to compare.
2517  * @returns 0 if the strings match.
2518  */
2519  int EquationSystem::NoCaseStringCompare(
2520  const string & s1,
2521  const string& s2)
2522  {
2523  //if (s1.size() < s2.size()) return -1;
2524  //if (s1.size() > s2.size()) return 1;
2525 
2526  string::const_iterator it1=s1.begin();
2527  string::const_iterator it2=s2.begin();
2528 
2529  // Stop when either string's end has been reached
2530  while ( (it1!=s1.end()) && (it2!=s2.end()) )
2531  {
2532  if(::toupper(*it1) != ::toupper(*it2)) //letters differ?
2533  {
2534  // Return -1 to indicate smaller than, 1 otherwise
2535  return (::toupper(*it1) < ::toupper(*it2)) ? -1 : 1;
2536  }
2537 
2538  // Proceed to the next character in each string
2539  ++it1;
2540  ++it2;
2541  }
2542 
2543  size_t size1=s1.size();
2544  size_t size2=s2.size();// cache lengths
2545 
2546  // Return -1, 0 or 1 according to strings' lengths
2547  if (size1==size2)
2548  {
2549  return 0;
2550  }
2551 
2552  return (size1 < size2) ? -1 : 1;
2553  }
2554 
2555  Array<OneD, bool> EquationSystem::v_GetSystemSingularChecks()
2556  {
2557  return Array<OneD, bool>(m_session->GetVariables().size(), false);
2558  }
2559 
2560  void EquationSystem::v_GetFluxVector(
2561  const int i, Array<OneD,
2562  Array<OneD, NekDouble> > &physfield,
2563  Array<OneD, Array<OneD, NekDouble> > &flux)
2564  {
2565  ASSERTL0(false, "v_GetFluxVector: This function is not valid "
2566  "for the Base class");
2567  }
2568 
2569  void EquationSystem::v_GetFluxVector(
2570  const int i, const int j,
2571  Array<OneD, Array<OneD, NekDouble> > &physfield,
2572  Array<OneD, Array<OneD, NekDouble> > &flux)
2573  {
2574  ASSERTL0(false, "v_GetqFluxVector: This function is not valid "
2575  "for the Base class");
2576  }
2577 
2578  void EquationSystem::v_GetFluxVector(
2579  const int i, Array<OneD,
2580  Array<OneD, NekDouble> > &physfield,
2581  Array<OneD, Array<OneD, NekDouble> > &fluxX,
2582  Array<OneD, Array<OneD, NekDouble> > &fluxY)
2583  {
2584  ASSERTL0(false, "v_GetFluxVector: This function is not valid "
2585  "for the Base class");
2586  }
2587 
2588  void EquationSystem::v_NumericalFlux(
2589  Array<OneD, Array<OneD, NekDouble> > &physfield,
2590  Array<OneD, Array<OneD, NekDouble> > &numflux)
2591  {
2592  ASSERTL0(false, "v_NumericalFlux: This function is not valid "
2593  "for the Base class");
2594  }
2595 
2596  void EquationSystem::v_NumericalFlux(
2597  Array<OneD, Array<OneD, NekDouble> > &physfield,
2598  Array<OneD, Array<OneD, NekDouble> > &numfluxX,
2599  Array<OneD, Array<OneD, NekDouble> > &numfluxY )
2600  {
2601  ASSERTL0(false, "v_NumericalFlux: This function is not valid "
2602  "for the Base class");
2603  }
2604 
2605  void EquationSystem::v_NumFluxforScalar(
2606  const Array<OneD, Array<OneD, NekDouble> > &ufield,
2607  Array<OneD, Array<OneD, Array<OneD, NekDouble> > > &uflux)
2608  {
2609  ASSERTL0(false, "v_NumFluxforScalar: This function is not valid "
2610  "for the Base class");
2611  }
2612 
2613  void EquationSystem::v_NumFluxforVector(
2614  const Array<OneD, Array<OneD, NekDouble> > &ufield,
2615  Array<OneD, Array<OneD, Array<OneD, NekDouble> > > &qfield,
2616  Array<OneD, Array<OneD, NekDouble > > &qflux)
2617  {
2618  ASSERTL0(false, "v_NumFluxforVector: This function is not valid "
2619  "for the Base class");
2620  }
2621 
2622  MultiRegions::ExpListSharedPtr EquationSystem::v_GetPressure()
2623  {
2624  ASSERTL0(false, "This function is not valid for the Base class");
2625  MultiRegions::ExpListSharedPtr null;
2626  return null;
2627  }
2628 
2629  void EquationSystem::v_ExtraFldOutput(
2630  std::vector<Array<OneD, NekDouble> > &fieldcoeffs,
2631  std::vector<std::string> &variables)
2632  {
2633  }
2634 
2635  }
2636 }
Nektar::SolverUtils::EquationSystem::EvaluateFunctionPts
SOLVER_UTILS_EXPORT void EvaluateFunctionPts(std::string pFieldName, Array< OneD, NekDouble > &pArray, const std::string &pFunctionName, const NekDouble &pTime=0.0, const int domain=0)
Definition: EquationSystem.cpp:958
Nektar::SolverUtils::EquationSystem::v_GenerateSummary
virtual SOLVER_UTILS_EXPORT void v_GenerateSummary(SummaryList &l)
Virtual function for generating summary information.
Definition: EquationSystem.cpp:1678
Nektar::SolverUtils::EquationSystem::m_singleMode
bool m_singleMode
Flag to determine if single homogeneous mode is used.
Definition: EquationSystem.h:502
Nektar::SolverUtils::EquationSystem::NumFluxforScalar
SOLVER_UTILS_EXPORT void NumFluxforScalar(const Array< OneD, Array< OneD, NekDouble > > &ufield, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &uflux)
Definition: EquationSystem.h:982
Nektar::SolverUtils::EquationSystem::AdvectionNonConservativeForm
f SOLVER_UTILS_EXPORT void AdvectionNonConservativeForm(const Array< OneD, Array< OneD, NekDouble > > &V, const Array< OneD, const NekDouble > &u, Array< OneD, NekDouble > &outarray, Array< OneD, NekDouble > &wk=NullNekDouble1DArray)
Compute the non-conservative advection.
Definition: EquationSystem.cpp:1802
ExpList3DHomogeneous1D.h
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
Nektar::SolverUtils::EquationSystem::InitialiseBaseFlow
SOLVER_UTILS_EXPORT void InitialiseBaseFlow(Array< OneD, Array< OneD, NekDouble > > &base)
Perform initialisation of the base flow.
Definition: EquationSystem.cpp:1427
Nektar::FieldUtils::Interpolator
A class that contains algorithms for interpolation between pts fields, expansions and different meshe...
Definition: Interpolator.h:78
ExpList2D.h
Nektar::SolverUtils::EquationSystemFactory
LibUtilities::NekFactory< std::string, EquationSystem, const LibUtilities::SessionReaderSharedPtr & > EquationSystemFactory
Datatype of the NekFactory used to instantiate classes derived from the EquationSystem class...
Definition: EquationSystem.h:66
Nektar::SolverUtils::EquationSystem::DescribeFunction
SOLVER_UTILS_EXPORT std::string DescribeFunction(std::string pFieldName, const std::string &pFunctionName, const int domain)
Provide a description of a function for a given field name.
Definition: EquationSystem.cpp:1119
Nektar::SpatialDomains::MeshGraph::Read
static boost::shared_ptr< MeshGraph > Read(const LibUtilities::SessionReaderSharedPtr &pSession, DomainRangeShPtr &rng=NullDomainRangeShPtr)
Definition: MeshGraph.cpp:124
Nektar::MemoryManager::AllocateSharedPtr
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
Definition: NekMemoryManager.hpp:235
ContField3D.h
Nektar::SolverUtils::EquationSystem::LoadPts
SOLVER_UTILS_EXPORT void LoadPts(std::string funcFilename, std::string filename, Nektar::LibUtilities::PtsFieldSharedPtr &outPts)
Definition: EquationSystem.cpp:1058
Nektar::SolverUtils::EquationSystem::m_time
NekDouble m_time
Current time of simulation.
Definition: EquationSystem.h:480
Nektar::SolverUtils::EquationSystem::WeakAdvectionNonConservativeForm
SOLVER_UTILS_EXPORT void WeakAdvectionNonConservativeForm(const Array< OneD, Array< OneD, NekDouble > > &V, const Array< OneD, const NekDouble > &u, Array< OneD, NekDouble > &outarray, bool UseContCoeffs=false)
Compute the inner product .
Definition: EquationSystem.cpp:1769
Nektar::SolverUtils::EquationSystem::m_loadedFldFields
std::map< std::string, std::pair< std::string, loadedFldField > > m_loadedFldFields
Definition: EquationSystem.h:466
Nektar::LibUtilities::Import
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:279
Nektar::SolverUtils::EquationSystem::m_timestep
NekDouble m_timestep
Time step size.
Definition: EquationSystem.h:486
Nektar::MemoryManager
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
Definition: NekMemoryManager.hpp:102
Nektar::SolverUtils::EquationSystem::NoCaseStringCompare
SOLVER_UTILS_EXPORT int NoCaseStringCompare(const std::string &s1, const std::string &s2)
Perform a case-insensitive string comparison.
Definition: EquationSystem.cpp:2519
Nektar::SolverUtils::SummaryList
std::vector< std::pair< std::string, std::string > > SummaryList
Definition: Misc.h:47
Nektar::SolverUtils::EquationSystem::GetNvariables
SOLVER_UTILS_EXPORT int GetNvariables()
Definition: EquationSystem.h:865
Nektar::SolverUtils::EquationSystem::m_halfMode
bool m_halfMode
Flag to determine if half homogeneous mode is used.
Definition: EquationSystem.h:504
Nektar::SolverUtils::EquationSystem::m_LhomZ
NekDouble m_LhomZ
physical length in Z direction (if homogeneous)
Definition: EquationSystem.h:550
Nektar::SolverUtils::EquationSystem::v_L2Error
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.
Definition: EquationSystem.cpp:1170
ContField2D.h
Nektar::SolverUtils::EquationSystem::m_checkIfSystemSingular
Array< OneD, bool > m_checkIfSystemSingular
Flag to indicate if the fields should be checked for singularity.
Definition: EquationSystem.h:528
Nektar::SolverUtils::EquationSystem::m_expdim
int m_expdim
Expansion dimension.
Definition: EquationSystem.h:500
Nektar::SolverUtils::EquationSystem::Checkpoint_Output
SOLVER_UTILS_EXPORT void Checkpoint_Output(const int n)
Write checkpoint file of m_fields.
Definition: EquationSystem.cpp:2142
Nektar::SolverUtils::GetDiffusionFactory
DiffusionFactory & GetDiffusionFactory()
Definition: Diffusion.cpp:42
Vmath::Vvtvp
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
Definition: Vmath.cpp:442
Nektar::SolverUtils::EquationSystem::m_projectionType
enum MultiRegions::ProjectionType m_projectionType
Type of projection; e.g continuous or discontinuous.
Definition: EquationSystem.h:520
Nektar::LibUtilities::eModified_A
Principle Modified Functions .
Definition: BasisType.h:49
Nektar::SolverUtils::EquationSystem::m_loadedPtsFields
std::map< std::string, std::pair< std::string, LibUtilities::PtsFieldSharedPtr > > m_loadedPtsFields
pts fields we already read from disk: {funcFilename: (filename, ptsfield)}
Definition: EquationSystem.h:464
Nektar::SolverUtils::EquationSystem::ErrorExtraPoints
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].
Definition: EquationSystem.cpp:1279
std
STL namespace.
Nektar::SolverUtils::EquationSystem::WeakAdvectionGreensDivergenceForm
SOLVER_UTILS_EXPORT void WeakAdvectionGreensDivergenceForm(const Array< OneD, Array< OneD, NekDouble > > &F, Array< OneD, NekDouble > &outarray)
Compute the inner product .
Definition: EquationSystem.cpp:1728
Nektar::SolverUtils::EquationSystem::m_checktime
NekDouble m_checktime
Time between checkpoints.
Definition: EquationSystem.h:490
Nektar::SolverUtils::EquationSystem::m_useFFT
bool m_useFFT
Flag to determine if FFT is used for homogeneous transform.
Definition: EquationSystem.h:508
Nektar::SolverUtils::EquationSystem::m_LhomY
NekDouble m_LhomY
physical length in Y direction (if homogeneous)
Definition: EquationSystem.h:549
Nektar::MultiRegions::ContField2DSharedPtr
boost::shared_ptr< ContField2D > ContField2DSharedPtr
Definition: ContField2D.h:287
Nektar::SolverUtils::EquationSystem::m_specHP_dealiasing
bool m_specHP_dealiasing
Flag to determine if dealisising is usde for the Spectral/hp element discretisation.
Definition: EquationSystem.h:518
Nektar::SolverUtils::EquationSystem::GetNumExpModes
SOLVER_UTILS_EXPORT int GetNumExpModes()
Definition: EquationSystem.h:854
Nektar::LibUtilities::eFourier
Fourier Expansion .
Definition: BasisType.h:52
Nektar::SolverUtils::EquationSystem::v_ExtraFldOutput
virtual SOLVER_UTILS_EXPORT void v_ExtraFldOutput(std::vector< Array< OneD, NekDouble > > &fieldcoeffs, std::vector< std::string > &variables)
Definition: EquationSystem.cpp:2629
Nektar::SolverUtils::EquationSystem::PrintProgressbar
SOLVER_UTILS_EXPORT void PrintProgressbar(const int position, const int goal) const
Definition: EquationSystem.h:690
ExpList3DHomogeneous2D.h
Nektar::SolverUtils::EquationSystem::v_LinfError
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...
Definition: EquationSystem.cpp:1230
Nektar::SolverUtils::EquationSystem::m_npointsZ
int m_npointsZ
number of points in Z direction (if homogeneous)
Definition: EquationSystem.h:554
Nektar::SolverUtils::EquationSystem::~EquationSystem
virtual SOLVER_UTILS_EXPORT ~EquationSystem()
Destructor.
Definition: EquationSystem.cpp:674
Nektar::SolverUtils::EquationSystem::m_sessionName
std::string m_sessionName
Name of the session.
Definition: EquationSystem.h:478
Nektar::SolverUtils::EquationSystem::m_nchk
int m_nchk
Number of checkpoints written so far.
Definition: EquationSystem.h:492
Nektar::SolverUtils::EquationSystem::WeakAdvectionDivergenceForm
SOLVER_UTILS_EXPORT void WeakAdvectionDivergenceForm(const Array< OneD, Array< OneD, NekDouble > > &F, Array< OneD, NekDouble > &outarray)
Compute the inner product .
Definition: EquationSystem.cpp:1741
Nektar::LibUtilities::SessionReaderSharedPtr
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:51
Nektar::LibUtilities::FieldMetaDataMap
std::map< std::string, std::string > FieldMetaDataMap
Definition: FieldIO.h:54
Nektar::SolverUtils::EquationSystem::m_checksteps
int m_checksteps
Number of steps between checkpoints.
Definition: EquationSystem.h:496
Nektar::SolverUtils::EquationSystem::EvaluateFunctionExp
SOLVER_UTILS_EXPORT void EvaluateFunctionExp(std::string pFieldName, Array< OneD, NekDouble > &pArray, const std::string &pFunctionName, const NekDouble &pTime=0.0, const int domain=0)
Definition: EquationSystem.cpp:791
Nektar::LibUtilities::PtsFieldSharedPtr
boost::shared_ptr< PtsField > PtsFieldSharedPtr
Definition: PtsField.h:178
Nektar::SolverUtils::EquationSystem::m_comm
LibUtilities::CommSharedPtr m_comm
Communicator.
Definition: EquationSystem.h:456
Nektar::LibUtilities::eGaussRadauMAlpha1Beta0
Gauss Radau pinned at x=-1, .
Definition: PointsType.h:59
ContField3DHomogeneous2D.h
Nektar::SolverUtils::EquationSystem::m_fintime
NekDouble m_fintime
Finish time of the simulation.
Definition: EquationSystem.h:484
Nektar::SolverUtils::EquationSystem::v_InitObject
virtual SOLVER_UTILS_EXPORT void v_InitObject()
Initialisation object for EquationSystem.
Definition: EquationSystem.cpp:120
Nektar::SolverUtils::EquationSystem::m_tanbasis
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_tanbasis
2 x m_spacedim x nq
Definition: EquationSystem.h:526
Nektar::MultiRegions::eGlobal
Global coefficients.
Definition: MultiRegions.hpp:52
Nektar::SolverUtils::EquationSystem::m_steps
int m_steps
Number of steps to take.
Definition: EquationSystem.h:494
ContField3DHomogeneous1D.h
Nektar::SolverUtils::EquationSystem::m_traceNormals
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Array holding trace normals for DG simulations in the forwards direction.
Definition: EquationSystem.h:522
Nektar::SolverUtils::EquationSystem::m_HomoDirec
int m_HomoDirec
number of homogenous directions
Definition: EquationSystem.h:556
Nektar::SolverUtils::EquationSystem::FwdTransFields
SOLVER_UTILS_EXPORT void FwdTransFields()
Definition: EquationSystem.cpp:1708
Nektar::SolverUtils::EquationSystem::WeakDGDiffusion
SOLVER_UTILS_EXPORT void WeakDGDiffusion(const Array< OneD, Array< OneD, NekDouble > > &InField, Array< OneD, Array< OneD, NekDouble > > &OutField, bool NumericalFluxIncludesNormal=true, bool InFieldIsInPhysSpace=false)
Calculate weak DG Diffusion in the LDG form.
Definition: EquationSystem.cpp:1980
Nektar::SolverUtils::EquationSystem::EvaluateFunctionFld
SOLVER_UTILS_EXPORT void EvaluateFunctionFld(std::string pFieldName, Array< OneD, NekDouble > &pArray, const std::string &pFunctionName, const NekDouble &pTime=0.0, const int domain=0)
Definition: EquationSystem.cpp:824
Nektar::LibUtilities::eFourierEvenlySpaced
1D Evenly-spaced points using Fourier Fit
Definition: PointsType.h:66
Nektar::SolverUtils::EquationSystem::m_multipleModes
bool m_multipleModes
Flag to determine if use multiple homogenenous modes are used.
Definition: EquationSystem.h:506
Nektar::LibUtilities::eFourierHalfModeRe
Fourier Modified expansions with just the real part of the first mode .
Definition: BasisType.h:59
Nektar::SolverUtils::EquationSystem::v_GetSystemSingularChecks
virtual SOLVER_UTILS_EXPORT Array< OneD, bool > v_GetSystemSingularChecks()
Definition: EquationSystem.cpp:2555
Nektar::SolverUtils::EquationSystem::NumericalFlux
SOLVER_UTILS_EXPORT void NumericalFlux(Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &numflux)
Definition: EquationSystem.h:969
Nektar::SolverUtils::EquationSystem::SetUpBaseFields
SOLVER_UTILS_EXPORT void SetUpBaseFields(SpatialDomains::MeshGraphSharedPtr &mesh)
Definition: EquationSystem.cpp:1440
Nektar::MultiRegions::ContField2D
This class is the abstraction of a global continuous two- dimensional spectral/hp element expansion w...
Definition: ContField2D.h:56
Nektar::SolverUtils::EquationSystem::m_npointsY
int m_npointsY
number of points in Y direction (if homogeneous)
Definition: EquationSystem.h:553
Nektar::MultiRegions::ContField1D
Abstraction of a global continuous one-dimensional spectral/hp element expansion which approximates t...
Definition: ContField1D.h:56
Nektar::DNekScalBlkMat
NekMatrix< NekMatrix< NekMatrix< NekDouble, StandardMatrixTag >, ScaledMatrixTag >, BlockMatrixTag > DNekScalBlkMat
Definition: NekTypeDefs.hpp:66
Nektar::SolverUtils::EquationSystem::v_TransPhysToCoeff
virtual SOLVER_UTILS_EXPORT void v_TransPhysToCoeff()
Virtual function for transformation to coefficient space.
Definition: EquationSystem.cpp:1671
Nektar::SolverUtils::EquationSystem::v_GetFluxVector
virtual SOLVER_UTILS_EXPORT void v_GetFluxVector(const int i, Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &flux)
Definition: EquationSystem.cpp:2560
Nektar::LibUtilities::eModified_B
Principle Modified Functions .
Definition: BasisType.h:50
Nektar::SolverUtils::AddSummaryItem
void AddSummaryItem(SummaryList &l, const std::string &name, const std::string &value)
Adds a summary item to the summary info list.
Definition: Misc.cpp:50
Nektar::SolverUtils::EquationSystem::v_Output
virtual SOLVER_UTILS_EXPORT void v_Output(void)
Definition: EquationSystem.cpp:1688
Nektar::SolverUtils::EquationSystem::m_boundaryConditions
SpatialDomains::BoundaryConditionsSharedPtr m_boundaryConditions
Pointer to boundary conditions object.
Definition: EquationSystem.h:474
Nektar::MultiRegions::ExpListSharedPtr
boost::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
Definition: LocTraceToTraceMap.h:56
Nektar::SolverUtils::EquationSystem::m_spacedim
int m_spacedim
Spatial dimension (>= expansion dim).
Definition: EquationSystem.h:498
Nektar::MultiRegions::DisContField1D
This class is the abstraction of a global discontinuous two- dimensional spectral/hp element expansio...
Definition: DisContField1D.h:55
Nektar::SolverUtils::EquationSystem::m_fieldMetaDataMap
LibUtilities::FieldMetaDataMap m_fieldMetaDataMap
Map to identify relevant solver info to dump in output fields.
Definition: EquationSystem.h:530
Nektar::SolverUtils::GetAdvectionFactory
AdvectionFactory & GetAdvectionFactory()
Gets the factory for initialising advection objects.
Definition: Advection.cpp:46
Nektar::LibUtilities::PointsKey
Defines a specification for a set of points.
Definition: Points.h:58
Vmath::Neg
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:396
Nektar::LibUtilities::FieldIOSharedPtr
boost::shared_ptr< FieldIO > FieldIOSharedPtr
Definition: FieldIO.h:309
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:44
Nektar::SolverUtils::EquationSystem::v_SetInitialConditions
virtual SOLVER_UTILS_EXPORT void v_SetInitialConditions(NekDouble initialtime=0.0, bool dumpInitialConditions=true, const int domain=0)
Definition: EquationSystem.cpp:1351
Nektar::SolverUtils::EquationSystem::GetSession
SOLVER_UTILS_EXPORT LibUtilities::SessionReaderSharedPtr GetSession()
Get Session name.
Definition: EquationSystem.h:131
Nektar::SolverUtils::EquationSystem::GetFluxVector
SOLVER_UTILS_EXPORT void GetFluxVector(const int i, Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &flux)
Definition: EquationSystem.h:947
Nektar::MultiRegions::ExpList2DSharedPtr
boost::shared_ptr< ExpList2D > ExpList2DSharedPtr
Shared pointer to an ExpList2D object.
Definition: ExpList2D.h:49
Nektar::SolverUtils::EquationSystem::EvaluateFunction
SOLVER_UTILS_EXPORT void EvaluateFunction(Array< OneD, Array< OneD, NekDouble > > &pArray, std::string pFunctionName, const NekDouble pTime=0.0, const int domain=0)
Evaluates a function as specified in the session file.
Definition: EquationSystem.cpp:688
Nektar::GlobalMapping::MappingSharedPtr
GLOBAL_MAPPING_EXPORT typedef boost::shared_ptr< Mapping > MappingSharedPtr
A shared pointer to a Mapping object.
Definition: Mapping.h:51
Nektar::LibUtilities::PtsIO::Import
void Import(const string &inFile, PtsFieldSharedPtr &ptsField, FieldMetaDataMap &fieldmetadatamap=NullFieldMetaDataMap)
Import a pts field from file.
Definition: PtsIO.cpp:69
Nektar::LibUtilities::EquationSharedPtr
boost::shared_ptr< Equation > EquationSharedPtr
Definition: SessionReader.h:100
Nektar::SolverUtils::EquationSystem::v_DoInitialise
virtual SOLVER_UTILS_EXPORT void v_DoInitialise()
Virtual function for initialisation implementation.
Definition: EquationSystem.cpp:1421
Nektar::SolverUtils::EquationSystem::m_homogen_dealiasing
bool m_homogen_dealiasing
Flag to determine if dealiasing is used for homogeneous simulations.
Definition: EquationSystem.h:513
Nektar::SolverUtils::GetEquationSystemFactory
EquationSystemFactory & GetEquationSystemFactory()
Definition: EquationSystem.cpp:83
MatrixKey.h
Nektar::LibUtilities::eFourierHalfModeIm
Fourier Modified expansions with just the imaginary part of the first mode .
Definition: BasisType.h:60
Nektar::SolverUtils::EquationSystem::ImportFldToMultiDomains
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.
Definition: EquationSystem.cpp:2312
Nektar::SolverUtils::EquationSystem::SetBoundaryConditions
SOLVER_UTILS_EXPORT void SetBoundaryConditions(NekDouble time)
Evaluates the boundary conditions at the given time.
Definition: EquationSystem.cpp:1152
ExpList3D.h
Nektar::MultiRegions::DirCartesianMap
MultiRegions::Direction const DirCartesianMap[]
Definition: ExpList.h:86
Nektar::LibUtilities::FieldIO::CreateDefault
static boost::shared_ptr< FieldIO > CreateDefault(const LibUtilities::SessionReaderSharedPtr session)
Returns an object for the default FieldIO method.
Definition: FieldIO.cpp:181
Nektar::SolverUtils::EquationSystem::WriteFld
SOLVER_UTILS_EXPORT void WriteFld(const std::string &outname)
Write field data to the given filename.
Definition: EquationSystem.cpp:2180
Nektar::LibUtilities::FieldIO::CreateForFile
static boost::shared_ptr< FieldIO > CreateForFile(const LibUtilities::SessionReaderSharedPtr session, const std::string &filename)
Construct a FieldIO object for a given input filename.
Definition: FieldIO.cpp:212
Nektar::SolverUtils::EquationSystem::GetNpoints
SOLVER_UTILS_EXPORT int GetNpoints()
Definition: EquationSystem.h:910
Nektar::SolverUtils::EquationSystem::m_fields
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
Definition: EquationSystem.h:468
Nektar::SolverUtils::EquationSystem::m_session
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
Definition: EquationSystem.h:458
Nektar::LibUtilities::eFourierSingleMode
Fourier ModifiedExpansion with just the first mode .
Definition: BasisType.h:58
EquationSystem.h
Nektar::SolverUtils::EquationSystem::m_base
Array< OneD, MultiRegions::ExpListSharedPtr > m_base
Base fields.
Definition: EquationSystem.h:470
Equation.h
Mapping.h
Nektar::SolverUtils::EquationSystem::GetTraceNpoints
SOLVER_UTILS_EXPORT int GetTraceNpoints()
Definition: EquationSystem.h:880
Nektar::SolverUtils::EquationSystem::GetNcoeffs
SOLVER_UTILS_EXPORT int GetNcoeffs()
Definition: EquationSystem.h:844
Nektar::LibUtilities::eFourierSingleModeSpaced
1D Non Evenly-spaced points for Single Mode analysis
Definition: PointsType.h:67
Nektar::SolverUtils::EquationSystem::m_fld
LibUtilities::FieldIOSharedPtr m_fld
Field input/output.
Definition: EquationSystem.h:460
ContField1D.h
Nektar::SolverUtils::EquationSystem::v_NumFluxforVector
virtual SOLVER_UTILS_EXPORT void v_NumFluxforVector(const Array< OneD, Array< OneD, NekDouble > > &ufield, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &qfield, Array< OneD, Array< OneD, NekDouble > > &qflux)
Definition: EquationSystem.cpp:2613
Nektar::SolverUtils::EquationSystem::v_NegatedOp
virtual SOLVER_UTILS_EXPORT bool v_NegatedOp()
Virtual function to identify if operator is negated in DoSolve.
Definition: EquationSystem.cpp:1655
Nektar::SolverUtils::EquationSystem::v_GetPressure
virtual SOLVER_UTILS_EXPORT MultiRegions::ExpListSharedPtr v_GetPressure(void)
Definition: EquationSystem.cpp:2622
Nektar::SolverUtils::EquationSystem::NumFluxforVector
SOLVER_UTILS_EXPORT void NumFluxforVector(const Array< OneD, Array< OneD, NekDouble > > &ufield, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &qfield, Array< OneD, Array< OneD, NekDouble > > &qflux)
Definition: EquationSystem.h:989
Nektar::SolverUtils::EquationSystem::ZeroPhysFields
SOLVER_UTILS_EXPORT void ZeroPhysFields()
Definition: EquationSystem.cpp:1696
Nektar::SolverUtils::EquationSystem::SessionSummary
SOLVER_UTILS_EXPORT void SessionSummary(SummaryList &vSummary)
Write out a session summary.
Definition: EquationSystem.cpp:2428
Nektar::SolverUtils::EquationSystem::m_graph
SpatialDomains::MeshGraphSharedPtr m_graph
Pointer to graph defining mesh.
Definition: EquationSystem.h:476
Nektar::SolverUtils::EquationSystem::v_NumFluxforScalar
virtual SOLVER_UTILS_EXPORT void v_NumFluxforScalar(const Array< OneD, Array< OneD, NekDouble > > &ufield, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &uflux)
Definition: EquationSystem.cpp:2605
AdvectionSystem.h
Nektar::SolverUtils::EquationSystem::WeakDGAdvection
SOLVER_UTILS_EXPORT void WeakDGAdvection(const Array< OneD, Array< OneD, NekDouble > > &InField, Array< OneD, Array< OneD, NekDouble > > &OutField, bool NumericalFluxIncludesNormal=true, bool InFieldIsInPhysSpace=false, int nvariables=0)
Calculate the weak discontinuous Galerkin advection.
Definition: EquationSystem.cpp:1866
Nektar::LocalRegions::MatrixKey::opLess
Used to lookup the create function in NekManager.
Definition: MatrixKey.h:68
Nektar::SolverUtils::EquationSystem::Checkpoint_BaseFlow
SOLVER_UTILS_EXPORT void Checkpoint_BaseFlow(const int n)
Write base flow file of m_fields.
Definition: EquationSystem.cpp:2168
Nektar::SolverUtils::EquationSystem::ImportFldBase
SOLVER_UTILS_EXPORT void ImportFldBase(std::string pInfile, SpatialDomains::MeshGraphSharedPtr pGraph)
Definition: EquationSystem.cpp:1607
Nektar::SolverUtils::EquationSystem::v_DoSolve
virtual SOLVER_UTILS_EXPORT void v_DoSolve()
Virtual function for solve implementation.
Definition: EquationSystem.cpp:1645
Nektar::SolverUtils::EquationSystem::m_interpolators
std::map< std::string, FieldUtils::Interpolator > m_interpolators
Map of interpolator objects.
Definition: EquationSystem.h:462
Nektar::MultiRegions::ContField3DSharedPtr
boost::shared_ptr< ContField3D > ContField3DSharedPtr
Definition: ContField3D.h:206
Vmath::Zero
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:373
Nektar::SolverUtils::EquationSystem::m_NumQuadPointsError
int m_NumQuadPointsError
Number of Quadrature points used to work out the error.
Definition: EquationSystem.h:533
ASSERTL1
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:228
Nektar::SpatialDomains::MeshGraphSharedPtr
boost::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:442
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1061
Nektar::LibUtilities::NullFieldMetaDataMap
static FieldMetaDataMap NullFieldMetaDataMap
Definition: FieldIO.h:55
Diffusion.h
Nektar::DNekScalMat
NekMatrix< NekMatrix< NekDouble, StandardMatrixTag >, ScaledMatrixTag > DNekScalMat
Definition: NekMatrixFwd.hpp:72
Nektar::GlobalMapping::Mapping::Load
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:265
Nektar::SolverUtils::EquationSystem::v_EvaluateExactSolution
virtual SOLVER_UTILS_EXPORT void v_EvaluateExactSolution(unsigned int field, Array< OneD, NekDouble > &outfield, const NekDouble time)
Definition: EquationSystem.cpp:1402
Nektar::SolverUtils::EquationSystem::v_TransCoeffToPhys
virtual SOLVER_UTILS_EXPORT void v_TransCoeffToPhys()
Virtual function for transformation to physical space.
Definition: EquationSystem.cpp:1663
Nektar::LibUtilities::BasisKey
Describes the specification for a Basis.
Definition: Basis.h:50
Nektar::SolverUtils::EquationSystem::v_NumericalFlux
virtual SOLVER_UTILS_EXPORT void v_NumericalFlux(Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &numflux)
Definition: EquationSystem.cpp:2588
Nektar::SolverUtils::EquationSystem::m_HomogeneousType
enum HomogeneousType m_HomogeneousType
Definition: EquationSystem.h:546
Nektar::LibUtilities::eGaussLobattoLegendre
1D Gauss-Lobatto-Legendre quadrature points
Definition: PointsType.h:52
Vmath::Vadd
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.cpp:299
Vmath::Vmul
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:183
Nektar::LibUtilities::NekFactory
Provides a generic Factory class.
Definition: NekFactory.hpp:116
Nektar::SolverUtils::EquationSystem::ImportFld
SOLVER_UTILS_EXPORT void ImportFld(const std::string &infile, Array< OneD, MultiRegions::ExpListSharedPtr > &pFields)
Input field data from the given file.
Definition: EquationSystem.cpp:2270
Generated on Wed Mar 8 2017 12:28:54 for Nektar++ by   doxygen 1.8.8