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