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