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  *
1463  */
1464  void EquationSystem::v_TransCoeffToPhys()
1465  {
1466 
1467  }
1468 
1469  /**
1470  *
1471  */
1472  void EquationSystem::v_TransPhysToCoeff()
1473  {
1474 
1475  }
1476 
1477 
1478  /// Virtual function for generating summary information.
1479  void EquationSystem::v_GenerateSummary(SummaryList& l)
1480  {
1481  SessionSummary(l);
1482  }
1483 
1484  /**
1485  * Write the field data to file. The file is named according to the session
1486  * name with the extension .fld appended.
1487  */
1488  void EquationSystem::v_Output(void)
1489  {
1490  WriteFld(m_sessionName + ".fld");
1491  }
1492 
1493  /**
1494  * Zero the physical fields.
1495  */
1496  void EquationSystem::ZeroPhysFields(void)
1497  {
1498  for (int i = 0; i < m_fields.num_elements(); i++)
1499  {
1500  Vmath::Zero(m_fields[i]->GetNpoints(),
1501  m_fields[i]->UpdatePhys(),1);
1502  }
1503  }
1504 
1505  /**
1506  * FwdTrans the m_fields members
1507  */
1508  void EquationSystem::FwdTransFields(void)
1509  {
1510  for (int i = 0; i < m_fields.num_elements(); i++)
1511  {
1512  m_fields[i]->FwdTrans(m_fields[i]->GetPhys(),
1513  m_fields[i]->UpdateCoeffs());
1514  m_fields[i]->SetPhysState(false);
1515  }
1516  }
1517 
1518  /**
1519  * Computes the weak Green form of advection terms (without boundary
1520  * integral), i.e. \f$ (\nabla \phi \cdot F) \f$ where for example
1521  * \f$ F=uV \f$.
1522  * @param F Fields.
1523  * @param outarray Storage for result.
1524  *
1525  * \note Assuming all fields are of the same expansion and order so that
1526  * we can use the parameters of m_fields[0].
1527  */
1528  void EquationSystem::WeakAdvectionGreensDivergenceForm(
1529  const Array<OneD, Array<OneD, NekDouble> > &F,
1530  Array<OneD, NekDouble> &outarray)
1531  {
1532  // Use dimension of velocity vector to dictate dimension of operation
1533  int ndim = F.num_elements();
1534  int nCoeffs = m_fields[0]->GetNcoeffs();
1535 
1536  Array<OneD, NekDouble> iprod(nCoeffs);
1537  Vmath::Zero(nCoeffs, outarray, 1);
1538 
1539  for (int i = 0; i < ndim; ++i)
1540  {
1541  m_fields[0]->IProductWRTDerivBase(i, F[i], iprod);
1542  Vmath::Vadd(nCoeffs, iprod, 1, outarray, 1, outarray, 1);
1543  }
1544  }
1545 
1546  /**
1547  * Calculate Inner product of the divergence advection form
1548  * \f$(\phi, \nabla \cdot F)\f$, where for example \f$ F = uV \f$.
1549  * @param F Fields.
1550  * @param outarray Storage for result.
1551  */
1552  void EquationSystem::WeakAdvectionDivergenceForm(
1553  const Array<OneD, Array<OneD, NekDouble> > &F,
1554  Array<OneD, NekDouble> &outarray)
1555  {
1556  // Use dimension of Velocity vector to dictate dimension of operation
1557  int ndim = F.num_elements();
1558  int nPointsTot = m_fields[0]->GetNpoints();
1559  Array<OneD, NekDouble> tmp(nPointsTot);
1560  Array<OneD, NekDouble> div(nPointsTot, 0.0);
1561 
1562  // Evaluate the divergence
1563  for (int i = 0; i < ndim; ++i)
1564  {
1565  //m_fields[0]->PhysDeriv(i,F[i],tmp);
1566  m_fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[i],F[i],tmp);
1567  Vmath::Vadd(nPointsTot, tmp, 1, div, 1, div, 1);
1568  }
1569 
1570  m_fields[0]->IProductWRTBase(div, outarray);
1571  }
1572 
1573  /**
1574  * Calculate Inner product of the divergence advection form
1575  * \f$ (\phi, V\cdot \nabla u) \f$
1576  * @param V Fields.
1577  * @param u Fields.
1578  * @param outarray Storage for result.
1579  */
1580  void EquationSystem::WeakAdvectionNonConservativeForm(
1581  const Array<OneD, Array<OneD, NekDouble> > &V,
1582  const Array<OneD, const NekDouble> &u,
1583  Array<OneD, NekDouble> &outarray,
1584  bool UseContCoeffs)
1585  {
1586  // use dimension of Velocity vector to dictate dimension of operation
1587  int ndim = V.num_elements();
1588 
1589  int nPointsTot = m_fields[0]->GetNpoints();
1590  Array<OneD, NekDouble> tmp(nPointsTot);
1591  Array<OneD, NekDouble> wk(ndim * nPointsTot, 0.0);
1592 
1593  AdvectionNonConservativeForm(V, u, tmp, wk);
1594 
1595  if (UseContCoeffs)
1596  {
1597  m_fields[0]->IProductWRTBase(tmp, outarray,
1598  MultiRegions::eGlobal);
1599  }
1600  else
1601  {
1602  m_fields[0]->IProductWRTBase_IterPerExp(tmp, outarray);
1603  }
1604  }
1605 
1606  /**
1607  * Calculate the inner product \f$ V\cdot \nabla u \f$
1608  * @param V Fields.
1609  * @param u Fields.
1610  * @param outarray Storage for result.
1611  * @param wk Workspace.
1612  */
1613  void EquationSystem::AdvectionNonConservativeForm(
1614  const Array<OneD, Array<OneD, NekDouble> > &V,
1615  const Array<OneD, const NekDouble> &u,
1616  Array<OneD, NekDouble> &outarray,
1617  Array<OneD, NekDouble> &wk)
1618  {
1619  // Use dimension of Velocity vector to dictate dimension of operation
1620  int ndim = V.num_elements();
1621  //int ndim = m_expdim;
1622 
1623  // ToDo: here we should add a check that V has right dimension
1624 
1625  int nPointsTot = m_fields[0]->GetNpoints();
1626  Array<OneD, NekDouble> grad0,grad1,grad2;
1627 
1628  // Check to see if wk space is defined
1629  if (wk.num_elements())
1630  {
1631  grad0 = wk;
1632  }
1633  else
1634  {
1635  grad0 = Array<OneD, NekDouble> (nPointsTot);
1636  }
1637 
1638  // Evaluate V\cdot Grad(u)
1639  switch(ndim)
1640  {
1641  case 1:
1642  m_fields[0]->PhysDeriv(u,grad0);
1643  Vmath::Vmul(nPointsTot, grad0, 1, V[0], 1, outarray,1);
1644  break;
1645  case 2:
1646  grad1 = Array<OneD, NekDouble> (nPointsTot);
1647  m_fields[0]->PhysDeriv(u, grad0, grad1);
1648  Vmath::Vmul (nPointsTot, grad0, 1, V[0], 1, outarray, 1);
1649  Vmath::Vvtvp(nPointsTot, grad1, 1, V[1], 1,
1650  outarray, 1, outarray, 1);
1651  break;
1652  case 3:
1653  grad1 = Array<OneD, NekDouble> (nPointsTot);
1654  grad2 = Array<OneD, NekDouble> (nPointsTot);
1655  m_fields[0]->PhysDeriv(u,grad0,grad1,grad2);
1656  Vmath::Vmul (nPointsTot, grad0, 1, V[0], 1, outarray, 1);
1657  Vmath::Vvtvp(nPointsTot, grad1, 1, V[1], 1,
1658  outarray, 1, outarray, 1);
1659  Vmath::Vvtvp(nPointsTot, grad2, 1, V[2], 1,
1660  outarray, 1, outarray, 1);
1661  break;
1662  default:
1663  ASSERTL0(false,"dimension unknown");
1664  }
1665  }
1666 
1667  /**
1668  * @brief Calculate weak DG advection in the form \f$ \langle\phi,
1669  * \hat{F}\cdot n\rangle - (\nabla \phi \cdot F) \f$
1670  *
1671  * @param InField Fields.
1672  * @param OutField Storage for result.
1673  * @param NumericalFluxIncludesNormal Default: true.
1674  * @param InFieldIsPhysSpace Default: false.
1675  * @param nvariables Number of fields.
1676  */
1677  void EquationSystem::WeakDGAdvection(
1678  const Array<OneD, Array<OneD, NekDouble> >& InField,
1679  Array<OneD, Array<OneD, NekDouble> >& OutField,
1680  bool NumericalFluxIncludesNormal,
1681  bool InFieldIsInPhysSpace,
1682  int nvariables)
1683  {
1684  int i;
1685  int nVelDim = m_expdim;
1686  int nPointsTot = GetNpoints();
1687  int ncoeffs = GetNcoeffs();
1688  int nTracePointsTot = GetTraceNpoints();
1689 
1690  if (!nvariables)
1691  {
1692  nvariables = m_fields.num_elements();
1693  }
1694 
1695  Array<OneD, Array<OneD, NekDouble> > fluxvector(nVelDim);
1696  Array<OneD, Array<OneD, NekDouble> > physfield (nvariables);
1697 
1698  for(i = 0; i < nVelDim; ++i)
1699  {
1700  fluxvector[i] = Array<OneD, NekDouble>(nPointsTot);
1701  }
1702 
1703  // Get the variables in physical space
1704  // already in physical space
1705  if (InFieldIsInPhysSpace == true)
1706  {
1707  for (i = 0; i < nvariables; ++i)
1708  {
1709  physfield[i] = InField[i];
1710  }
1711  }
1712  // otherwise do a backward transformation
1713  else
1714  {
1715  for(i = 0; i < nvariables; ++i)
1716  {
1717  // Could make this point to m_fields[i]->UpdatePhys();
1718  physfield[i] = Array<OneD, NekDouble>(nPointsTot);
1719  m_fields[i]->BwdTrans(InField[i],physfield[i]);
1720  }
1721  }
1722 
1723  // Get the advection part (without numerical flux)
1724  for (i = 0; i < nvariables; ++i)
1725  {
1726  // Get the ith component of the flux vector in (physical space)
1727  GetFluxVector(i, physfield, fluxvector);
1728 
1729  // Calculate the i^th value of (\grad_i \phi, F)
1730  WeakAdvectionGreensDivergenceForm(fluxvector,OutField[i]);
1731  }
1732 
1733  // Get the numerical flux and add to the modal coeffs
1734  // if the NumericalFluxs function already includes the
1735  // normal in the output
1736  if (NumericalFluxIncludesNormal == true)
1737  {
1738  Array<OneD, Array<OneD, NekDouble> > numflux (nvariables);
1739 
1740  for (i = 0; i < nvariables; ++i)
1741  {
1742  numflux[i] = Array<OneD, NekDouble>(nTracePointsTot);
1743  }
1744 
1745  // Evaluate numerical flux in physical space which may in
1746  // general couple all component of vectors
1747  NumericalFlux(physfield, numflux);
1748 
1749  // Evaulate <\phi, \hat{F}\cdot n> - OutField[i]
1750  for (i = 0; i < nvariables; ++i)
1751  {
1752  Vmath::Neg(ncoeffs,OutField[i],1);
1753  m_fields[i]->AddTraceIntegral(numflux[i],OutField[i]);
1754  m_fields[i]->SetPhysState(false);
1755  }
1756  }
1757  // if the NumericalFlux function does not include the
1758  // normal in the output
1759  else
1760  {
1761  Array<OneD, Array<OneD, NekDouble> > numfluxX (nvariables);
1762  Array<OneD, Array<OneD, NekDouble> > numfluxY (nvariables);
1763 
1764  for (i = 0; i < nvariables; ++i)
1765  {
1766  numfluxX[i] = Array<OneD, NekDouble>(nTracePointsTot);
1767  numfluxY[i] = Array<OneD, NekDouble>(nTracePointsTot);
1768  }
1769 
1770  // Evaluate numerical flux in physical space which may in
1771  // general couple all component of vectors
1772  NumericalFlux(physfield, numfluxX, numfluxY);
1773 
1774  // Evaulate <\phi, \hat{F}\cdot n> - OutField[i]
1775  for(i = 0; i < nvariables; ++i)
1776  {
1777  Vmath::Neg(ncoeffs,OutField[i],1);
1778  m_fields[i]->AddTraceIntegral(numfluxX[i], numfluxY[i],
1779  OutField[i]);
1780  m_fields[i]->SetPhysState(false);
1781  }
1782  }
1783  }
1784 
1785  /**
1786  * Calculate weak DG Diffusion in the LDG form
1787  * \f$ \langle\psi, \hat{u}\cdot n\rangle
1788  * - \langle\nabla\psi \cdot u\rangle
1789  * \langle\phi, \hat{q}\cdot n\rangle - (\nabla \phi \cdot q) \rangle \f$
1790  */
1791  void EquationSystem::WeakDGDiffusion(
1792  const Array<OneD, Array<OneD, NekDouble> >& InField,
1793  Array<OneD, Array<OneD, NekDouble> >& OutField,
1794  bool NumericalFluxIncludesNormal,
1795  bool InFieldIsInPhysSpace)
1796  {
1797  int i, j, k;
1798  int nPointsTot = GetNpoints();
1799  int ncoeffs = GetNcoeffs();
1800  int nTracePointsTot = GetTraceNpoints();
1801  int nvariables = m_fields.num_elements();
1802  int nqvar = 2;
1803 
1804  Array<OneD, NekDouble> qcoeffs (ncoeffs);
1805  Array<OneD, NekDouble> temp (ncoeffs);
1806 
1807  Array<OneD, Array<OneD, NekDouble> > fluxvector (m_spacedim);
1808  Array<OneD, Array<OneD, NekDouble> > ufield (nvariables);
1809 
1810  Array<OneD, Array<OneD, Array<OneD, NekDouble> > > flux (nqvar);
1811  Array<OneD, Array<OneD, Array<OneD, NekDouble> > > qfield (nqvar);
1812 
1813  for (j = 0; j < nqvar; ++j)
1814  {
1815  qfield[j] = Array<OneD, Array<OneD, NekDouble> >(nqvar);
1816  flux[j] = Array<OneD, Array<OneD, NekDouble> >(nqvar);
1817 
1818  for (i = 0; i< nvariables; ++i)
1819  {
1820  ufield[i] = Array<OneD, NekDouble>(nPointsTot, 0.0);
1821  qfield[j][i] = Array<OneD, NekDouble>(nPointsTot, 0.0);
1822  flux[j][i] = Array<OneD, NekDouble>(nTracePointsTot, 0.0);
1823  }
1824  }
1825 
1826  for (k = 0; k < m_spacedim; ++k)
1827  {
1828  fluxvector[k] = Array<OneD, NekDouble>(nPointsTot, 0.0);
1829  }
1830 
1831  // Get the variables in physical space already in physical space
1832  if (InFieldIsInPhysSpace == true)
1833  {
1834  for (i = 0; i < nvariables; ++i)
1835  {
1836  ufield[i] = InField[i];
1837  }
1838  }
1839  // Otherwise do a backward transformation
1840  else
1841  {
1842  for (i = 0; i < nvariables; ++i)
1843  {
1844  // Could make this point to m_fields[i]->UpdatePhys();
1845  ufield[i] = Array<OneD, NekDouble>(nPointsTot);
1846  m_fields[i]->BwdTrans(InField[i],ufield[i]);
1847  }
1848  }
1849 
1850  // ##########################################################
1851  // Compute q_{\eta} and q_{\xi} from su
1852  // Obtain Numerical Fluxes
1853  // ##########################################################
1854  NumFluxforScalar(ufield, flux);
1855 
1856  for (j = 0; j < nqvar; ++j)
1857  {
1858  for (i = 0; i < nvariables; ++i)
1859  {
1860  // Get the ith component of the flux vector in
1861  // (physical space) fluxvector = m_tanbasis * u,
1862  // where m_tanbasis = 2 by m_spacedim by nPointsTot
1863  if (m_tanbasis.num_elements())
1864  {
1865  for (k = 0; k < m_spacedim; ++k)
1866  {
1867  Vmath::Vmul(nPointsTot, m_tanbasis[j][k], 1,
1868  ufield[i], 1, fluxvector[k], 1);
1869  }
1870  }
1871  else
1872  {
1873  GetFluxVector(i, j, ufield, fluxvector);
1874  }
1875 
1876  // Calculate the i^th value of (\grad_i \phi, F)
1877  WeakAdvectionGreensDivergenceForm(fluxvector, qcoeffs);
1878 
1879  Vmath::Neg(ncoeffs,qcoeffs,1);
1880  m_fields[i]->AddTraceIntegral(flux[j][i], qcoeffs);
1881  m_fields[i]->SetPhysState(false);
1882 
1883  // Add weighted mass matrix = M ( \nabla \cdot Tanbasis )
1884 // if(m_gradtan.num_elements())
1885 // {
1886 // MultiRegions::GlobalMatrixKey key(StdRegions::eMass,
1887 // m_gradtan[j]);
1888 // m_fields[i]->MultiRegions::ExpList::GeneralMatrixOp(key,
1889 // InField[i], temp);
1890 // Vmath::Svtvp(ncoeffs, -1.0, temp, 1, qcoeffs, 1,
1891 // qcoeffs, 1);
1892 // }
1893 
1894  //Multiply by the inverse of mass matrix
1895  m_fields[i]->MultiplyByElmtInvMass(qcoeffs, qcoeffs);
1896 
1897  // Back to physical space
1898  m_fields[i]->BwdTrans(qcoeffs, qfield[j][i]);
1899  }
1900  }
1901 
1902 
1903  // ##########################################################
1904  // Compute u from q_{\eta} and q_{\xi}
1905  // ##########################################################
1906 
1907  // Obtain Numerical Fluxes
1908  NumFluxforVector(ufield, qfield, flux[0]);
1909 
1910  for (i = 0; i < nvariables; ++i)
1911  {
1912  // L = L(tan_eta) q_eta + L(tan_xi) q_xi
1913  OutField[i] = Array<OneD, NekDouble>(ncoeffs, 0.0);
1914  temp = Array<OneD, NekDouble>(ncoeffs, 0.0);
1915 
1916  if (m_tanbasis.num_elements())
1917  {
1918  for (j = 0; j < nqvar; ++j)
1919  {
1920  for (k = 0; k < m_spacedim; ++k)
1921  {
1922  Vmath::Vmul(nPointsTot, m_tanbasis[j][k], 1,
1923  qfield[j][i], 1, fluxvector[k], 1);
1924  }
1925 
1926  WeakAdvectionGreensDivergenceForm(fluxvector, temp);
1927  Vmath::Vadd(ncoeffs, temp, 1, OutField[i], 1,
1928  OutField[i], 1);
1929  }
1930  }
1931  else
1932  {
1933  for (k = 0; k < m_spacedim; ++k)
1934  {
1935  Vmath::Vcopy(nPointsTot, qfield[k][i], 1,
1936  fluxvector[k], 1);
1937  }
1938 
1939  WeakAdvectionGreensDivergenceForm(fluxvector, OutField[i]);
1940  }
1941 
1942  // Evaulate <\phi, \hat{F}\cdot n> - OutField[i]
1943  Vmath::Neg(ncoeffs,OutField[i],1);
1944  m_fields[i]->AddTraceIntegral(flux[0][i], OutField[i]);
1945  m_fields[i]->SetPhysState(false);
1946  }
1947  }
1948 
1949  /**
1950  * Write the n-th checkpoint file.
1951  * @param n The index of the checkpoint file.
1952  */
1953  void EquationSystem::Checkpoint_Output(const int n)
1954  {
1955  std::string outname = m_sessionName + "_" +
1956  boost::lexical_cast<std::string>(n);
1957 
1958  WriteFld(outname + ".chk");
1959  }
1960 
1961  /**
1962  * Write the n-th checkpoint file.
1963  * @param n The index of the checkpoint file.
1964  */
1965  void EquationSystem::Checkpoint_Output(
1966  const int n,
1967  MultiRegions::ExpListSharedPtr &field,
1968  std::vector<Array<OneD, NekDouble> > &fieldcoeffs,
1969  std::vector<std::string> &variables)
1970  {
1971  char chkout[16] = "";
1972  sprintf(chkout, "%d", n);
1973  std::string outname = m_sessionName + "_" + chkout + ".chk";
1974  WriteFld(outname, field, fieldcoeffs, variables);
1975  }
1976 
1977  /**
1978  * Write the n-th base flow into a .chk file
1979  * @param n The index of the base flow file.
1980  */
1981  void EquationSystem::Checkpoint_BaseFlow(const int n)
1982  {
1983  std::string outname = m_sessionName + "_BaseFlow_" +
1984  boost::lexical_cast<std::string>(n);
1985 
1986  WriteFld(outname + ".chk");
1987  }
1988 
1989  /**
1990  * Writes the field data to a file with the given filename.
1991  * @param outname Filename to write to.
1992  */
1993  void EquationSystem::WriteFld(const std::string &outname)
1994  {
1995  std::vector<Array<OneD, NekDouble> > fieldcoeffs(
1996  m_fields.num_elements());
1997  std::vector<std::string> variables(m_fields.num_elements());
1998 
1999  for (int i = 0; i < m_fields.num_elements(); ++i)
2000  {
2001  if (m_fields[i]->GetNcoeffs() == m_fields[0]->GetNcoeffs())
2002  {
2003  fieldcoeffs[i] = m_fields[i]->UpdateCoeffs();
2004  }
2005  else
2006  {
2007  fieldcoeffs[i] = Array<OneD,NekDouble>(m_fields[0]->
2008  GetNcoeffs());
2009  m_fields[0]->ExtractCoeffsToCoeffs(m_fields[i],
2010  m_fields[i]->GetCoeffs(),
2011  fieldcoeffs[i]);
2012  }
2013  variables[i] = m_boundaryConditions->GetVariable(i);
2014  }
2015 
2016  v_ExtraFldOutput(fieldcoeffs, variables);
2017 
2018  WriteFld(outname, m_fields[0], fieldcoeffs, variables);
2019  }
2020 
2021 
2022 
2023  /**
2024  * Writes the field data to a file with the given filename.
2025  * @param outname Filename to write to.
2026  * @param field ExpList on which data is based.
2027  * @param fieldcoeffs An array of array of expansion coefficients.
2028  * @param variables An array of variable names.
2029  */
2030  void EquationSystem::WriteFld(
2031  const std::string &outname,
2032  MultiRegions::ExpListSharedPtr &field,
2033  std::vector<Array<OneD, NekDouble> > &fieldcoeffs,
2034  std::vector<std::string> &variables)
2035  {
2036  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef
2037  = field->GetFieldDefinitions();
2038  std::vector<std::vector<NekDouble> > FieldData(FieldDef.size());
2039 
2040  // Copy Data into FieldData and set variable
2041  for(int j = 0; j < fieldcoeffs.size(); ++j)
2042  {
2043  for(int i = 0; i < FieldDef.size(); ++i)
2044  {
2045  // Could do a search here to find correct variable
2046  FieldDef[i]->m_fields.push_back(variables[j]);
2047  field->AppendFieldData(FieldDef[i], FieldData[i],
2048  fieldcoeffs[j]);
2049  }
2050  }
2051 
2052  // Update time in field info if required
2053  if(m_fieldMetaDataMap.find("Time") != m_fieldMetaDataMap.end())
2054  {
2055  m_fieldMetaDataMap["Time"] = boost::lexical_cast<std::string>(m_time);
2056  }
2057 
2058  m_fld->Write(outname, FieldDef, FieldData, m_fieldMetaDataMap);
2059  }
2060 
2061 
2062  /**
2063  * Import field from infile and load into \a m_fields. This routine will
2064  * also perform a \a BwdTrans to ensure data is in both the physical and
2065  * coefficient storage.
2066  * @param infile Filename to read.
2067  */
2068  void EquationSystem::ImportFld(
2069  const std::string &infile,
2070  Array<OneD, MultiRegions::ExpListSharedPtr> &pFields)
2071  {
2072  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef;
2073  std::vector<std::vector<NekDouble> > FieldData;
2074 
2075  m_fld->Import(infile,FieldDef,FieldData);
2076 
2077  // Copy FieldData into m_fields
2078  for(int j = 0; j < pFields.num_elements(); ++j)
2079  {
2080  Vmath::Zero(pFields[j]->GetNcoeffs(),
2081  pFields[j]->UpdateCoeffs(),1);
2082 
2083  for(int i = 0; i < FieldDef.size(); ++i)
2084  {
2085  ASSERTL1(FieldDef[i]->m_fields[j] ==
2086  m_session->GetVariable(j),
2087  std::string("Order of ") + infile
2088  + std::string(" data and that defined in "
2089  "m_boundaryconditions differs"));
2090 
2091  pFields[j]->ExtractDataToCoeffs(FieldDef[i], FieldData[i],
2092  FieldDef[i]->m_fields[j],
2093  pFields[j]->UpdateCoeffs());
2094  }
2095  pFields[j]->BwdTrans(pFields[j]->GetCoeffs(),
2096  pFields[j]->UpdatePhys());
2097  }
2098  }
2099 
2100 
2101 
2102  /**
2103  * Import field from infile and load into \a m_fields. This routine will
2104  * also perform a \a BwdTrans to ensure data is in both the physical and
2105  * coefficient storage.
2106  * @param infile Filename to read.
2107  * If optionan \a ndomains is specified it assumes we loop over nodmains for each nvariables.
2108  */
2109  void EquationSystem::ImportFldToMultiDomains(
2110  const std::string &infile,
2111  Array<OneD, MultiRegions::ExpListSharedPtr> &pFields,
2112  const int ndomains)
2113  {
2114  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef;
2115  std::vector<std::vector<NekDouble> > FieldData;
2116 
2117  LibUtilities::Import(infile,FieldDef,FieldData);
2118 
2119  int nvariables = GetNvariables();
2120 
2121  ASSERTL0(ndomains*nvariables == pFields.num_elements(),"Number of fields does not match the number of variables and domains");
2122 
2123  // Copy FieldData into m_fields
2124  for(int j = 0; j < ndomains; ++j)
2125  {
2126  for(int i = 0; i < nvariables; ++i)
2127  {
2128  Vmath::Zero(pFields[j*nvariables+i]->GetNcoeffs(),pFields[j*nvariables+i]->UpdateCoeffs(),1);
2129 
2130  for(int n = 0; n < FieldDef.size(); ++n)
2131  {
2132  ASSERTL1(FieldDef[n]->m_fields[i] == m_session->GetVariable(i),
2133  std::string("Order of ") + infile
2134  + std::string(" data and that defined in "
2135  "m_boundaryconditions differs"));
2136 
2137  pFields[j*nvariables+i]->ExtractDataToCoeffs(FieldDef[n], FieldData[n],
2138  FieldDef[n]->m_fields[i],
2139  pFields[j*nvariables+i]->UpdateCoeffs());
2140  }
2141  pFields[j*nvariables+i]->BwdTrans(pFields[j*nvariables+i]->GetCoeffs(),
2142  pFields[j*nvariables+i]->UpdatePhys());
2143  }
2144  }
2145  }
2146 
2147  /**
2148  * Import field from infile and load into \a pField. This routine will
2149  * also perform a \a BwdTrans to ensure data is in both the physical and
2150  * coefficient storage.
2151  */
2152  void EquationSystem::ImportFld(
2153  const std::string &infile,
2154  MultiRegions::ExpListSharedPtr &pField,
2155  std::string &pFieldName)
2156  {
2157  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef;
2158  std::vector<std::vector<NekDouble> > FieldData;
2159 
2160  m_fld->Import(infile,FieldDef,FieldData);
2161  int idx = -1;
2162 
2163  Vmath::Zero(pField->GetNcoeffs(),pField->UpdateCoeffs(),1);
2164 
2165  for(int i = 0; i < FieldDef.size(); ++i)
2166  {
2167  // find the index of the required field in the file.
2168  for(int j = 0; j < FieldData.size(); ++j)
2169  {
2170  if (FieldDef[i]->m_fields[j] == pFieldName)
2171  {
2172  idx = j;
2173  }
2174  }
2175  ASSERTL1(idx >= 0, "Field " + pFieldName + " not found.");
2176 
2177  pField->ExtractDataToCoeffs(FieldDef[i], FieldData[i],
2178  FieldDef[i]->m_fields[idx],
2179  pField->UpdateCoeffs());
2180  }
2181  pField->BwdTrans(pField->GetCoeffs(), pField->UpdatePhys());
2182  }
2183 
2184  /**
2185  * Import field from infile and load into the array \a coeffs.
2186  *
2187  * @param infile Filename to read.
2188  * @param fieldStr an array of string identifying fields to be imported
2189  * @param coeffs and array of array of coefficients to store imported data
2190  */
2191  void EquationSystem::ImportFld(
2192  const std::string &infile,
2193  std::vector< std::string> &fieldStr,
2194  Array<OneD, Array<OneD, NekDouble> > &coeffs)
2195  {
2196 
2197  ASSERTL0(fieldStr.size() <= coeffs.num_elements(),
2198  "length of fieldstr should be the same as pFields");
2199 
2200  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef;
2201  std::vector<std::vector<NekDouble> > FieldData;
2202 
2203  m_fld->Import(infile,FieldDef,FieldData);
2204 
2205  // Copy FieldData into m_fields
2206  for(int j = 0; j < fieldStr.size(); ++j)
2207  {
2208  Vmath::Zero(coeffs[j].num_elements(),coeffs[j],1);
2209  for(int i = 0; i < FieldDef.size(); ++i)
2210  {
2211  m_fields[0]->ExtractDataToCoeffs(FieldDef[i], FieldData[i],
2212  fieldStr[j], coeffs[j]);
2213  }
2214  }
2215  }
2216 
2217  /**
2218  * Write out a summary of the session data.
2219  * @param out Output stream to write data to.
2220  */
2221  void EquationSystem::SessionSummary(SummaryList& s)
2222  {
2223  AddSummaryItem(s, "EquationType", m_session->GetSolverInfo("EQTYPE"));
2224  AddSummaryItem(s, "Session Name", m_sessionName);
2225  AddSummaryItem(s, "Spatial Dim.", m_spacedim);
2226  AddSummaryItem(s, "Max SEM Exp. Order", m_fields[0]->EvalBasisNumModesMax());
2227  if(m_HomogeneousType == eHomogeneous1D)
2228  {
2229  AddSummaryItem(s, "Quasi-3D", "Homogeneous in z-direction");
2230  AddSummaryItem(s, "Expansion Dim.", m_expdim + 1);
2231  AddSummaryItem(s, "Num. Hom. Modes (z)", m_npointsZ);
2232  AddSummaryItem(s, "Hom. length (LZ)", "m_LhomZ");
2233  AddSummaryItem(s, "FFT Type", m_useFFT ? "FFTW" : "MVM");
2234  AddSummaryItem(s, "Selected Mode", m_MultipleModes
2235  ? boost::lexical_cast<string>(m_NumMode) : "ALL");
2236  }
2237  else if(m_HomogeneousType == eHomogeneous2D)
2238  {
2239  AddSummaryItem(s, "Quasi-3D", "Homogeneous in yz-plane");
2240  AddSummaryItem(s, "Expansion Dim.", m_expdim + 2);
2241  AddSummaryItem(s, "Num. Hom. Modes (y)", m_npointsY);
2242  AddSummaryItem(s, "Num. Hom. Modes (z)", m_npointsZ);
2243  AddSummaryItem(s, "Hom. length (LY)", "m_LhomY");
2244  AddSummaryItem(s, "Hom. length (LZ)", "m_LhomZ");
2245  AddSummaryItem(s, "FFT Type", m_useFFT ? "FFTW" : "MVM");
2246  }
2247  else
2248  {
2249  AddSummaryItem(s, "Expansion Dim.", m_expdim);
2250  }
2251 
2252  if (m_session->DefinesSolverInfo("UpwindType"))
2253  {
2254  AddSummaryItem(s, "Riemann Solver",
2255  m_session->GetSolverInfo("UpwindType"));
2256  }
2257 
2258  if (m_session->DefinesSolverInfo("AdvectionType"))
2259  {
2260  std::string AdvectionType;
2261  AdvectionType = m_session->GetSolverInfo("AdvectionType");
2262  AddSummaryItem(s, "Advection Type", GetAdvectionFactory().
2263  GetClassDescription(AdvectionType));
2264  }
2265 
2266  if (m_projectionType == MultiRegions::eGalerkin)
2267  {
2268  AddSummaryItem(s, "Projection Type", "Continuous Galerkin");
2269  }
2270  else if (m_projectionType == MultiRegions::eDiscontinuous)
2271  {
2272  AddSummaryItem(s, "Projection Type", "Discontinuous Galerkin");
2273  }
2274  else if (m_projectionType == MultiRegions::eMixed_CG_Discontinuous)
2275  {
2276  AddSummaryItem(s, "Projection Type",
2277  "Mixed Continuous Galerkin and Discontinuous");
2278  }
2279 
2280  if (m_session->DefinesSolverInfo("DiffusionType"))
2281  {
2282  std::string DiffusionType;
2283  DiffusionType = m_session->GetSolverInfo("DiffusionType");
2284  AddSummaryItem(s, "Diffusion Type", GetDiffusionFactory().
2285  GetClassDescription(DiffusionType));
2286  }
2287  }
2288 
2289  /**
2290  * Performs a case-insensitive string comparison (from web).
2291  * @param s1 First string to compare.
2292  * @param s2 Second string to compare.
2293  * @returns 0 if the strings match.
2294  */
2295  int EquationSystem::NoCaseStringCompare(
2296  const string & s1,
2297  const string& s2)
2298  {
2299  //if (s1.size() < s2.size()) return -1;
2300  //if (s1.size() > s2.size()) return 1;
2301 
2302  string::const_iterator it1=s1.begin();
2303  string::const_iterator it2=s2.begin();
2304 
2305  // Stop when either string's end has been reached
2306  while ( (it1!=s1.end()) && (it2!=s2.end()) )
2307  {
2308  if(::toupper(*it1) != ::toupper(*it2)) //letters differ?
2309  {
2310  // Return -1 to indicate smaller than, 1 otherwise
2311  return (::toupper(*it1) < ::toupper(*it2)) ? -1 : 1;
2312  }
2313 
2314  // Proceed to the next character in each string
2315  ++it1;
2316  ++it2;
2317  }
2318 
2319  size_t size1=s1.size();
2320  size_t size2=s2.size();// cache lengths
2321 
2322  // Return -1, 0 or 1 according to strings' lengths
2323  if (size1==size2)
2324  {
2325  return 0;
2326  }
2327 
2328  return (size1 < size2) ? -1 : 1;
2329  }
2330 
2331  Array<OneD, bool> EquationSystem::v_GetSystemSingularChecks()
2332  {
2333  return Array<OneD, bool>(m_session->GetVariables().size(), false);
2334  }
2335 
2336  void EquationSystem::v_GetFluxVector(
2337  const int i, Array<OneD,
2338  Array<OneD, NekDouble> > &physfield,
2339  Array<OneD, Array<OneD, NekDouble> > &flux)
2340  {
2341  ASSERTL0(false, "v_GetFluxVector: This function is not valid "
2342  "for the Base class");
2343  }
2344 
2345  void EquationSystem::v_GetFluxVector(
2346  const int i, const int j,
2347  Array<OneD, Array<OneD, NekDouble> > &physfield,
2348  Array<OneD, Array<OneD, NekDouble> > &flux)
2349  {
2350  ASSERTL0(false, "v_GetqFluxVector: This function is not valid "
2351  "for the Base class");
2352  }
2353 
2354  void EquationSystem::v_GetFluxVector(
2355  const int i, Array<OneD,
2356  Array<OneD, NekDouble> > &physfield,
2357  Array<OneD, Array<OneD, NekDouble> > &fluxX,
2358  Array<OneD, Array<OneD, NekDouble> > &fluxY)
2359  {
2360  ASSERTL0(false, "v_GetFluxVector: This function is not valid "
2361  "for the Base class");
2362  }
2363 
2364  void EquationSystem::v_NumericalFlux(
2365  Array<OneD, Array<OneD, NekDouble> > &physfield,
2366  Array<OneD, Array<OneD, NekDouble> > &numflux)
2367  {
2368  ASSERTL0(false, "v_NumericalFlux: This function is not valid "
2369  "for the Base class");
2370  }
2371 
2372  void EquationSystem::v_NumericalFlux(
2373  Array<OneD, Array<OneD, NekDouble> > &physfield,
2374  Array<OneD, Array<OneD, NekDouble> > &numfluxX,
2375  Array<OneD, Array<OneD, NekDouble> > &numfluxY )
2376  {
2377  ASSERTL0(false, "v_NumericalFlux: This function is not valid "
2378  "for the Base class");
2379  }
2380 
2381  void EquationSystem::v_NumFluxforScalar(
2382  const Array<OneD, Array<OneD, NekDouble> > &ufield,
2383  Array<OneD, Array<OneD, Array<OneD, NekDouble> > > &uflux)
2384  {
2385  ASSERTL0(false, "v_NumFluxforScalar: This function is not valid "
2386  "for the Base class");
2387  }
2388 
2389  void EquationSystem::v_NumFluxforVector(
2390  const Array<OneD, Array<OneD, NekDouble> > &ufield,
2391  Array<OneD, Array<OneD, Array<OneD, NekDouble> > > &qfield,
2392  Array<OneD, Array<OneD, NekDouble > > &qflux)
2393  {
2394  ASSERTL0(false, "v_NumFluxforVector: This function is not valid "
2395  "for the Base class");
2396  }
2397 
2398  MultiRegions::ExpListSharedPtr EquationSystem::v_GetPressure()
2399  {
2400  ASSERTL0(false, "This function is not valid for the Base class");
2401  MultiRegions::ExpListSharedPtr null;
2402  return null;
2403  }
2404 
2405  void EquationSystem::v_ExtraFldOutput(
2406  std::vector<Array<OneD, NekDouble> > &fieldcoeffs,
2407  std::vector<std::string> &variables)
2408  {
2409  }
2410 
2411  }
2412 }
Nektar::SolverUtils::EquationSystem::v_GenerateSummary
virtual SOLVER_UTILS_EXPORT void v_GenerateSummary(SummaryList &l)
Virtual function for generating summary information.
Definition: EquationSystem.cpp:1479
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:906
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:1613
ExpList3DHomogeneous1D.h
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:135
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
Nektar::SolverUtils::EquationSystem::m_SingleMode
bool m_SingleMode
Flag to determine if single homogeneous mode is used.
Definition: EquationSystem.h:450
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:432
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:1580
Nektar::SolverUtils::EquationSystem::m_timestep
NekDouble m_timestep
Time step size.
Definition: EquationSystem.h:436
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:789
Nektar::SolverUtils::EquationSystem::m_LhomZ
NekDouble m_LhomZ
physical length in Z direction (if homogeneous)
Definition: EquationSystem.h:498
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:476
Nektar::SolverUtils::EquationSystem::m_expdim
int m_expdim
Expansion dimension.
Definition: EquationSystem.h:448
Nektar::SolverUtils::EquationSystem::Checkpoint_Output
SOLVER_UTILS_EXPORT void Checkpoint_Output(const int n)
Write checkpoint file of m_fields.
Definition: EquationSystem.cpp:1953
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:468
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:418
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:1528
Nektar::SolverUtils::EquationSystem::m_checktime
NekDouble m_checktime
Time between checkpoints.
Definition: EquationSystem.h:440
Nektar::SolverUtils::EquationSystem::m_useFFT
bool m_useFFT
Flag to determine if FFT is used for homogeneous transform.
Definition: EquationSystem.h:456
Nektar::SolverUtils::EquationSystem::m_LhomY
NekDouble m_LhomY
physical length in Y direction (if homogeneous)
Definition: EquationSystem.h:497
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:466
Nektar::SolverUtils::EquationSystem::GetNumExpModes
SOLVER_UTILS_EXPORT int GetNumExpModes()
Definition: EquationSystem.h:778
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:2405
Nektar::SolverUtils::EquationSystem::PrintProgressbar
SOLVER_UTILS_EXPORT void PrintProgressbar(const int position, const int goal) const
Definition: EquationSystem.h:614
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:502
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:430
Nektar::SolverUtils::EquationSystem::m_NumMode
int m_NumMode
Mode to use in case of single mode analysis.
Definition: EquationSystem.h:506
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:1552
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:444
Nektar::LibUtilities::PtsFieldSharedPtr
boost::shared_ptr< PtsField > PtsFieldSharedPtr
Definition: PtsField.h:247
Nektar::SolverUtils::EquationSystem::m_comm
LibUtilities::CommSharedPtr m_comm
Communicator.
Definition: EquationSystem.h:410
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:434
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:474
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:442
Nektar::SolverUtils::EquationSystem::NoCaseStringCompare
SOLVER_UTILS_EXPORT int NoCaseStringCompare(const string &s1, const string &s2)
Perform a case-insensitive string comparison.
Definition: EquationSystem.cpp:2295
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:470
Nektar::SolverUtils::EquationSystem::m_HomoDirec
int m_HomoDirec
number of homogenous directions
Definition: EquationSystem.h:504
Nektar::SolverUtils::EquationSystem::FwdTransFields
SOLVER_UTILS_EXPORT void FwdTransFields()
Definition: EquationSystem.cpp:1508
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:1791
Nektar::LibUtilities::eFourierEvenlySpaced
1D Evenly-spaced points using Fourier Fit
Definition: PointsType.h:64
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:2331
Nektar::SolverUtils::EquationSystem::NumericalFlux
SOLVER_UTILS_EXPORT void NumericalFlux(Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &numflux)
Definition: EquationSystem.h:893
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:501
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:1472
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:2336
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:1488
Nektar::SolverUtils::EquationSystem::m_boundaryConditions
SpatialDomains::BoundaryConditionsSharedPtr m_boundaryConditions
Pointer to boundary conditions object.
Definition: EquationSystem.h:426
Nektar::SolverUtils::EquationSystem::m_HalfMode
bool m_HalfMode
Flag to determine if half homogeneous mode is used.
Definition: EquationSystem.h:452
Nektar::SolverUtils::EquationSystem::m_spacedim
int m_spacedim
Spatial dimension (>= expansion dim).
Definition: EquationSystem.h:446
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:478
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:871
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:461
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:2109
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:1993
Nektar::SolverUtils::EquationSystem::GetNpoints
SOLVER_UTILS_EXPORT int GetNpoints()
Definition: EquationSystem.h:834
Nektar::SolverUtils::EquationSystem::m_fields
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
Definition: EquationSystem.h:420
Nektar::SolverUtils::EquationSystem::m_session
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
Definition: EquationSystem.h:412
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:422
Equation.h
Nektar::SolverUtils::EquationSystem::GetTraceNpoints
SOLVER_UTILS_EXPORT int GetTraceNpoints()
Definition: EquationSystem.h:804
Nektar::SolverUtils::EquationSystem::GetNcoeffs
SOLVER_UTILS_EXPORT int GetNcoeffs()
Definition: EquationSystem.h:768
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:414
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:2389
Nektar::SolverUtils::EquationSystem::v_GetPressure
virtual SOLVER_UTILS_EXPORT MultiRegions::ExpListSharedPtr v_GetPressure(void)
Definition: EquationSystem.cpp:2398
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:913
Nektar::SolverUtils::EquationSystem::ZeroPhysFields
SOLVER_UTILS_EXPORT void ZeroPhysFields()
Definition: EquationSystem.cpp:1496
Nektar::SolverUtils::EquationSystem::SessionSummary
SOLVER_UTILS_EXPORT void SessionSummary(SummaryList &vSummary)
Write out a session summary.
Definition: EquationSystem.cpp:2221
Nektar::SolverUtils::EquationSystem::m_graph
SpatialDomains::MeshGraphSharedPtr m_graph
Pointer to graph defining mesh.
Definition: EquationSystem.h:428
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:2381
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:1677
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:1981
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:481
ASSERTL1
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:165
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:1016
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::m_MultipleModes
bool m_MultipleModes
Flag to determine if use multiple homogenenous modes are used.
Definition: EquationSystem.h:454
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:1464
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:416
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:2364
Nektar::SolverUtils::EquationSystem::m_HomogeneousType
enum HomogeneousType m_HomogeneousType
Definition: EquationSystem.h:494
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:2068
Generated on Tue Dec 8 2015 23:41:28 for Nektar++ by   doxygen 1.8.9.1