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