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