Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
CompressibleFlowSystem.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File CompressibleFlowSystem.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: Auxiliary functions for the compressible flow system
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 #include <iostream>
36 #include <iomanip>
37 
38 #include <CompressibleFlowSolver/EquationSystems/CompressibleFlowSystem.h>
39 #include <LocalRegions/TriExp.h>
40 #include <LocalRegions/QuadExp.h>
41 #include <LocalRegions/HexExp.h>
42 #include <MultiRegions/ExpList.h>
43 #include <LibUtilities/Foundations/InterpCoeff.h>
44 #include <MultiRegions/AssemblyMap/AssemblyMapDG.h>
45 
46 using namespace std;
47 
48 namespace Nektar
49 {
50  string CompressibleFlowSystem::className =
51  SolverUtils::GetEquationSystemFactory().RegisterCreatorFunction(
52  "CompressibleFlowSystem",
53  CompressibleFlowSystem::create,
54  "Auxiliary functions for the compressible flow system.");
55 
56  CompressibleFlowSystem::CompressibleFlowSystem(
57  const LibUtilities::SessionReaderSharedPtr& pSession)
58  : UnsteadySystem(pSession)
59  {
60  }
61 
62  /**
63  * @brief Initialization object for CompressibleFlowSystem class.
64  */
65  void CompressibleFlowSystem::v_InitObject()
66  {
67  UnsteadySystem::v_InitObject();
68 
69  ASSERTL0(m_session->DefinesSolverInfo("UPWINDTYPE"),
70  "No UPWINDTYPE defined in session.");
71 
72  // Do not forwards transform initial condition
73  m_homoInitialFwd = false;
74 
75  // Set up locations of velocity vector.
76  m_vecLocs = Array<OneD, Array<OneD, NekDouble> >(1);
77  m_vecLocs[0] = Array<OneD, NekDouble>(m_spacedim);
78  for (int i = 0; i < m_spacedim; ++i)
79  {
80  m_vecLocs[0][i] = 1 + i;
81  }
82 
83  // Get gamma parameter from session file.
84  ASSERTL0(m_session->DefinesParameter("Gamma"),
85  "Compressible flow sessions must define a Gamma parameter.");
86  m_session->LoadParameter("Gamma", m_gamma, 1.4);
87 
88  // Get E0 parameter from session file.
89  ASSERTL0(m_session->DefinesParameter("pInf"),
90  "Compressible flow sessions must define a pInf parameter.");
91  m_session->LoadParameter("pInf", m_pInf, 101325);
92 
93  // Get rhoInf parameter from session file.
94  ASSERTL0(m_session->DefinesParameter("rhoInf"),
95  "Compressible flow sessions must define a rhoInf parameter.");
96  m_session->LoadParameter("rhoInf", m_rhoInf, 1.225);
97 
98  // Get uInf parameter from session file.
99  ASSERTL0(m_session->DefinesParameter("uInf"),
100  "Compressible flow sessions must define a uInf parameter.");
101  m_session->LoadParameter("uInf", m_uInf, 0.1);
102 
103  m_UInf = m_uInf;
104 
105  // Get vInf parameter from session file.
106  if (m_spacedim == 2 || m_spacedim == 3)
107  {
108  ASSERTL0(m_session->DefinesParameter("vInf"),
109  "Compressible flow sessions must define a vInf parameter"
110  "for 2D/3D problems.");
111  m_session->LoadParameter("vInf", m_vInf, 0.0);
112  m_UInf = sqrt(m_uInf*m_uInf + m_vInf*m_vInf);
113  }
114 
115  // Get wInf parameter from session file.
116  if (m_spacedim == 3)
117  {
118  ASSERTL0(m_session->DefinesParameter("wInf"),
119  "Compressible flow sessions must define a wInf parameter"
120  "for 3D problems.");
121  m_session->LoadParameter("wInf", m_wInf, 0.0);
122  m_UInf = sqrt(m_uInf*m_uInf + m_vInf*m_vInf + m_wInf*m_wInf);
123  }
124 
125  m_session->LoadParameter ("GasConstant", m_gasConstant, 287.058);
126  m_session->LoadParameter ("Twall", m_Twall, 300.15);
127  m_session->LoadSolverInfo("ViscosityType", m_ViscosityType, "Constant");
128  m_session->LoadParameter ("mu", m_mu, 1.78e-05);
129  m_session->LoadParameter ("Skappa", m_Skappa, -2.048);
130  m_session->LoadParameter ("Kappa", m_Kappa, 0.0);
131  m_session->LoadParameter ("mu0", m_mu0, 1.0);
132  m_session->LoadParameter ("FL", m_FacL, 0.0);
133  m_session->LoadParameter ("FH", m_FacH, 0.0);
134  m_session->LoadParameter ("hFactor", m_hFactor, 1.0);
135  m_session->LoadParameter ("epsMax", m_eps_max, 1.0);
136  m_session->LoadParameter ("C1", m_C1, 3.0);
137  m_session->LoadParameter ("C2", m_C2, 5.0);
138  m_session->LoadSolverInfo("ShockCaptureType",
139  m_shockCaptureType, "Off");
140  m_session->LoadParameter ("thermalConductivity",
141  m_thermalConductivity, 0.0257);
142 
143  m_EqTypeStr = m_session->GetSolverInfo("EQTYPE");
144 
145  m_Cp = m_gamma / (m_gamma - 1.0) * m_gasConstant;
146  m_Prandtl = m_Cp * m_mu / m_thermalConductivity;
147 
148  m_session->LoadParameter("amplitude", m_amplitude, 0.001);
149  m_session->LoadParameter("omega", m_omega, 1.0);
150  m_session->LoadParameter("SteadyStateTol", m_steadyStateTol, 0.0);
151 
152  // Forcing terms for the sponge region
153  m_forcing = SolverUtils::Forcing::Load(m_session, m_fields,
154  m_fields.num_elements());
155 
156  // Loop over Boundary Regions for PressureOutflowFileBC
157  int nvariables = m_fields.num_elements();
158  Array<OneD, Array<OneD, NekDouble> > tmpStorage(nvariables);
159  for (int n = 0; n < m_fields[0]->GetBndConditions().num_elements(); ++n)
160  {
161  // PressureOutflowFile Boundary Condition
162  if (m_fields[0]->GetBndConditions()[n]->GetUserDefined() ==
163  "PressureOutflowFile")
164  {
165  int numBCPts = m_fields[0]->
166  GetBndCondExpansions()[n]->GetNpoints();
167  m_pressureStorage = Array<OneD, NekDouble>(numBCPts, 0.0);
168  for (int i = 0; i < nvariables; ++i)
169  {
170  tmpStorage[i] = Array<OneD, NekDouble>(numBCPts, 0.0);
171 
172  Vmath::Vcopy(
173  numBCPts,
174  m_fields[i]->GetBndCondExpansions()[n]->GetPhys(), 1,
175  tmpStorage[i], 1);
176  }
177  GetPressure(tmpStorage, m_pressureStorage);
178  }
179  }
180 
181  // Loop over Boundary Regions for PressureInflowFileBC
182  m_fieldStorage = Array<OneD, Array<OneD, NekDouble> > (nvariables);
183  for (int n = 0; n < m_fields[0]->GetBndConditions().num_elements(); ++n)
184  {
185  // PressureInflowFile Boundary Condition
186  if (m_fields[0]->GetBndConditions()[n]->GetUserDefined() ==
187  "PressureInflowFile")
188  {
189  int numBCPts = m_fields[0]->
190  GetBndCondExpansions()[n]->GetNpoints();
191  for (int i = 0; i < nvariables; ++i)
192  {
193  m_fieldStorage[i] = Array<OneD, NekDouble>(numBCPts, 0.0);
194  Vmath::Vcopy(
195  numBCPts,
196  m_fields[i]->GetBndCondExpansions()[n]->GetPhys(), 1,
197  m_fieldStorage[i], 1);
198  }
199  }
200  }
201 
202  // Type of advection class to be used
203  switch(m_projectionType)
204  {
205  // Continuous field
206  case MultiRegions::eGalerkin:
207  {
208  ASSERTL0(false, "Continuous field not supported.");
209  break;
210  }
211  // Discontinuous field
212  case MultiRegions::eDiscontinuous:
213  {
214  string advName, diffName, riemName;
215 
216  // Setting up advection and diffusion operators
217  m_session->LoadSolverInfo("AdvectionType", advName, "WeakDG");
218  m_session->LoadSolverInfo("DiffusionType", diffName, "LDGNS");
219 
220  m_advection = SolverUtils::GetAdvectionFactory()
221  .CreateInstance(advName, advName);
222  m_diffusion = SolverUtils::GetDiffusionFactory()
223  .CreateInstance(diffName, diffName);
224 
225  if (m_specHP_dealiasing)
226  {
227  m_advection->SetFluxVector(&CompressibleFlowSystem::
228  GetFluxVectorDeAlias, this);
229  m_diffusion->SetFluxVectorNS(
230  &CompressibleFlowSystem::GetViscousFluxVectorDeAlias,
231  this);
232  }
233  else
234  {
235  m_advection->SetFluxVector (&CompressibleFlowSystem::
236  GetFluxVector, this);
237  m_diffusion->SetFluxVectorNS(&CompressibleFlowSystem::
238  GetViscousFluxVector, this);
239  }
240 
241  if (m_shockCaptureType=="Smooth" && m_EqTypeStr=="EulerADCFE")
242  {
243  m_advection->SetFluxVector(&CompressibleFlowSystem::
244  GetFluxVector, this);
245 
246  m_diffusion->SetArtificialDiffusionVector(
247  &CompressibleFlowSystem::GetSmoothArtificialViscosity, this);
248  }
249  if (m_shockCaptureType=="NonSmooth" && m_EqTypeStr=="EulerADCFE")
250  {
251  m_advection->SetFluxVector(&CompressibleFlowSystem::
252  GetFluxVector, this);
253 
254  m_diffusion->SetArtificialDiffusionVector(
255  &CompressibleFlowSystem::GetArtificialDynamicViscosity, this);
256  }
257 
258  // Setting up Riemann solver for advection operator
259  m_session->LoadSolverInfo("UpwindType", riemName, "Average");
260 
261  m_riemannSolver = SolverUtils::GetRiemannSolverFactory()
262  .CreateInstance(riemName);
263 
264  // Setting up upwind solver for diffusion operator
265  m_riemannSolverLDG = SolverUtils::GetRiemannSolverFactory()
266  .CreateInstance("UpwindLDG");
267 
268  // Setting up parameters for advection operator Riemann solver
269  m_riemannSolver->SetParam (
270  "gamma", &CompressibleFlowSystem::GetGamma, this);
271  m_riemannSolver->SetAuxVec(
272  "vecLocs", &CompressibleFlowSystem::GetVecLocs, this);
273  m_riemannSolver->SetVector(
274  "N", &CompressibleFlowSystem::GetNormals, this);
275 
276  // Setting up parameters for diffusion operator Riemann solver
277  m_riemannSolverLDG->SetParam (
278  "gamma", &CompressibleFlowSystem::GetGamma, this);
279  m_riemannSolverLDG->SetVector(
280  "vecLocs", &CompressibleFlowSystem::GetVecLocs, this);
281  m_riemannSolverLDG->SetVector(
282  "N", &CompressibleFlowSystem::GetNormals, this);
283 
284  // Concluding initialisation of advection / diffusion operators
285  m_advection->SetRiemannSolver (m_riemannSolver);
286  m_diffusion->SetRiemannSolver (m_riemannSolverLDG);
287  m_advection->InitObject (m_session, m_fields);
288  m_diffusion->InitObject (m_session, m_fields);
289  break;
290  }
291  default:
292  {
293  ASSERTL0(false, "Unsupported projection type.");
294  break;
295  }
296  }
297  }
298 
299  /**
300  * @brief Destructor for CompressibleFlowSystem class.
301  */
302  CompressibleFlowSystem::~CompressibleFlowSystem()
303  {
304 
305  }
306 
307  /**
308  * @brief Print out a summary with some relevant information.
309  */
310  void CompressibleFlowSystem::v_GenerateSummary(SolverUtils::SummaryList& s)
311  {
312  UnsteadySystem::v_GenerateSummary(s);
313  }
314 
315  /**
316  * @brief Set boundary conditions which can be:
317  * a) Wall and Symmerty BCs implemented at CompressibleFlowSystem level
318  * since they are compressible solver specific;
319  * b) Time dependent BCs.
320  *
321  * @param inarray: fields array.
322  * @param time: time.
323  */
324  void CompressibleFlowSystem::SetCommonBC(
325  const std::string &userDefStr,
326  const int n,
327  const NekDouble time,
328  int &cnt,
329  Array<OneD, Array<OneD, NekDouble> > &Fwd,
330  Array<OneD, Array<OneD, NekDouble> > &inarray)
331  {
332  std::string varName;
333  int nvariables = m_fields.num_elements();
334 
335  if(!userDefStr.empty())
336  {
337  if(boost::iequals(userDefStr,"Wall"))
338  {
339  WallBC(n, cnt, Fwd, inarray);
340  }
341  else if(boost::iequals(userDefStr,"WallViscous") ||
342  boost::iequals(userDefStr,"WallAdiabatic"))
343  {
344  // Wall Boundary Condition
345  WallViscousBC(n, cnt, Fwd, inarray);
346  }
347  else if(boost::iequals(userDefStr,"Symmetry"))
348  {
349  // Symmetric Boundary Condition
350  SymmetryBC(n, cnt, Fwd, inarray);
351  }
352  else if(boost::iequals(userDefStr,"RiemannInvariant"))
353  {
354  // Riemann invariant characteristic Boundary Condition
355  RiemannInvariantBC(n, cnt, Fwd, inarray);
356  }
357  else if(boost::iequals(userDefStr,"PressureOutflowNonReflective"))
358  {
359  // Pressure outflow non-reflective Boundary Condition
360  PressureOutflowNonReflectiveBC(n, cnt, Fwd, inarray);
361  }
362  else if(boost::iequals(userDefStr,"PressureOutflow"))
363  {
364  // Pressure outflow Boundary Condition
365  PressureOutflowBC(n, cnt, Fwd, inarray);
366  }
367  else if(boost::iequals(userDefStr,"PressureOutflowFile"))
368  {
369  // Pressure outflow Boundary Condition from file
370  PressureOutflowFileBC(n, cnt, Fwd, inarray);
371  }
372  else if(boost::iequals(userDefStr,"PressureInflowFile"))
373  {
374  // Pressure inflow Boundary Condition from file
375  PressureInflowFileBC(n, cnt, Fwd, inarray);
376  }
377  else if(boost::iequals(userDefStr,"ExtrapOrder0"))
378  {
379  // Extrapolation of the data at the boundaries
380  ExtrapOrder0BC(n, cnt, Fwd, inarray);
381  }
382  else if(boost::iequals(userDefStr,"TimeDependent"))
383  {
384  for (int i = 0; i < nvariables; ++i)
385  {
386  varName = m_session->GetVariable(i);
387  m_fields[i]->EvaluateBoundaryConditions(time, varName);
388  }
389  }
390  else
391  {
392  string errmsg = "Unrecognised boundary condition: ";
393  errmsg += userDefStr;
394  ASSERTL0(false,errmsg.c_str());
395  }
396  }
397  }
398 
399  /**
400  * @brief Wall boundary conditions for compressible flow problems.
401  */
402  void CompressibleFlowSystem::WallBC(
403  int bcRegion,
404  int cnt,
405  Array<OneD, Array<OneD, NekDouble> > &Fwd,
406  Array<OneD, Array<OneD, NekDouble> > &physarray)
407  {
408  int i;
409  int nVariables = physarray.num_elements();
410 
411  const Array<OneD, const int> &traceBndMap
412  = m_fields[0]->GetTraceBndMap();
413 
414  // Adjust the physical values of the trace to take
415  // user defined boundaries into account
416  int e, id1, id2, nBCEdgePts, eMax;
417 
418  eMax = m_fields[0]->GetBndCondExpansions()[bcRegion]->GetExpSize();
419 
420  for (e = 0; e < eMax; ++e)
421  {
422  nBCEdgePts = m_fields[0]->GetBndCondExpansions()[bcRegion]->
423  GetExp(e)->GetTotPoints();
424  id1 = m_fields[0]->GetBndCondExpansions()[bcRegion]->
425  GetPhys_Offset(e);
426  id2 = m_fields[0]->GetTrace()->GetPhys_Offset(traceBndMap[cnt+e]);
427 
428  // Boundary condition for epsilon term.
429  if (nVariables == m_spacedim+3)
430  {
431  NekDouble factor = 1.0;
432  NekDouble factor2 = 1.0;
433 
434  Array<OneD, NekDouble > tmp2(nBCEdgePts, 0.0);
435  Vmath::Smul(nBCEdgePts,
436  factor,
437  &Fwd[nVariables-1][id2], 1,
438  &tmp2[0], 1);
439 
440  Vmath::Vsub(nBCEdgePts,
441  &Fwd[nVariables-1][id2], 1,
442  &tmp2[0], 1,
443  &Fwd[nVariables-1][id2], 1);
444 
445  Vmath::Smul(nBCEdgePts,
446  factor2,
447  &Fwd[nVariables-1][id2], 1,
448  &Fwd[nVariables-1][id2], 1);
449 
450  }
451 
452  // For 2D/3D, define: v* = v - 2(v.n)n
453  Array<OneD, NekDouble> tmp(nBCEdgePts, 0.0);
454 
455  // Calculate (v.n)
456  for (i = 0; i < m_spacedim; ++i)
457  {
458  Vmath::Vvtvp(nBCEdgePts,
459  &Fwd[1+i][id2], 1,
460  &m_traceNormals[i][id2], 1,
461  &tmp[0], 1,
462  &tmp[0], 1);
463  }
464 
465  // Calculate 2.0(v.n)
466  Vmath::Smul(nBCEdgePts, -2.0, &tmp[0], 1, &tmp[0], 1);
467 
468  // Calculate v* = v - 2.0(v.n)n
469  for (i = 0; i < m_spacedim; ++i)
470  {
471  Vmath::Vvtvp(nBCEdgePts,
472  &tmp[0], 1,
473  &m_traceNormals[i][id2], 1,
474  &Fwd[1+i][id2], 1,
475  &Fwd[1+i][id2], 1);
476  }
477 
478  // Copy boundary adjusted values into the boundary expansion
479  for (i = 0; i < nVariables; ++i)
480  {
481  Vmath::Vcopy(nBCEdgePts, &Fwd[i][id2], 1,
482  &(m_fields[i]->GetBndCondExpansions()[bcRegion]->
483  UpdatePhys())[id1], 1);
484  }
485  }
486  }
487 
488  /**
489  * @brief Wall boundary conditions for viscous compressible flow problems.
490  */
491  void CompressibleFlowSystem::WallViscousBC(
492  int bcRegion,
493  int cnt,
494  Array<OneD, Array<OneD, NekDouble> > &Fwd,
495  Array<OneD, Array<OneD, NekDouble> > &physarray)
496  {
497  int i;
498  int nVariables = physarray.num_elements();
499 
500  const Array<OneD, const int> &traceBndMap
501  = m_fields[0]->GetTraceBndMap();
502 
503  // Take into account that for PDE based shock capturing, eps = 0 at the
504  // wall. Adjust the physical values of the trace to take user defined
505  // boundaries into account
506  int e, id1, id2, nBCEdgePts, eMax;
507 
508  eMax = m_fields[0]->GetBndCondExpansions()[bcRegion]->GetExpSize();
509 
510  for (e = 0; e < eMax; ++e)
511  {
512  nBCEdgePts = m_fields[0]->GetBndCondExpansions()[bcRegion]->
513  GetExp(e)->GetTotPoints();
514  id1 = m_fields[0]->GetBndCondExpansions()[bcRegion]->
515  GetPhys_Offset(e);
516  id2 = m_fields[0]->GetTrace()->GetPhys_Offset(traceBndMap[cnt+e]);
517 
518  if (nVariables == m_spacedim+3)
519  {
520  NekDouble factor = 0.0;
521  NekDouble factor2 = 1.0;
522 
523  Array<OneD, NekDouble > tmp2(nBCEdgePts, 0.0);
524  Vmath::Smul(nBCEdgePts,
525  factor,
526  &Fwd[nVariables-1][id2], 1,
527  &tmp2[0], 1);
528 
529  Vmath::Vsub(nBCEdgePts,
530  &Fwd[nVariables-1][id2], 1,
531  &tmp2[0], 1,
532  &Fwd[nVariables-1][id2], 1);
533 
534  Vmath::Smul(nBCEdgePts,
535  factor2,
536  &Fwd[nVariables-1][id2], 1,
537  &Fwd[nVariables-1][id2], 1);
538  }
539 
540  for (i = 0; i < m_spacedim; i++)
541  {
542  Vmath::Neg(nBCEdgePts, &Fwd[i+1][id2], 1);
543  }
544 
545  // Copy boundary adjusted values into the boundary expansion
546  for (i = 0; i < nVariables; ++i)
547  {
548  Vmath::Vcopy(nBCEdgePts, &Fwd[i][id2], 1,
549  &(m_fields[i]->GetBndCondExpansions()[bcRegion]->
550  UpdatePhys())[id1], 1);
551  }
552  }
553  }
554 
555  /**
556  * @brief Symmetry boundary conditions for compressible flow problems.
557  */
558  void CompressibleFlowSystem::SymmetryBC(
559  int bcRegion,
560  int cnt,
561  Array<OneD, Array<OneD, NekDouble> > &Fwd,
562  Array<OneD, Array<OneD, NekDouble> > &physarray)
563  {
564  int i;
565  int nVariables = physarray.num_elements();
566 
567  const Array<OneD, const int> &traceBndMap
568  = m_fields[0]->GetTraceBndMap();
569 
570  // Take into account that for PDE based shock capturing, eps = 0 at the
571  // wall.
572  int e, id1, id2, nBCEdgePts, eMax;
573 
574  eMax = m_fields[0]->GetBndCondExpansions()[bcRegion]->GetExpSize();
575 
576  for (e = 0; e < eMax; ++e)
577  {
578  nBCEdgePts = m_fields[0]->GetBndCondExpansions()[bcRegion]->
579  GetExp(e)->GetTotPoints();
580  id1 = m_fields[0]->GetBndCondExpansions()[bcRegion]->
581  GetPhys_Offset(e);
582  id2 = m_fields[0]->GetTrace()->GetPhys_Offset(traceBndMap[cnt+e]);
583 
584  if (nVariables == m_spacedim+3)
585  {
586  NekDouble factor = 0.0;
587  NekDouble factor2 = 1.0;
588 
589  Array<OneD, NekDouble > tmp2(nBCEdgePts, 0.0);
590  Vmath::Smul(nBCEdgePts,
591  factor,
592  &Fwd[nVariables-1][id2], 1,
593  &tmp2[0], 1);
594 
595  Vmath::Vsub(nBCEdgePts,
596  &Fwd[nVariables-1][id2], 1,
597  &tmp2[0], 1,
598  &Fwd[nVariables-1][id2], 1);
599 
600  Vmath::Smul(nBCEdgePts,
601  factor2,
602  &Fwd[nVariables-1][id2], 1,
603  &Fwd[nVariables-1][id2], 1);
604  }
605 
606  // For 2D/3D, define: v* = v - 2(v.n)n
607  Array<OneD, NekDouble> tmp(nBCEdgePts, 0.0);
608 
609  // Calculate (v.n)
610  for (i = 0; i < m_spacedim; ++i)
611  {
612  Vmath::Vvtvp(nBCEdgePts,
613  &Fwd[1+i][id2], 1,
614  &m_traceNormals[i][id2], 1,
615  &tmp[0], 1,
616  &tmp[0], 1);
617  }
618 
619  // Calculate 2.0(v.n)
620  Vmath::Smul(nBCEdgePts, -2.0, &tmp[0], 1, &tmp[0], 1);
621 
622  // Calculate v* = v - 2.0(v.n)n
623  for (i = 0; i < m_spacedim; ++i)
624  {
625  Vmath::Vvtvp(nBCEdgePts,
626  &tmp[0], 1,
627  &m_traceNormals[i][id2], 1,
628  &Fwd[1+i][id2], 1,
629  &Fwd[1+i][id2], 1);
630  }
631 
632  // Copy boundary adjusted values into the boundary expansion
633  for (i = 0; i < nVariables; ++i)
634  {
635  Vmath::Vcopy(nBCEdgePts, &Fwd[i][id2], 1,
636  &(m_fields[i]->GetBndCondExpansions()[bcRegion]->
637  UpdatePhys())[id1], 1);
638  }
639  }
640  }
641 
642  /**
643  * @brief Outflow characteristic boundary conditions for compressible
644  * flow problems.
645  */
646  void CompressibleFlowSystem::RiemannInvariantBC(
647  int bcRegion,
648  int cnt,
649  Array<OneD, Array<OneD, NekDouble> > &Fwd,
650  Array<OneD, Array<OneD, NekDouble> > &physarray)
651  {
652  int i, j;
653  int nTracePts = GetTraceTotPoints();
654  int nDimensions = m_spacedim;
655 
656  const Array<OneD, const int> &traceBndMap
657  = m_fields[0]->GetTraceBndMap();
658 
659  NekDouble gamma = m_gamma;
660  NekDouble gammaInv = 1.0 / gamma;
661  NekDouble gammaMinusOne = gamma - 1.0;
662  NekDouble gammaMinusOneInv = 1.0 / gammaMinusOne;
663 
664  Array<OneD, NekDouble> tmp1 (nTracePts, 0.0);
665  Array<OneD, NekDouble> tmp2 (nTracePts, 0.0);
666  Array<OneD, NekDouble> VnInf(nTracePts, 0.0);
667  Array<OneD, NekDouble> velInf(nDimensions, 0.0);
668 
669  // Computing the normal velocity for characteristics coming
670  // from outside the computational domain
671  velInf[0] = m_uInf;
672  Vmath::Smul(nTracePts, m_uInf, m_traceNormals[0], 1, VnInf, 1);
673  if (nDimensions == 2 || nDimensions == 3)
674  {
675  velInf[1] = m_vInf;
676  Vmath::Smul(nTracePts, m_vInf, m_traceNormals[1], 1, tmp1, 1);
677  Vmath::Vadd(nTracePts, VnInf, 1, tmp1, 1, VnInf, 1);
678  }
679  if (nDimensions == 3)
680  {
681  velInf[2] = m_wInf;
682  Vmath::Smul(nTracePts, m_wInf, m_traceNormals[2], 1, tmp2, 1);
683  Vmath::Vadd(nTracePts, VnInf, 1, tmp2, 1, VnInf, 1);
684  }
685 
686  // Computing the normal velocity for characteristics coming
687  // from inside the computational domain
688  Array<OneD, NekDouble > Vn (nTracePts, 0.0);
689  Array<OneD, NekDouble > Vel(nTracePts, 0.0);
690  for (i = 0; i < nDimensions; ++i)
691  {
692  Vmath::Vdiv(nTracePts, Fwd[i+1], 1, Fwd[0], 1, Vel, 1);
693  Vmath::Vvtvp(nTracePts, m_traceNormals[i], 1, Vel, 1, Vn, 1, Vn, 1);
694  }
695 
696  // Computing the absolute value of the velocity in order to compute the
697  // Mach number to decide whether supersonic or subsonic
698  Array<OneD, NekDouble > absVel(nTracePts, 0.0);
699  for (i = 0; i < nDimensions; ++i)
700  {
701  Vmath::Vdiv(nTracePts, Fwd[i+1], 1, Fwd[0], 1, tmp1, 1);
702  Vmath::Vmul(nTracePts, tmp1, 1, tmp1, 1, tmp1, 1);
703  Vmath::Vadd(nTracePts, tmp1, 1, absVel, 1, absVel, 1);
704  }
705  Vmath::Vsqrt(nTracePts, absVel, 1, absVel, 1);
706 
707  // Get speed of sound
708  Array<OneD, NekDouble > soundSpeed(nTracePts);
709  Array<OneD, NekDouble > pressure (nTracePts);
710 
711  for (i = 0; i < nTracePts; i++)
712  {
713  if (m_spacedim == 1)
714  {
715  pressure[i] = (gammaMinusOne) * (Fwd[2][i] -
716  0.5 * (Fwd[1][i] * Fwd[1][i] / Fwd[0][i]));
717  }
718  else if (m_spacedim == 2)
719  {
720  pressure[i] = (gammaMinusOne) * (Fwd[3][i] -
721  0.5 * (Fwd[1][i] * Fwd[1][i] / Fwd[0][i] +
722  Fwd[2][i] * Fwd[2][i] / Fwd[0][i]));
723  }
724  else
725  {
726  pressure[i] = (gammaMinusOne) * (Fwd[4][i] -
727  0.5 * (Fwd[1][i] * Fwd[1][i] / Fwd[0][i] +
728  Fwd[2][i] * Fwd[2][i] / Fwd[0][i] +
729  Fwd[3][i] * Fwd[3][i] / Fwd[0][i]));
730  }
731 
732  soundSpeed[i] = sqrt(gamma * pressure[i] / Fwd[0][i]);
733  }
734 
735  // Get Mach
736  Array<OneD, NekDouble > Mach(nTracePts, 0.0);
737  Vmath::Vdiv(nTracePts, Vn, 1, soundSpeed, 1, Mach, 1);
738  Vmath::Vabs(nTracePts, Mach, 1, Mach, 1);
739 
740  // Auxiliary variables
741  int eMax;
742  int e, id1, id2, nBCEdgePts, pnt;
743  NekDouble cPlus, rPlus, cMinus, rMinus, VDBC, VNBC;
744  Array<OneD, NekDouble> velBC(nDimensions, 0.0);
745  Array<OneD, NekDouble> rhoVelBC(nDimensions, 0.0);
746  NekDouble rhoBC, EBC, cBC, sBC, pBC;
747 
748  eMax = m_fields[0]->GetBndCondExpansions()[bcRegion]->GetExpSize();
749 
750  // Loop on bcRegions
751  for (e = 0; e < eMax; ++e)
752  {
753  nBCEdgePts = m_fields[0]->GetBndCondExpansions()[bcRegion]->
754  GetExp(e)->GetTotPoints();
755 
756  id1 = m_fields[0]->GetBndCondExpansions()[bcRegion]->
757  GetPhys_Offset(e);
758  id2 = m_fields[0]->GetTrace()->GetPhys_Offset(traceBndMap[cnt+e]);
759 
760  // Loop on the points of the bcRegion
761  for (i = 0; i < nBCEdgePts; i++)
762  {
763  pnt = id2+i;
764 
765  // Impose inflow Riemann invariant
766  if (Vn[pnt] <= 0.0)
767  {
768  // Subsonic flows
769  if (Mach[pnt] < 1.00)
770  {
771  // + Characteristic from inside
772  cPlus = sqrt(gamma * pressure[pnt] / Fwd[0][pnt]);
773  rPlus = Vn[pnt] + 2.0 * cPlus * gammaMinusOneInv;
774 
775  // - Characteristic from boundary
776  cMinus = sqrt(gamma * m_pInf / m_rhoInf);
777  rMinus = VnInf[pnt] - 2.0 * cMinus * gammaMinusOneInv;
778  }
779  else
780  {
781  // + Characteristic from inside
782  cPlus = sqrt(gamma * m_pInf / m_rhoInf);
783  rPlus = VnInf[pnt] + 2.0 * cPlus * gammaMinusOneInv;
784 
785  // + Characteristic from inside
786  cMinus = sqrt(gamma * m_pInf / m_rhoInf);
787  rMinus = VnInf[pnt] - 2.0 * cPlus * gammaMinusOneInv;
788  }
789 
790  // Riemann boundary variables
791  VNBC = 0.5 * (rPlus + rMinus);
792  cBC = 0.25 * gammaMinusOne * (rPlus - rMinus);
793  VDBC = VNBC - VnInf[pnt];
794 
795  // Thermodynamic boundary variables
796  sBC = m_pInf / (pow(m_rhoInf, gamma));
797  rhoBC = pow((cBC * cBC) / (gamma * sBC), gammaMinusOneInv);
798  pBC = rhoBC * cBC * cBC * gammaInv;
799 
800  // Kinetic energy initialiasation
801  NekDouble EkBC = 0.0;
802 
803  // Boundary velocities
804  for ( j = 0; j < nDimensions; ++j)
805  {
806  velBC[j] = velInf[j] + VDBC * m_traceNormals[j][pnt];
807  rhoVelBC[j] = rhoBC * velBC[j];
808  EkBC += 0.5 * rhoBC * velBC[j]*velBC[j];
809  }
810 
811  // Boundary energy
812  EBC = pBC * gammaMinusOneInv + EkBC;
813 
814  // Imposing Riemann Invariant boundary conditions
815  (m_fields[0]->GetBndCondExpansions()[bcRegion]->
816  UpdatePhys())[id1+i] = rhoBC;
817  for (j = 0; j < nDimensions; ++j)
818  {
819  (m_fields[j+1]->GetBndCondExpansions()[bcRegion]->
820  UpdatePhys())[id1+i] = rhoVelBC[j];
821  }
822  (m_fields[nDimensions+1]->GetBndCondExpansions()[bcRegion]->
823  UpdatePhys())[id1+i] = EBC;
824 
825  }
826  else // Impose outflow Riemann invariant
827  {
828  // Subsonic flows
829  if (Mach[pnt] < 1.00)
830  {
831  // + Characteristic from inside
832  cPlus = sqrt(gamma * pressure[pnt] / Fwd[0][pnt]);
833  rPlus = Vn[pnt] + 2.0 * cPlus * gammaMinusOneInv;
834 
835  // - Characteristic from boundary
836  cMinus = sqrt(gamma * m_pInf / m_rhoInf);
837  rMinus = VnInf[pnt] - 2.0 * cMinus * gammaMinusOneInv;
838  }
839  else
840  {
841  // + Characteristic from inside
842  cPlus = sqrt(gamma * pressure[pnt] / Fwd[0][pnt]);
843  rPlus = Vn[pnt] + 2.0 * cPlus * gammaMinusOneInv;
844 
845  // + Characteristic from inside
846  cMinus = sqrt(gamma * pressure[pnt] / Fwd[0][pnt]);
847  rMinus = Vn[pnt] - 2.0 * cPlus * gammaMinusOneInv;
848  }
849 
850  // Riemann boundary variables
851  VNBC = 0.5 * (rPlus + rMinus);
852  cBC = 0.25 * gammaMinusOne * (rPlus - rMinus);
853  VDBC = VNBC - Vn[pnt];
854 
855  // Thermodynamic boundary variables
856  sBC = pressure[pnt] / (pow(Fwd[0][pnt], gamma));
857  rhoBC = pow((cBC * cBC) / (gamma * sBC), gammaMinusOneInv);
858  pBC = rhoBC * cBC * cBC * gammaInv;
859 
860  // Kinetic energy initialiasation
861  NekDouble EkBC = 0.0;
862 
863  // Boundary velocities
864  for ( j = 0; j < nDimensions; ++j)
865  {
866  velBC[j] = Fwd[j+1][pnt] / Fwd[0][pnt] +
867  VDBC * m_traceNormals[j][pnt];
868  rhoVelBC[j] = rhoBC * velBC[j];
869  EkBC += 0.5 * rhoBC * velBC[j]*velBC[j];
870  }
871 
872  // Boundary energy
873  EBC = pBC * gammaMinusOneInv + EkBC;
874 
875  // Imposing Riemann Invariant boundary conditions
876  (m_fields[0]->GetBndCondExpansions()[bcRegion]->
877  UpdatePhys())[id1+i] = rhoBC;
878  for (j = 0; j < nDimensions; ++j)
879  {
880  (m_fields[j+1]->GetBndCondExpansions()[bcRegion]->
881  UpdatePhys())[id1+i] = rhoVelBC[j];
882  }
883  (m_fields[nDimensions+1]->GetBndCondExpansions()[bcRegion]->
884  UpdatePhys())[id1+i] = EBC;
885  }
886  }
887  }
888  }
889 
890  /**
891  * @brief Pressure outflow non-reflective
892  * boundary conditions for compressible flow
893  * problems.
894  */
895  void CompressibleFlowSystem::PressureOutflowNonReflectiveBC(
896  int bcRegion,
897  int cnt,
898  Array<OneD, Array<OneD, NekDouble> > &Fwd,
899  Array<OneD, Array<OneD, NekDouble> > &physarray)
900  {
901  int i, j;
902  int nTracePts = GetTraceTotPoints();
903  int nVariables = physarray.num_elements();
904  int nDimensions = m_spacedim;
905 
906  const Array<OneD, const int> &traceBndMap
907  = m_fields[0]->GetTraceBndMap();
908 
909  NekDouble gamma = m_gamma;
910  NekDouble gammaMinusOne = gamma - 1.0;
911  NekDouble gammaMinusOneInv = 1.0 / gammaMinusOne;
912 
913  Array<OneD, NekDouble> tmp1 (nTracePts, 0.0);
914  Array<OneD, NekDouble> tmp2 (nTracePts, 0.0);
915  Array<OneD, NekDouble> VnInf(nTracePts, 0.0);
916  Array<OneD, NekDouble> velInf(nDimensions, 0.0);
917 
918  // Computing the normal velocity for characteristics coming
919  // from outside the computational domain
920  velInf[0] = m_uInf;
921  Vmath::Smul(nTracePts, m_uInf, m_traceNormals[0], 1, VnInf, 1);
922  if (nDimensions == 2 || nDimensions == 3)
923  {
924  velInf[1] = m_vInf;
925  Vmath::Smul(nTracePts, m_vInf, m_traceNormals[1], 1, tmp1, 1);
926  Vmath::Vadd(nTracePts, VnInf, 1, tmp1, 1, VnInf, 1);
927  }
928  if (nDimensions == 3)
929  {
930  velInf[2] = m_wInf;
931  Vmath::Smul(nTracePts, m_wInf, m_traceNormals[2], 1, tmp2, 1);
932  Vmath::Vadd(nTracePts, VnInf, 1, tmp2, 1, VnInf, 1);
933  }
934 
935  // Computing the normal velocity for characteristics coming
936  // from inside the computational domain
937  Array<OneD, NekDouble > Vn (nTracePts, 0.0);
938  Array<OneD, NekDouble > Vel(nTracePts, 0.0);
939  for (i = 0; i < nDimensions; ++i)
940  {
941  Vmath::Vdiv(nTracePts, Fwd[i+1], 1, Fwd[0], 1, Vel, 1);
942  Vmath::Vvtvp(nTracePts, m_traceNormals[i], 1, Vel, 1, Vn, 1, Vn, 1);
943  }
944 
945  // Computing the absolute value of the velocity in order to compute the
946  // Mach number to decide whether supersonic or subsonic
947  Array<OneD, NekDouble > absVel(nTracePts, 0.0);
948  for (i = 0; i < nDimensions; ++i)
949  {
950  Vmath::Vdiv(nTracePts, Fwd[i+1], 1, Fwd[0], 1, tmp1, 1);
951  Vmath::Vmul(nTracePts, tmp1, 1, tmp1, 1, tmp1, 1);
952  Vmath::Vadd(nTracePts, tmp1, 1, absVel, 1, absVel, 1);
953  }
954  Vmath::Vsqrt(nTracePts, absVel, 1, absVel, 1);
955 
956  // Get speed of sound
957  Array<OneD, NekDouble > soundSpeed(nTracePts);
958  Array<OneD, NekDouble > pressure (nTracePts);
959 
960  for (i = 0; i < nTracePts; i++)
961  {
962  if (m_spacedim == 1)
963  {
964  pressure[i] = (gammaMinusOne) * (Fwd[2][i] -
965  0.5 * (Fwd[1][i] * Fwd[1][i] / Fwd[0][i]));
966  }
967  else if (m_spacedim == 2)
968  {
969  pressure[i] = (gammaMinusOne) * (Fwd[3][i] -
970  0.5 * (Fwd[1][i] * Fwd[1][i] / Fwd[0][i] +
971  Fwd[2][i] * Fwd[2][i] / Fwd[0][i]));
972  }
973  else
974  {
975  pressure[i] = (gammaMinusOne) * (Fwd[4][i] -
976  0.5 * (Fwd[1][i] * Fwd[1][i] / Fwd[0][i] +
977  Fwd[2][i] * Fwd[2][i] / Fwd[0][i] +
978  Fwd[3][i] * Fwd[3][i] / Fwd[0][i]));
979  }
980 
981  soundSpeed[i] = sqrt(gamma * pressure[i] / Fwd[0][i]);
982  }
983 
984  // Get Mach
985  Array<OneD, NekDouble > Mach(nTracePts, 0.0);
986  Vmath::Vdiv(nTracePts, Vn, 1, soundSpeed, 1, Mach, 1);
987  Vmath::Vabs(nTracePts, Mach, 1, Mach, 1);
988 
989  // Auxiliary variables
990  int e, id1, id2, npts, pnt;
991  NekDouble rhoeb;
992 
993  // Loop on the bcRegions
994  for (e = 0; e < m_fields[0]->GetBndCondExpansions()[bcRegion]->
995  GetExpSize(); ++e)
996  {
997  npts = m_fields[0]->GetBndCondExpansions()[bcRegion]->
998  GetExp(e)->GetTotPoints();
999  id1 = m_fields[0]->GetBndCondExpansions()[bcRegion]->
1000  GetPhys_Offset(e) ;
1001  id2 = m_fields[0]->GetTrace()->GetPhys_Offset(traceBndMap[cnt+e]);
1002 
1003  // Loop on points of bcRegion 'e'
1004  for (i = 0; i < npts; i++)
1005  {
1006  pnt = id2+i;
1007 
1008  // Subsonic flows
1009  if (Mach[pnt] < 0.99)
1010  {
1011  // Kinetic energy calculation
1012  NekDouble Ek = 0.0;
1013  for (j = 1; j < nVariables-1; ++j)
1014  {
1015  Ek += 0.5 * (Fwd[j][pnt] * Fwd[j][pnt]) / Fwd[0][pnt];
1016  }
1017 
1018  rhoeb = m_pInf * gammaMinusOneInv + Ek;
1019 
1020  // Partial extrapolation for subsonic cases
1021  for (j = 0; j < nVariables-1; ++j)
1022  {
1023  (m_fields[j]->GetBndCondExpansions()[bcRegion]->
1024  UpdatePhys())[id1+i] = Fwd[j][pnt];
1025  }
1026 
1027  (m_fields[nVariables-1]->GetBndCondExpansions()[bcRegion]->
1028  UpdatePhys())[id1+i] = 2.0 * rhoeb - Fwd[nVariables-1][pnt];
1029  }
1030  // Supersonic flows
1031  else
1032  {
1033  for (j = 0; j < nVariables; ++j)
1034  {
1035  // Extrapolation for supersonic cases
1036  (m_fields[j]->GetBndCondExpansions()[bcRegion]->
1037  UpdatePhys())[id1+i] = Fwd[j][pnt];
1038  }
1039  }
1040  }
1041  }
1042  }
1043 
1044  /**
1045  * @brief Pressure outflow boundary conditions for compressible flow
1046  * problems.
1047  */
1048  void CompressibleFlowSystem::PressureOutflowBC(
1049  int bcRegion,
1050  int cnt,
1051  Array<OneD, Array<OneD, NekDouble> > &Fwd,
1052  Array<OneD, Array<OneD, NekDouble> > &physarray)
1053  {
1054  int i, j;
1055  int nTracePts = GetTraceTotPoints();
1056  int nVariables = physarray.num_elements();
1057  int nDimensions = m_spacedim;
1058 
1059  const Array<OneD, const int> &traceBndMap
1060  = m_fields[0]->GetTraceBndMap();
1061 
1062  NekDouble gamma = m_gamma;
1063  NekDouble gammaMinusOne = gamma - 1.0;
1064  NekDouble gammaMinusOneInv = 1.0 / gammaMinusOne;
1065 
1066  Array<OneD, NekDouble> tmp1 (nTracePts, 0.0);
1067  Array<OneD, NekDouble> tmp2 (nTracePts, 0.0);
1068  Array<OneD, NekDouble> VnInf(nTracePts, 0.0);
1069  Array<OneD, NekDouble> velInf(nDimensions, 0.0);
1070 
1071  // Computing the normal velocity for characteristics coming
1072  // from outside the computational domain
1073  velInf[0] = m_uInf;
1074  Vmath::Smul(nTracePts, m_uInf, m_traceNormals[0], 1, VnInf, 1);
1075  if (nDimensions == 2 || nDimensions == 3)
1076  {
1077  velInf[1] = m_vInf;
1078  Vmath::Smul(nTracePts, m_vInf, m_traceNormals[1], 1, tmp1, 1);
1079  Vmath::Vadd(nTracePts, VnInf, 1, tmp1, 1, VnInf, 1);
1080  }
1081  if (nDimensions == 3)
1082  {
1083  velInf[2] = m_wInf;
1084  Vmath::Smul(nTracePts, m_wInf, m_traceNormals[2], 1, tmp2, 1);
1085  Vmath::Vadd(nTracePts, VnInf, 1, tmp2, 1, VnInf, 1);
1086  }
1087 
1088  // Computing the normal velocity for characteristics coming
1089  // from inside the computational domain
1090  Array<OneD, NekDouble > Vn (nTracePts, 0.0);
1091  Array<OneD, NekDouble > Vel(nTracePts, 0.0);
1092  for (i = 0; i < nDimensions; ++i)
1093  {
1094  Vmath::Vdiv(nTracePts, Fwd[i+1], 1, Fwd[0], 1, Vel, 1);
1095  Vmath::Vvtvp(nTracePts, m_traceNormals[i], 1, Vel, 1, Vn, 1, Vn, 1);
1096  }
1097 
1098  // Computing the absolute value of the velocity in order to compute the
1099  // Mach number to decide whether supersonic or subsonic
1100  Array<OneD, NekDouble > absVel(nTracePts, 0.0);
1101  for (i = 0; i < nDimensions; ++i)
1102  {
1103  Vmath::Vdiv(nTracePts, Fwd[i+1], 1, Fwd[0], 1, tmp1, 1);
1104  Vmath::Vmul(nTracePts, tmp1, 1, tmp1, 1, tmp1, 1);
1105  Vmath::Vadd(nTracePts, tmp1, 1, absVel, 1, absVel, 1);
1106  }
1107  Vmath::Vsqrt(nTracePts, absVel, 1, absVel, 1);
1108 
1109  // Get speed of sound
1110  Array<OneD, NekDouble > soundSpeed(nTracePts);
1111  Array<OneD, NekDouble > pressure (nTracePts);
1112 
1113  for (i = 0; i < nTracePts; i++)
1114  {
1115  if (m_spacedim == 1)
1116  {
1117  pressure[i] = (gammaMinusOne) * (Fwd[2][i] -
1118  0.5 * (Fwd[1][i] * Fwd[1][i] / Fwd[0][i]));
1119  }
1120  else if (m_spacedim == 2)
1121  {
1122  pressure[i] = (gammaMinusOne) * (Fwd[3][i] -
1123  0.5 * (Fwd[1][i] * Fwd[1][i] / Fwd[0][i] +
1124  Fwd[2][i] * Fwd[2][i] / Fwd[0][i]));
1125  }
1126  else
1127  {
1128  pressure[i] = (gammaMinusOne) * (Fwd[4][i] -
1129  0.5 * (Fwd[1][i] * Fwd[1][i] / Fwd[0][i] +
1130  Fwd[2][i] * Fwd[2][i] / Fwd[0][i] +
1131  Fwd[3][i] * Fwd[3][i] / Fwd[0][i]));
1132  }
1133 
1134  soundSpeed[i] = sqrt(gamma * pressure[i] / Fwd[0][i]);
1135  }
1136 
1137  // Get Mach
1138  Array<OneD, NekDouble > Mach(nTracePts, 0.0);
1139  Vmath::Vdiv(nTracePts, Vn, 1, soundSpeed, 1, Mach, 1);
1140  Vmath::Vabs(nTracePts, Mach, 1, Mach, 1);
1141 
1142  // Auxiliary variables
1143  int e, id1, id2, npts, pnt;
1144  NekDouble rhoeb;
1145 
1146  // Loop on the bcRegions
1147  for (e = 0; e < m_fields[0]->GetBndCondExpansions()[bcRegion]->
1148  GetExpSize(); ++e)
1149  {
1150  npts = m_fields[0]->GetBndCondExpansions()[bcRegion]->
1151  GetExp(e)->GetTotPoints();
1152  id1 = m_fields[0]->GetBndCondExpansions()[bcRegion]->
1153  GetPhys_Offset(e) ;
1154  id2 = m_fields[0]->GetTrace()->GetPhys_Offset(traceBndMap[cnt+e]);
1155 
1156  // Loop on points of bcRegion 'e'
1157  for (i = 0; i < npts; i++)
1158  {
1159  pnt = id2+i;
1160 
1161  // Subsonic flows
1162  if (Mach[pnt] < 0.99)
1163  {
1164  // Kinetic energy calculation
1165  NekDouble Ek = 0.0;
1166  for (j = 1; j < nVariables-1; ++j)
1167  {
1168  Ek += 0.5 * (Fwd[j][pnt] * Fwd[j][pnt]) / Fwd[0][pnt];
1169  }
1170 
1171  rhoeb = m_pInf * gammaMinusOneInv + Ek;
1172 
1173  // Partial extrapolation for subsonic cases
1174  for (j = 0; j < nVariables-1; ++j)
1175  {
1176  (m_fields[j]->GetBndCondExpansions()[bcRegion]->
1177  UpdatePhys())[id1+i] = Fwd[j][pnt];
1178  }
1179 
1180  (m_fields[nVariables-1]->GetBndCondExpansions()[bcRegion]->
1181  UpdatePhys())[id1+i] = rhoeb;
1182  }
1183  // Supersonic flows
1184  else
1185  {
1186  for (j = 0; j < nVariables; ++j)
1187  {
1188  // Extrapolation for supersonic cases
1189  (m_fields[j]->GetBndCondExpansions()[bcRegion]->
1190  UpdatePhys())[id1+i] = Fwd[j][pnt];
1191  }
1192  }
1193  }
1194  }
1195  }
1196 
1197 
1198  /**
1199  * @brief Pressure outflow boundary conditions for compressible flow
1200  * problems.
1201  */
1202  void CompressibleFlowSystem::PressureOutflowFileBC(
1203  int bcRegion,
1204  int cnt,
1205  Array<OneD, Array<OneD, NekDouble> > &Fwd,
1206  Array<OneD, Array<OneD, NekDouble> > &physarray)
1207  {
1208  int i, j;
1209  int nTracePts = GetTraceTotPoints();
1210  int nVariables = physarray.num_elements();
1211  int nDimensions = m_spacedim;
1212 
1213  const Array<OneD, const int> &traceBndMap
1214  = m_fields[0]->GetTraceBndMap();
1215 
1216  NekDouble gamma = m_gamma;
1217  NekDouble gammaMinusOne = gamma - 1.0;
1218  NekDouble gammaMinusOneInv = 1.0 / gammaMinusOne;
1219 
1220  Array<OneD, NekDouble> tmp1 (nTracePts, 0.0);
1221  Array<OneD, NekDouble> tmp2 (nTracePts, 0.0);
1222  Array<OneD, NekDouble> VnInf(nTracePts, 0.0);
1223  Array<OneD, NekDouble> velInf(nDimensions, 0.0);
1224 
1225  // Computing the normal velocity for characteristics coming
1226  // from outside the computational domain
1227  velInf[0] = m_uInf;
1228  Vmath::Smul(nTracePts, m_uInf, m_traceNormals[0], 1, VnInf, 1);
1229  if (nDimensions == 2 || nDimensions == 3)
1230  {
1231  velInf[1] = m_vInf;
1232  Vmath::Smul(nTracePts, m_vInf, m_traceNormals[1], 1, tmp1, 1);
1233  Vmath::Vadd(nTracePts, VnInf, 1, tmp1, 1, VnInf, 1);
1234  }
1235  if (nDimensions == 3)
1236  {
1237  velInf[2] = m_wInf;
1238  Vmath::Smul(nTracePts, m_wInf, m_traceNormals[2], 1, tmp2, 1);
1239  Vmath::Vadd(nTracePts, VnInf, 1, tmp2, 1, VnInf, 1);
1240  }
1241 
1242  // Computing the normal velocity for characteristics coming
1243  // from inside the computational domain
1244  Array<OneD, NekDouble > Vn (nTracePts, 0.0);
1245  Array<OneD, NekDouble > Vel(nTracePts, 0.0);
1246  for (i = 0; i < nDimensions; ++i)
1247  {
1248  Vmath::Vdiv(nTracePts, Fwd[i+1], 1, Fwd[0], 1, Vel, 1);
1249  Vmath::Vvtvp(nTracePts, m_traceNormals[i], 1, Vel, 1, Vn, 1, Vn, 1);
1250  }
1251 
1252  // Computing the absolute value of the velocity in order to compute the
1253  // Mach number to decide whether supersonic or subsonic
1254  Array<OneD, NekDouble > absVel(nTracePts, 0.0);
1255  for (i = 0; i < nDimensions; ++i)
1256  {
1257  Vmath::Vdiv(nTracePts, Fwd[i+1], 1, Fwd[0], 1, tmp1, 1);
1258  Vmath::Vmul(nTracePts, tmp1, 1, tmp1, 1, tmp1, 1);
1259  Vmath::Vadd(nTracePts, tmp1, 1, absVel, 1, absVel, 1);
1260  }
1261  Vmath::Vsqrt(nTracePts, absVel, 1, absVel, 1);
1262 
1263  // Get speed of sound
1264  Array<OneD, NekDouble > soundSpeed (nTracePts, 0.0);
1265  Array<OneD, NekDouble > pressure (nTracePts, 0.0);
1266 
1267  for (i = 0; i < nTracePts; i++)
1268  {
1269  if (m_spacedim == 1)
1270  {
1271  pressure[i] = (gammaMinusOne) * (Fwd[2][i] -
1272  0.5 * (Fwd[1][i] * Fwd[1][i] / Fwd[0][i]));
1273  }
1274  else if (m_spacedim == 2)
1275  {
1276  pressure[i] = (gammaMinusOne) * (Fwd[3][i] -
1277  0.5 * (Fwd[1][i] * Fwd[1][i] / Fwd[0][i] +
1278  Fwd[2][i] * Fwd[2][i] / Fwd[0][i]));
1279  }
1280  else
1281  {
1282  pressure[i] = (gammaMinusOne) * (Fwd[4][i] -
1283  0.5 * (Fwd[1][i] * Fwd[1][i] / Fwd[0][i] +
1284  Fwd[2][i] * Fwd[2][i] / Fwd[0][i] +
1285  Fwd[3][i] * Fwd[3][i] / Fwd[0][i]));
1286  }
1287 
1288  soundSpeed[i] = sqrt(gamma * pressure[i] / Fwd[0][i]);
1289  }
1290 
1291  // Get Mach
1292  Array<OneD, NekDouble > Mach(nTracePts, 0.0);
1293  Vmath::Vdiv(nTracePts, Vn, 1, soundSpeed, 1, Mach, 1);
1294  Vmath::Vabs(nTracePts, Mach, 1, Mach, 1);
1295 
1296  // Auxiliary variables
1297  int e, id1, id2, npts, pnt;
1298  NekDouble rhoeb;
1299 
1300  // Loop on the bcRegions
1301  for (e = 0; e < m_fields[0]->GetBndCondExpansions()[bcRegion]->
1302  GetExpSize(); ++e)
1303  {
1304  npts = m_fields[0]->GetBndCondExpansions()[bcRegion]->
1305  GetExp(e)->GetTotPoints();
1306  id1 = m_fields[0]->GetBndCondExpansions()[bcRegion]->
1307  GetPhys_Offset(e);
1308  id2 = m_fields[0]->GetTrace()->GetPhys_Offset(traceBndMap[cnt+e]);
1309 
1310  // Loop on points of bcRegion 'e'
1311  for (i = 0; i < npts; i++)
1312  {
1313  pnt = id2+i;
1314 
1315  // Subsonic flows
1316  if (Mach[pnt] < 0.99)
1317  {
1318  // Kinetic energy calculation
1319  NekDouble Ek = 0.0;
1320  for (j = 1; j < nVariables-1; ++j)
1321  {
1322  Ek += 0.5 * (Fwd[j][pnt] * Fwd[j][pnt]) / Fwd[0][pnt];
1323  }
1324 
1325  rhoeb = m_pressureStorage[id1+i] * gammaMinusOneInv + Ek;
1326 
1327  // Partial extrapolation for subsonic cases
1328  for (j = 0; j < nVariables-1; ++j)
1329  {
1330  (m_fields[j]->GetBndCondExpansions()[bcRegion]->
1331  UpdatePhys())[id1+i] = Fwd[j][pnt];
1332  }
1333 
1334  (m_fields[nVariables-1]->GetBndCondExpansions()[bcRegion]->
1335  UpdatePhys())[id1+i] = rhoeb;
1336  }
1337  // Supersonic flows
1338  else
1339  {
1340  for (j = 0; j < nVariables; ++j)
1341  {
1342  // Extrapolation for supersonic cases
1343  (m_fields[j]->GetBndCondExpansions()[bcRegion]->
1344  UpdatePhys())[id1+i] = Fwd[j][pnt];
1345  }
1346  }
1347  }
1348  }
1349  }
1350 
1351 
1352  /**
1353  * @brief Pressure inflow boundary conditions for compressible flow
1354  * problems where either the density and the velocities are assigned from a
1355  * file or the full state is assigned from a file (depending on the problem
1356  * type, either subsonic or supersonic).
1357  */
1358  void CompressibleFlowSystem::PressureInflowFileBC(
1359  int bcRegion,
1360  int cnt,
1361  Array<OneD, Array<OneD, NekDouble> > &Fwd,
1362  Array<OneD, Array<OneD, NekDouble> > &physarray)
1363  {
1364  int i, j;
1365  int nTracePts = GetTraceTotPoints();
1366  int nVariables = physarray.num_elements();
1367  int nDimensions = m_spacedim;
1368 
1369  const Array<OneD, const int> &traceBndMap
1370  = m_fields[0]->GetTraceBndMap();
1371 
1372  NekDouble gamma = m_gamma;
1373  NekDouble gammaMinusOne = gamma - 1.0;
1374  NekDouble gammaMinusOneInv = 1.0 / gammaMinusOne;
1375 
1376  Array<OneD, NekDouble> tmp1 (nTracePts, 0.0);
1377  Array<OneD, NekDouble> tmp2 (nTracePts, 0.0);
1378  Array<OneD, NekDouble> VnInf(nTracePts, 0.0);
1379  Array<OneD, NekDouble> velInf(nDimensions, 0.0);
1380 
1381  // Computing the normal velocity for characteristics coming
1382  // from outside the computational domain
1383  velInf[0] = m_uInf;
1384  Vmath::Smul(nTracePts, m_uInf, m_traceNormals[0], 1, VnInf, 1);
1385  if (nDimensions == 2 || nDimensions == 3)
1386  {
1387  velInf[1] = m_vInf;
1388  Vmath::Smul(nTracePts, m_vInf, m_traceNormals[1], 1, tmp1, 1);
1389  Vmath::Vadd(nTracePts, VnInf, 1, tmp1, 1, VnInf, 1);
1390  }
1391  if (nDimensions == 3)
1392  {
1393  velInf[2] = m_wInf;
1394  Vmath::Smul(nTracePts, m_wInf, m_traceNormals[2], 1, tmp2, 1);
1395  Vmath::Vadd(nTracePts, VnInf, 1, tmp2, 1, VnInf, 1);
1396  }
1397 
1398  // Computing the normal velocity for characteristics coming
1399  // from inside the computational domain
1400  Array<OneD, NekDouble > Vn (nTracePts, 0.0);
1401  Array<OneD, NekDouble > Vel(nTracePts, 0.0);
1402  for (i = 0; i < nDimensions; ++i)
1403  {
1404  Vmath::Vdiv(nTracePts, Fwd[i+1], 1, Fwd[0], 1, Vel, 1);
1405  Vmath::Vvtvp(nTracePts, m_traceNormals[i], 1, Vel, 1, Vn, 1, Vn, 1);
1406  }
1407 
1408  // Computing the absolute value of the velocity in order to compute the
1409  // Mach number to decide whether supersonic or subsonic
1410  Array<OneD, NekDouble > absVel(nTracePts, 0.0);
1411  for (i = 0; i < nDimensions; ++i)
1412  {
1413  Vmath::Vdiv(nTracePts, Fwd[i+1], 1, Fwd[0], 1, tmp1, 1);
1414  Vmath::Vmul(nTracePts, tmp1, 1, tmp1, 1, tmp1, 1);
1415  Vmath::Vadd(nTracePts, tmp1, 1, absVel, 1, absVel, 1);
1416  }
1417  Vmath::Vsqrt(nTracePts, absVel, 1, absVel, 1);
1418 
1419  // Get speed of sound
1420  Array<OneD, NekDouble > soundSpeed (nTracePts, 0.0);
1421  Array<OneD, NekDouble > pressure (nTracePts, 0.0);
1422 
1423  for (i = 0; i < nTracePts; i++)
1424  {
1425  if (m_spacedim == 1)
1426  {
1427  pressure[i] = (gammaMinusOne) * (Fwd[2][i] -
1428  0.5 * (Fwd[1][i] * Fwd[1][i] / Fwd[0][i]));
1429  }
1430  else if (m_spacedim == 2)
1431  {
1432  pressure[i] = (gammaMinusOne) * (Fwd[3][i] -
1433  0.5 * (Fwd[1][i] * Fwd[1][i] / Fwd[0][i] +
1434  Fwd[2][i] * Fwd[2][i] / Fwd[0][i]));
1435  }
1436  else
1437  {
1438  pressure[i] = (gammaMinusOne) * (Fwd[4][i] -
1439  0.5 * (Fwd[1][i] * Fwd[1][i] / Fwd[0][i] +
1440  Fwd[2][i] * Fwd[2][i] / Fwd[0][i] +
1441  Fwd[3][i] * Fwd[3][i] / Fwd[0][i]));
1442  }
1443 
1444  soundSpeed[i] = sqrt(gamma * pressure[i] / Fwd[0][i]);
1445  }
1446 
1447  // Get Mach
1448  Array<OneD, NekDouble > Mach(nTracePts, 0.0);
1449  Vmath::Vdiv(nTracePts, Vn, 1, soundSpeed, 1, Mach, 1);
1450  Vmath::Vabs(nTracePts, Mach, 1, Mach, 1);
1451 
1452  // Auxiliary variables
1453  int e, id1, id2, npts, pnt;
1454  NekDouble rhoeb;
1455 
1456  // Loop on the bcRegions
1457  for (e = 0; e < m_fields[0]->GetBndCondExpansions()[bcRegion]->
1458  GetExpSize(); ++e)
1459  {
1460  npts = m_fields[0]->GetBndCondExpansions()[bcRegion]->
1461  GetExp(e)->GetTotPoints();
1462  id1 = m_fields[0]->GetBndCondExpansions()[bcRegion]->
1463  GetPhys_Offset(e);
1464  id2 = m_fields[0]->GetTrace()->GetPhys_Offset(traceBndMap[cnt+e]);
1465 
1466  // Loop on points of bcRegion 'e'
1467  for (i = 0; i < npts; i++)
1468  {
1469  pnt = id2+i;
1470 
1471  // Subsonic flows
1472  if (Mach[pnt] < 0.99)
1473  {
1474  // Partial extrapolation for subsonic cases
1475  for (j = 0; j < nVariables-1; ++j)
1476  {
1477  (m_fields[j]->GetBndCondExpansions()[bcRegion]->
1478  UpdatePhys())[id1+i] = m_fieldStorage[j][id1+i];
1479  }
1480 
1481  // Kinetic energy calculation
1482  NekDouble Ek = 0.0;
1483  for (j = 1; j < nVariables-1; ++j)
1484  {
1485  Ek += 0.5 * (m_fieldStorage[j][id1+i] *
1486  m_fieldStorage[j][id1+i]) /
1487  m_fieldStorage[0][id1+i];
1488  }
1489  rhoeb = gammaMinusOneInv * pressure[pnt] + Ek;
1490 
1491  (m_fields[nVariables-1]->GetBndCondExpansions()[bcRegion]->
1492  UpdatePhys())[id1+i] = rhoeb;
1493  }
1494  // Supersonic flows
1495  else
1496  {
1497  for (j = 0; j < nVariables; ++j)
1498  {
1499  // Extrapolation for supersonic cases
1500  (m_fields[j]->GetBndCondExpansions()[bcRegion]->
1501  UpdatePhys())[id1+i] = Fwd[j][pnt];
1502  }
1503  }
1504  }
1505  }
1506  }
1507 
1508 
1509  /**
1510  * @brief Extrapolation of order 0 for all the variables such that,
1511  * at the boundaries, a trivial Riemann problem is solved.
1512  */
1513  void CompressibleFlowSystem::ExtrapOrder0BC(
1514  int bcRegion,
1515  int cnt,
1516  Array<OneD, Array<OneD, NekDouble> > &Fwd,
1517  Array<OneD, Array<OneD, NekDouble> > &physarray)
1518  {
1519  int i, j;
1520  int e, pnt;
1521  int id1, id2, nBCEdgePts;
1522  int nVariables = physarray.num_elements();
1523  int nDimensions = m_spacedim;
1524 
1525  const Array<OneD, const int> &traceBndMap
1526  = m_fields[0]->GetTraceBndMap();
1527 
1528  int eMax;
1529 
1530  eMax = m_fields[0]->GetBndCondExpansions()[bcRegion]->GetExpSize();
1531 
1532  // Loop on bcRegions
1533  for (e = 0; e < eMax; ++e)
1534  {
1535  nBCEdgePts = m_fields[0]->GetBndCondExpansions()[bcRegion]->
1536  GetExp(e)->GetTotPoints();
1537  id1 = m_fields[0]->GetBndCondExpansions()[bcRegion]->
1538  GetPhys_Offset(e) ;
1539  id2 = m_fields[0]->GetTrace()->GetPhys_Offset(traceBndMap[cnt+e]);
1540 
1541  // Loop on points of bcRegion 'e'
1542  for (i = 0; i < nBCEdgePts; i++)
1543  {
1544  pnt = id2+i;
1545 
1546  // Setting up bcs for density
1547  (m_fields[0]->GetBndCondExpansions()[bcRegion]->
1548  UpdatePhys())[id1+i] = Fwd[0][pnt];
1549 
1550  // Setting up bcs for velocities
1551  for (j = 1; j <=nDimensions; ++j)
1552  {
1553  (m_fields[j]->GetBndCondExpansions()[bcRegion]->
1554  UpdatePhys())[id1+i] = Fwd[j][pnt];
1555  }
1556 
1557  // Setting up bcs for energy
1558  (m_fields[nVariables-1]->GetBndCondExpansions()[bcRegion]->
1559  UpdatePhys())[id1+i] = Fwd[nVariables-1][pnt];
1560  }
1561  }
1562  }
1563 
1564  /**
1565  * @brief Return the flux vector for the compressible Euler equations.
1566  *
1567  * @param i Component of the flux vector to calculate.
1568  * @param physfield Fields.
1569  * @param flux Resulting flux.
1570  */
1571  void CompressibleFlowSystem::GetFluxVector(
1572  const Array<OneD, Array<OneD, NekDouble> > &physfield,
1573  Array<OneD, Array<OneD, Array<OneD, NekDouble> > > &flux)
1574  {
1575  int i, j;
1576  int nq = physfield[0].num_elements();
1577  int nVariables = m_fields.num_elements();
1578 
1579  Array<OneD, NekDouble> pressure(nq);
1580  Array<OneD, Array<OneD, NekDouble> > velocity(m_spacedim);
1581 
1582  // Flux vector for the rho equation
1583  for (i = 0; i < m_spacedim; ++i)
1584  {
1585  velocity[i] = Array<OneD, NekDouble>(nq);
1586  Vmath::Vcopy(nq, physfield[i+1], 1, flux[0][i], 1);
1587  }
1588 
1589  GetVelocityVector(physfield, velocity);
1590  GetPressure (physfield, velocity, pressure);
1591 
1592  // Flux vector for the velocity fields
1593  for (i = 0; i < m_spacedim; ++i)
1594  {
1595  for (j = 0; j < m_spacedim; ++j)
1596  {
1597  Vmath::Vmul(nq, velocity[j], 1, physfield[i+1], 1,
1598  flux[i+1][j], 1);
1599  }
1600 
1601  // Add pressure to appropriate field
1602  Vmath::Vadd(nq, flux[i+1][i], 1, pressure, 1, flux[i+1][i], 1);
1603  }
1604 
1605  // Flux vector for energy.
1606  Vmath::Vadd(nq, physfield[m_spacedim+1], 1, pressure, 1,
1607  pressure, 1);
1608 
1609  for (j = 0; j < m_spacedim; ++j)
1610  {
1611  Vmath::Vmul(nq, velocity[j], 1, pressure, 1,
1612  flux[m_spacedim+1][j], 1);
1613  }
1614 
1615  // For the smooth viscosity model
1616  if (nVariables == m_spacedim+3)
1617  {
1618  // Add a zero row for the advective fluxes
1619  for (j = 0; j < m_spacedim; ++j)
1620  {
1621  Vmath::Zero(nq, flux[m_spacedim+2][j], 1);
1622  }
1623  }
1624  }
1625 
1626  /**
1627  * @brief Return the flux vector for the compressible Euler equations
1628  * by using the de-aliasing technique.
1629  *
1630  * @param i Component of the flux vector to calculate.
1631  * @param physfield Fields.
1632  * @param flux Resulting flux.
1633  */
1634  void CompressibleFlowSystem::GetFluxVectorDeAlias(
1635  const Array<OneD, Array<OneD, NekDouble> > &physfield,
1636  Array<OneD, Array<OneD, Array<OneD, NekDouble> > > &flux)
1637  {
1638  int i, j;
1639  int nq = physfield[0].num_elements();
1640  int nVariables = m_fields.num_elements();
1641 
1642  // Factor to rescale 1d points in dealiasing
1643  NekDouble OneDptscale = 2;
1644  nq = m_fields[0]->Get1DScaledTotPoints(OneDptscale);
1645 
1646  Array<OneD, NekDouble> pressure(nq);
1647  Array<OneD, Array<OneD, NekDouble> > velocity(m_spacedim);
1648 
1649  Array<OneD, Array<OneD, NekDouble> > physfield_interp(nVariables);
1650  Array<OneD, Array<OneD, Array<OneD, NekDouble> > > flux_interp(
1651  nVariables);
1652 
1653  for (i = 0; i < nVariables; ++ i)
1654  {
1655  physfield_interp[i] = Array<OneD, NekDouble>(nq);
1656  flux_interp[i] = Array<OneD, Array<OneD, NekDouble> >(m_spacedim);
1657  m_fields[0]->PhysInterp1DScaled(
1658  OneDptscale, physfield[i], physfield_interp[i]);
1659 
1660  for (j = 0; j < m_spacedim; ++j)
1661  {
1662  flux_interp[i][j] = Array<OneD, NekDouble>(nq);
1663  }
1664  }
1665 
1666  // Flux vector for the rho equation
1667  for (i = 0; i < m_spacedim; ++i)
1668  {
1669  velocity[i] = Array<OneD, NekDouble>(nq);
1670 
1671  // Galerkin project solution back to original space
1672  m_fields[0]->PhysGalerkinProjection1DScaled(
1673  OneDptscale, physfield_interp[i+1], flux[0][i]);
1674  }
1675 
1676  GetVelocityVector(physfield_interp, velocity);
1677  GetPressure (physfield_interp, velocity, pressure);
1678 
1679  // Evaluation of flux vector for the velocity fields
1680  for (i = 0; i < m_spacedim; ++i)
1681  {
1682  for (j = 0; j < m_spacedim; ++j)
1683  {
1684  Vmath::Vmul(nq, velocity[j], 1, physfield_interp[i+1], 1,
1685  flux_interp[i+1][j], 1);
1686  }
1687 
1688  // Add pressure to appropriate field
1689  Vmath::Vadd(nq, flux_interp[i+1][i], 1, pressure,1,
1690  flux_interp[i+1][i], 1);
1691  }
1692 
1693  // Galerkin project solution back to origianl space
1694  for (i = 0; i < m_spacedim; ++i)
1695  {
1696  for (j = 0; j < m_spacedim; ++j)
1697  {
1698  m_fields[0]->PhysGalerkinProjection1DScaled(
1699  OneDptscale, flux_interp[i+1][j], flux[i+1][j]);
1700  }
1701  }
1702 
1703  // Evaluation of flux vector for energy
1704  Vmath::Vadd(nq, physfield_interp[m_spacedim+1], 1, pressure, 1,
1705  pressure, 1);
1706 
1707  for (j = 0; j < m_spacedim; ++j)
1708  {
1709  Vmath::Vmul(nq, velocity[j], 1, pressure, 1,
1710  flux_interp[m_spacedim+1][j], 1);
1711 
1712  // Galerkin project solution back to origianl space
1713  m_fields[0]->PhysGalerkinProjection1DScaled(
1714  OneDptscale,
1715  flux_interp[m_spacedim+1][j],
1716  flux[m_spacedim+1][j]);
1717  }
1718  }
1719 
1720 
1721  /**
1722  * @brief Return the flux vector for the LDG diffusion problem.
1723  * \todo Complete the viscous flux vector
1724  */
1725  void CompressibleFlowSystem::GetViscousFluxVector(
1726  const Array<OneD, Array<OneD, NekDouble> > &physfield,
1727  Array<OneD, Array<OneD, Array<OneD, NekDouble> > > &derivativesO1,
1728  Array<OneD, Array<OneD, Array<OneD, NekDouble> > > &viscousTensor)
1729  {
1730  int j, k;
1731  int nVariables = m_fields.num_elements();
1732  int nPts = physfield[0].num_elements();
1733 
1734  // Stokes hypotesis
1735  const NekDouble lambda = -2.0/3.0;
1736 
1737  // Auxiliary variables
1738  Array<OneD, NekDouble > mu (nPts, 0.0);
1739  Array<OneD, NekDouble > thermalConductivity(nPts, 0.0);
1740  Array<OneD, NekDouble > mu2 (nPts, 0.0);
1741  Array<OneD, NekDouble > divVel (nPts, 0.0);
1742 
1743  // Variable viscosity through the Sutherland's law
1744  if (m_ViscosityType == "Variable")
1745  {
1746  GetDynamicViscosity(physfield[nVariables-2], mu);
1747  NekDouble tRa = m_Cp / m_Prandtl;
1748  Vmath::Smul(nPts, tRa, &mu[0], 1, &thermalConductivity[0], 1);
1749  }
1750  else
1751  {
1752  Vmath::Fill(nPts, m_mu, &mu[0], 1);
1753  Vmath::Fill(nPts, m_thermalConductivity,
1754  &thermalConductivity[0], 1);
1755  }
1756 
1757  // Computing diagonal terms of viscous stress tensor
1758  Array<OneD, Array<OneD, NekDouble> > tmp(m_spacedim);
1759  Array<OneD, Array<OneD, NekDouble> > Sgg(m_spacedim);
1760 
1761  // mu2 = 2 * mu
1762  Vmath::Smul(nPts, 2.0, &mu[0], 1, &mu2[0], 1);
1763 
1764  // Velocity divergence
1765  for (j = 0; j < m_spacedim; ++j)
1766  {
1767  Vmath::Vadd(nPts, &divVel[0], 1, &derivativesO1[j][j][0], 1,
1768  &divVel[0], 1);
1769  }
1770 
1771  // Velocity divergence scaled by lambda * mu
1772  Vmath::Smul(nPts, lambda, &divVel[0], 1, &divVel[0], 1);
1773  Vmath::Vmul(nPts, &mu[0], 1, &divVel[0], 1, &divVel[0], 1);
1774 
1775  // Diagonal terms of viscous stress tensor (Sxx, Syy, Szz)
1776  // Sjj = 2 * mu * du_j/dx_j - (2 / 3) * mu * sum_j(du_j/dx_j)
1777  for (j = 0; j < m_spacedim; ++j)
1778  {
1779  tmp[j] = Array<OneD, NekDouble>(nPts, 0.0);
1780  Sgg[j] = Array<OneD, NekDouble>(nPts, 0.0);
1781 
1782  Vmath::Vmul(nPts, &mu2[0], 1, &derivativesO1[j][j][0], 1,
1783  &tmp[j][0], 1);
1784 
1785  Vmath::Vadd(nPts, &tmp[j][0], 1, &divVel[0], 1, &Sgg[j][0], 1);
1786  }
1787 
1788  // Extra diagonal terms of viscous stress tensor (Sxy, Sxz, Syz)
1789  // Note: they exist for 2D and 3D problems only
1790  Array<OneD, NekDouble > Sxy(nPts, 0.0);
1791  Array<OneD, NekDouble > Sxz(nPts, 0.0);
1792  Array<OneD, NekDouble > Syz(nPts, 0.0);
1793 
1794  if (m_spacedim == 2)
1795  {
1796  // Sxy = (du/dy + dv/dx)
1797  Vmath::Vadd(nPts, &derivativesO1[0][1][0], 1,
1798  &derivativesO1[1][0][0], 1, &Sxy[0], 1);
1799 
1800  // Sxy = mu * (du/dy + dv/dx)
1801  Vmath::Vmul(nPts, &mu[0], 1, &Sxy[0], 1, &Sxy[0], 1);
1802  }
1803  else if (m_spacedim == 3)
1804  {
1805  // Sxy = (du/dy + dv/dx)
1806  Vmath::Vadd(nPts, &derivativesO1[0][1][0], 1,
1807  &derivativesO1[1][0][0], 1, &Sxy[0], 1);
1808 
1809  // Sxz = (du/dz + dw/dx)
1810  Vmath::Vadd(nPts, &derivativesO1[0][2][0], 1,
1811  &derivativesO1[2][0][0], 1, &Sxz[0], 1);
1812 
1813  // Syz = (dv/dz + dw/dy)
1814  Vmath::Vadd(nPts, &derivativesO1[1][2][0], 1,
1815  &derivativesO1[2][1][0], 1, &Syz[0], 1);
1816 
1817  // Sxy = mu * (du/dy + dv/dx)
1818  Vmath::Vmul(nPts, &mu[0], 1, &Sxy[0], 1, &Sxy[0], 1);
1819 
1820  // Sxz = mu * (du/dy + dv/dx)
1821  Vmath::Vmul(nPts, &mu[0], 1, &Sxz[0], 1, &Sxz[0], 1);
1822 
1823  // Syz = mu * (du/dy + dv/dx)
1824  Vmath::Vmul(nPts, &mu[0], 1, &Syz[0], 1, &Syz[0], 1);
1825  }
1826 
1827  // Energy-related terms
1828  Array<OneD, NekDouble > STx(nPts, 0.0);
1829  Array<OneD, NekDouble > STy(nPts, 0.0);
1830  Array<OneD, NekDouble > STz(nPts, 0.0);
1831 
1832  // Building the viscous flux vector
1833 
1834  // Viscous flux vector for the rho equation
1835  for (k = 0; k < m_spacedim; ++k)
1836  {
1837  Vmath::Zero(nPts, viscousTensor[k][0], 1);
1838  }
1839 
1840  if (m_spacedim == 1)
1841  {
1842  Array<OneD, NekDouble > tmp1(nPts, 0.0);
1843 
1844  // u * Sxx
1845  Vmath::Vmul(nPts, &physfield[0][0], 1, &Sgg[0][0], 1, &STx[0], 1);
1846 
1847  // k * dT/dx
1848  Vmath::Vmul(nPts, &thermalConductivity[0], 1,
1849  &derivativesO1[0][1][0], 1, &tmp1[0], 1);
1850 
1851  // STx = u * Sxx + (K / mu) * dT/dx
1852  Vmath::Vadd(nPts, &STx[0], 1, &tmp1[0], 1, &STx[0], 1);
1853  }
1854  else if (m_spacedim == 2)
1855  {
1856  Array<OneD, NekDouble > tmp1(nPts, 0.0);
1857  Array<OneD, NekDouble > tmp2(nPts, 0.0);
1858 
1859  // Computation of STx
1860 
1861  // u * Sxx
1862  Vmath::Vmul(nPts, &physfield[0][0], 1, &Sgg[0][0], 1, &STx[0], 1);
1863 
1864  // v * Sxy
1865  Vmath::Vmul(nPts, &physfield[1][0], 1, &Sxy[0], 1, &tmp1[0], 1);
1866 
1867  // k * dT/dx
1868  Vmath::Vmul(nPts, &thermalConductivity[0], 1,
1869  &derivativesO1[0][2][0], 1, &tmp2[0], 1);
1870 
1871  // STx = u * Sxx + v * Sxy + K * dT/dx
1872  Vmath::Vadd(nPts, &STx[0], 1, &tmp1[0], 1, &STx[0], 1);
1873  Vmath::Vadd(nPts, &STx[0], 1, &tmp2[0], 1, &STx[0], 1);
1874 
1875  // Computation of STy
1876 
1877  // Re-initialise temporary arrays
1878  Vmath::Zero(nPts, &tmp1[0], 1);
1879  Vmath::Zero(nPts, &tmp2[0], 1);
1880 
1881  // v * Syy
1882  Vmath::Vmul(nPts, &physfield[1][0], 1, &Sgg[1][0], 1, &STy[0], 1);
1883 
1884  // u * Sxy
1885  Vmath::Vmul(nPts, &physfield[0][0], 1, &Sxy[0], 1, &tmp1[0], 1);
1886 
1887  // k * dT/dy
1888  Vmath::Vmul(nPts, &thermalConductivity[0], 1,
1889  &derivativesO1[1][2][0], 1, &tmp2[0], 1);
1890 
1891  // STy = v * Syy + u * Sxy + K * dT/dy
1892  Vmath::Vadd(nPts, &STy[0], 1, &tmp1[0], 1, &STy[0], 1);
1893  Vmath::Vadd(nPts, &STy[0], 1, &tmp2[0], 1, &STy[0], 1);
1894  }
1895  else if (m_spacedim == 3)
1896  {
1897  Array<OneD, NekDouble > tmp1(nPts, 0.0);
1898  Array<OneD, NekDouble > tmp2(nPts, 0.0);
1899  Array<OneD, NekDouble > tmp3(nPts, 0.0);
1900 
1901  // Computation of STx
1902 
1903  // u * Sxx
1904  Vmath::Vmul(nPts, &physfield[0][0], 1, &Sgg[0][0], 1, &STx[0], 1);
1905 
1906  // v * Sxy
1907  Vmath::Vmul(nPts, &physfield[1][0], 1, &Sxy[0], 1, &tmp1[0], 1);
1908 
1909  // v * Sxz
1910  Vmath::Vmul(nPts, &physfield[2][0], 1, &Sxz[0], 1, &tmp2[0], 1);
1911 
1912  // k * dT/dx
1913  Vmath::Vmul(nPts, &thermalConductivity[0], 1,
1914  &derivativesO1[0][3][0], 1, &tmp3[0], 1);
1915 
1916  // STx = u * Sxx + v * Sxy + w * Sxz + (K / mu) * dT/dx
1917  Vmath::Vadd(nPts, &STx[0], 1, &tmp1[0], 1, &STx[0], 1);
1918  Vmath::Vadd(nPts, &STx[0], 1, &tmp2[0], 1, &STx[0], 1);
1919  Vmath::Vadd(nPts, &STx[0], 1, &tmp3[0], 1, &STx[0], 1);
1920 
1921  // Computation of STy
1922 
1923  // Re-initialise temporary arrays
1924  Vmath::Zero(nPts, &tmp1[0], 1);
1925  Vmath::Zero(nPts, &tmp2[0], 1);
1926  Vmath::Zero(nPts, &tmp3[0], 1);
1927 
1928  // v * Syy
1929  Vmath::Vmul(nPts, &physfield[1][0], 1, &Sgg[1][0], 1, &STy[0], 1);
1930 
1931  // u * Sxy
1932  Vmath::Vmul(nPts, &physfield[0][0], 1, &Sxy[0], 1, &tmp1[0], 1);
1933 
1934  // w * Syz
1935  Vmath::Vmul(nPts, &physfield[2][0], 1, &Syz[0], 1, &tmp2[0], 1);
1936 
1937  // k * dT/dy
1938  Vmath::Vmul(nPts, &thermalConductivity[0], 1,
1939  &derivativesO1[1][3][0], 1, &tmp3[0], 1);
1940 
1941  // STy = v * Syy + u * Sxy + w * Syz + K * dT/dy
1942  Vmath::Vadd(nPts, &STy[0], 1, &tmp1[0], 1, &STy[0], 1);
1943  Vmath::Vadd(nPts, &STy[0], 1, &tmp2[0], 1, &STy[0], 1);
1944  Vmath::Vadd(nPts, &STy[0], 1, &tmp3[0], 1, &STy[0], 1);
1945 
1946  // Computation of STz
1947 
1948  // Re-initialise temporary arrays
1949  Vmath::Zero(nPts, &tmp1[0], 1);
1950  Vmath::Zero(nPts, &tmp2[0], 1);
1951  Vmath::Zero(nPts, &tmp3[0], 1);
1952 
1953  // w * Szz
1954  Vmath::Vmul(nPts, &physfield[2][0], 1, &Sgg[2][0], 1, &STz[0], 1);
1955 
1956  // u * Sxz
1957  Vmath::Vmul(nPts, &physfield[0][0], 1, &Sxz[0], 1, &tmp1[0], 1);
1958 
1959  // v * Syz
1960  Vmath::Vmul(nPts, &physfield[1][0], 1, &Syz[0], 1, &tmp2[0], 1);
1961 
1962  // k * dT/dz
1963  Vmath::Vmul(nPts, &thermalConductivity[0], 1,
1964  &derivativesO1[2][3][0], 1, &tmp3[0], 1);
1965 
1966  // STz = w * Szz + u * Sxz + v * Syz + K * dT/dz
1967  Vmath::Vadd(nPts, &STz[0], 1, &tmp1[0], 1, &STz[0], 1);
1968  Vmath::Vadd(nPts, &STz[0], 1, &tmp2[0], 1, &STz[0], 1);
1969  Vmath::Vadd(nPts, &STz[0], 1, &tmp3[0], 1, &STz[0], 1);
1970  }
1971 
1972  switch (m_spacedim)
1973  {
1974  case 1:
1975  {
1976  // f_11v = f_rho = 0
1977  Vmath::Zero(nPts, &viscousTensor[0][0][0], 1);
1978 
1979  // f_21v = f_rhou
1980  Vmath::Vcopy(nPts, &Sgg[0][0], 1, &viscousTensor[0][1][0], 1);
1981 
1982  // f_31v = f_E
1983  Vmath::Vcopy(nPts, &STx[0], 1, &viscousTensor[0][2][0], 1);
1984  break;
1985  }
1986  case 2:
1987  {
1988  // f_11v = f_rho1 = 0
1989  Vmath::Zero(nPts, &viscousTensor[0][0][0], 1);
1990  // f_12v = f_rho2 = 0
1991  Vmath::Zero(nPts, &viscousTensor[1][0][0], 1);
1992 
1993  // f_21v = f_rhou1
1994  Vmath::Vcopy(nPts, &Sgg[0][0], 1, &viscousTensor[0][1][0], 1);
1995  // f_22v = f_rhou2
1996  Vmath::Vcopy(nPts, &Sxy[0], 1, &viscousTensor[1][1][0], 1);
1997 
1998  // f_31v = f_rhov1
1999  Vmath::Vcopy(nPts, &Sxy[0], 1, &viscousTensor[0][2][0], 1);
2000  // f_32v = f_rhov2
2001  Vmath::Vcopy(nPts, &Sgg[1][0], 1, &viscousTensor[1][2][0], 1);
2002 
2003  // f_41v = f_E1
2004  Vmath::Vcopy(nPts, &STx[0], 1, &viscousTensor[0][3][0], 1);
2005  // f_42v = f_E2
2006  Vmath::Vcopy(nPts, &STy[0], 1, &viscousTensor[1][3][0], 1);
2007  break;
2008  }
2009  case 3:
2010  {
2011  // f_11v = f_rho1 = 0
2012  Vmath::Zero(nPts, &viscousTensor[0][0][0], 1);
2013  // f_12v = f_rho2 = 0
2014  Vmath::Zero(nPts, &viscousTensor[1][0][0], 1);
2015  // f_13v = f_rho3 = 0
2016  Vmath::Zero(nPts, &viscousTensor[2][0][0], 1);
2017 
2018  // f_21v = f_rhou1
2019  Vmath::Vcopy(nPts, &Sgg[0][0], 1, &viscousTensor[0][1][0], 1);
2020  // f_22v = f_rhou2
2021  Vmath::Vcopy(nPts, &Sxy[0], 1, &viscousTensor[1][1][0], 1);
2022  // f_23v = f_rhou3
2023  Vmath::Vcopy(nPts, &Sxz[0], 1, &viscousTensor[2][1][0], 1);
2024 
2025  // f_31v = f_rhov1
2026  Vmath::Vcopy(nPts, &Sxy[0], 1, &viscousTensor[0][2][0], 1);
2027  // f_32v = f_rhov2
2028  Vmath::Vcopy(nPts, &Sgg[1][0], 1, &viscousTensor[1][2][0], 1);
2029  // f_33v = f_rhov3
2030  Vmath::Vcopy(nPts, &Syz[0], 1, &viscousTensor[2][2][0], 1);
2031 
2032  // f_31v = f_rhow1
2033  Vmath::Vcopy(nPts, &Sxz[0], 1, &viscousTensor[0][3][0], 1);
2034  // f_32v = f_rhow2
2035  Vmath::Vcopy(nPts, &Syz[0], 1, &viscousTensor[1][3][0], 1);
2036  // f_33v = f_rhow3
2037  Vmath::Vcopy(nPts, &Sgg[2][0], 1, &viscousTensor[2][3][0], 1);
2038 
2039  // f_41v = f_E1
2040  Vmath::Vcopy(nPts, &STx[0], 1, &viscousTensor[0][4][0], 1);
2041  // f_42v = f_E2
2042  Vmath::Vcopy(nPts, &STy[0], 1, &viscousTensor[1][4][0], 1);
2043  // f_43v = f_E3
2044  Vmath::Vcopy(nPts, &STz[0], 1, &viscousTensor[2][4][0], 1);
2045  break;
2046  }
2047  default:
2048  {
2049  ASSERTL0(false, "Illegal expansion dimension");
2050  }
2051  }
2052  }
2053 
2054  /**
2055  * @brief Return the flux vector for the LDG diffusion problem.
2056  * \todo Complete the viscous flux vector
2057  */
2058  void CompressibleFlowSystem::GetViscousFluxVectorDeAlias(
2059  const Array<OneD, Array<OneD, NekDouble> > &physfield,
2060  Array<OneD, Array<OneD, Array<OneD, NekDouble> > > &derivativesO1,
2061  Array<OneD, Array<OneD, Array<OneD, NekDouble> > > &viscousTensor)
2062  {
2063 #if 0
2064  int i, j, k;
2065  int nVariables = m_fields.num_elements();
2066  int nPts = physfield[0].num_elements();
2067 
2068  int variables_phys = physfield.num_elements();
2069 
2070  // Factor to rescale 1d points in dealiasing.
2071  NekDouble OneDptscale = 2;
2072 
2073  // Get number of points to dealias a cubic non-linearity
2074  nPts = m_fields[0]->Get1DScaledTotPoints(OneDptscale);
2075 
2076  int nVariables_aux = derivativesO1[0].num_elements();
2077 
2078  Array<OneD, Array<OneD, NekDouble> > physfield_interp(variables_phys);
2079  Array<OneD, Array<OneD, Array<OneD, NekDouble> > > derivativesO1_interp(
2080  m_spacedim);
2081  Array<OneD, Array<OneD, Array<OneD, NekDouble> > >viscousTensor_interp(
2082  m_spacedim);
2083 
2084  for (i = 0; i < m_spacedim; ++ i)
2085  {
2086  viscousTensor_interp[i] = Array<OneD, Array<OneD, NekDouble> >(
2087  nVariables);
2088  for (j = 0; j < nVariables; ++j)
2089  {
2090  viscousTensor_interp[i][j] = Array<OneD, NekDouble>(nPts);
2091  }
2092  }
2093 
2094  // Stokes hypotesis
2095  NekDouble lambda = -2.0/3.0;
2096 
2097  // Auxiliary variables
2098  Array<OneD, NekDouble > mu (nPts, 0.0);
2099  Array<OneD, NekDouble > mu2 (nPts, 0.0);
2100  Array<OneD, NekDouble > divVel (nPts, 0.0);
2101  Array<OneD, NekDouble > pressure (nPts, 0.0);
2102  Array<OneD, NekDouble > temperature(nPts, 0.0);
2103 
2104  for (i = 0; i < nVariables; ++i)
2105  {
2106  m_fields[0]->PhysInterp1DScaled(
2107  OneDptscale, physfield[i], fields_interp[i]);
2108  }
2109 
2110  for (i = 0; i < variables_phys; ++i)
2111  {
2112  physfield_interp[i] = Array<OneD, NekDouble> (nPts);
2113 
2114  // Interpolation to higher space
2115  m_fields[0]->PhysInterp1DScaled(OneDptscale, physfield[i],
2116  physfield_interp[i]);
2117  }
2118 
2119  for (i = 0; i < m_spacedim; ++i)
2120  {
2121  derivativesO1_interp[i] = Array<OneD, Array<OneD, NekDouble> >(
2122  nVariables_aux);
2123  for (j = 0; j < nVariables_aux; ++j)
2124  {
2125  derivativesO1_interp[i][j] = Array<OneD, NekDouble>(nPts);
2126  m_fields[0]->PhysInterp1DScaled(
2127  OneDptscale, derivativesO1[i][j],
2128  derivativesO1_interp[i][j]);
2129  }
2130  }
2131 
2132  // Thermodynamic related quantities
2133  GetPressure(fields_interp, pressure);
2134  GetTemperature(fields_interp, pressure, temperature);
2135 
2136  // Variable viscosity through the Sutherland's law
2137  if (m_ViscosityType == "Variable")
2138  {
2139  GetDynamicViscosity(fields_interp[variables_phys-1], mu);
2140  }
2141  else
2142  {
2143  Vmath::Sadd(nPts, m_mu, &mu[0], 1, &mu[0], 1);
2144  }
2145 
2146  // Computing diagonal terms of viscous stress tensor
2147  Array<OneD, Array<OneD, NekDouble> > tmp(m_spacedim);
2148  Array<OneD, Array<OneD, NekDouble> > Sgg(m_spacedim);
2149 
2150  // mu2 = 2 * mu
2151  Vmath::Smul(nPts, 2.0, &mu[0], 1, &mu2[0], 1);
2152 
2153  // Velocity divergence
2154  for (j = 0; j < m_spacedim; ++j)
2155  {
2156  Vmath::Vadd(nPts, &divVel[0], 1, &derivativesO1_interp[j][j][0], 1,
2157  &divVel[0], 1);
2158  }
2159 
2160  // Velocity divergence scaled by lambda * mu
2161  Vmath::Smul(nPts, lambda, &divVel[0], 1, &divVel[0], 1);
2162  Vmath::Vmul(nPts, &mu[0], 1, &divVel[0], 1, &divVel[0], 1);
2163 
2164  // Digonal terms of viscous stress tensor (Sxx, Syy, Szz)
2165  // Sjj = 2 * mu * du_j/dx_j - (2 / 3) * mu * sum_j(du_j/dx_j)
2166  for (j = 0; j < m_spacedim; ++j)
2167  {
2168  tmp[j] = Array<OneD, NekDouble>(nPts, 0.0);
2169  Sgg[j] = Array<OneD, NekDouble>(nPts, 0.0);
2170 
2171  Vmath::Vmul(nPts, &mu2[0], 1, &derivativesO1_interp[j][j][0], 1,
2172  &tmp[j][0], 1);
2173 
2174  Vmath::Vadd(nPts, &tmp[j][0], 1, &divVel[0], 1, &Sgg[j][0], 1);
2175  }
2176 
2177  // Extra diagonal terms of viscous stress tensor (Sxy, Sxz, Syz)
2178  // Note: they exist for 2D and 3D problems only
2179  Array<OneD, NekDouble > Sxy(nPts, 0.0);
2180  Array<OneD, NekDouble > Sxz(nPts, 0.0);
2181  Array<OneD, NekDouble > Syz(nPts, 0.0);
2182 
2183  if (m_spacedim == 2)
2184  {
2185  // Sxy = (du/dy + dv/dx)
2186  Vmath::Vadd(nPts, &derivativesO1_interp[0][1][0], 1,
2187  &derivativesO1_interp[1][0][0], 1, &Sxy[0], 1);
2188 
2189  // Sxy = mu * (du/dy + dv/dx)
2190  Vmath::Vmul(nPts, &mu[0], 1, &Sxy[0], 1, &Sxy[0], 1);
2191  }
2192  else if (m_spacedim == 3)
2193  {
2194  // Sxy = (du/dy + dv/dx)
2195  Vmath::Vadd(nPts, &derivativesO1_interp[0][1][0], 1,
2196  &derivativesO1_interp[1][0][0], 1, &Sxy[0], 1);
2197 
2198  // Sxz = (du/dz + dw/dx)
2199  Vmath::Vadd(nPts, &derivativesO1_interp[0][2][0], 1,
2200  &derivativesO1_interp[2][0][0], 1, &Sxz[0], 1);
2201 
2202  // Syz = (dv/dz + dw/dy)
2203  Vmath::Vadd(nPts, &derivativesO1_interp[1][2][0], 1,
2204  &derivativesO1_interp[2][1][0], 1, &Syz[0], 1);
2205 
2206  // Sxy = mu * (du/dy + dv/dx)
2207  Vmath::Vmul(nPts, &mu[0], 1, &Sxy[0], 1, &Sxy[0], 1);
2208 
2209  // Sxz = mu * (du/dy + dv/dx)
2210  Vmath::Vmul(nPts, &mu[0], 1, &Sxz[0], 1, &Sxz[0], 1);
2211 
2212  // Syz = mu * (du/dy + dv/dx)
2213  Vmath::Vmul(nPts, &mu[0], 1, &Syz[0], 1, &Syz[0], 1);
2214  }
2215 
2216  // Energy-related terms
2217  Array<OneD, NekDouble > STx(nPts, 0.0);
2218  Array<OneD, NekDouble > STy(nPts, 0.0);
2219  Array<OneD, NekDouble > STz(nPts, 0.0);
2220  // Building the viscous flux vector
2221  if (i == 0)
2222  {
2223  // Viscous flux vector for the rho equation
2224  for (k = 0; k < m_spacedim; ++k)
2225  {
2226  Vmath::Zero(nPts, viscousTensor_interp[k][i], 1);
2227  }
2228  }
2229 
2230  if (m_spacedim == 1)
2231  {
2232  Array<OneD, NekDouble > tmp1(nPts, 0.0);
2233 
2234  // u * Sxx
2235  Vmath::Vmul(nPts, &physfield_interp[0][0], 1,
2236  &Sgg[0][0], 1, &STx[0], 1);
2237 
2238  // k * dT/dx
2239  Vmath::Smul(nPts, m_thermalConductivity,
2240  &derivativesO1_interp[0][1][0], 1, &tmp1[0], 1);
2241 
2242  // STx = u * Sxx + (K / mu) * dT/dx
2243  Vmath::Vadd(nPts, &STx[0], 1, &tmp1[0], 1, &STx[0], 1);
2244  }
2245  else if (m_spacedim == 2)
2246  {
2247  Array<OneD, NekDouble > tmp1(nPts, 0.0);
2248  Array<OneD, NekDouble > tmp2(nPts, 0.0);
2249 
2250  // Computation of STx
2251 
2252  // u * Sxx
2253  Vmath::Vmul(nPts, &physfield_interp[0][0], 1,
2254  &Sgg[0][0], 1, &STx[0], 1);
2255 
2256  // v * Sxy
2257  Vmath::Vmul(nPts, &physfield_interp[1][0], 1,
2258  &Sxy[0], 1, &tmp1[0], 1);
2259 
2260  // k * dT/dx
2261  Vmath::Smul(nPts, m_thermalConductivity,
2262  &derivativesO1_interp[0][2][0], 1, &tmp2[0], 1);
2263 
2264  // STx = u * Sxx + v * Sxy + K * dT/dx
2265  Vmath::Vadd(nPts, &STx[0], 1, &tmp1[0], 1, &STx[0], 1);
2266  Vmath::Vadd(nPts, &STx[0], 1, &tmp2[0], 1, &STx[0], 1);
2267 
2268  // Computation of STy
2269 
2270  // Re-initialise temporary arrays
2271  Vmath::Zero(nPts, &tmp1[0], 1);
2272  Vmath::Zero(nPts, &tmp2[0], 1);
2273 
2274  // v * Syy
2275  Vmath::Vmul(nPts, &physfield_interp[1][0], 1,
2276  &Sgg[1][0], 1, &STy[0], 1);
2277 
2278  // u * Sxy
2279  Vmath::Vmul(nPts, &physfield_interp[0][0], 1,
2280  &Sxy[0], 1, &tmp1[0], 1);
2281 
2282  // k * dT/dy
2283  Vmath::Smul(nPts, m_thermalConductivity,
2284  &derivativesO1_interp[1][2][0], 1, &tmp2[0], 1);
2285 
2286  // STy = v * Syy + u * Sxy + K * dT/dy
2287  Vmath::Vadd(nPts, &STy[0], 1, &tmp1[0], 1, &STy[0], 1);
2288  Vmath::Vadd(nPts, &STy[0], 1, &tmp2[0], 1, &STy[0], 1);
2289  }
2290  else if (m_spacedim == 3)
2291  {
2292  Array<OneD, NekDouble > tmp1(nPts, 0.0);
2293  Array<OneD, NekDouble > tmp2(nPts, 0.0);
2294  Array<OneD, NekDouble > tmp3(nPts, 0.0);
2295 
2296  // Computation of STx
2297 
2298  // u * Sxx
2299  Vmath::Vmul(nPts, &physfield_interp[0][0], 1,
2300  &Sgg[0][0], 1, &STx[0], 1);
2301 
2302  // v * Sxy
2303  Vmath::Vmul(nPts, &physfield_interp[1][0], 1,
2304  &Sxy[0], 1, &tmp1[0], 1);
2305 
2306  // v * Sxy
2307  Vmath::Vmul(nPts, &physfield_interp[2][0], 1,
2308  &Sxz[0], 1, &tmp2[0], 1);
2309 
2310  // k * dT/dx
2311  Vmath::Smul(nPts, m_thermalConductivity,
2312  &derivativesO1_interp[0][3][0], 1, &tmp3[0], 1);
2313 
2314  // STx = u * Sxx + v * Sxy + w * Sxz + (K / mu) * dT/dx
2315  Vmath::Vadd(nPts, &STx[0], 1, &tmp1[0], 1, &STx[0], 1);
2316  Vmath::Vadd(nPts, &STx[0], 1, &tmp2[0], 1, &STx[0], 1);
2317  Vmath::Vadd(nPts, &STx[0], 1, &tmp3[0], 1, &STx[0], 1);
2318 
2319  // Computation of STy
2320 
2321  // Re-initialise temporary arrays
2322  Vmath::Zero(nPts, &tmp1[0], 1);
2323  Vmath::Zero(nPts, &tmp2[0], 1);
2324  Vmath::Zero(nPts, &tmp3[0], 1);
2325 
2326  // v * Syy
2327  Vmath::Vmul(nPts, &physfield_interp[1][0], 1,
2328  &Sgg[1][0], 1, &STy[0], 1);
2329 
2330  // u * Sxy
2331  Vmath::Vmul(nPts, &physfield_interp[0][0], 1,
2332  &Sxy[0], 1, &tmp1[0], 1);
2333 
2334  // w * Syz
2335  Vmath::Vmul(nPts, &physfield_interp[2][0], 1,
2336  &Syz[0], 1, &tmp2[0], 1);
2337 
2338  // k * dT/dy
2339  Vmath::Smul(nPts, m_thermalConductivity,
2340  &derivativesO1_interp[1][3][0], 1, &tmp3[0], 1);
2341 
2342  // STy = v * Syy + u * Sxy + w * Syz + K * dT/dy
2343  Vmath::Vadd(nPts, &STy[0], 1, &tmp1[0], 1, &STy[0], 1);
2344  Vmath::Vadd(nPts, &STy[0], 1, &tmp2[0], 1, &STy[0], 1);
2345  Vmath::Vadd(nPts, &STy[0], 1, &tmp3[0], 1, &STy[0], 1);
2346 
2347  // Computation of STz
2348 
2349  // Re-initialise temporary arrays
2350  Vmath::Zero(nPts, &tmp1[0], 1);
2351  Vmath::Zero(nPts, &tmp2[0], 1);
2352  Vmath::Zero(nPts, &tmp3[0], 1);
2353 
2354  // w * Szz
2355  Vmath::Vmul(nPts, &physfield_interp[2][0], 1,
2356  &Sgg[2][0], 1, &STz[0], 1);
2357 
2358  // u * Sxz
2359  Vmath::Vmul(nPts, &physfield_interp[0][0], 1,
2360  &Sxz[0], 1, &tmp1[0], 1);
2361 
2362  // v * Syz
2363  Vmath::Vmul(nPts, &physfield_interp[1][0], 1,
2364  &Syz[0], 1, &tmp2[0], 1);
2365 
2366  // k * dT/dz
2367  Vmath::Smul(nPts, m_thermalConductivity,
2368  &derivativesO1_interp[2][3][0], 1, &tmp3[0], 1);
2369 
2370  // STz = w * Szz + u * Sxz + v * Syz + K * dT/dz
2371  Vmath::Vadd(nPts, &STz[0], 1, &tmp1[0], 1, &STz[0], 1);
2372  Vmath::Vadd(nPts, &STz[0], 1, &tmp2[0], 1, &STz[0], 1);
2373  Vmath::Vadd(nPts, &STz[0], 1, &tmp3[0], 1, &STz[0], 1);
2374  }
2375 
2376  switch (m_spacedim)
2377  {
2378  case 1:
2379  {
2380 
2381  int nq = physfield[0].num_elements();
2382  // f_11v = f_rho = 0
2383  Vmath::Zero(nq, &viscousTensor_interp[0][0][0], 1);
2384 
2385  // f_21v = f_rhou
2386  Vmath::Vcopy(nPts, &Sgg[0][0], 1, &viscousTensor_interp[0][1][0], 1);
2387 
2388  // f_31v = f_E
2389  Vmath::Vcopy(nPts, &STx[0], 1, &viscousTensor_interp[0][2][0], 1);
2390  break;
2391  }
2392  case 2:
2393  {
2394  int nq = physfield[0].num_elements();
2395  // f_11v = f_rho1 = 0
2396  Vmath::Zero(nq, &viscousTensor_interp[0][0][0], 1);
2397  // f_12v = f_rho2 = 0
2398  Vmath::Zero(nq, &viscousTensor_interp[1][0][0], 1);
2399 
2400  // f_21v = f_rhou1
2401  Vmath::Vcopy(nPts, &Sgg[0][0], 1, &viscousTensor_interp[0][1][0], 1);
2402  // f_22v = f_rhou2
2403  Vmath::Vcopy(nPts, &Sxy[0], 1, &viscousTensor_interp[1][1][0], 1);
2404 
2405  // f_31v = f_rhov1
2406  Vmath::Vcopy(nPts, &Sxy[0], 1, &viscousTensor_interp[0][2][0], 1);
2407  // f_32v = f_rhov2
2408  Vmath::Vcopy(nPts, &Sgg[1][0], 1, &viscousTensor_interp[1][2][0], 1);
2409 
2410  // f_41v = f_E1
2411  Vmath::Vcopy(nPts, &STx[0], 1, &viscousTensor_interp[0][3][0], 1);
2412  // f_42v = f_E2
2413  Vmath::Vcopy(nPts, &STy[0], 1, &viscousTensor_interp[1][3][0], 1);
2414  break;
2415  }
2416  case 3:
2417  {
2418  int nq = physfield[0].num_elements();
2419  // f_11v = f_rho1 = 0
2420  Vmath::Zero(nq, &viscousTensor_interp[0][0][0], 1);
2421  // f_12v = f_rho2 = 0
2422  Vmath::Zero(nq, &viscousTensor_interp[1][0][0], 1);
2423  // f_13v = f_rho3 = 0
2424  Vmath::Zero(nq, &viscousTensor_interp[2][0][0], 1);
2425 
2426  // f_21v = f_rhou1
2427  Vmath::Vcopy(nPts, &Sgg[0][0], 1, &viscousTensor_interp[0][1][0], 1);
2428  // f_22v = f_rhou2
2429  Vmath::Vcopy(nPts, &Sxy[0], 1, &viscousTensor_interp[1][1][0], 1);
2430  // f_23v = f_rhou3
2431  Vmath::Vcopy(nPts, &Sxz[0], 1, &viscousTensor_interp[2][1][0], 1);
2432 
2433  // f_31v = f_rhov1
2434  Vmath::Vcopy(nPts, &Sxy[0], 1, &viscousTensor_interp[0][2][0], 1);
2435  // f_32v = f_rhov2
2436  Vmath::Vcopy(nPts, &Sgg[1][0], 1, &viscousTensor_interp[1][2][0], 1);
2437  // f_33v = f_rhov3
2438  Vmath::Vcopy(nPts, &Syz[0], 1, &viscousTensor_interp[2][2][0], 1);
2439 
2440  // f_31v = f_rhow1
2441  Vmath::Vcopy(nPts, &Sxz[0], 1, &viscousTensor_interp[0][3][0], 1);
2442  // f_32v = f_rhow2
2443  Vmath::Vcopy(nPts, &Syz[0], 1, &viscousTensor_interp[1][3][0], 1);
2444  // f_33v = f_rhow3
2445  Vmath::Vcopy(nPts, &Sgg[2][0], 1, &viscousTensor_interp[2][3][0], 1);
2446 
2447  // f_41v = f_E1
2448  Vmath::Vcopy(nPts, &STx[0], 1, &viscousTensor_interp[0][4][0], 1);
2449  // f_42v = f_E2
2450  Vmath::Vcopy(nPts, &STy[0], 1, &viscousTensor_interp[1][4][0], 1);
2451  // f_43v = f_E3
2452  Vmath::Vcopy(nPts, &STz[0], 1, &viscousTensor_interp[2][4][0], 1);
2453  break;
2454  }
2455  default:
2456  {
2457  ASSERTL0(false, "Illegal expansion dimension");
2458  }
2459  }
2460 
2461  for (i = 0; i < m_spacedim; ++i)
2462  {
2463  for (j = 1; j < nVariables; ++j)
2464  {
2465  // Galerkin project solution back to origianl space
2466  m_fields[0]->PhysGalerkinProjection1DScaled(
2467  OneDptscale,
2468  viscousTensor_interp[i][j],
2469  viscousTensor[i][j]);
2470  }
2471  }
2472 #endif
2473 }
2474 
2475  /**
2476  * @brief Calculate the pressure field \f$ p =
2477  * (\gamma-1)(E-\frac{1}{2}\rho\| \mathbf{v} \|^2) \f$ assuming an ideal
2478  * gas law.
2479  *
2480  * @param physfield Input momentum.
2481  * @param pressure Computed pressure field.
2482  */
2483  void CompressibleFlowSystem::GetPressure(
2484  const Array<OneD, const Array<OneD, NekDouble> > &physfield,
2485  Array<OneD, NekDouble> &pressure)
2486  {
2487  int nBCEdgePts = physfield[0].num_elements();
2488  NekDouble alpha = -0.5;
2489 
2490  // Calculate ||rho v||^2
2491  Vmath::Vmul(nBCEdgePts, physfield[1], 1, physfield[1], 1, pressure, 1);
2492  for (int i = 1; i < m_spacedim; ++i)
2493  {
2494  Vmath::Vvtvp(nBCEdgePts, physfield[1+i], 1, physfield[1+i], 1,
2495  pressure, 1, pressure, 1);
2496  }
2497  // Divide by rho to get rho*||v||^2
2498  Vmath::Vdiv (nBCEdgePts, pressure, 1, physfield[0], 1, pressure, 1);
2499  // pressure <- E - 0.5*pressure
2500  Vmath::Svtvp(nBCEdgePts, alpha,
2501  pressure, 1, physfield[m_spacedim+1], 1, pressure, 1);
2502  // Multiply by (gamma-1)
2503  Vmath::Smul (nBCEdgePts, m_gamma-1, pressure, 1, pressure, 1);
2504  }
2505 
2506  /**
2507  * @brief Compute the enthalpy term \f$ H = E + p/rho \$.
2508  */
2509  void CompressibleFlowSystem::GetEnthalpy(
2510  const Array<OneD, const Array<OneD, NekDouble> > &physfield,
2511  Array<OneD, NekDouble> &pressure,
2512  Array<OneD, NekDouble> &enthalpy)
2513  {
2514  int npts = m_fields[0]->GetTotPoints();
2515  Array<OneD, NekDouble> tmp(npts, 0.0);
2516 
2517  // Calculate E = rhoE/rho
2518  Vmath::Vdiv(npts, physfield[m_spacedim+1], 1, physfield[0], 1, tmp, 1);
2519  // Calculate p/rho
2520  Vmath::Vdiv(npts, pressure, 1, physfield[0], 1, enthalpy, 1);
2521  // Calculate H = E + p/rho
2522  Vmath::Vadd(npts, tmp, 1, enthalpy, 1, enthalpy, 1);
2523  }
2524 
2525  /**
2526  * @brief Calculate the pressure field \f$ p =
2527  * (\gamma-1)(E-\frac{1}{2}\rho\| \mathbf{v} \|^2) \f$ assuming an ideal
2528  * gas law.
2529  *
2530  * This is a slightly optimised way to calculate the pressure field which
2531  * avoids division by the density field if the velocity field has already
2532  * been calculated.
2533  *
2534  * @param physfield Input momentum.
2535  * @param velocity Velocity vector.
2536  * @param pressure Computed pressure field.
2537  */
2538  void CompressibleFlowSystem::GetPressure(
2539  const Array<OneD, const Array<OneD, NekDouble> > &physfield,
2540  const Array<OneD, const Array<OneD, NekDouble> > &velocity,
2541  Array<OneD, NekDouble> &pressure)
2542  {
2543  int nBCEdgePts = physfield[0].num_elements();
2544  NekDouble alpha = -0.5;
2545 
2546  // Calculate ||\rho v||^2.
2547  Vmath::Vmul (nBCEdgePts, velocity[0], 1, physfield[1], 1, pressure, 1);
2548  for (int i = 1; i < m_spacedim; ++i)
2549  {
2550  Vmath::Vvtvp(nBCEdgePts, velocity[i], 1, physfield[1+i], 1,
2551  pressure, 1, pressure, 1);
2552  }
2553  // pressure <- E - 0.5*pressure
2554  Vmath::Svtvp(nBCEdgePts, alpha,
2555  pressure, 1, physfield[m_spacedim+1], 1, pressure, 1);
2556  // Multiply by (gamma-1).
2557  Vmath::Smul (nBCEdgePts, m_gamma-1, pressure, 1, pressure, 1);
2558  }
2559 
2560  /**
2561  * @brief Compute the velocity field \f$ \mathbf{v} \f$ given the momentum
2562  * \f$ \rho\mathbf{v} \f$.
2563  *
2564  * @param physfield Momentum field.
2565  * @param velocity Velocity field.
2566  */
2567  void CompressibleFlowSystem::GetVelocityVector(
2568  const Array<OneD, Array<OneD, NekDouble> > &physfield,
2569  Array<OneD, Array<OneD, NekDouble> > &velocity)
2570  {
2571  const int nBCEdgePts = physfield[0].num_elements();
2572 
2573  for (int i = 0; i < m_spacedim; ++i)
2574  {
2575  Vmath::Vdiv(nBCEdgePts, physfield[1+i], 1, physfield[0], 1,
2576  velocity[i], 1);
2577  }
2578  }
2579 
2580  /**
2581  * @brief Compute the temperature \f$ T = p/\rho R \f$.
2582  *
2583  * @param physfield Input physical field.
2584  * @param pressure The pressure field corresponding to physfield.
2585  * @param temperature The resulting temperature \f$ T \f$.
2586  */
2587  void CompressibleFlowSystem::GetTemperature(
2588  const Array<OneD, const Array<OneD, NekDouble> > &physfield,
2589  Array<OneD, NekDouble > &pressure,
2590  Array<OneD, NekDouble > &temperature)
2591  {
2592  const int nq = physfield[0].num_elements();
2593 
2594  Vmath::Vdiv(nq, pressure, 1, physfield[0], 1, temperature, 1);
2595  Vmath::Smul(nq, 1.0/m_gasConstant, temperature, 1, temperature, 1);
2596  }
2597 
2598  /**
2599  * @brief Compute the sound speed \f$ c = sqrt(\gamma p/\rho) \f$.
2600  *
2601  * @param physfield Input physical field.
2602  * @param pressure The pressure field corresponding to physfield.
2603  * @param soundspeed The resulting sound speed \f$ c \f$.
2604  */
2605  void CompressibleFlowSystem::GetSoundSpeed(
2606  const Array<OneD, Array<OneD, NekDouble> > &physfield,
2607  Array<OneD, NekDouble > &pressure,
2608  Array<OneD, NekDouble > &soundspeed)
2609  {
2610  const int nq = m_fields[0]->GetTotPoints();
2611  Vmath::Vdiv (nq, pressure, 1, physfield[0], 1, soundspeed, 1);
2612  Vmath::Smul (nq, m_gamma, soundspeed, 1, soundspeed, 1);
2613  Vmath::Vsqrt(nq, soundspeed, 1, soundspeed, 1);
2614  }
2615 
2616  /**
2617  * @brief Compute the mach number \f$ M = \| \mathbf{v} \|^2 / c \f$.
2618  *
2619  * @param physfield Input physical field.
2620  * @param soundfield The speed of sound corresponding to physfield.
2621  * @param mach The resulting mach number \f$ M \f$.
2622  */
2623  void CompressibleFlowSystem::GetMach(
2624  Array<OneD, Array<OneD, NekDouble> > &physfield,
2625  Array<OneD, NekDouble > &soundspeed,
2626  Array<OneD, NekDouble > &mach)
2627  {
2628  const int nq = m_fields[0]->GetTotPoints();
2629 
2630  Vmath::Vmul(nq, physfield[1], 1, physfield[1], 1, mach, 1);
2631 
2632  for (int i = 1; i < m_spacedim; ++i)
2633  {
2634  Vmath::Vvtvp(nq, physfield[1+i], 1, physfield[1+i], 1,
2635  mach, 1, mach, 1);
2636  }
2637 
2638  Vmath::Vdiv(nq, mach, 1, physfield[0], 1, mach, 1);
2639  Vmath::Vdiv(nq, mach, 1, physfield[0], 1, mach, 1);
2640  Vmath::Vsqrt(nq, mach, 1, mach, 1);
2641 
2642  Vmath::Vdiv(nq, mach, 1, soundspeed, 1, mach, 1);
2643  }
2644 
2645  /**
2646  * @brief Compute the dynamic viscosity using the Sutherland's law
2647  * \f$ \mu = \mu_star * (T / T_star)^3/2 * (T_star + 110) / (T + 110) \f$,
2648  * where: \mu_star = 1.7894 * 10^-5 Kg / (m * s)
2649  * T_star = 288.15 K
2650  *
2651  * @param physfield Input physical field.
2652  * @param mu The resulting dynamic viscosity.
2653  */
2654  void CompressibleFlowSystem::GetDynamicViscosity(
2655  const Array<OneD, const NekDouble> &temperature,
2656  Array<OneD, NekDouble> &mu)
2657  {
2658  const int nPts = temperature.num_elements();
2659  NekDouble mu_star = m_mu;
2660  NekDouble T_star = m_pInf / (m_rhoInf * m_gasConstant);
2661  NekDouble ratio;
2662 
2663  for (int i = 0; i < nPts; ++i)
2664  {
2665  ratio = temperature[i] / T_star;
2666  mu[i] = mu_star * ratio * sqrt(ratio) *
2667  (T_star + 110.0) / (temperature[i] + 110.0);
2668  }
2669  }
2670 
2671  /**
2672  * @brief Perform post-integration checks, presently just to check steady
2673  * state behaviour.
2674  */
2675  bool CompressibleFlowSystem::v_PostIntegrate(int step)
2676  {
2677  if (m_steadyStateTol > 0.0)
2678  {
2679  bool doOutput = step % m_infosteps == 0;
2680  if (CalcSteadyState(doOutput))
2681  {
2682  if (m_comm->GetRank() == 0)
2683  {
2684  cout << "Reached Steady State to tolerance "
2685  << m_steadyStateTol << endl;
2686  }
2687  return true;
2688  }
2689  }
2690  return false;
2691  }
2692 
2693  /**
2694  * @brief Calculate whether the system has reached a steady state by
2695  * observing residuals to a user-defined tolerance.
2696  */
2697  bool CompressibleFlowSystem::CalcSteadyState(bool output)
2698  {
2699  const int nPoints = GetTotPoints();
2700  const int nFields = m_fields.num_elements();
2701 
2702  // Holds L2 errors.
2703  Array<OneD, NekDouble> L2 (nFields);
2704  Array<OneD, NekDouble> residual(nFields);
2705 
2706  for (int i = 0; i < nFields; ++i)
2707  {
2708  Array<OneD, NekDouble> diff(nPoints);
2709 
2710  Vmath::Vsub(nPoints, m_fields[i]->GetPhys(), 1, m_un[i], 1, diff, 1);
2711  Vmath::Vmul(nPoints, diff, 1, diff, 1, diff, 1);
2712  residual[i] = Vmath::Vsum(nPoints, diff, 1);
2713  }
2714 
2715  m_comm->AllReduce(residual, LibUtilities::ReduceSum);
2716 
2717  // L2 error
2718  L2[0] = sqrt(residual[0]) / m_rhoInf;
2719 
2720  for (int i = 1; i < nFields-1; ++i)
2721  {
2722  L2[i] = sqrt(residual[i]) / m_UInf / m_rhoInf;
2723  }
2724 
2725  NekDouble Einf = m_pInf / (m_gamma-1.0) + 0.5 * m_rhoInf * m_UInf;
2726  L2[nFields-1] = sqrt(residual[nFields-1]) / Einf;
2727 
2728  if (m_comm->GetRank() == 0 && output)
2729  {
2730  // Output time
2731  m_errFile << setprecision(8) << setw(17) << scientific << m_time;
2732 
2733  // Output residuals
2734  for (int i = 0; i < nFields; ++i)
2735  {
2736  m_errFile << setprecision(11) << setw(22) << scientific
2737  << L2[i];
2738  }
2739 
2740  m_errFile << endl;
2741  }
2742 
2743  // Calculate maximum L2 error
2744  NekDouble maxL2 = Vmath::Vmax(nFields, L2, 1);
2745 
2746  if (m_session->DefinesCmdLineArgument("verbose") &&
2747  m_comm->GetRank() == 0 && output)
2748  {
2749  cout << "-- Maximum L^2 residual: " << maxL2 << endl;
2750  }
2751 
2752  if (maxL2 <= m_steadyStateTol)
2753  {
2754  return true;
2755  }
2756 
2757  for (int i = 0; i < m_fields.num_elements(); ++i)
2758  {
2759  Vmath::Vcopy(nPoints, m_fields[i]->GetPhys(), 1, m_un[i], 1);
2760  }
2761 
2762  return false;
2763  }
2764 
2765  /**
2766  * @brief Calculate entropy.
2767  */
2768  void CompressibleFlowSystem::GetEntropy(
2769  const Array<OneD, const Array<OneD, NekDouble> > &physfield,
2770  const Array<OneD, const NekDouble> &pressure,
2771  const Array<OneD, const NekDouble> &temperature,
2772  Array<OneD, NekDouble> &entropy)
2773  {
2774  NekDouble entropy_sum = 0.0;
2775  const int npts = m_fields[0]->GetTotPoints();
2776  const NekDouble temp_inf = m_pInf/(m_rhoInf*m_gasConstant);;
2777  Array<OneD, NekDouble> L2entropy(npts, 0.0);
2778 
2779  for (int i = 0; i < npts; ++i)
2780  {
2781  entropy[i] = m_gamma / (m_gamma - 1.0) * m_gasConstant *
2782  log(temperature[i]/temp_inf) - m_gasConstant *
2783  log(pressure[i] / m_pInf);
2784  }
2785 
2786  Vmath::Vmul(npts, entropy, 1, entropy, 1, L2entropy, 1);
2787 
2788  entropy_sum = Vmath::Vsum(npts, L2entropy, 1);
2789 
2790  entropy_sum = sqrt(entropy_sum);
2791 
2792  std::ofstream m_file( "L2entropy.txt", std::ios_base::app);
2793 
2794  m_file << setprecision(16) << scientific << entropy_sum << endl;
2795  //m_file << Vmath::Vmax(entropy.num_elements(),entropy,1) << endl;
2796 
2797  m_file.close();
2798  }
2799 
2800  /**
2801  * @brief Calculate the maximum timestep subject to CFL restrictions.
2802  */
2803  NekDouble CompressibleFlowSystem::v_GetTimeStep(
2804  const Array<OneD, const Array<OneD, NekDouble> > &inarray)
2805  {
2806  int n;
2807  int nElements = m_fields[0]->GetExpSize();
2808  const Array<OneD, int> ExpOrder = GetNumExpModesPerExp();
2809 
2810  Array<OneD, NekDouble> tstep (nElements, 0.0);
2811  Array<OneD, NekDouble> stdVelocity(nElements);
2812 
2813  // Get standard velocity to compute the time-step limit
2814  GetStdVelocity(inarray, stdVelocity);
2815 
2816  // Factors to compute the time-step limit
2817  NekDouble minLength = 0.0;
2818  NekDouble alpha = MaxTimeStepEstimator();
2819  NekDouble cLambda = 0.2; // Spencer book-317
2820 
2821  // Loop over elements to compute the time-step limit for each element
2822  for(n = 0; n < nElements; ++n)
2823  {
2824  int npoints = m_fields[0]->GetExp(n)->GetTotPoints();
2825  Array<OneD, NekDouble> one2D(npoints, 1.0);
2826  NekDouble Area = m_fields[0]->GetExp(n)->Integral(one2D);
2827 
2828  minLength = sqrt(Area);
2829  if (m_fields[0]->GetExp(n)->as<LocalRegions::TriExp>())
2830  {
2831  minLength *= 2.0;
2832  }
2833 
2834  tstep[n] = m_cflSafetyFactor * alpha * minLength
2835  / (stdVelocity[n] * cLambda
2836  * (ExpOrder[n] - 1) * (ExpOrder[n] - 1));
2837  }
2838 
2839  // Get the minimum time-step limit and return the time-step
2840  NekDouble TimeStep = Vmath::Vmin(nElements, tstep, 1);
2841  m_comm->AllReduce(TimeStep, LibUtilities::ReduceMin);
2842  return TimeStep;
2843  }
2844 
2845  /**
2846  * @brief Set up logic for residual calculation.
2847  */
2848  void CompressibleFlowSystem::v_SetInitialConditions(
2849  NekDouble initialtime,
2850  bool dumpInitialConditions,
2851  const int domain)
2852  {
2853  if (m_session->DefinesParameter("SteadyStateTol"))
2854  {
2855  const int nPoints = m_fields[0]->GetTotPoints();
2856  m_un = Array<OneD, Array<OneD, NekDouble> > (
2857  m_fields.num_elements());
2858 
2859  for (int i = 0; i < m_fields.num_elements(); ++i)
2860  {
2861  cout << Vmath::Vmax(nPoints, m_fields[i]->GetPhys(), 1) << endl;
2862  m_un[i] = Array<OneD, NekDouble>(nPoints);
2863  Vmath::Vcopy(nPoints, m_fields[i]->GetPhys(), 1, m_un[i], 1);
2864  }
2865 
2866  if (m_comm->GetRank() == 0)
2867  {
2868  std::string fName = m_session->GetSessionName() +
2869  std::string(".res");
2870  m_errFile.open(fName.c_str());
2871  m_errFile << "# "
2872  << setw(15) << left << "Time"
2873  << setw(22) << left << "rho";
2874 
2875  std::string velFields[3] = {"u", "v", "w"};
2876 
2877  for (int i = 0; i < m_fields.num_elements()-2; ++i)
2878  {
2879  m_errFile << setw(22) << "rho"+velFields[i];
2880  }
2881 
2882  m_errFile << setw(22) << left << "E" << endl;
2883  }
2884  }
2885  }
2886 
2887 
2888  /**
2889  * @brief Compute the advection velocity in the standard space
2890  * for each element of the expansion.
2891  *
2892  * @param inarray Momentum field.
2893  * @param stdV Standard velocity field.
2894  */
2895  void CompressibleFlowSystem::GetStdVelocity(
2896  const Array<OneD, const Array<OneD, NekDouble> > &inarray,
2897  Array<OneD, NekDouble> &stdV)
2898  {
2899  int nTotQuadPoints = GetTotPoints();
2900  int n_element = m_fields[0]->GetExpSize();
2901  int nBCEdgePts = 0;
2902 
2903  // Getting the velocity vector on the 2D normal space
2904  Array<OneD, Array<OneD, NekDouble> > velocity (m_spacedim);
2905  Array<OneD, Array<OneD, NekDouble> > stdVelocity(m_spacedim);
2906  Array<OneD, NekDouble> pressure (nTotQuadPoints);
2907  Array<OneD, NekDouble> soundspeed (nTotQuadPoints);
2908  LibUtilities::PointsKeyVector ptsKeys;
2909 
2910  // Zero output array
2911  Vmath::Zero(stdV.num_elements(), stdV, 1);
2912 
2913  for (int i = 0; i < m_spacedim; ++i)
2914  {
2915  velocity [i] = Array<OneD, NekDouble>(nTotQuadPoints);
2916  stdVelocity[i] = Array<OneD, NekDouble>(nTotQuadPoints, 0.0);
2917  }
2918 
2919  GetVelocityVector(inarray, velocity);
2920  GetPressure (inarray, velocity, pressure);
2921  GetSoundSpeed (inarray, pressure, soundspeed);
2922 
2923  for(int el = 0; el < n_element; ++el)
2924  {
2925  ptsKeys = m_fields[0]->GetExp(el)->GetPointsKeys();
2926 
2927  // Possible bug: not multiply by jacobian??
2928  const SpatialDomains::GeomFactorsSharedPtr metricInfo =
2929  m_fields[0]->GetExp(el)->GetGeom()->GetMetricInfo();
2930  const Array<TwoD, const NekDouble> &gmat =
2931  m_fields[0]->GetExp(el)->GetGeom()->GetMetricInfo()
2932  ->GetDerivFactors(ptsKeys);
2933 
2934  int nq = m_fields[0]->GetExp(el)->GetTotPoints();
2935 
2936  if(metricInfo->GetGtype() == SpatialDomains::eDeformed)
2937  {
2938  // d xi/ dx = gmat = 1/J * d x/d xi
2939  for (int i = 0; i < m_spacedim; ++i)
2940  {
2941  Vmath::Vmul(nq, gmat[i], 1, velocity[0], 1,
2942  stdVelocity[i], 1);
2943  for (int j = 1; j < m_spacedim; ++j)
2944  {
2945  Vmath::Vvtvp(nq, gmat[m_spacedim*j+i], 1, velocity[j],
2946  1, stdVelocity[i], 1, stdVelocity[i], 1);
2947  }
2948  }
2949  }
2950  else
2951  {
2952  for (int i = 0; i < m_spacedim; ++i)
2953  {
2954  Vmath::Smul(nq, gmat[i][0], velocity[0], 1,
2955  stdVelocity[i], 1);
2956  for (int j = 1; j < m_spacedim; ++j)
2957  {
2958  Vmath::Svtvp(nq, gmat[m_spacedim*j+i][0], velocity[j],
2959  1, stdVelocity[i], 1, stdVelocity[i], 1);
2960  }
2961  }
2962  }
2963 
2964  for (int i = 0; i < nq; ++i)
2965  {
2966  NekDouble pntVelocity = 0.0;
2967  for (int j = 0; j < m_spacedim; ++j)
2968  {
2969  pntVelocity += stdVelocity[j][i]*stdVelocity[j][i];
2970  }
2971  pntVelocity = sqrt(pntVelocity) + soundspeed[nBCEdgePts];
2972  if (pntVelocity > stdV[el])
2973  {
2974  stdV[el] = pntVelocity;
2975  }
2976  nBCEdgePts++;
2977  }
2978  }
2979  }
2980 
2981  /**
2982  * @brief Set the denominator to compute the time step when a cfl
2983  * control is employed. This function is no longer used but is still
2984  * here for being utilised in the future.
2985  *
2986  * @param n Order of expansion element by element.
2987  */
2988  NekDouble CompressibleFlowSystem::GetStabilityLimit(int n)
2989  {
2990  ASSERTL0(n <= 20, "Illegal modes dimension for CFL calculation "
2991  "(P has to be less then 20)");
2992 
2993  NekDouble CFLDG[21] = { 2.0000, 6.0000, 11.8424, 19.1569,
2994  27.8419, 37.8247, 49.0518, 61.4815,
2995  75.0797, 89.8181, 105.6700, 122.6200,
2996  140.6400, 159.7300, 179.8500, 201.0100,
2997  223.1800, 246.3600, 270.5300, 295.6900,
2998  321.8300}; //CFLDG 1D [0-20]
2999  NekDouble CFL = 0.0;
3000 
3001  if (m_projectionType == MultiRegions::eDiscontinuous)
3002  {
3003  CFL = CFLDG[n];
3004  }
3005  else
3006  {
3007  ASSERTL0(false, "Continuous Galerkin stability coefficients "
3008  "not introduced yet.");
3009  }
3010 
3011  return CFL;
3012  }
3013 
3014  /**
3015  * @brief Compute the vector of denominators to compute the time step
3016  * when a cfl control is employed. This function is no longer used but
3017  * is still here for being utilised in the future.
3018  *
3019  * @param ExpOrder Order of expansion element by element.
3020  */
3021  Array<OneD, NekDouble> CompressibleFlowSystem::GetStabilityLimitVector(
3022  const Array<OneD,int> &ExpOrder)
3023  {
3024  int i;
3025  Array<OneD,NekDouble> returnval(m_fields[0]->GetExpSize(), 0.0);
3026  for (i =0; i<m_fields[0]->GetExpSize(); i++)
3027  {
3028  returnval[i] = GetStabilityLimit(ExpOrder[i]);
3029  }
3030  return returnval;
3031  }
3032 
3033  void CompressibleFlowSystem::GetSensor(
3034  const Array<OneD, const Array<OneD, NekDouble> > &physarray,
3035  Array<OneD, NekDouble> &Sensor,
3036  Array<OneD, NekDouble> &SensorKappa)
3037  {
3038  int e, NumModesElement, nQuadPointsElement;
3039  int nTotQuadPoints = GetTotPoints();
3040  int nElements = m_fields[0]->GetExpSize();
3041 
3042  // Find solution (SolP) at p = P;
3043  // The input array (physarray) is the solution at p = P;
3044 
3045  Array<OneD,int> ExpOrderElement = GetNumExpModesPerExp();
3046 
3047  Array<OneD, NekDouble> SolP (nTotQuadPoints, 0.0);
3048  Array<OneD, NekDouble> SolPmOne(nTotQuadPoints, 0.0);
3049  Array<OneD, NekDouble> SolNorm (nTotQuadPoints, 0.0);
3050 
3051  Vmath::Vcopy(nTotQuadPoints, physarray[0], 1, SolP, 1);
3052 
3053  int CoeffsCount = 0;
3054 
3055  for (e = 0; e < nElements; e++)
3056  {
3057  NumModesElement = ExpOrderElement[e];
3058  int nQuadPointsElement = m_fields[0]->GetExp(e)->GetTotPoints();
3059  int nCoeffsElement = m_fields[0]->GetExp(e)->GetNcoeffs();
3060  int numCutOff = NumModesElement - 1;
3061 
3062  // Set-up of the Orthogonal basis for a Quadrilateral element which
3063  // is needed to obtain thesolution at P = p - 1;
3064 
3065  Array<OneD, NekDouble> SolPElementPhys (nQuadPointsElement, 0.0);
3066  Array<OneD, NekDouble> SolPElementCoeffs(nCoeffsElement, 0.0);
3067 
3068  Array<OneD, NekDouble> SolPmOneElementPhys(nQuadPointsElement, 0.0);
3069  Array<OneD, NekDouble> SolPmOneElementCoeffs(nCoeffsElement, 0.0);
3070 
3071  // create vector the save the solution points per element at P = p;
3072 
3073  for (int i = 0; i < nQuadPointsElement; i++)
3074  {
3075  SolPElementPhys[i] = SolP[CoeffsCount+i];
3076  }
3077 
3078  m_fields[0]->GetExp(e)->FwdTrans(SolPElementPhys,
3079  SolPElementCoeffs);
3080 
3081  // ReduceOrderCoeffs reduces the polynomial order of the solution
3082  // that is represented by the coeffs given as an inarray. This is
3083  // done by projecting the higher order solution onto the orthogonal
3084  // basis and padding the higher order coefficients with zeros.
3085 
3086  m_fields[0]->GetExp(e)->ReduceOrderCoeffs(numCutOff,
3087  SolPElementCoeffs,
3088  SolPmOneElementCoeffs);
3089 
3090  m_fields[0]->GetExp(e)->BwdTrans(SolPmOneElementCoeffs,
3091  SolPmOneElementPhys);
3092 
3093  for (int i = 0; i < nQuadPointsElement; i++)
3094  {
3095  SolPmOne[CoeffsCount+i] = SolPmOneElementPhys[i];
3096  }
3097 
3098  NekDouble SolPmeanNumerator = 0.0;
3099  NekDouble SolPmeanDenumerator = 0.0;
3100 
3101  // Determining the norm of the numerator of the Sensor
3102 
3103  Vmath::Vsub(nQuadPointsElement,
3104  SolPElementPhys, 1,
3105  SolPmOneElementPhys, 1,
3106  SolNorm, 1);
3107 
3108  Vmath::Vmul(nQuadPointsElement,
3109  SolNorm, 1,
3110  SolNorm, 1,
3111  SolNorm, 1);
3112 
3113  for (int i = 0; i < nQuadPointsElement; i++)
3114  {
3115  SolPmeanNumerator += SolNorm[i];
3116  SolPmeanDenumerator += SolPElementPhys[i];
3117  }
3118 
3119  for (int i = 0; i < nQuadPointsElement; ++i)
3120  {
3121  Sensor[CoeffsCount+i] =
3122  sqrt(SolPmeanNumerator / nQuadPointsElement) /
3123  sqrt(SolPmeanDenumerator / nQuadPointsElement);
3124 
3125  Sensor[CoeffsCount+i] = log10(Sensor[CoeffsCount+i]);
3126  }
3127  CoeffsCount += nQuadPointsElement;
3128  }
3129 
3130  CoeffsCount = 0.0;
3131 
3132  for (e = 0; e < nElements; e++)
3133  {
3134  NumModesElement = ExpOrderElement[e];
3135  NekDouble ThetaS = m_mu0;
3136  NekDouble Phi0 = m_Skappa;
3137  NekDouble DeltaPhi = m_Kappa;
3138  nQuadPointsElement = m_fields[0]->GetExp(e)->GetTotPoints();
3139 
3140  for (int i = 0; i < nQuadPointsElement; i++)
3141  {
3142  if (Sensor[CoeffsCount+i] <= (Phi0 - DeltaPhi))
3143  {
3144  SensorKappa[CoeffsCount+i] = 0;
3145  }
3146  else if(Sensor[CoeffsCount+i] >= (Phi0 + DeltaPhi))
3147  {
3148  SensorKappa[CoeffsCount+i] = ThetaS;
3149  }
3150  else if(abs(Sensor[CoeffsCount+i]-Phi0) < DeltaPhi)
3151  {
3152  SensorKappa[CoeffsCount+i] =
3153  ThetaS / 2 * (1 + sin(M_PI * (Sensor[CoeffsCount+i] -
3154  Phi0) / (2 * DeltaPhi)));
3155  }
3156  }
3157 
3158  CoeffsCount += nQuadPointsElement;
3159  }
3160 
3161  }
3162 
3163  void CompressibleFlowSystem::GetForcingTerm(
3164  const Array<OneD, const Array<OneD, NekDouble> > &inarray,
3165  Array<OneD, Array<OneD, NekDouble> > outarrayForcing)
3166  {
3167  const int nPts = m_fields[0]->GetTotPoints();
3168  const int nvariables = m_fields.num_elements();
3169  const int nElements = m_fields[0]->GetExpSize();
3170 
3171  Array<OneD, NekDouble> Sensor(nPts, 0.0);
3172  Array<OneD, NekDouble> SensorKappa(nPts, 0.0);
3173  Array <OneD, NekDouble > Lambda(nPts, 0.0);
3174  Array <OneD, NekDouble > Tau(nPts, 1.0);
3175  Array <OneD, NekDouble > soundspeed(nPts, 0.0);
3176  Array <OneD, NekDouble > pressure(nPts, 0.0);
3177  Array <OneD, NekDouble > temperature(nPts, 0.0);
3178  Array <OneD, NekDouble > absVelocity(nPts, 0.0);
3179  Array <OneD, NekDouble > hel(nPts, 0.0);
3180  Array <OneD, NekDouble > h_minmin(m_spacedim, 0.0);
3181 
3182  Array<OneD,int> pOrderElmt = GetNumExpModesPerExp();
3183  Array<OneD, NekDouble> pOrder (nPts, 0.0);
3184 
3185  // Thermodynamic related quantities
3186  GetPressure(inarray, pressure);
3187  GetTemperature(inarray, pressure, temperature);
3188  GetSoundSpeed(inarray, pressure, soundspeed);
3189  GetAbsoluteVelocity(inarray, absVelocity);
3190  GetSensor(inarray, Sensor, SensorKappa);
3191 
3192  // Determine the maximum wavespeed
3193  Vmath::Vadd(nPts, absVelocity, 1, soundspeed, 1, Lambda, 1);
3194 
3195  NekDouble LambdaMax = Vmath::Vmax(nPts, Lambda, 1);
3196 
3197  int PointCount = 0;
3198 
3199  for (int e = 0; e < nElements; e++)
3200  {
3201  int nQuadPointsElement = m_fields[0]->GetExp(e)->GetTotPoints();
3202 
3203  for (int n = 0; n < nQuadPointsElement; n++)
3204  {
3205  pOrder[n + PointCount] = pOrderElmt[e];
3206 
3207  // order 1.0e-06
3208  Tau[n + PointCount] =
3209  1.0 / (m_C1*pOrder[n + PointCount]*LambdaMax);
3210 
3211  outarrayForcing[nvariables-1][n + PointCount] =
3212  1 / Tau[n + PointCount] * (m_hFactor * LambdaMax /
3213  pOrder[n + PointCount] *
3214  SensorKappa[n + PointCount] -
3215  inarray[nvariables-1][n + PointCount]);
3216  }
3217  PointCount += nQuadPointsElement;
3218  }
3219  }
3220 
3221  void CompressibleFlowSystem::GetElementDimensions(
3222  Array<OneD, Array<OneD, NekDouble> > &outarray,
3223  Array<OneD, NekDouble > &hmin)
3224  {
3225  // So far, this function is only implemented for quads
3226  const int nElements = m_fields[0]->GetExpSize();
3227 
3228  SpatialDomains::HexGeomSharedPtr ElHexGeom;
3229 
3230  NekDouble hx = 0.0;
3231  NekDouble hy = 0.0;
3232 
3233  for (int e = 0; e < nElements; e++)
3234  {
3235  NekDouble nedges = m_fields[0]->GetExp(e)->GetNedges();
3236  Array <OneD, NekDouble> L1(nedges, 0.0);
3237 
3238  for (int j = 0; j < nedges; ++j)
3239  {
3240 
3241  NekDouble x0 = 0.0;
3242  NekDouble y0 = 0.0;
3243  NekDouble z0 = 0.0;
3244 
3245  NekDouble x1 = 0.0;
3246  NekDouble y1 = 0.0;
3247  NekDouble z1 = 0.0;
3248 
3249  if (boost::dynamic_pointer_cast<SpatialDomains::QuadGeom>(
3250  m_fields[0]->GetExp(e)->GetGeom()))
3251  {
3252  SpatialDomains::QuadGeomSharedPtr ElQuadGeom =
3253  boost::dynamic_pointer_cast<SpatialDomains::QuadGeom>(
3254  m_fields[0]->GetExp(e)->GetGeom());
3255 
3256  ElQuadGeom->GetEdge(j)->GetVertex(0)->GetCoords(x0, y0, z0);
3257  ElQuadGeom->GetEdge(j)->GetVertex(1)->GetCoords(x1, y1, z1);
3258 
3259  L1[j] = sqrt(pow((x0-x1),2)+pow((y0-y1),2)+pow((z0-z1),2));
3260  }
3261  else
3262  {
3263  ASSERTL0(false, "GetElementDimensions() is only "
3264  "implemented for quadrilateral elements");
3265  }
3266  }
3267  // determine the minimum length in x and y direction
3268  // still have to find a better estimate when dealing
3269  // with unstructured meshes
3270  if(boost::dynamic_pointer_cast<SpatialDomains::QuadGeom>(
3271  m_fields[0]->GetExp(e)->GetGeom()))
3272  {
3273  hx = min(L1[0], L1[2]);
3274  hy = min(L1[1], L1[3]);
3275 
3276  outarray[0][e] = hx;
3277  outarray[1][e] = hy;
3278  }
3279  }
3280 
3281  hmin[0] = Vmath::Vmin(outarray[0].num_elements(), outarray[0], 1);
3282  hmin[1] = Vmath::Vmin(outarray[1].num_elements(), outarray[1], 1);
3283  }
3284 
3285  void CompressibleFlowSystem::GetAbsoluteVelocity(
3286  const Array<OneD, const Array<OneD, NekDouble> > &inarray,
3287  Array<OneD, NekDouble> &Vtot)
3288  {
3289  int nTotQuadPoints = GetTotPoints();
3290 
3291  // Getting the velocity vector on the 2D normal space
3292  Array<OneD, Array<OneD, NekDouble> > velocity (m_spacedim);
3293 
3294  Vmath::Zero(Vtot.num_elements(), Vtot, 1);
3295 
3296  for (int i = 0; i < m_spacedim; ++i)
3297  {
3298  velocity[i] = Array<OneD, NekDouble>(nTotQuadPoints);
3299  }
3300 
3301  GetVelocityVector(inarray, velocity);
3302 
3303  for (int i = 0; i < m_spacedim; ++i)
3304  {
3305  Vmath::Vvtvp(nTotQuadPoints,
3306  velocity[i], 1,
3307  velocity[i], 1,
3308  Vtot, 1,
3309  Vtot, 1);
3310  }
3311 
3312  Vmath::Vsqrt(Vtot.num_elements(),Vtot,1,Vtot,1);
3313  }
3314 
3315  void CompressibleFlowSystem::GetSmoothArtificialViscosity(
3316  const Array<OneD, Array<OneD, NekDouble> > &physfield,
3317  Array<OneD, NekDouble > &eps_bar)
3318  {
3319  int nvariables = physfield.num_elements();
3320  int nPts = m_fields[0]->GetTotPoints();
3321 
3322  Array<OneD, NekDouble > pressure (nPts, 0.0);
3323  Array<OneD, NekDouble > temperature(nPts, 0.0);
3324  Array<OneD, NekDouble > sensor (nPts, 0.0);
3325  Array<OneD, NekDouble > SensorKappa(nPts, 0.0);
3326  Array<OneD, NekDouble > absVelocity(nPts, 0.0);
3327  Array<OneD, NekDouble > soundspeed (nPts, 0.0);
3328  Array<OneD, NekDouble > Lambda (nPts, 0.0);
3329  Array<OneD, NekDouble > mu_var (nPts, 0.0);
3330  Array<OneD, NekDouble > h_minmin (m_spacedim, 0.0);
3331  Vmath::Zero(nPts, eps_bar, 1);
3332 
3333  // Thermodynamic related quantities
3334  GetPressure(physfield, pressure);
3335  GetTemperature(physfield, pressure, temperature);
3336  GetSoundSpeed(physfield, pressure, soundspeed);
3337  GetAbsoluteVelocity(physfield, absVelocity);
3338  GetSensor(physfield, sensor, SensorKappa);
3339 
3340  // Determine the maximum wavespeed
3341  Vmath::Vadd(nPts, absVelocity, 1, soundspeed, 1, Lambda, 1);
3342 
3343  // Determine hbar = hx_i/h
3344  Array<OneD,int> pOrderElmt = GetNumExpModesPerExp();
3345 
3346  NekDouble ThetaH = m_FacH;
3347  NekDouble ThetaL = m_FacL;
3348 
3349  NekDouble Phi0 = (ThetaH+ThetaL)/2;
3350  NekDouble DeltaPhi = ThetaH-Phi0;
3351 
3352  Vmath::Zero(eps_bar.num_elements(), eps_bar, 1);
3353 
3354  /*Vmath::Smul(eps_bar.num_elements(),
3355  m_eps_max,
3356  &physfield[nvariables-1][0], 1,
3357  &eps_bar[0], 1);*/
3358 
3359  for (int e = 0; e < eps_bar.num_elements(); e++)
3360  {
3361  if (physfield[nvariables-1][e] <= (Phi0 - DeltaPhi))
3362  {
3363  eps_bar[e] = 0;
3364  }
3365  else if(physfield[nvariables-1][e] >= (Phi0 + DeltaPhi))
3366  {
3367  eps_bar[e] = m_mu0;
3368  }
3369  else if(abs(physfield[nvariables-1][e]-Phi0) < DeltaPhi)
3370  {
3371  eps_bar[e] = m_mu0/2*(1+sin(M_PI*
3372  (physfield[nvariables-1][e]-Phi0)/(2*DeltaPhi)));
3373  }
3374  }
3375 
3376  }
3377 
3378  void CompressibleFlowSystem::GetArtificialDynamicViscosity(
3379  const Array<OneD, Array<OneD, NekDouble> > &physfield,
3380  Array<OneD, NekDouble > &mu_var)
3381  {
3382  const int nElements = m_fields[0]->GetExpSize();
3383  int PointCount = 0;
3384  int nTotQuadPoints = GetTotPoints();
3385 
3386  Array<OneD, NekDouble> S_e (nTotQuadPoints, 0.0);
3387  Array<OneD, NekDouble> se (nTotQuadPoints, 0.0);
3388  Array<OneD, NekDouble> Sensor (nTotQuadPoints, 0.0);
3389  Array<OneD, NekDouble> SensorKappa(nTotQuadPoints, 0.0);
3390  Array<OneD, NekDouble> absVelocity(nTotQuadPoints, 0.0);
3391  Array<OneD, NekDouble> soundspeed (nTotQuadPoints, 0.0);
3392  Array<OneD, NekDouble> pressure (nTotQuadPoints, 0.0);
3393 
3394  GetAbsoluteVelocity(physfield, absVelocity);
3395  GetPressure (physfield, pressure);
3396  GetSoundSpeed (physfield, pressure, soundspeed);
3397  GetSensor (physfield, Sensor, SensorKappa);
3398 
3399  Array<OneD, int> pOrderElmt = GetNumExpModesPerExp();
3400  Array<OneD, NekDouble> Lambda(nTotQuadPoints, 1.0);
3401  Vmath::Vadd(nTotQuadPoints, absVelocity, 1, soundspeed, 1, Lambda, 1);
3402 
3403  for (int e = 0; e < nElements; e++)
3404  {
3405  // Threshold value specified in C. Biottos thesis. Based on a 1D
3406  // shock tube problem S_k = log10(1/p^4). See G.E. Barter and
3407  // D.L. Darmofal. Shock Capturing with PDE-based artificial
3408  // diffusion for DGFEM: Part 1 Formulation, Journal of Computational
3409  // Physics 229 (2010) 1810-1827 for further reference
3410 
3411  // Adjustable depending on the coarsness of the mesh. Might want to
3412  // move this variable into the session file
3413 
3414  int nQuadPointsElement = m_fields[0]->GetExp(e)->GetTotPoints();
3415  Array <OneD, NekDouble> one2D(nQuadPointsElement, 1.0);
3416 
3417  for (int n = 0; n < nQuadPointsElement; n++)
3418  {
3419  NekDouble mu_0 = m_mu0;
3420 
3421  if (Sensor[n+PointCount] < (m_Skappa-m_Kappa))
3422  {
3423  mu_var[n+PointCount] = 0;
3424  }
3425  else if (Sensor[n+PointCount] >= (m_Skappa-m_Kappa)
3426  && Sensor[n+PointCount] <= (m_Skappa+m_Kappa))
3427  {
3428  mu_var[n+PointCount] = mu_0 * (0.5 * (1 + sin(
3429  M_PI * (Sensor[n+PointCount] -
3430  m_Skappa - m_Kappa) /
3431  (2*m_Kappa))));
3432  }
3433  else if (Sensor[n+PointCount] > (m_Skappa+m_Kappa))
3434  {
3435  mu_var[n+PointCount] = mu_0;
3436  }
3437  }
3438 
3439  PointCount += nQuadPointsElement;
3440  }
3441  }
3442 
3443  void CompressibleFlowSystem::SetVarPOrderElmt(
3444  const Array<OneD, const Array<OneD, NekDouble> > &physfield,
3445  Array<OneD, NekDouble > &PolyOrder)
3446  {
3447  int e;
3448  NekDouble s_ds, s_sm, s_fl;
3449 
3450  int nElements = m_fields[0]->GetExpSize();
3451  int npts = m_fields[0]->GetTotPoints();
3452 
3453  Array<OneD, NekDouble > Sensor (npts, 0.0);
3454  Array<OneD, NekDouble > SensorKappa(npts, 0.0);
3455  Array<OneD, NekDouble > se (npts, 0.0);
3456 
3457  GetSensor(physfield, Sensor, SensorKappa);
3458 
3459  int nQuadPointsElement = 0;
3460  int npCount = 0;
3461  int MinOrder = 2;
3462  int MaxOrder = 12;
3463  int MinOrderShock = 4;
3464 
3465  std::ofstream m_file( "VariablePComposites.txt", std::ios_base::app);
3466  for (int e = 0; e < nElements; e++)
3467  {
3468  m_file << "<C ID=\"" << e+1 << "\"> Q[" << e << "] </C>"<< endl;
3469  }
3470  m_file.close();
3471 
3472  std::ofstream m_file2( "VariablePExpansions.txt", std::ios_base::app);
3473 
3474  for (e = 0; e < nElements; e++)
3475  {
3476  nQuadPointsElement = m_fields[0]->GetExp(e)->GetTotPoints();
3477 
3478  // Define thresholds
3479  // Ideally, these threshold values could be given
3480  // in the Session File
3481  s_ds = -5.0;
3482  //s_ds = s_0*log10(PolyOrder[e]);
3483  s_sm = -6;
3484  s_fl = -7;
3485 
3486 
3487  for (int i = 0; i < nQuadPointsElement; i++)
3488  {
3489  se[npCount + i] = (Sensor[npCount + i]);
3490 
3491  if (se[npCount + i] > s_ds)
3492  {
3493  if (PolyOrder[npCount + i] > MinOrderShock)
3494  {
3495  PolyOrder[npCount + i] = PolyOrder[npCount + i] - 1;
3496  }
3497  else if(PolyOrder[e] < MinOrderShock)
3498  {
3499  PolyOrder[npCount + i] = PolyOrder[npCount + i] + 1;
3500  }
3501 
3502  }
3503  else if (se[npCount + i] > s_sm && se[npCount + i] < s_ds)
3504  {
3505  if (PolyOrder[npCount + i] < MaxOrder)
3506  {
3507  PolyOrder[npCount + i] = PolyOrder[npCount + i] + 2;
3508  }
3509  }
3510  else if (se[npCount + i] > s_fl && se[npCount + i] < s_sm)
3511  {
3512  PolyOrder[npCount + i] = PolyOrder[npCount + i] + 1;
3513  }
3514  else if (se[npCount + i] < s_fl)
3515  {
3516  if (PolyOrder[npCount + i] > MinOrder)
3517  {
3518  PolyOrder[npCount + i] = PolyOrder[npCount + i];
3519  }
3520  }
3521  }
3522  m_file2 << "<E COMPOSITE= \"C[" << e+1
3523  << "]\" NUMMODES=\"" << PolyOrder[npCount + 1]
3524  << "\" TYPE=\"MODIFIED\" FIELDS=\"rho,rhou,rhov,rhow,E\" />"
3525  << endl;
3526  npCount += nQuadPointsElement;
3527  }
3528 
3529  m_file2.close();
3530  }
3531 
3532  void CompressibleFlowSystem::v_ExtraFldOutput(
3533  std::vector<Array<OneD, NekDouble> > &fieldcoeffs,
3534  std::vector<std::string> &variables)
3535  {
3536  const int nPhys = m_fields[0]->GetNpoints();
3537  const int nCoeffs = m_fields[0]->GetNcoeffs();
3538  Array<OneD, Array<OneD, NekDouble> > tmp(m_fields.num_elements());
3539 
3540  for (int i = 0; i < m_fields.num_elements(); ++i)
3541  {
3542  tmp[i] = m_fields[i]->GetPhys();
3543  }
3544 
3545  Array<OneD, NekDouble> pressure(nPhys), soundspeed(nPhys), mach(nPhys);
3546  Array<OneD, NekDouble> sensor(nPhys), SensorKappa(nPhys), smooth(nPhys);
3547 
3548  GetPressure (tmp, pressure);
3549  GetSoundSpeed(tmp, pressure, soundspeed);
3550  GetMach (tmp, soundspeed, mach);
3551  GetSensor (tmp, sensor, SensorKappa);
3552  GetSmoothArtificialViscosity (tmp, smooth);
3553 
3554  Array<OneD, NekDouble> pFwd(nCoeffs), sFwd(nCoeffs), mFwd(nCoeffs);
3555  Array<OneD, NekDouble> sensFwd(nCoeffs), smoothFwd(nCoeffs);
3556 
3557  m_fields[0]->FwdTrans(pressure, pFwd);
3558  m_fields[0]->FwdTrans(soundspeed, sFwd);
3559  m_fields[0]->FwdTrans(mach, mFwd);
3560  m_fields[0]->FwdTrans(sensor, sensFwd);
3561  m_fields[0]->FwdTrans(smooth, smoothFwd);
3562 
3563  variables.push_back ("p");
3564  variables.push_back ("a");
3565  variables.push_back ("Mach");
3566  variables.push_back ("Sensor");
3567  variables.push_back ("SmoothVisc");
3568  fieldcoeffs.push_back(pFwd);
3569  fieldcoeffs.push_back(sFwd);
3570  fieldcoeffs.push_back(mFwd);
3571  fieldcoeffs.push_back(sensFwd);
3572  fieldcoeffs.push_back(smoothFwd);
3573  }
3574 }
Nektar::CompressibleFlowSystem::GetMach
void GetMach(Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, NekDouble > &soundspeed, Array< OneD, NekDouble > &mach)
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:188
Tau
Definition: Tau.hpp:39
Nektar::CompressibleFlowSystem::v_InitObject
virtual void v_InitObject()
Initialization object for CompressibleFlowSystem class.
Definition: CompressibleFlowSystem.cpp:65
Nektar::CompressibleFlowSystem::m_fieldStorage
Array< OneD, Array< OneD, NekDouble > > m_fieldStorage
Definition: CompressibleFlowSystem.h:132
Nektar::LibUtilities::PointsKeyVector
std::vector< PointsKey > PointsKeyVector
Definition: Points.h:220
Nektar::CompressibleFlowSystem::~CompressibleFlowSystem
virtual ~CompressibleFlowSystem()
Destructor for CompressibleFlowSystem class.
Definition: CompressibleFlowSystem.cpp:302
Nektar::CompressibleFlowSystem::GetSoundSpeed
void GetSoundSpeed(const Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, NekDouble > &pressure, Array< OneD, NekDouble > &soundspeed)
Nektar::SolverUtils::UnsteadySystem::m_homoInitialFwd
bool m_homoInitialFwd
Flag to determine if simulation should start in homogeneous forward transformed state.
Definition: UnsteadySystem.h:80
Nektar::LibUtilities::NekFactory::CreateInstance
tBaseSharedPtr CreateInstance(tKey idKey BOOST_PP_COMMA_IF(MAX_PARAM) BOOST_PP_ENUM_BINARY_PARAMS(MAX_PARAM, tParam, x))
Create an instance of the class referred to by idKey.
Definition: NekFactory.hpp:162
Nektar::CompressibleFlowSystem::WallBC
void WallBC(int bcRegion, int cnt, Array< OneD, Array< OneD, NekDouble > > &Fwd, Array< OneD, Array< OneD, NekDouble > > &physarray)
Wall boundary conditions for compressible flow problems.
Definition: CompressibleFlowSystem.cpp:402
Nektar::SolverUtils::EquationSystem::m_time
NekDouble m_time
Current time of simulation.
Definition: EquationSystem.h:472
Nektar::CompressibleFlowSystem::GetViscousFluxVector
void GetViscousFluxVector(const Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &derivatives, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &viscousTensor)
Return the flux vector for the LDG diffusion problem.
Definition: CompressibleFlowSystem.cpp:1725
Nektar::CompressibleFlowSystem::GetStdVelocity
void GetStdVelocity(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, NekDouble > &stdV)
Vmath::Vsqrt
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
Definition: Vmath.cpp:394
Nektar::CompressibleFlowSystem::m_riemannSolverLDG
SolverUtils::RiemannSolverSharedPtr m_riemannSolverLDG
Definition: CompressibleFlowSystem.h:80
Vmath::Vmax
T Vmax(int n, const T *x, const int incx)
Return the maximum element in x – called vmax to avoid conflict with max.
Definition: Vmath.cpp:765
Nektar::CompressibleFlowSystem::v_PostIntegrate
virtual bool v_PostIntegrate(int step)
Nektar::SolverUtils::SummaryList
std::vector< std::pair< std::string, std::string > > SummaryList
Definition: Misc.h:47
Vmath::Vmin
T Vmin(int n, const T *x, const int incx)
Return the minimum element in x - called vmin to avoid conflict with min.
Definition: Vmath.cpp:857
Nektar::SolverUtils::EquationSystem::GetNumExpModesPerExp
SOLVER_UTILS_EXPORT const Array< OneD, int > GetNumExpModesPerExp()
Definition: EquationSystem.h:826
Vmath::Fill
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition: Vmath.cpp:46
Nektar::SpatialDomains::QuadGeomSharedPtr
boost::shared_ptr< QuadGeom > QuadGeomSharedPtr
Definition: QuadGeom.h:54
Nektar::CompressibleFlowSystem::GetStabilityLimit
NekDouble GetStabilityLimit(int n)
Function to calculate the stability limit for DG/CG.
Nektar::SolverUtils::GetDiffusionFactory
DiffusionFactory & GetDiffusionFactory()
Definition: Diffusion.cpp:42
Vmath::Svtvp
void Svtvp(int n, const T alpha, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
svtvp (scalar times vector plus vector): z = alpha*x + y
Definition: Vmath.cpp:471
Vmath::Vvtvp
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
Definition: Vmath.cpp:428
Nektar::SolverUtils::EquationSystem::m_projectionType
enum MultiRegions::ProjectionType m_projectionType
Type of projection; e.g continuous or discontinuous.
Definition: EquationSystem.h:514
Nektar::CompressibleFlowSystem::SymmetryBC
void SymmetryBC(int bcRegion, int cnt, Array< OneD, Array< OneD, NekDouble > > &Fwd, Array< OneD, Array< OneD, NekDouble > > &physarray)
Symmetry boundary conditions for compressible flow problems.
Definition: CompressibleFlowSystem.cpp:558
Nektar::CompressibleFlowSystem::RiemannInvariantBC
void RiemannInvariantBC(int bcRegion, int cnt, Array< OneD, Array< OneD, NekDouble > > &Fwd, Array< OneD, Array< OneD, NekDouble > > &physarray)
Outflow characteristic boundary conditions for compressible flow problems.
Definition: CompressibleFlowSystem.cpp:646
std
STL namespace.
Nektar::SolverUtils::EquationSystem::m_specHP_dealiasing
bool m_specHP_dealiasing
Flag to determine if dealisising is usde for the Spectral/hp element discretisation.
Definition: EquationSystem.h:512
Nektar::CompressibleFlowSystem::GetTemperature
void GetTemperature(const Array< OneD, const Array< OneD, NekDouble > > &physfield, Array< OneD, NekDouble > &pressure, Array< OneD, NekDouble > &temperature)
Nektar::CompressibleFlowSystem::GetArtificialDynamicViscosity
void GetArtificialDynamicViscosity(const Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, NekDouble > &mu_var)
Nektar::SolverUtils::Forcing::Load
static SOLVER_UTILS_EXPORT std::vector< ForcingSharedPtr > Load(const LibUtilities::SessionReaderSharedPtr &pSession, const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const unsigned int &pNumForcingFields=0)
Definition: Forcing.cpp:86
Nektar::CompressibleFlowSystem::m_vecLocs
Array< OneD, Array< OneD, NekDouble > > m_vecLocs
Definition: CompressibleFlowSystem.h:83
Nektar::CompressibleFlowSystem::m_un
Array< OneD, Array< OneD, NekDouble > > m_un
Definition: CompressibleFlowSystem.h:135
Vmath::Vdiv
void Vdiv(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x/y.
Definition: Vmath.cpp:227
Nektar::CompressibleFlowSystem::GetViscousFluxVectorDeAlias
void GetViscousFluxVectorDeAlias(const Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &derivatives, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &viscousTensor)
Return the flux vector for the LDG diffusion problem.
Definition: CompressibleFlowSystem.cpp:2058
Nektar::SolverUtils::UnsteadySystem::MaxTimeStepEstimator
SOLVER_UTILS_EXPORT NekDouble MaxTimeStepEstimator()
Get the maximum timestep estimator for cfl control.
Definition: UnsteadySystem.cpp:133
Nektar::SpatialDomains::HexGeomSharedPtr
boost::shared_ptr< HexGeom > HexGeomSharedPtr
Definition: HexGeom.h:110
Nektar::LibUtilities::SessionReaderSharedPtr
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:51
Vmath::Vabs
void Vabs(int n, const T *x, const int incx, T *y, const int incy)
vabs: y = |x|
Definition: Vmath.cpp:410
CompressibleFlowSystem.h
Nektar::SolverUtils::EquationSystem::m_comm
LibUtilities::CommSharedPtr m_comm
Communicator.
Definition: EquationSystem.h:450
Nektar::SolverUtils::EquationSystem::GetTotPoints
SOLVER_UTILS_EXPORT int GetTotPoints()
Definition: EquationSystem.h:866
Nektar::CompressibleFlowSystem::PressureOutflowBC
void PressureOutflowBC(int bcRegion, int cnt, Array< OneD, Array< OneD, NekDouble > > &Fwd, Array< OneD, Array< OneD, NekDouble > > &physarray)
Pressure outflow boundary conditions for compressible flow problems.
Definition: CompressibleFlowSystem.cpp:1048
Nektar::CompressibleFlowSystem::PressureOutflowFileBC
void PressureOutflowFileBC(int bcRegion, int cnt, Array< OneD, Array< OneD, NekDouble > > &Fwd, Array< OneD, Array< OneD, NekDouble > > &physarray)
Pressure outflow boundary conditions for compressible flow problems.
Definition: CompressibleFlowSystem.cpp:1202
Nektar::CompressibleFlowSystem::PressureOutflowNonReflectiveBC
void PressureOutflowNonReflectiveBC(int bcRegion, int cnt, Array< OneD, Array< OneD, NekDouble > > &Fwd, Array< OneD, Array< OneD, NekDouble > > &physarray)
Pressure outflow non-reflective boundary conditions for compressible flow problems.
Definition: CompressibleFlowSystem.cpp:895
Nektar::CompressibleFlowSystem::v_GenerateSummary
virtual void v_GenerateSummary(SolverUtils::SummaryList &s)
Print a summary of time stepping parameters.
Definition: CompressibleFlowSystem.cpp:310
Nektar::SolverUtils::EquationSystem::m_traceNormals
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Array holding trace normals for DG simulations in the forwards direction.
Definition: EquationSystem.h:516
Nektar::CompressibleFlowSystem::SetVarPOrderElmt
void SetVarPOrderElmt(const Array< OneD, const Array< OneD, NekDouble > > &physfield, Array< OneD, NekDouble > &PolyOrder)
Nektar::CompressibleFlowSystem::GetVelocityVector
void GetVelocityVector(const Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &velocity)
Nektar::SolverUtils::UnsteadySystem::v_GenerateSummary
virtual SOLVER_UTILS_EXPORT void v_GenerateSummary(SummaryList &s)
Print a summary of time stepping parameters.
Definition: UnsteadySystem.cpp:460
Nektar::CompressibleFlowSystem::WallViscousBC
void WallViscousBC(int bcRegion, int cnt, Array< OneD, Array< OneD, NekDouble > > &Fwd, Array< OneD, Array< OneD, NekDouble > > &physarray)
Wall boundary conditions for viscous compressible flow problems.
Definition: CompressibleFlowSystem.cpp:491
Nektar::CompressibleFlowSystem::GetVecLocs
const Array< OneD, const Array< OneD, NekDouble > > & GetVecLocs()
Definition: CompressibleFlowSystem.h:294
TriExp.h
Vmath::Smul
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
Definition: Vmath.cpp:199
ExpList.h
Nektar::CompressibleFlowSystem::GetFluxVectorDeAlias
void GetFluxVectorDeAlias(const Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &flux)
Return the flux vector for the compressible Euler equations by using the de-aliasing technique...
Definition: CompressibleFlowSystem.cpp:1634
Nektar::CompressibleFlowSystem::GetEntropy
void GetEntropy(const Array< OneD, const Array< OneD, NekDouble > > &physfield, const Array< OneD, const NekDouble > &pressure, const Array< OneD, const NekDouble > &temperature, Array< OneD, NekDouble > &entropy)
Nektar::SolverUtils::UnsteadySystem
Base class for unsteady solvers.
Definition: UnsteadySystem.h:48
Nektar::SolverUtils::GetRiemannSolverFactory
RiemannSolverFactory & GetRiemannSolverFactory()
Definition: RiemannSolver.cpp:63
HexExp.h
Nektar::CompressibleFlowSystem::PressureInflowFileBC
void PressureInflowFileBC(int bcRegion, int cnt, Array< OneD, Array< OneD, NekDouble > > &Fwd, Array< OneD, Array< OneD, NekDouble > > &physarray)
Pressure inflow boundary conditions for compressible flow problems where either the density and the v...
Definition: CompressibleFlowSystem.cpp:1358
npts
static std::string npts
Definition: InputFld.cpp:43
Nektar::SolverUtils::EquationSystem::m_spacedim
int m_spacedim
Spatial dimension (>= expansion dim).
Definition: EquationSystem.h:492
Nektar::SolverUtils::GetAdvectionFactory
AdvectionFactory & GetAdvectionFactory()
Gets the factory for initialising advection objects.
Definition: Advection.cpp:46
Vmath::Neg
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:382
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:44
Nektar::SolverUtils::UnsteadySystem::v_InitObject
virtual SOLVER_UTILS_EXPORT void v_InitObject()
Init object for UnsteadySystem class.
Definition: UnsteadySystem.cpp:78
QuadExp.h
Nektar::SolverUtils::UnsteadySystem::m_cflSafetyFactor
NekDouble m_cflSafetyFactor
CFL safety factor (comprise between 0 to 1).
Definition: UnsteadySystem.h:59
Nektar::CompressibleFlowSystem::v_ExtraFldOutput
virtual void v_ExtraFldOutput(std::vector< Array< OneD, NekDouble > > &fieldcoeffs, std::vector< std::string > &variables)
Nektar::CompressibleFlowSystem::CalcSteadyState
bool CalcSteadyState(bool output)
Vmath::Sadd
void Sadd(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Add vector y = alpha + x.
Definition: Vmath.cpp:301
Nektar::CompressibleFlowSystem::m_pressureStorage
Array< OneD, NekDouble > m_pressureStorage
Definition: CompressibleFlowSystem.h:129
Nektar::SolverUtils::GetEquationSystemFactory
EquationSystemFactory & GetEquationSystemFactory()
Definition: EquationSystem.cpp:83
Nektar::CompressibleFlowSystem::GetStabilityLimitVector
Array< OneD, NekDouble > GetStabilityLimitVector(const Array< OneD, int > &ExpOrder)
Function to calculate the stability limit for DG/CG (a vector of them).
Nektar::CompressibleFlowSystem::GetAbsoluteVelocity
void GetAbsoluteVelocity(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, NekDouble > &Vtot)
Nektar::CompressibleFlowSystem::m_advection
SolverUtils::AdvectionSharedPtr m_advection
Definition: CompressibleFlowSystem.h:81
Vmath::Vsub
void Vsub(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Subtract vector z = x-y.
Definition: Vmath.cpp:329
Nektar::CompressibleFlowSystem::m_diffusion
SolverUtils::DiffusionSharedPtr m_diffusion
Definition: CompressibleFlowSystem.h:82
Nektar::SolverUtils::EquationSystem::GetTraceTotPoints
SOLVER_UTILS_EXPORT int GetTraceTotPoints()
Definition: EquationSystem.h:841
Nektar::SolverUtils::EquationSystem::GetPhys_Offset
SOLVER_UTILS_EXPORT int GetPhys_Offset(int n)
Definition: EquationSystem.h:856
Nektar::CompressibleFlowSystem::GetNormals
const Array< OneD, const Array< OneD, NekDouble > > & GetNormals()
Definition: CompressibleFlowSystem.h:299
Nektar::SpatialDomains::GeomFactorsSharedPtr
boost::shared_ptr< GeomFactors > GeomFactorsSharedPtr
Pointer to a GeomFactors object.
Definition: GeomFactors.h:62
Nektar::SolverUtils::EquationSystem::m_fields
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
Definition: EquationSystem.h:460
Nektar::CompressibleFlowSystem::v_SetInitialConditions
virtual void v_SetInitialConditions(NekDouble initialtime=0.0, bool dumpInitialConditions=true, const int domain=0)
Nektar::CompressibleFlowSystem::GetDynamicViscosity
void GetDynamicViscosity(const Array< OneD, const NekDouble > &temperature, Array< OneD, NekDouble > &mu)
Nektar::CompressibleFlowSystem::v_GetTimeStep
virtual NekDouble v_GetTimeStep(const Array< OneD, const Array< OneD, NekDouble > > &inarray)
Return the timestep to be used for the next step in the time-marching loop.
Nektar::SolverUtils::EquationSystem::GetExpSize
SOLVER_UTILS_EXPORT int GetExpSize()
Definition: EquationSystem.h:851
Nektar::CompressibleFlowSystem::GetEnthalpy
void GetEnthalpy(const Array< OneD, const Array< OneD, NekDouble > > &physfield, Array< OneD, NekDouble > &pressure, Array< OneD, NekDouble > &enthalpy)
Nektar::SolverUtils::EquationSystem::m_session
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
Definition: EquationSystem.h:452
Nektar::CompressibleFlowSystem::GetSmoothArtificialViscosity
void GetSmoothArtificialViscosity(const Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, NekDouble > &eps_bar)
Nektar::CompressibleFlowSystem::SetCommonBC
void SetCommonBC(const std::string &userDefStr, const int n, const NekDouble time, int &cnt, Array< OneD, Array< OneD, NekDouble > > &Fwd, Array< OneD, Array< OneD, NekDouble > > &inarray)
Set boundary conditions which can be: a) Wall and Symmerty BCs implemented at CompressibleFlowSystem ...
Definition: CompressibleFlowSystem.cpp:324
Nektar::CompressibleFlowSystem::GetForcingTerm
void GetForcingTerm(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > outarrayForcing)
Nektar::SolverUtils::EquationSystem::GetPressure
SOLVER_UTILS_EXPORT MultiRegions::ExpListSharedPtr GetPressure()
Get pressure field if available.
Definition: EquationSystem.h:748
Nektar::SolverUtils::UnsteadySystem::m_infosteps
int m_infosteps
Number of time steps between outputting status information.
Definition: UnsteadySystem.h:63
Nektar::CompressibleFlowSystem::GetSensor
void GetSensor(const Array< OneD, const Array< OneD, NekDouble > > &physarray, Array< OneD, NekDouble > &Sensor, Array< OneD, NekDouble > &SensorKappa)
InterpCoeff.h
AssemblyMapDG.h
Nektar::CompressibleFlowSystem::m_riemannSolver
SolverUtils::RiemannSolverSharedPtr m_riemannSolver
Definition: CompressibleFlowSystem.h:79
Nektar::CompressibleFlowSystem::GetFluxVector
void GetFluxVector(const Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &flux)
Return the flux vector for the compressible Euler equations.
Definition: CompressibleFlowSystem.cpp:1571
Nektar::CompressibleFlowSystem::ExtrapOrder0BC
void ExtrapOrder0BC(int bcRegion, int cnt, Array< OneD, Array< OneD, NekDouble > > &Fwd, Array< OneD, Array< OneD, NekDouble > > &physarray)
Extrapolation of order 0 for all the variables such that, at the boundaries, a trivial Riemann proble...
Definition: CompressibleFlowSystem.cpp:1513
Vmath::Vsum
T Vsum(int n, const T *x, const int incx)
Subtract return sum(x)
Definition: Vmath.cpp:723
Vmath::Zero
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:359
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1047
Nektar::SpatialDomains::eDeformed
Geometry is curved or has non-constant factors.
Definition: SpatialDomains.hpp:57
Nektar::CompressibleFlowSystem::GetElementDimensions
void GetElementDimensions(Array< OneD, Array< OneD, NekDouble > > &outarray, Array< OneD, NekDouble > &hmin)
Nektar::CompressibleFlowSystem::m_forcing
std::vector< SolverUtils::ForcingSharedPtr > m_forcing
Definition: CompressibleFlowSystem.h:122
Vmath::Vadd
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.cpp:285
Vmath::Vmul
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.cpp:169
Nektar::LibUtilities::NekFactory::RegisterCreatorFunction
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215
Generated on Fri Sep 9 2016 12:03:13 for Nektar++ by   doxygen 1.8.8