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