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 <iomanip>
36 
37 #include <CompressibleFlowSolver/EquationSystems/CompressibleFlowSystem.h>
38 #include <LocalRegions/TriExp.h>
39 #include <LocalRegions/QuadExp.h>
40 #include <LocalRegions/HexExp.h>
41 #include <MultiRegions/ExpList.h>
42 #include <LibUtilities/Foundations/InterpCoeff.h>
43 #include <MultiRegions/AssemblyMap/AssemblyMapDG.h>
44 
45 namespace Nektar
46 {
47  string CompressibleFlowSystem::className =
48  SolverUtils::GetEquationSystemFactory().RegisterCreatorFunction(
49  "CompressibleFlowSystem",
50  CompressibleFlowSystem::create,
51  "Auxiliary functions for the compressible flow system.");
52 
53  CompressibleFlowSystem::CompressibleFlowSystem(
54  const LibUtilities::SessionReaderSharedPtr& pSession)
55  : UnsteadySystem(pSession)
56  {
57  }
58 
59  /**
60  * @brief Initialization object for CompressibleFlowSystem class.
61  */
62  void CompressibleFlowSystem::v_InitObject()
63  {
64  UnsteadySystem::v_InitObject();
65 
66  ASSERTL0(m_session->DefinesSolverInfo("UPWINDTYPE"),
67  "No UPWINDTYPE defined in session.");
68 
69  // Do not forwards transform initial condition
70  m_homoInitialFwd = false;
71 
72  // Set up locations of velocity vector.
73  m_vecLocs = Array<OneD, Array<OneD, NekDouble> >(1);
74  m_vecLocs[0] = Array<OneD, NekDouble>(m_spacedim);
75  for (int i = 0; i < m_spacedim; ++i)
76  {
77  m_vecLocs[0][i] = 1 + i;
78  }
79 
80  // Get gamma parameter from session file.
81  ASSERTL0(m_session->DefinesParameter("Gamma"),
82  "Compressible flow sessions must define a Gamma parameter.");
83  m_session->LoadParameter("Gamma", m_gamma, 1.4);
84 
85  // Get E0 parameter from session file.
86  ASSERTL0(m_session->DefinesParameter("pInf"),
87  "Compressible flow sessions must define a pInf parameter.");
88  m_session->LoadParameter("pInf", m_pInf, 101325);
89 
90  // Get rhoInf parameter from session file.
91  ASSERTL0(m_session->DefinesParameter("rhoInf"),
92  "Compressible flow sessions must define a rhoInf parameter.");
93  m_session->LoadParameter("rhoInf", m_rhoInf, 1.225);
94 
95  // Get uInf parameter from session file.
96  ASSERTL0(m_session->DefinesParameter("uInf"),
97  "Compressible flow sessions must define a uInf parameter.");
98  m_session->LoadParameter("uInf", m_uInf, 0.1);
99 
100  // Get vInf parameter from session file.
101  if (m_spacedim == 2 || m_spacedim == 3)
102  {
103  ASSERTL0(m_session->DefinesParameter("vInf"),
104  "Compressible flow sessions must define a vInf parameter"
105  "for 2D/3D problems.");
106  m_session->LoadParameter("vInf", m_vInf, 0.0);
107  }
108 
109  // Get wInf parameter from session file.
110  if (m_spacedim == 3)
111  {
112  ASSERTL0(m_session->DefinesParameter("wInf"),
113  "Compressible flow sessions must define a wInf parameter"
114  "for 3D problems.");
115  m_session->LoadParameter("wInf", m_wInf, 0.0);
116  }
117 
118  m_session->LoadParameter ("GasConstant", m_gasConstant, 287.058);
119  m_session->LoadParameter ("Twall", m_Twall, 300.15);
120  m_session->LoadSolverInfo("ViscosityType", m_ViscosityType, "Constant");
121  m_session->LoadParameter ("mu", m_mu, 1.78e-05);
122  m_session->LoadParameter ("Skappa", m_Skappa, -2.048);
123  m_session->LoadParameter ("Kappa", m_Kappa, 0.0);
124  m_session->LoadParameter ("mu0", m_mu0, 1.0);
125  m_session->LoadParameter ("FL", m_FacL, 0.0);
126  m_session->LoadParameter ("FH", m_FacH, 0.0);
127  m_session->LoadParameter ("hFactor", m_hFactor, 1.0);
128  m_session->LoadParameter ("epsMax", m_eps_max, 1.0);
129  m_session->LoadParameter ("C1", m_C1, 3.0);
130  m_session->LoadParameter ("C2", m_C2, 5.0);
131  m_session->LoadSolverInfo("ShockCaptureType",
132  m_shockCaptureType, "Off");
133  m_session->LoadParameter ("thermalConductivity",
134  m_thermalConductivity, 0.0257);
135 
136  m_EqTypeStr = m_session->GetSolverInfo("EQTYPE");
137 
138  m_Cp = m_gamma / (m_gamma - 1.0) * m_gasConstant;
139  m_Prandtl = m_Cp * m_mu / m_thermalConductivity;
140 
141  // Type of advection class to be used
142  switch(m_projectionType)
143  {
144  // Continuous field
145  case MultiRegions::eGalerkin:
146  {
147  ASSERTL0(false, "Continuous field not supported.");
148  break;
149  }
150  // Discontinuous field
151  case MultiRegions::eDiscontinuous:
152  {
153  string advName, diffName, riemName;
154 
155  // Setting up advection and diffusion operators
156  m_session->LoadSolverInfo("AdvectionType", advName, "WeakDG");
157  m_session->LoadSolverInfo("DiffusionType", diffName, "LDGNS");
158 
159  m_advection = SolverUtils::GetAdvectionFactory()
160  .CreateInstance(advName, advName);
161  m_diffusion = SolverUtils::GetDiffusionFactory()
162  .CreateInstance(diffName, diffName);
163 
164  if (m_specHP_dealiasing)
165  {
166  m_advection->SetFluxVector(&CompressibleFlowSystem::
167  GetFluxVectorDeAlias, this);
168  m_diffusion->SetFluxVectorNS(
169  &CompressibleFlowSystem::GetViscousFluxVectorDeAlias,
170  this);
171  }
172  else
173  {
174  m_advection->SetFluxVector (&CompressibleFlowSystem::
175  GetFluxVector, this);
176  m_diffusion->SetFluxVectorNS(&CompressibleFlowSystem::
177  GetViscousFluxVector, this);
178  }
179 
180  if (m_shockCaptureType=="Smooth" && m_EqTypeStr=="EulerADCFE")
181  {
182  m_advection->SetFluxVector(&CompressibleFlowSystem::
183  GetFluxVector, this);
184 
185  m_diffusion->SetArtificialDiffusionVector(
186  &CompressibleFlowSystem::GetSmoothArtificialViscosity, this);
187  }
188  if (m_shockCaptureType=="NonSmooth" && m_EqTypeStr=="EulerADCFE")
189  {
190  m_advection->SetFluxVector(&CompressibleFlowSystem::
191  GetFluxVector, this);
192 
193  m_diffusion->SetArtificialDiffusionVector(
194  &CompressibleFlowSystem::GetArtificialDynamicViscosity, this);
195  }
196 
197  // Setting up Riemann solver for advection operator
198  m_session->LoadSolverInfo("UpwindType", riemName, "Average");
199 
200  m_riemannSolver = SolverUtils::GetRiemannSolverFactory()
201  .CreateInstance(riemName);
202 
203  // Setting up upwind solver for diffusion operator
204  m_riemannSolverLDG = SolverUtils::GetRiemannSolverFactory()
205  .CreateInstance("UpwindLDG");
206 
207  // Setting up parameters for advection operator Riemann solver
208  m_riemannSolver->SetParam (
209  "gamma", &CompressibleFlowSystem::GetGamma, this);
210  m_riemannSolver->SetAuxVec(
211  "vecLocs", &CompressibleFlowSystem::GetVecLocs, this);
212  m_riemannSolver->SetVector(
213  "N", &CompressibleFlowSystem::GetNormals, this);
214 
215  // Setting up parameters for diffusion operator Riemann solver
216  m_riemannSolverLDG->SetParam (
217  "gamma", &CompressibleFlowSystem::GetGamma, this);
218  m_riemannSolverLDG->SetVector(
219  "vecLocs", &CompressibleFlowSystem::GetVecLocs, this);
220  m_riemannSolverLDG->SetVector(
221  "N", &CompressibleFlowSystem::GetNormals, this);
222 
223  // Concluding initialisation of advection / diffusion operators
224  m_advection->SetRiemannSolver (m_riemannSolver);
225  m_diffusion->SetRiemannSolver (m_riemannSolverLDG);
226  m_advection->InitObject (m_session, m_fields);
227  m_diffusion->InitObject (m_session, m_fields);
228  break;
229  }
230  default:
231  {
232  ASSERTL0(false, "Unsupported projection type.");
233  break;
234  }
235  }
236  }
237 
238  /**
239  * @brief Destructor for CompressibleFlowSystem class.
240  */
241  CompressibleFlowSystem::~CompressibleFlowSystem()
242  {
243 
244  }
245 
246  /**
247  * @brief Print out a summary with some relevant information.
248  */
249  void CompressibleFlowSystem::v_GenerateSummary(SolverUtils::SummaryList& s)
250  {
251  UnsteadySystem::v_GenerateSummary(s);
252  }
253 
254  /**
255  * @brief Wall boundary conditions for compressible flow problems.
256  */
257  void CompressibleFlowSystem::WallBC(
258  int bcRegion,
259  int cnt,
260  Array<OneD, Array<OneD, NekDouble> > &physarray)
261  {
262  int i;
263  int nTracePts = GetTraceTotPoints();
264  int nVariables = physarray.num_elements();
265 
266  const Array<OneD, const int> &traceBndMap
267  = m_fields[0]->GetTraceBndMap();
268 
269  // Get physical values of the forward trace
270  Array<OneD, Array<OneD, NekDouble> > Fwd(nVariables);
271  for (i = 0; i < nVariables; ++i)
272  {
273  Fwd[i] = Array<OneD, NekDouble>(nTracePts);
274  m_fields[i]->ExtractTracePhys(physarray[i], Fwd[i]);
275  }
276 
277  // Adjust the physical values of the trace to take
278  // user defined boundaries into account
279  int e, id1, id2, nBCEdgePts, eMax;
280 
281  eMax = m_fields[0]->GetBndCondExpansions()[bcRegion]->GetExpSize();
282 
283  for (e = 0; e < eMax; ++e)
284  {
285  nBCEdgePts = m_fields[0]->GetBndCondExpansions()[bcRegion]->
286  GetExp(e)->GetTotPoints();
287  id1 = m_fields[0]->GetBndCondExpansions()[bcRegion]->
288  GetPhys_Offset(e);
289  id2 = m_fields[0]->GetTrace()->GetPhys_Offset(traceBndMap[cnt+e]);
290 
291  // Boundary condition for epsilon term.
292  if (nVariables == m_spacedim+3)
293  {
294  NekDouble factor = 1.0;
295  NekDouble factor2 = 1.0;
296 
297  Array<OneD, NekDouble > tmp2(nBCEdgePts, 0.0);
298  Vmath::Smul(nBCEdgePts,
299  factor,
300  &Fwd[nVariables-1][id2], 1,
301  &tmp2[0], 1);
302 
303  Vmath::Vsub(nBCEdgePts,
304  &Fwd[nVariables-1][id2], 1,
305  &tmp2[0], 1,
306  &Fwd[nVariables-1][id2], 1);
307 
308  Vmath::Smul(nBCEdgePts,
309  factor2,
310  &Fwd[nVariables-1][id2], 1,
311  &Fwd[nVariables-1][id2], 1);
312 
313  }
314 
315  // For 2D/3D, define: v* = v - 2(v.n)n
316  Array<OneD, NekDouble> tmp(nBCEdgePts, 0.0);
317 
318  // Calculate (v.n)
319  for (i = 0; i < m_spacedim; ++i)
320  {
321  Vmath::Vvtvp(nBCEdgePts,
322  &Fwd[1+i][id2], 1,
323  &m_traceNormals[i][id2], 1,
324  &tmp[0], 1,
325  &tmp[0], 1);
326  }
327 
328  // Calculate 2.0(v.n)
329  Vmath::Smul(nBCEdgePts, -2.0, &tmp[0], 1, &tmp[0], 1);
330 
331  // Calculate v* = v - 2.0(v.n)n
332  for (i = 0; i < m_spacedim; ++i)
333  {
334  Vmath::Vvtvp(nBCEdgePts,
335  &tmp[0], 1,
336  &m_traceNormals[i][id2], 1,
337  &Fwd[1+i][id2], 1,
338  &Fwd[1+i][id2], 1);
339  }
340 
341  // Copy boundary adjusted values into the boundary expansion
342  for (i = 0; i < nVariables; ++i)
343  {
344  Vmath::Vcopy(nBCEdgePts, &Fwd[i][id2], 1,
345  &(m_fields[i]->GetBndCondExpansions()[bcRegion]->
346  UpdatePhys())[id1], 1);
347  }
348  }
349  }
350 
351  /**
352  * @brief Wall boundary conditions for viscous compressible flow problems.
353  */
354  void CompressibleFlowSystem::WallViscousBC(
355  int bcRegion,
356  int cnt,
357  Array<OneD, Array<OneD, NekDouble> > &physarray)
358  {
359  int i;
360  int nTracePts = GetTraceTotPoints();
361  int nVariables = physarray.num_elements();
362 
363  const Array<OneD, const int> &traceBndMap
364  = m_fields[0]->GetTraceBndMap();
365 
366  // Get physical values of the forward trace
367  Array<OneD, Array<OneD, NekDouble> > Fwd(nVariables);
368  for (i = 0; i < nVariables; ++i)
369  {
370  Fwd[i] = Array<OneD, NekDouble>(nTracePts);
371  m_fields[i]->ExtractTracePhys(physarray[i], Fwd[i]);
372  }
373 
374  // Take into account that for PDE based shock capturing, eps = 0 at the
375  // wall. Adjust the physical values of the trace to take user defined
376  // boundaries into account
377  int e, id1, id2, nBCEdgePts, eMax;
378 
379  eMax = m_fields[0]->GetBndCondExpansions()[bcRegion]->GetExpSize();
380 
381  for (e = 0; e < eMax; ++e)
382  {
383  nBCEdgePts = m_fields[0]->GetBndCondExpansions()[bcRegion]->
384  GetExp(e)->GetTotPoints();
385  id1 = m_fields[0]->GetBndCondExpansions()[bcRegion]->
386  GetPhys_Offset(e);
387  id2 = m_fields[0]->GetTrace()->GetPhys_Offset(traceBndMap[cnt+e]);
388 
389  if (nVariables == m_spacedim+3)
390  {
391  NekDouble factor = 0.0;
392  NekDouble factor2 = 1.0;
393 
394  Array<OneD, NekDouble > tmp2(nBCEdgePts, 0.0);
395  Vmath::Smul(nBCEdgePts,
396  factor,
397  &Fwd[nVariables-1][id2], 1,
398  &tmp2[0], 1);
399 
400  Vmath::Vsub(nBCEdgePts,
401  &Fwd[nVariables-1][id2], 1,
402  &tmp2[0], 1,
403  &Fwd[nVariables-1][id2], 1);
404 
405  Vmath::Smul(nBCEdgePts,
406  factor2,
407  &Fwd[nVariables-1][id2], 1,
408  &Fwd[nVariables-1][id2], 1);
409  }
410 
411  for (i = 0; i < m_spacedim; i++)
412  {
413  Vmath::Neg(nBCEdgePts, &Fwd[i+1][id2], 1);
414  }
415 
416  // Copy boundary adjusted values into the boundary expansion
417  for (i = 0; i < nVariables; ++i)
418  {
419  Vmath::Vcopy(nBCEdgePts, &Fwd[i][id2], 1,
420  &(m_fields[i]->GetBndCondExpansions()[bcRegion]->
421  UpdatePhys())[id1], 1);
422  }
423  }
424  }
425 
426  /**
427  * @brief Symmetry boundary conditions for compressible flow problems.
428  */
429  void CompressibleFlowSystem::SymmetryBC(
430  int bcRegion,
431  int cnt,
432  Array<OneD, Array<OneD, NekDouble> > &physarray)
433  {
434  int i;
435  int nTracePts = GetTraceTotPoints();
436  int nVariables = physarray.num_elements();
437 
438  const Array<OneD, const int> &traceBndMap
439  = m_fields[0]->GetTraceBndMap();
440 
441  // Get physical values of the forward trace (from exp to phys)
442  Array<OneD, Array<OneD, NekDouble> > Fwd(nVariables);
443  for (i = 0; i < nVariables; ++i)
444  {
445  Fwd[i] = Array<OneD, NekDouble>(nTracePts);
446  m_fields[i]->ExtractTracePhys(physarray[i], Fwd[i]);
447  }
448 
449  // Take into account that for PDE based shock capturing, eps = 0 at the
450  // wall.
451  int e, id1, id2, nBCEdgePts, eMax;
452 
453  eMax = m_fields[0]->GetBndCondExpansions()[bcRegion]->GetExpSize();
454 
455  for (e = 0; e < eMax; ++e)
456  {
457  nBCEdgePts = m_fields[0]->GetBndCondExpansions()[bcRegion]->
458  GetExp(e)->GetTotPoints();
459  id1 = m_fields[0]->GetBndCondExpansions()[bcRegion]->
460  GetPhys_Offset(e);
461  id2 = m_fields[0]->GetTrace()->GetPhys_Offset(traceBndMap[cnt+e]);
462 
463  if (nVariables == m_spacedim+3)
464  {
465  NekDouble factor = 0.0;
466  NekDouble factor2 = 1.0;
467 
468  Array<OneD, NekDouble > tmp2(nBCEdgePts, 0.0);
469  Vmath::Smul(nBCEdgePts,
470  factor,
471  &Fwd[nVariables-1][id2], 1,
472  &tmp2[0], 1);
473 
474  Vmath::Vsub(nBCEdgePts,
475  &Fwd[nVariables-1][id2], 1,
476  &tmp2[0], 1,
477  &Fwd[nVariables-1][id2], 1);
478 
479  Vmath::Smul(nBCEdgePts,
480  factor2,
481  &Fwd[nVariables-1][id2], 1,
482  &Fwd[nVariables-1][id2], 1);
483  }
484 
485  // For 2D/3D, define: v* = v - 2(v.n)n
486  Array<OneD, NekDouble> tmp(nBCEdgePts, 0.0);
487 
488  // Calculate (v.n)
489  for (i = 0; i < m_spacedim; ++i)
490  {
491  Vmath::Vvtvp(nBCEdgePts,
492  &Fwd[1+i][id2], 1,
493  &m_traceNormals[i][id2], 1,
494  &tmp[0], 1,
495  &tmp[0], 1);
496  }
497 
498  // Calculate 2.0(v.n)
499  Vmath::Smul(nBCEdgePts, -2.0, &tmp[0], 1, &tmp[0], 1);
500 
501  // Calculate v* = v - 2.0(v.n)n
502  for (i = 0; i < m_spacedim; ++i)
503  {
504  Vmath::Vvtvp(nBCEdgePts,
505  &tmp[0], 1,
506  &m_traceNormals[i][id2], 1,
507  &Fwd[1+i][id2], 1,
508  &Fwd[1+i][id2], 1);
509  }
510 
511  // Copy boundary adjusted values into the boundary expansion
512  for (i = 0; i < nVariables; ++i)
513  {
514  Vmath::Vcopy(nBCEdgePts, &Fwd[i][id2], 1,
515  &(m_fields[i]->GetBndCondExpansions()[bcRegion]->
516  UpdatePhys())[id1], 1);
517  }
518  }
519  }
520 
521  /**
522  * @brief Outflow characteristic boundary conditions for compressible
523  * flow problems.
524  */
525  void CompressibleFlowSystem::RiemannInvariantBC(
526  int bcRegion,
527  int cnt,
528  Array<OneD, Array<OneD, NekDouble> > &physarray)
529  {
530  int i, j;
531  int nTracePts = GetTraceTotPoints();
532  int nVariables = physarray.num_elements();
533  int nDimensions = m_spacedim;
534 
535  const Array<OneD, const int> &traceBndMap
536  = m_fields[0]->GetTraceBndMap();
537 
538  NekDouble gamma = m_gamma;
539  NekDouble gammaInv = 1.0 / gamma;
540  NekDouble gammaMinusOne = gamma - 1.0;
541  NekDouble gammaMinusOneInv = 1.0 / gammaMinusOne;
542 
543  Array<OneD, NekDouble> tmp1 (nTracePts, 0.0);
544  Array<OneD, NekDouble> tmp2 (nTracePts, 0.0);
545  Array<OneD, NekDouble> VnInf(nTracePts, 0.0);
546  Array<OneD, NekDouble> velInf(nDimensions, 0.0);
547 
548  // Computing the normal velocity for characteristics coming
549  // from outside the computational domain
550  velInf[0] = m_uInf;
551  Vmath::Smul(nTracePts, m_uInf, m_traceNormals[0], 1, VnInf, 1);
552  if (nDimensions == 2 || nDimensions == 3)
553  {
554  velInf[1] = m_vInf;
555  Vmath::Smul(nTracePts, m_vInf, m_traceNormals[1], 1, tmp1, 1);
556  Vmath::Vadd(nTracePts, VnInf, 1, tmp1, 1, VnInf, 1);
557  }
558  if (nDimensions == 3)
559  {
560  velInf[2] = m_wInf;
561  Vmath::Smul(nTracePts, m_wInf, m_traceNormals[2], 1, tmp2, 1);
562  Vmath::Vadd(nTracePts, VnInf, 1, tmp2, 1, VnInf, 1);
563  }
564 
565  // Get physical values of the forward trace
566  Array<OneD, Array<OneD, NekDouble> > Fwd(nVariables);
567  for (i = 0; i < nVariables; ++i)
568  {
569  Fwd[i] = Array<OneD, NekDouble>(nTracePts);
570  m_fields[i]->ExtractTracePhys(physarray[i], Fwd[i]);
571  }
572 
573  // Computing the normal velocity for characteristics coming
574  // from inside the computational domain
575  Array<OneD, NekDouble > Vn (nTracePts, 0.0);
576  Array<OneD, NekDouble > Vel(nTracePts, 0.0);
577  for (i = 0; i < nDimensions; ++i)
578  {
579  Vmath::Vdiv(nTracePts, Fwd[i+1], 1, Fwd[0], 1, Vel, 1);
580  Vmath::Vvtvp(nTracePts, m_traceNormals[i], 1, Vel, 1, Vn, 1, Vn, 1);
581  }
582 
583  // Computing the absolute value of the velocity in order to compute the
584  // Mach number to decide whether supersonic or subsonic
585  Array<OneD, NekDouble > absVel(nTracePts, 0.0);
586  for (i = 0; i < nDimensions; ++i)
587  {
588  Vmath::Vdiv(nTracePts, Fwd[i+1], 1, Fwd[0], 1, tmp1, 1);
589  Vmath::Vmul(nTracePts, tmp1, 1, tmp1, 1, tmp1, 1);
590  Vmath::Vadd(nTracePts, tmp1, 1, absVel, 1, absVel, 1);
591  }
592  Vmath::Vsqrt(nTracePts, absVel, 1, absVel, 1);
593 
594  // Get speed of sound
595  Array<OneD, NekDouble > soundSpeed(nTracePts);
596  Array<OneD, NekDouble > pressure (nTracePts);
597 
598  for (i = 0; i < nTracePts; i++)
599  {
600  if (m_spacedim == 1)
601  {
602  pressure[i] = (gammaMinusOne) * (Fwd[2][i] -
603  0.5 * (Fwd[1][i] * Fwd[1][i] / Fwd[0][i]));
604  }
605  else if (m_spacedim == 2)
606  {
607  pressure[i] = (gammaMinusOne) * (Fwd[3][i] -
608  0.5 * (Fwd[1][i] * Fwd[1][i] / Fwd[0][i] +
609  Fwd[2][i] * Fwd[2][i] / Fwd[0][i]));
610  }
611  else
612  {
613  pressure[i] = (gammaMinusOne) * (Fwd[4][i] -
614  0.5 * (Fwd[1][i] * Fwd[1][i] / Fwd[0][i] +
615  Fwd[2][i] * Fwd[2][i] / Fwd[0][i] +
616  Fwd[3][i] * Fwd[3][i] / Fwd[0][i]));
617  }
618 
619  soundSpeed[i] = sqrt(gamma * pressure[i] / Fwd[0][i]);
620  }
621 
622  // Get Mach
623  Array<OneD, NekDouble > Mach(nTracePts, 0.0);
624  Vmath::Vdiv(nTracePts, Vn, 1, soundSpeed, 1, Mach, 1);
625  Vmath::Vabs(nTracePts, Mach, 1, Mach, 1);
626 
627  // Auxiliary variables
628  int eMax;
629  int e, id1, id2, nBCEdgePts, pnt;
630  NekDouble cPlus, rPlus, cMinus, rMinus, VDBC, VNBC;
631  Array<OneD, NekDouble> velBC(nDimensions, 0.0);
632  Array<OneD, NekDouble> rhoVelBC(nDimensions, 0.0);
633  NekDouble rhoBC, EBC, cBC, sBC, pBC;
634 
635  eMax = m_fields[0]->GetBndCondExpansions()[bcRegion]->GetExpSize();
636 
637  // Loop on bcRegions
638  for (e = 0; e < eMax; ++e)
639  {
640  nBCEdgePts = m_fields[0]->GetBndCondExpansions()[bcRegion]->
641  GetExp(e)->GetNumPoints(0);
642 
643  id1 = m_fields[0]->GetBndCondExpansions()[bcRegion]->
644  GetPhys_Offset(e);
645  id2 = m_fields[0]->GetTrace()->GetPhys_Offset(traceBndMap[cnt+e]);
646 
647  // Loop on the points of the bcRegion
648  for (i = 0; i < nBCEdgePts; i++)
649  {
650  pnt = id2+i;
651 
652  // Impose inflow Riemann invariant
653  if (Vn[pnt] <= 0.0)
654  {
655  // Subsonic flows
656  if (Mach[pnt] < 1.00)
657  {
658  // + Characteristic from inside
659  cPlus = sqrt(gamma * pressure[pnt] / Fwd[0][pnt]);
660  rPlus = Vn[pnt] + 2.0 * cPlus * gammaMinusOneInv;
661 
662  // - Characteristic from boundary
663  cMinus = sqrt(gamma * m_pInf / m_rhoInf);
664  rMinus = VnInf[pnt] - 2.0 * cMinus * gammaMinusOneInv;
665  }
666  else
667  {
668  // + Characteristic from inside
669  cPlus = sqrt(gamma * m_pInf / m_rhoInf);
670  rPlus = VnInf[pnt] + 2.0 * cPlus * gammaMinusOneInv;
671 
672  // + Characteristic from inside
673  cMinus = sqrt(gamma * m_pInf / m_rhoInf);
674  rMinus = VnInf[pnt] - 2.0 * cPlus * gammaMinusOneInv;
675  }
676 
677  // Riemann boundary variables
678  VNBC = 0.5 * (rPlus + rMinus);
679  cBC = 0.25 * gammaMinusOne * (rPlus - rMinus);
680  VDBC = VNBC - VnInf[pnt];
681 
682  // Thermodynamic boundary variables
683  sBC = m_pInf / (pow(m_rhoInf, gamma));
684  rhoBC = pow((cBC * cBC) / (gamma * sBC), gammaMinusOneInv);
685  pBC = rhoBC * cBC * cBC * gammaInv;
686 
687  // Kinetic energy initialiasation
688  NekDouble EkBC = 0.0;
689 
690  // Boundary velocities
691  for ( j = 0; j < nDimensions; ++j)
692  {
693  velBC[j] = velInf[j] + VDBC * m_traceNormals[j][pnt];
694  rhoVelBC[j] = rhoBC * velBC[j];
695  EkBC += 0.5 * rhoBC * velBC[j]*velBC[j];
696  }
697 
698  // Boundary energy
699  EBC = pBC * gammaMinusOneInv + EkBC;
700 
701  // Imposing Riemann Invariant boundary conditions
702  (m_fields[0]->GetBndCondExpansions()[bcRegion]->
703  UpdatePhys())[id1+i] = rhoBC;
704  for (j = 0; j < nDimensions; ++j)
705  {
706  (m_fields[j+1]->GetBndCondExpansions()[bcRegion]->
707  UpdatePhys())[id1+i] = rhoVelBC[j];
708  }
709  (m_fields[nDimensions+1]->GetBndCondExpansions()[bcRegion]->
710  UpdatePhys())[id1+i] = EBC;
711 
712  }
713  else // Impose outflow Riemann invariant
714  {
715  // Subsonic flows
716  if (Mach[pnt] < 1.00)
717  {
718  // + Characteristic from inside
719  cPlus = sqrt(gamma * pressure[pnt] / Fwd[0][pnt]);
720  rPlus = Vn[pnt] + 2.0 * cPlus * gammaMinusOneInv;
721 
722  // - Characteristic from boundary
723  cMinus = sqrt(gamma * m_pInf / m_rhoInf);
724  rMinus = VnInf[pnt] - 2.0 * cMinus * gammaMinusOneInv;
725  }
726  else
727  {
728  // + Characteristic from inside
729  cPlus = sqrt(gamma * pressure[pnt] / Fwd[0][pnt]);
730  rPlus = Vn[pnt] + 2.0 * cPlus * gammaMinusOneInv;
731 
732  // + Characteristic from inside
733  cMinus = sqrt(gamma * pressure[pnt] / Fwd[0][pnt]);
734  rMinus = Vn[pnt] - 2.0 * cPlus * gammaMinusOneInv;
735  }
736 
737  // Riemann boundary variables
738  VNBC = 0.5 * (rPlus + rMinus);
739  cBC = 0.25 * gammaMinusOne * (rPlus - rMinus);
740  VDBC = VNBC - Vn[pnt];
741 
742  // Thermodynamic boundary variables
743  sBC = pressure[pnt] / (pow(Fwd[0][pnt], gamma));
744  rhoBC = pow((cBC * cBC) / (gamma * sBC), gammaMinusOneInv);
745  pBC = rhoBC * cBC * cBC * gammaInv;
746 
747  // Kinetic energy initialiasation
748  NekDouble EkBC = 0.0;
749 
750  // Boundary velocities
751  for ( j = 0; j < nDimensions; ++j)
752  {
753  velBC[j] = Fwd[j+1][pnt] / Fwd[0][pnt] +
754  VDBC * m_traceNormals[j][pnt];
755  rhoVelBC[j] = rhoBC * velBC[j];
756  EkBC += 0.5 * rhoBC * velBC[j]*velBC[j];
757  }
758 
759  // Boundary energy
760  EBC = pBC * gammaMinusOneInv + EkBC;
761 
762  // Imposing Riemann Invariant boundary conditions
763  (m_fields[0]->GetBndCondExpansions()[bcRegion]->
764  UpdatePhys())[id1+i] = rhoBC;
765  for (j = 0; j < nDimensions; ++j)
766  {
767  (m_fields[j+1]->GetBndCondExpansions()[bcRegion]->
768  UpdatePhys())[id1+i] = rhoVelBC[j];
769  }
770  (m_fields[nDimensions+1]->GetBndCondExpansions()[bcRegion]->
771  UpdatePhys())[id1+i] = EBC;
772  }
773  }
774  }
775  }
776 
777  /**
778  * @brief Extrapolation of order 0 for all the variables such that,
779  * at the boundaries, a trivial Riemann problem is solved.
780  */
781  void CompressibleFlowSystem::ExtrapOrder0BC(
782  int bcRegion,
783  int cnt,
784  Array<OneD, Array<OneD, NekDouble> > &physarray)
785  {
786  int i, j;
787  int e, pnt;
788  int id1, id2, nBCEdgePts;
789  int nTracePts = GetTraceTotPoints();
790  int nVariables = physarray.num_elements();
791  int nDimensions = m_spacedim;
792 
793  const Array<OneD, const int> &traceBndMap
794  = m_fields[0]->GetTraceBndMap();
795 
796  // Get physical values of the forward trace
797  Array<OneD, Array<OneD, NekDouble> > Fwd(nVariables);
798  for (i = 0; i < nVariables; ++i)
799  {
800  Fwd[i] = Array<OneD, NekDouble>(nTracePts);
801  m_fields[i]->ExtractTracePhys(physarray[i], Fwd[i]);
802  }
803 
804  int eMax;
805 
806  eMax = m_fields[0]->GetBndCondExpansions()[bcRegion]->GetExpSize();
807 
808  // Loop on bcRegions
809  for (e = 0; e < eMax; ++e)
810  {
811  nBCEdgePts = m_fields[0]->GetBndCondExpansions()[bcRegion]->
812  GetExp(e)->GetNumPoints(0);
813  id1 = m_fields[0]->GetBndCondExpansions()[bcRegion]->
814  GetPhys_Offset(e) ;
815  id2 = m_fields[0]->GetTrace()->GetPhys_Offset(traceBndMap[cnt+e]);
816 
817  // Loop on points of bcRegion 'e'
818  for (i = 0; i < nBCEdgePts; i++)
819  {
820  pnt = id2+i;
821 
822  // Setting up bcs for density
823  (m_fields[0]->GetBndCondExpansions()[bcRegion]->
824  UpdatePhys())[id1+i] = Fwd[0][pnt];
825 
826  // Setting up bcs for velocities
827  for (j = 1; j <=nDimensions; ++j)
828  {
829  (m_fields[j]->GetBndCondExpansions()[bcRegion]->
830  UpdatePhys())[id1+i] = Fwd[j][pnt];
831  }
832 
833  // Setting up bcs for energy
834  (m_fields[nVariables-1]->GetBndCondExpansions()[bcRegion]->
835  UpdatePhys())[id1+i] = Fwd[nVariables-1][pnt];
836  }
837  }
838  }
839 
840  /**
841  * @brief Return the flux vector for the compressible Euler equations.
842  *
843  * @param i Component of the flux vector to calculate.
844  * @param physfield Fields.
845  * @param flux Resulting flux.
846  */
847  void CompressibleFlowSystem::GetFluxVector(
848  const Array<OneD, Array<OneD, NekDouble> > &physfield,
849  Array<OneD, Array<OneD, Array<OneD, NekDouble> > > &flux)
850  {
851  int i, j;
852  int nq = physfield[0].num_elements();
853  int nVariables = m_fields.num_elements();
854 
855  Array<OneD, NekDouble> pressure(nq);
856  Array<OneD, Array<OneD, NekDouble> > velocity(m_spacedim);
857 
858  // Flux vector for the rho equation
859  for (i = 0; i < m_spacedim; ++i)
860  {
861  velocity[i] = Array<OneD, NekDouble>(nq);
862  Vmath::Vcopy(nq, physfield[i+1], 1, flux[0][i], 1);
863  }
864 
865  GetVelocityVector(physfield, velocity);
866  GetPressure (physfield, velocity, pressure);
867 
868  // Flux vector for the velocity fields
869  for (i = 0; i < m_spacedim; ++i)
870  {
871  for (j = 0; j < m_spacedim; ++j)
872  {
873  Vmath::Vmul(nq, velocity[j], 1, physfield[i+1], 1,
874  flux[i+1][j], 1);
875  }
876 
877  // Add pressure to appropriate field
878  Vmath::Vadd(nq, flux[i+1][i], 1, pressure, 1, flux[i+1][i], 1);
879  }
880 
881  // Flux vector for energy.
882  Vmath::Vadd(nq, physfield[m_spacedim+1], 1, pressure, 1,
883  pressure, 1);
884 
885  for (j = 0; j < m_spacedim; ++j)
886  {
887  Vmath::Vmul(nq, velocity[j], 1, pressure, 1,
888  flux[m_spacedim+1][j], 1);
889  }
890 
891  // For the smooth viscosity model
892  if (nVariables == m_spacedim+3)
893  {
894  // Add a zero row for the advective fluxes
895  for (j = 0; j < m_spacedim; ++j)
896  {
897  Vmath::Zero(nq, flux[m_spacedim+2][j], 1);
898  }
899  }
900  }
901 
902  /**
903  * @brief Return the flux vector for the compressible Euler equations
904  * by using the de-aliasing technique.
905  *
906  * @param i Component of the flux vector to calculate.
907  * @param physfield Fields.
908  * @param flux Resulting flux.
909  */
910  void CompressibleFlowSystem::GetFluxVectorDeAlias(
911  const Array<OneD, Array<OneD, NekDouble> > &physfield,
912  Array<OneD, Array<OneD, Array<OneD, NekDouble> > > &flux)
913  {
914  int i, j;
915  int nq = physfield[0].num_elements();
916  int nVariables = m_fields.num_elements();
917 
918  // Factor to rescale 1d points in dealiasing
919  NekDouble OneDptscale = 2;
920  nq = m_fields[0]->Get1DScaledTotPoints(OneDptscale);
921 
922  Array<OneD, NekDouble> pressure(nq);
923  Array<OneD, Array<OneD, NekDouble> > velocity(m_spacedim);
924 
925  Array<OneD, Array<OneD, NekDouble> > physfield_interp(nVariables);
926  Array<OneD, Array<OneD, Array<OneD, NekDouble> > > flux_interp(
927  nVariables);
928 
929  for (i = 0; i < nVariables; ++ i)
930  {
931  physfield_interp[i] = Array<OneD, NekDouble>(nq);
932  flux_interp[i] = Array<OneD, Array<OneD, NekDouble> >(m_spacedim);
933  m_fields[0]->PhysInterp1DScaled(
934  OneDptscale, physfield[i], physfield_interp[i]);
935 
936  for (j = 0; j < m_spacedim; ++j)
937  {
938  flux_interp[i][j] = Array<OneD, NekDouble>(nq);
939  }
940  }
941 
942  // Flux vector for the rho equation
943  for (i = 0; i < m_spacedim; ++i)
944  {
945  velocity[i] = Array<OneD, NekDouble>(nq);
946 
947  // Galerkin project solution back to original space
948  m_fields[0]->PhysGalerkinProjection1DScaled(
949  OneDptscale, physfield_interp[i+1], flux[0][i]);
950  }
951 
952  GetVelocityVector(physfield_interp, velocity);
953  GetPressure (physfield_interp, velocity, pressure);
954 
955  // Evaluation of flux vector for the velocity fields
956  for (i = 0; i < m_spacedim; ++i)
957  {
958  for (j = 0; j < m_spacedim; ++j)
959  {
960  Vmath::Vmul(nq, velocity[j], 1, physfield_interp[i+1], 1,
961  flux_interp[i+1][j], 1);
962  }
963 
964  // Add pressure to appropriate field
965  Vmath::Vadd(nq, flux_interp[i+1][i], 1, pressure,1,
966  flux_interp[i+1][i], 1);
967  }
968 
969  // Galerkin project solution back to origianl space
970  for (i = 0; i < m_spacedim; ++i)
971  {
972  for (j = 0; j < m_spacedim; ++j)
973  {
974  m_fields[0]->PhysGalerkinProjection1DScaled(
975  OneDptscale, flux_interp[i+1][j], flux[i+1][j]);
976  }
977  }
978 
979  // Evaluation of flux vector for energy
980  Vmath::Vadd(nq, physfield_interp[m_spacedim+1], 1, pressure, 1,
981  pressure, 1);
982 
983  for (j = 0; j < m_spacedim; ++j)
984  {
985  Vmath::Vmul(nq, velocity[j], 1, pressure, 1,
986  flux_interp[m_spacedim+1][j], 1);
987 
988  // Galerkin project solution back to origianl space
989  m_fields[0]->PhysGalerkinProjection1DScaled(
990  OneDptscale,
991  flux_interp[m_spacedim+1][j],
992  flux[m_spacedim+1][j]);
993  }
994  }
995 
996 
997  /**
998  * @brief Return the flux vector for the LDG diffusion problem.
999  * \todo Complete the viscous flux vector
1000  */
1001  void CompressibleFlowSystem::GetViscousFluxVector(
1002  const Array<OneD, Array<OneD, NekDouble> > &physfield,
1003  Array<OneD, Array<OneD, Array<OneD, NekDouble> > > &derivativesO1,
1004  Array<OneD, Array<OneD, Array<OneD, NekDouble> > > &viscousTensor)
1005  {
1006  int j, k;
1007  int nVariables = m_fields.num_elements();
1008  int nPts = physfield[0].num_elements();
1009 
1010  // Stokes hypotesis
1011  const NekDouble lambda = -2.0/3.0;
1012 
1013  // Auxiliary variables
1014  Array<OneD, NekDouble > mu (nPts, 0.0);
1015  Array<OneD, NekDouble > thermalConductivity(nPts, 0.0);
1016  Array<OneD, NekDouble > mu2 (nPts, 0.0);
1017  Array<OneD, NekDouble > divVel (nPts, 0.0);
1018 
1019  // Variable viscosity through the Sutherland's law
1020  if (m_ViscosityType == "Variable")
1021  {
1022  GetDynamicViscosity(physfield[nVariables-2], mu);
1023  NekDouble tRa = m_Cp / m_Prandtl;
1024  Vmath::Smul(nPts, tRa, &mu[0], 1, &thermalConductivity[0], 1);
1025  }
1026  else
1027  {
1028  Vmath::Fill(nPts, m_mu, &mu[0], 1);
1029  Vmath::Fill(nPts, m_thermalConductivity,
1030  &thermalConductivity[0], 1);
1031  }
1032 
1033  // Computing diagonal terms of viscous stress tensor
1034  Array<OneD, Array<OneD, NekDouble> > tmp(m_spacedim);
1035  Array<OneD, Array<OneD, NekDouble> > Sgg(m_spacedim);
1036 
1037  // mu2 = 2 * mu
1038  Vmath::Smul(nPts, 2.0, &mu[0], 1, &mu2[0], 1);
1039 
1040  // Velocity divergence
1041  for (j = 0; j < m_spacedim; ++j)
1042  {
1043  Vmath::Vadd(nPts, &divVel[0], 1, &derivativesO1[j][j][0], 1,
1044  &divVel[0], 1);
1045  }
1046 
1047  // Velocity divergence scaled by lambda * mu
1048  Vmath::Smul(nPts, lambda, &divVel[0], 1, &divVel[0], 1);
1049  Vmath::Vmul(nPts, &mu[0], 1, &divVel[0], 1, &divVel[0], 1);
1050 
1051  // Diagonal terms of viscous stress tensor (Sxx, Syy, Szz)
1052  // Sjj = 2 * mu * du_j/dx_j - (2 / 3) * mu * sum_j(du_j/dx_j)
1053  for (j = 0; j < m_spacedim; ++j)
1054  {
1055  tmp[j] = Array<OneD, NekDouble>(nPts, 0.0);
1056  Sgg[j] = Array<OneD, NekDouble>(nPts, 0.0);
1057 
1058  Vmath::Vmul(nPts, &mu2[0], 1, &derivativesO1[j][j][0], 1,
1059  &tmp[j][0], 1);
1060 
1061  Vmath::Vadd(nPts, &tmp[j][0], 1, &divVel[0], 1, &Sgg[j][0], 1);
1062  }
1063 
1064  // Extra diagonal terms of viscous stress tensor (Sxy, Sxz, Syz)
1065  // Note: they exist for 2D and 3D problems only
1066  Array<OneD, NekDouble > Sxy(nPts, 0.0);
1067  Array<OneD, NekDouble > Sxz(nPts, 0.0);
1068  Array<OneD, NekDouble > Syz(nPts, 0.0);
1069 
1070  if (m_spacedim == 2)
1071  {
1072  // Sxy = (du/dy + dv/dx)
1073  Vmath::Vadd(nPts, &derivativesO1[0][1][0], 1,
1074  &derivativesO1[1][0][0], 1, &Sxy[0], 1);
1075 
1076  // Sxy = mu * (du/dy + dv/dx)
1077  Vmath::Vmul(nPts, &mu[0], 1, &Sxy[0], 1, &Sxy[0], 1);
1078  }
1079  else if (m_spacedim == 3)
1080  {
1081  // Sxy = (du/dy + dv/dx)
1082  Vmath::Vadd(nPts, &derivativesO1[0][1][0], 1,
1083  &derivativesO1[1][0][0], 1, &Sxy[0], 1);
1084 
1085  // Sxz = (du/dz + dw/dx)
1086  Vmath::Vadd(nPts, &derivativesO1[0][2][0], 1,
1087  &derivativesO1[2][0][0], 1, &Sxz[0], 1);
1088 
1089  // Syz = (dv/dz + dw/dy)
1090  Vmath::Vadd(nPts, &derivativesO1[1][2][0], 1,
1091  &derivativesO1[2][1][0], 1, &Syz[0], 1);
1092 
1093  // Sxy = mu * (du/dy + dv/dx)
1094  Vmath::Vmul(nPts, &mu[0], 1, &Sxy[0], 1, &Sxy[0], 1);
1095 
1096  // Sxz = mu * (du/dy + dv/dx)
1097  Vmath::Vmul(nPts, &mu[0], 1, &Sxz[0], 1, &Sxz[0], 1);
1098 
1099  // Syz = mu * (du/dy + dv/dx)
1100  Vmath::Vmul(nPts, &mu[0], 1, &Syz[0], 1, &Syz[0], 1);
1101  }
1102 
1103  // Energy-related terms
1104  Array<OneD, NekDouble > STx(nPts, 0.0);
1105  Array<OneD, NekDouble > STy(nPts, 0.0);
1106  Array<OneD, NekDouble > STz(nPts, 0.0);
1107 
1108  // Building the viscous flux vector
1109 
1110  // Viscous flux vector for the rho equation
1111  for (k = 0; k < m_spacedim; ++k)
1112  {
1113  Vmath::Zero(nPts, viscousTensor[k][0], 1);
1114  }
1115 
1116  if (m_spacedim == 1)
1117  {
1118  Array<OneD, NekDouble > tmp1(nPts, 0.0);
1119 
1120  // u * Sxx
1121  Vmath::Vmul(nPts, &physfield[0][0], 1, &Sgg[0][0], 1, &STx[0], 1);
1122 
1123  // k * dT/dx
1124  Vmath::Vmul(nPts, &thermalConductivity[0], 1,
1125  &derivativesO1[0][1][0], 1, &tmp1[0], 1);
1126 
1127  // STx = u * Sxx + (K / mu) * dT/dx
1128  Vmath::Vadd(nPts, &STx[0], 1, &tmp1[0], 1, &STx[0], 1);
1129  }
1130  else if (m_spacedim == 2)
1131  {
1132  Array<OneD, NekDouble > tmp1(nPts, 0.0);
1133  Array<OneD, NekDouble > tmp2(nPts, 0.0);
1134 
1135  // Computation of STx
1136 
1137  // u * Sxx
1138  Vmath::Vmul(nPts, &physfield[0][0], 1, &Sgg[0][0], 1, &STx[0], 1);
1139 
1140  // v * Sxy
1141  Vmath::Vmul(nPts, &physfield[1][0], 1, &Sxy[0], 1, &tmp1[0], 1);
1142 
1143  // k * dT/dx
1144  Vmath::Vmul(nPts, &thermalConductivity[0], 1,
1145  &derivativesO1[0][2][0], 1, &tmp2[0], 1);
1146 
1147  // STx = u * Sxx + v * Sxy + K * dT/dx
1148  Vmath::Vadd(nPts, &STx[0], 1, &tmp1[0], 1, &STx[0], 1);
1149  Vmath::Vadd(nPts, &STx[0], 1, &tmp2[0], 1, &STx[0], 1);
1150 
1151  // Computation of STy
1152 
1153  // Re-initialise temporary arrays
1154  Vmath::Zero(nPts, &tmp1[0], 1);
1155  Vmath::Zero(nPts, &tmp2[0], 1);
1156 
1157  // v * Syy
1158  Vmath::Vmul(nPts, &physfield[1][0], 1, &Sgg[1][0], 1, &STy[0], 1);
1159 
1160  // u * Sxy
1161  Vmath::Vmul(nPts, &physfield[0][0], 1, &Sxy[0], 1, &tmp1[0], 1);
1162 
1163  // k * dT/dy
1164  Vmath::Vmul(nPts, &thermalConductivity[0], 1,
1165  &derivativesO1[1][2][0], 1, &tmp2[0], 1);
1166 
1167  // STy = v * Syy + u * Sxy + K * dT/dy
1168  Vmath::Vadd(nPts, &STy[0], 1, &tmp1[0], 1, &STy[0], 1);
1169  Vmath::Vadd(nPts, &STy[0], 1, &tmp2[0], 1, &STy[0], 1);
1170  }
1171  else if (m_spacedim == 3)
1172  {
1173  Array<OneD, NekDouble > tmp1(nPts, 0.0);
1174  Array<OneD, NekDouble > tmp2(nPts, 0.0);
1175  Array<OneD, NekDouble > tmp3(nPts, 0.0);
1176 
1177  // Computation of STx
1178 
1179  // u * Sxx
1180  Vmath::Vmul(nPts, &physfield[0][0], 1, &Sgg[0][0], 1, &STx[0], 1);
1181 
1182  // v * Sxy
1183  Vmath::Vmul(nPts, &physfield[1][0], 1, &Sxy[0], 1, &tmp1[0], 1);
1184 
1185  // v * Sxz
1186  Vmath::Vmul(nPts, &physfield[2][0], 1, &Sxz[0], 1, &tmp2[0], 1);
1187 
1188  // k * dT/dx
1189  Vmath::Vmul(nPts, &thermalConductivity[0], 1,
1190  &derivativesO1[0][3][0], 1, &tmp3[0], 1);
1191 
1192  // STx = u * Sxx + v * Sxy + w * Sxz + (K / mu) * dT/dx
1193  Vmath::Vadd(nPts, &STx[0], 1, &tmp1[0], 1, &STx[0], 1);
1194  Vmath::Vadd(nPts, &STx[0], 1, &tmp2[0], 1, &STx[0], 1);
1195  Vmath::Vadd(nPts, &STx[0], 1, &tmp3[0], 1, &STx[0], 1);
1196 
1197  // Computation of STy
1198 
1199  // Re-initialise temporary arrays
1200  Vmath::Zero(nPts, &tmp1[0], 1);
1201  Vmath::Zero(nPts, &tmp2[0], 1);
1202  Vmath::Zero(nPts, &tmp3[0], 1);
1203 
1204  // v * Syy
1205  Vmath::Vmul(nPts, &physfield[1][0], 1, &Sgg[1][0], 1, &STy[0], 1);
1206 
1207  // u * Sxy
1208  Vmath::Vmul(nPts, &physfield[0][0], 1, &Sxy[0], 1, &tmp1[0], 1);
1209 
1210  // w * Syz
1211  Vmath::Vmul(nPts, &physfield[2][0], 1, &Syz[0], 1, &tmp2[0], 1);
1212 
1213  // k * dT/dy
1214  Vmath::Vmul(nPts, &thermalConductivity[0], 1,
1215  &derivativesO1[1][3][0], 1, &tmp3[0], 1);
1216 
1217  // STy = v * Syy + u * Sxy + w * Syz + K * dT/dy
1218  Vmath::Vadd(nPts, &STy[0], 1, &tmp1[0], 1, &STy[0], 1);
1219  Vmath::Vadd(nPts, &STy[0], 1, &tmp2[0], 1, &STy[0], 1);
1220  Vmath::Vadd(nPts, &STy[0], 1, &tmp3[0], 1, &STy[0], 1);
1221 
1222  // Computation of STz
1223 
1224  // Re-initialise temporary arrays
1225  Vmath::Zero(nPts, &tmp1[0], 1);
1226  Vmath::Zero(nPts, &tmp2[0], 1);
1227  Vmath::Zero(nPts, &tmp3[0], 1);
1228 
1229  // w * Szz
1230  Vmath::Vmul(nPts, &physfield[2][0], 1, &Sgg[2][0], 1, &STz[0], 1);
1231 
1232  // u * Sxz
1233  Vmath::Vmul(nPts, &physfield[0][0], 1, &Sxz[0], 1, &tmp1[0], 1);
1234 
1235  // v * Syz
1236  Vmath::Vmul(nPts, &physfield[1][0], 1, &Syz[0], 1, &tmp2[0], 1);
1237 
1238  // k * dT/dz
1239  Vmath::Vmul(nPts, &thermalConductivity[0], 1,
1240  &derivativesO1[2][3][0], 1, &tmp3[0], 1);
1241 
1242  // STz = w * Szz + u * Sxz + v * Syz + K * dT/dz
1243  Vmath::Vadd(nPts, &STz[0], 1, &tmp1[0], 1, &STz[0], 1);
1244  Vmath::Vadd(nPts, &STz[0], 1, &tmp2[0], 1, &STz[0], 1);
1245  Vmath::Vadd(nPts, &STz[0], 1, &tmp3[0], 1, &STz[0], 1);
1246  }
1247 
1248  switch (m_spacedim)
1249  {
1250  case 1:
1251  {
1252  // f_11v = f_rho = 0
1253  Vmath::Zero(nPts, &viscousTensor[0][0][0], 1);
1254 
1255  // f_21v = f_rhou
1256  Vmath::Vcopy(nPts, &Sgg[0][0], 1, &viscousTensor[0][1][0], 1);
1257 
1258  // f_31v = f_E
1259  Vmath::Vcopy(nPts, &STx[0], 1, &viscousTensor[0][2][0], 1);
1260  break;
1261  }
1262  case 2:
1263  {
1264  // f_11v = f_rho1 = 0
1265  Vmath::Zero(nPts, &viscousTensor[0][0][0], 1);
1266  // f_12v = f_rho2 = 0
1267  Vmath::Zero(nPts, &viscousTensor[1][0][0], 1);
1268 
1269  // f_21v = f_rhou1
1270  Vmath::Vcopy(nPts, &Sgg[0][0], 1, &viscousTensor[0][1][0], 1);
1271  // f_22v = f_rhou2
1272  Vmath::Vcopy(nPts, &Sxy[0], 1, &viscousTensor[1][1][0], 1);
1273 
1274  // f_31v = f_rhov1
1275  Vmath::Vcopy(nPts, &Sxy[0], 1, &viscousTensor[0][2][0], 1);
1276  // f_32v = f_rhov2
1277  Vmath::Vcopy(nPts, &Sgg[1][0], 1, &viscousTensor[1][2][0], 1);
1278 
1279  // f_41v = f_E1
1280  Vmath::Vcopy(nPts, &STx[0], 1, &viscousTensor[0][3][0], 1);
1281  // f_42v = f_E2
1282  Vmath::Vcopy(nPts, &STy[0], 1, &viscousTensor[1][3][0], 1);
1283  break;
1284  }
1285  case 3:
1286  {
1287  // f_11v = f_rho1 = 0
1288  Vmath::Zero(nPts, &viscousTensor[0][0][0], 1);
1289  // f_12v = f_rho2 = 0
1290  Vmath::Zero(nPts, &viscousTensor[1][0][0], 1);
1291  // f_13v = f_rho3 = 0
1292  Vmath::Zero(nPts, &viscousTensor[2][0][0], 1);
1293 
1294  // f_21v = f_rhou1
1295  Vmath::Vcopy(nPts, &Sgg[0][0], 1, &viscousTensor[0][1][0], 1);
1296  // f_22v = f_rhou2
1297  Vmath::Vcopy(nPts, &Sxy[0], 1, &viscousTensor[1][1][0], 1);
1298  // f_23v = f_rhou3
1299  Vmath::Vcopy(nPts, &Sxz[0], 1, &viscousTensor[2][1][0], 1);
1300 
1301  // f_31v = f_rhov1
1302  Vmath::Vcopy(nPts, &Sxy[0], 1, &viscousTensor[0][2][0], 1);
1303  // f_32v = f_rhov2
1304  Vmath::Vcopy(nPts, &Sgg[1][0], 1, &viscousTensor[1][2][0], 1);
1305  // f_33v = f_rhov3
1306  Vmath::Vcopy(nPts, &Syz[0], 1, &viscousTensor[2][2][0], 1);
1307 
1308  // f_31v = f_rhow1
1309  Vmath::Vcopy(nPts, &Sxz[0], 1, &viscousTensor[0][3][0], 1);
1310  // f_32v = f_rhow2
1311  Vmath::Vcopy(nPts, &Syz[0], 1, &viscousTensor[1][3][0], 1);
1312  // f_33v = f_rhow3
1313  Vmath::Vcopy(nPts, &Sgg[2][0], 1, &viscousTensor[2][3][0], 1);
1314 
1315  // f_41v = f_E1
1316  Vmath::Vcopy(nPts, &STx[0], 1, &viscousTensor[0][4][0], 1);
1317  // f_42v = f_E2
1318  Vmath::Vcopy(nPts, &STy[0], 1, &viscousTensor[1][4][0], 1);
1319  // f_43v = f_E3
1320  Vmath::Vcopy(nPts, &STz[0], 1, &viscousTensor[2][4][0], 1);
1321  break;
1322  }
1323  default:
1324  {
1325  ASSERTL0(false, "Illegal expansion dimension");
1326  }
1327  }
1328  }
1329 
1330  /**
1331  * @brief Return the flux vector for the LDG diffusion problem.
1332  * \todo Complete the viscous flux vector
1333  */
1334  void CompressibleFlowSystem::GetViscousFluxVectorDeAlias(
1335  const Array<OneD, Array<OneD, NekDouble> > &physfield,
1336  Array<OneD, Array<OneD, Array<OneD, NekDouble> > > &derivativesO1,
1337  Array<OneD, Array<OneD, Array<OneD, NekDouble> > > &viscousTensor)
1338  {
1339 #if 0
1340  int i, j, k;
1341  int nVariables = m_fields.num_elements();
1342  int nPts = physfield[0].num_elements();
1343 
1344  int variables_phys = physfield.num_elements();
1345 
1346  // Factor to rescale 1d points in dealiasing.
1347  NekDouble OneDptscale = 2;
1348 
1349  // Get number of points to dealias a cubic non-linearity
1350  nPts = m_fields[0]->Get1DScaledTotPoints(OneDptscale);
1351 
1352  int nVariables_aux = derivativesO1[0].num_elements();
1353 
1354  Array<OneD, Array<OneD, NekDouble> > physfield_interp(variables_phys);
1355  Array<OneD, Array<OneD, Array<OneD, NekDouble> > > derivativesO1_interp(
1356  m_spacedim);
1357  Array<OneD, Array<OneD, Array<OneD, NekDouble> > >viscousTensor_interp(
1358  m_spacedim);
1359 
1360  for (i = 0; i < m_spacedim; ++ i)
1361  {
1362  viscousTensor_interp[i] = Array<OneD, Array<OneD, NekDouble> >(
1363  nVariables);
1364  for (j = 0; j < nVariables; ++j)
1365  {
1366  viscousTensor_interp[i][j] = Array<OneD, NekDouble>(nPts);
1367  }
1368  }
1369 
1370  // Stokes hypotesis
1371  NekDouble lambda = -2.0/3.0;
1372 
1373  // Auxiliary variables
1374  Array<OneD, NekDouble > mu (nPts, 0.0);
1375  Array<OneD, NekDouble > mu2 (nPts, 0.0);
1376  Array<OneD, NekDouble > divVel (nPts, 0.0);
1377  Array<OneD, NekDouble > pressure (nPts, 0.0);
1378  Array<OneD, NekDouble > temperature(nPts, 0.0);
1379 
1380  for (i = 0; i < nVariables; ++i)
1381  {
1382  m_fields[0]->PhysInterp1DScaled(
1383  OneDptscale, physfield[i], fields_interp[i]);
1384  }
1385 
1386  for (i = 0; i < variables_phys; ++i)
1387  {
1388  physfield_interp[i] = Array<OneD, NekDouble> (nPts);
1389 
1390  // Interpolation to higher space
1391  m_fields[0]->PhysInterp1DScaled(OneDptscale, physfield[i],
1392  physfield_interp[i]);
1393  }
1394 
1395  for (i = 0; i < m_spacedim; ++i)
1396  {
1397  derivativesO1_interp[i] = Array<OneD, Array<OneD, NekDouble> >(
1398  nVariables_aux);
1399  for (j = 0; j < nVariables_aux; ++j)
1400  {
1401  derivativesO1_interp[i][j] = Array<OneD, NekDouble>(nPts);
1402  m_fields[0]->PhysInterp1DScaled(
1403  OneDptscale, derivativesO1[i][j],
1404  derivativesO1_interp[i][j]);
1405  }
1406  }
1407 
1408  // Thermodynamic related quantities
1409  GetPressure(fields_interp, pressure);
1410  GetTemperature(fields_interp, pressure, temperature);
1411 
1412  // Variable viscosity through the Sutherland's law
1413  if (m_ViscosityType == "Variable")
1414  {
1415  GetDynamicViscosity(fields_interp[variables_phys-1], mu);
1416  }
1417  else
1418  {
1419  Vmath::Sadd(nPts, m_mu, &mu[0], 1, &mu[0], 1);
1420  }
1421 
1422  // Computing diagonal terms of viscous stress tensor
1423  Array<OneD, Array<OneD, NekDouble> > tmp(m_spacedim);
1424  Array<OneD, Array<OneD, NekDouble> > Sgg(m_spacedim);
1425 
1426  // mu2 = 2 * mu
1427  Vmath::Smul(nPts, 2.0, &mu[0], 1, &mu2[0], 1);
1428 
1429  // Velocity divergence
1430  for (j = 0; j < m_spacedim; ++j)
1431  {
1432  Vmath::Vadd(nPts, &divVel[0], 1, &derivativesO1_interp[j][j][0], 1,
1433  &divVel[0], 1);
1434  }
1435 
1436  // Velocity divergence scaled by lambda * mu
1437  Vmath::Smul(nPts, lambda, &divVel[0], 1, &divVel[0], 1);
1438  Vmath::Vmul(nPts, &mu[0], 1, &divVel[0], 1, &divVel[0], 1);
1439 
1440  // Digonal terms of viscous stress tensor (Sxx, Syy, Szz)
1441  // Sjj = 2 * mu * du_j/dx_j - (2 / 3) * mu * sum_j(du_j/dx_j)
1442  for (j = 0; j < m_spacedim; ++j)
1443  {
1444  tmp[j] = Array<OneD, NekDouble>(nPts, 0.0);
1445  Sgg[j] = Array<OneD, NekDouble>(nPts, 0.0);
1446 
1447  Vmath::Vmul(nPts, &mu2[0], 1, &derivativesO1_interp[j][j][0], 1,
1448  &tmp[j][0], 1);
1449 
1450  Vmath::Vadd(nPts, &tmp[j][0], 1, &divVel[0], 1, &Sgg[j][0], 1);
1451  }
1452 
1453  // Extra diagonal terms of viscous stress tensor (Sxy, Sxz, Syz)
1454  // Note: they exist for 2D and 3D problems only
1455  Array<OneD, NekDouble > Sxy(nPts, 0.0);
1456  Array<OneD, NekDouble > Sxz(nPts, 0.0);
1457  Array<OneD, NekDouble > Syz(nPts, 0.0);
1458 
1459  if (m_spacedim == 2)
1460  {
1461  // Sxy = (du/dy + dv/dx)
1462  Vmath::Vadd(nPts, &derivativesO1_interp[0][1][0], 1,
1463  &derivativesO1_interp[1][0][0], 1, &Sxy[0], 1);
1464 
1465  // Sxy = mu * (du/dy + dv/dx)
1466  Vmath::Vmul(nPts, &mu[0], 1, &Sxy[0], 1, &Sxy[0], 1);
1467  }
1468  else if (m_spacedim == 3)
1469  {
1470  // Sxy = (du/dy + dv/dx)
1471  Vmath::Vadd(nPts, &derivativesO1_interp[0][1][0], 1,
1472  &derivativesO1_interp[1][0][0], 1, &Sxy[0], 1);
1473 
1474  // Sxz = (du/dz + dw/dx)
1475  Vmath::Vadd(nPts, &derivativesO1_interp[0][2][0], 1,
1476  &derivativesO1_interp[2][0][0], 1, &Sxz[0], 1);
1477 
1478  // Syz = (dv/dz + dw/dy)
1479  Vmath::Vadd(nPts, &derivativesO1_interp[1][2][0], 1,
1480  &derivativesO1_interp[2][1][0], 1, &Syz[0], 1);
1481 
1482  // Sxy = mu * (du/dy + dv/dx)
1483  Vmath::Vmul(nPts, &mu[0], 1, &Sxy[0], 1, &Sxy[0], 1);
1484 
1485  // Sxz = mu * (du/dy + dv/dx)
1486  Vmath::Vmul(nPts, &mu[0], 1, &Sxz[0], 1, &Sxz[0], 1);
1487 
1488  // Syz = mu * (du/dy + dv/dx)
1489  Vmath::Vmul(nPts, &mu[0], 1, &Syz[0], 1, &Syz[0], 1);
1490  }
1491 
1492  // Energy-related terms
1493  Array<OneD, NekDouble > STx(nPts, 0.0);
1494  Array<OneD, NekDouble > STy(nPts, 0.0);
1495  Array<OneD, NekDouble > STz(nPts, 0.0);
1496  // Building the viscous flux vector
1497  if (i == 0)
1498  {
1499  // Viscous flux vector for the rho equation
1500  for (k = 0; k < m_spacedim; ++k)
1501  {
1502  Vmath::Zero(nPts, viscousTensor_interp[k][i], 1);
1503  }
1504  }
1505 
1506  if (m_spacedim == 1)
1507  {
1508  Array<OneD, NekDouble > tmp1(nPts, 0.0);
1509 
1510  // u * Sxx
1511  Vmath::Vmul(nPts, &physfield_interp[0][0], 1,
1512  &Sgg[0][0], 1, &STx[0], 1);
1513 
1514  // k * dT/dx
1515  Vmath::Smul(nPts, m_thermalConductivity,
1516  &derivativesO1_interp[0][1][0], 1, &tmp1[0], 1);
1517 
1518  // STx = u * Sxx + (K / mu) * dT/dx
1519  Vmath::Vadd(nPts, &STx[0], 1, &tmp1[0], 1, &STx[0], 1);
1520  }
1521  else if (m_spacedim == 2)
1522  {
1523  Array<OneD, NekDouble > tmp1(nPts, 0.0);
1524  Array<OneD, NekDouble > tmp2(nPts, 0.0);
1525 
1526  // Computation of STx
1527 
1528  // u * Sxx
1529  Vmath::Vmul(nPts, &physfield_interp[0][0], 1,
1530  &Sgg[0][0], 1, &STx[0], 1);
1531 
1532  // v * Sxy
1533  Vmath::Vmul(nPts, &physfield_interp[1][0], 1,
1534  &Sxy[0], 1, &tmp1[0], 1);
1535 
1536  // k * dT/dx
1537  Vmath::Smul(nPts, m_thermalConductivity,
1538  &derivativesO1_interp[0][2][0], 1, &tmp2[0], 1);
1539 
1540  // STx = u * Sxx + v * Sxy + K * dT/dx
1541  Vmath::Vadd(nPts, &STx[0], 1, &tmp1[0], 1, &STx[0], 1);
1542  Vmath::Vadd(nPts, &STx[0], 1, &tmp2[0], 1, &STx[0], 1);
1543 
1544  // Computation of STy
1545 
1546  // Re-initialise temporary arrays
1547  Vmath::Zero(nPts, &tmp1[0], 1);
1548  Vmath::Zero(nPts, &tmp2[0], 1);
1549 
1550  // v * Syy
1551  Vmath::Vmul(nPts, &physfield_interp[1][0], 1,
1552  &Sgg[1][0], 1, &STy[0], 1);
1553 
1554  // u * Sxy
1555  Vmath::Vmul(nPts, &physfield_interp[0][0], 1,
1556  &Sxy[0], 1, &tmp1[0], 1);
1557 
1558  // k * dT/dy
1559  Vmath::Smul(nPts, m_thermalConductivity,
1560  &derivativesO1_interp[1][2][0], 1, &tmp2[0], 1);
1561 
1562  // STy = v * Syy + u * Sxy + K * dT/dy
1563  Vmath::Vadd(nPts, &STy[0], 1, &tmp1[0], 1, &STy[0], 1);
1564  Vmath::Vadd(nPts, &STy[0], 1, &tmp2[0], 1, &STy[0], 1);
1565  }
1566  else if (m_spacedim == 3)
1567  {
1568  Array<OneD, NekDouble > tmp1(nPts, 0.0);
1569  Array<OneD, NekDouble > tmp2(nPts, 0.0);
1570  Array<OneD, NekDouble > tmp3(nPts, 0.0);
1571 
1572  // Computation of STx
1573 
1574  // u * Sxx
1575  Vmath::Vmul(nPts, &physfield_interp[0][0], 1,
1576  &Sgg[0][0], 1, &STx[0], 1);
1577 
1578  // v * Sxy
1579  Vmath::Vmul(nPts, &physfield_interp[1][0], 1,
1580  &Sxy[0], 1, &tmp1[0], 1);
1581 
1582  // v * Sxy
1583  Vmath::Vmul(nPts, &physfield_interp[2][0], 1,
1584  &Sxz[0], 1, &tmp2[0], 1);
1585 
1586  // k * dT/dx
1587  Vmath::Smul(nPts, m_thermalConductivity,
1588  &derivativesO1_interp[0][3][0], 1, &tmp3[0], 1);
1589 
1590  // STx = u * Sxx + v * Sxy + w * Sxz + (K / mu) * dT/dx
1591  Vmath::Vadd(nPts, &STx[0], 1, &tmp1[0], 1, &STx[0], 1);
1592  Vmath::Vadd(nPts, &STx[0], 1, &tmp2[0], 1, &STx[0], 1);
1593  Vmath::Vadd(nPts, &STx[0], 1, &tmp3[0], 1, &STx[0], 1);
1594 
1595  // Computation of STy
1596 
1597  // Re-initialise temporary arrays
1598  Vmath::Zero(nPts, &tmp1[0], 1);
1599  Vmath::Zero(nPts, &tmp2[0], 1);
1600  Vmath::Zero(nPts, &tmp3[0], 1);
1601 
1602  // v * Syy
1603  Vmath::Vmul(nPts, &physfield_interp[1][0], 1,
1604  &Sgg[1][0], 1, &STy[0], 1);
1605 
1606  // u * Sxy
1607  Vmath::Vmul(nPts, &physfield_interp[0][0], 1,
1608  &Sxy[0], 1, &tmp1[0], 1);
1609 
1610  // w * Syz
1611  Vmath::Vmul(nPts, &physfield_interp[2][0], 1,
1612  &Syz[0], 1, &tmp2[0], 1);
1613 
1614  // k * dT/dy
1615  Vmath::Smul(nPts, m_thermalConductivity,
1616  &derivativesO1_interp[1][3][0], 1, &tmp3[0], 1);
1617 
1618  // STy = v * Syy + u * Sxy + w * Syz + K * dT/dy
1619  Vmath::Vadd(nPts, &STy[0], 1, &tmp1[0], 1, &STy[0], 1);
1620  Vmath::Vadd(nPts, &STy[0], 1, &tmp2[0], 1, &STy[0], 1);
1621  Vmath::Vadd(nPts, &STy[0], 1, &tmp3[0], 1, &STy[0], 1);
1622 
1623  // Computation of STz
1624 
1625  // Re-initialise temporary arrays
1626  Vmath::Zero(nPts, &tmp1[0], 1);
1627  Vmath::Zero(nPts, &tmp2[0], 1);
1628  Vmath::Zero(nPts, &tmp3[0], 1);
1629 
1630  // w * Szz
1631  Vmath::Vmul(nPts, &physfield_interp[2][0], 1,
1632  &Sgg[2][0], 1, &STz[0], 1);
1633 
1634  // u * Sxz
1635  Vmath::Vmul(nPts, &physfield_interp[0][0], 1,
1636  &Sxz[0], 1, &tmp1[0], 1);
1637 
1638  // v * Syz
1639  Vmath::Vmul(nPts, &physfield_interp[1][0], 1,
1640  &Syz[0], 1, &tmp2[0], 1);
1641 
1642  // k * dT/dz
1643  Vmath::Smul(nPts, m_thermalConductivity,
1644  &derivativesO1_interp[2][3][0], 1, &tmp3[0], 1);
1645 
1646  // STz = w * Szz + u * Sxz + v * Syz + K * dT/dz
1647  Vmath::Vadd(nPts, &STz[0], 1, &tmp1[0], 1, &STz[0], 1);
1648  Vmath::Vadd(nPts, &STz[0], 1, &tmp2[0], 1, &STz[0], 1);
1649  Vmath::Vadd(nPts, &STz[0], 1, &tmp3[0], 1, &STz[0], 1);
1650  }
1651 
1652  switch (m_spacedim)
1653  {
1654  case 1:
1655  {
1656 
1657  int nq = physfield[0].num_elements();
1658  // f_11v = f_rho = 0
1659  Vmath::Zero(nq, &viscousTensor_interp[0][0][0], 1);
1660 
1661  // f_21v = f_rhou
1662  Vmath::Vcopy(nPts, &Sgg[0][0], 1, &viscousTensor_interp[0][1][0], 1);
1663 
1664  // f_31v = f_E
1665  Vmath::Vcopy(nPts, &STx[0], 1, &viscousTensor_interp[0][2][0], 1);
1666  break;
1667  }
1668  case 2:
1669  {
1670  int nq = physfield[0].num_elements();
1671  // f_11v = f_rho1 = 0
1672  Vmath::Zero(nq, &viscousTensor_interp[0][0][0], 1);
1673  // f_12v = f_rho2 = 0
1674  Vmath::Zero(nq, &viscousTensor_interp[1][0][0], 1);
1675 
1676  // f_21v = f_rhou1
1677  Vmath::Vcopy(nPts, &Sgg[0][0], 1, &viscousTensor_interp[0][1][0], 1);
1678  // f_22v = f_rhou2
1679  Vmath::Vcopy(nPts, &Sxy[0], 1, &viscousTensor_interp[1][1][0], 1);
1680 
1681  // f_31v = f_rhov1
1682  Vmath::Vcopy(nPts, &Sxy[0], 1, &viscousTensor_interp[0][2][0], 1);
1683  // f_32v = f_rhov2
1684  Vmath::Vcopy(nPts, &Sgg[1][0], 1, &viscousTensor_interp[1][2][0], 1);
1685 
1686  // f_41v = f_E1
1687  Vmath::Vcopy(nPts, &STx[0], 1, &viscousTensor_interp[0][3][0], 1);
1688  // f_42v = f_E2
1689  Vmath::Vcopy(nPts, &STy[0], 1, &viscousTensor_interp[1][3][0], 1);
1690  break;
1691  }
1692  case 3:
1693  {
1694  int nq = physfield[0].num_elements();
1695  // f_11v = f_rho1 = 0
1696  Vmath::Zero(nq, &viscousTensor_interp[0][0][0], 1);
1697  // f_12v = f_rho2 = 0
1698  Vmath::Zero(nq, &viscousTensor_interp[1][0][0], 1);
1699  // f_13v = f_rho3 = 0
1700  Vmath::Zero(nq, &viscousTensor_interp[2][0][0], 1);
1701 
1702  // f_21v = f_rhou1
1703  Vmath::Vcopy(nPts, &Sgg[0][0], 1, &viscousTensor_interp[0][1][0], 1);
1704  // f_22v = f_rhou2
1705  Vmath::Vcopy(nPts, &Sxy[0], 1, &viscousTensor_interp[1][1][0], 1);
1706  // f_23v = f_rhou3
1707  Vmath::Vcopy(nPts, &Sxz[0], 1, &viscousTensor_interp[2][1][0], 1);
1708 
1709  // f_31v = f_rhov1
1710  Vmath::Vcopy(nPts, &Sxy[0], 1, &viscousTensor_interp[0][2][0], 1);
1711  // f_32v = f_rhov2
1712  Vmath::Vcopy(nPts, &Sgg[1][0], 1, &viscousTensor_interp[1][2][0], 1);
1713  // f_33v = f_rhov3
1714  Vmath::Vcopy(nPts, &Syz[0], 1, &viscousTensor_interp[2][2][0], 1);
1715 
1716  // f_31v = f_rhow1
1717  Vmath::Vcopy(nPts, &Sxz[0], 1, &viscousTensor_interp[0][3][0], 1);
1718  // f_32v = f_rhow2
1719  Vmath::Vcopy(nPts, &Syz[0], 1, &viscousTensor_interp[1][3][0], 1);
1720  // f_33v = f_rhow3
1721  Vmath::Vcopy(nPts, &Sgg[2][0], 1, &viscousTensor_interp[2][3][0], 1);
1722 
1723  // f_41v = f_E1
1724  Vmath::Vcopy(nPts, &STx[0], 1, &viscousTensor_interp[0][4][0], 1);
1725  // f_42v = f_E2
1726  Vmath::Vcopy(nPts, &STy[0], 1, &viscousTensor_interp[1][4][0], 1);
1727  // f_43v = f_E3
1728  Vmath::Vcopy(nPts, &STz[0], 1, &viscousTensor_interp[2][4][0], 1);
1729  break;
1730  }
1731  default:
1732  {
1733  ASSERTL0(false, "Illegal expansion dimension");
1734  }
1735  }
1736 
1737  for (i = 0; i < m_spacedim; ++i)
1738  {
1739  for (j = 1; j < nVariables; ++j)
1740  {
1741  // Galerkin project solution back to origianl space
1742  m_fields[0]->PhysGalerkinProjection1DScaled(
1743  OneDptscale,
1744  viscousTensor_interp[i][j],
1745  viscousTensor[i][j]);
1746  }
1747  }
1748 #endif
1749 }
1750 
1751  /**
1752  * @brief Calculate the pressure field \f$ p =
1753  * (\gamma-1)(E-\frac{1}{2}\rho\| \mathbf{v} \|^2) \f$ assuming an ideal
1754  * gas law.
1755  *
1756  * @param physfield Input momentum.
1757  * @param pressure Computed pressure field.
1758  */
1759  void CompressibleFlowSystem::GetPressure(
1760  const Array<OneD, const Array<OneD, NekDouble> > &physfield,
1761  Array<OneD, NekDouble> &pressure)
1762  {
1763  int nBCEdgePts = physfield[0].num_elements();
1764  NekDouble alpha = -0.5;
1765 
1766  // Calculate ||rho v||^2
1767  Vmath::Vmul(nBCEdgePts, physfield[1], 1, physfield[1], 1, pressure, 1);
1768  for (int i = 1; i < m_spacedim; ++i)
1769  {
1770  Vmath::Vvtvp(nBCEdgePts, physfield[1+i], 1, physfield[1+i], 1,
1771  pressure, 1, pressure, 1);
1772  }
1773  // Divide by rho to get rho*||v||^2
1774  Vmath::Vdiv (nBCEdgePts, pressure, 1, physfield[0], 1, pressure, 1);
1775  // pressure <- E - 0.5*pressure
1776  Vmath::Svtvp(nBCEdgePts, alpha,
1777  pressure, 1, physfield[m_spacedim+1], 1, pressure, 1);
1778  // Multiply by (gamma-1)
1779  Vmath::Smul (nBCEdgePts, m_gamma-1, pressure, 1, pressure, 1);
1780  }
1781 
1782  /**
1783  * @brief Compute the enthalpy term \f$ H = E + p/rho \$.
1784  */
1785  void CompressibleFlowSystem::GetEnthalpy(
1786  const Array<OneD, const Array<OneD, NekDouble> > &physfield,
1787  Array<OneD, NekDouble> &pressure,
1788  Array<OneD, NekDouble> &enthalpy)
1789  {
1790  int npts = m_fields[0]->GetTotPoints();
1791  Array<OneD, NekDouble> tmp(npts, 0.0);
1792 
1793  // Calculate E = rhoE/rho
1794  Vmath::Vdiv(npts, physfield[m_spacedim+1], 1, physfield[0], 1, tmp, 1);
1795  // Calculate p/rho
1796  Vmath::Vdiv(npts, pressure, 1, physfield[0], 1, enthalpy, 1);
1797  // Calculate H = E + p/rho
1798  Vmath::Vadd(npts, tmp, 1, enthalpy, 1, enthalpy, 1);
1799  }
1800 
1801  /**
1802  * @brief Calculate the pressure field \f$ p =
1803  * (\gamma-1)(E-\frac{1}{2}\rho\| \mathbf{v} \|^2) \f$ assuming an ideal
1804  * gas law.
1805  *
1806  * This is a slightly optimised way to calculate the pressure field which
1807  * avoids division by the density field if the velocity field has already
1808  * been calculated.
1809  *
1810  * @param physfield Input momentum.
1811  * @param velocity Velocity vector.
1812  * @param pressure Computed pressure field.
1813  */
1814  void CompressibleFlowSystem::GetPressure(
1815  const Array<OneD, const Array<OneD, NekDouble> > &physfield,
1816  const Array<OneD, const Array<OneD, NekDouble> > &velocity,
1817  Array<OneD, NekDouble> &pressure)
1818  {
1819  int nBCEdgePts = physfield[0].num_elements();
1820  NekDouble alpha = -0.5;
1821 
1822  // Calculate ||\rho v||^2.
1823  Vmath::Vmul (nBCEdgePts, velocity[0], 1, physfield[1], 1, pressure, 1);
1824  for (int i = 1; i < m_spacedim; ++i)
1825  {
1826  Vmath::Vvtvp(nBCEdgePts, velocity[i], 1, physfield[1+i], 1,
1827  pressure, 1, pressure, 1);
1828  }
1829  // pressure <- E - 0.5*pressure
1830  Vmath::Svtvp(nBCEdgePts, alpha,
1831  pressure, 1, physfield[m_spacedim+1], 1, pressure, 1);
1832  // Multiply by (gamma-1).
1833  Vmath::Smul (nBCEdgePts, m_gamma-1, pressure, 1, pressure, 1);
1834  }
1835 
1836  /**
1837  * @brief Compute the velocity field \f$ \mathbf{v} \f$ given the momentum
1838  * \f$ \rho\mathbf{v} \f$.
1839  *
1840  * @param physfield Momentum field.
1841  * @param velocity Velocity field.
1842  */
1843  void CompressibleFlowSystem::GetVelocityVector(
1844  const Array<OneD, Array<OneD, NekDouble> > &physfield,
1845  Array<OneD, Array<OneD, NekDouble> > &velocity)
1846  {
1847  const int nBCEdgePts = physfield[0].num_elements();
1848 
1849  for (int i = 0; i < m_spacedim; ++i)
1850  {
1851  Vmath::Vdiv(nBCEdgePts, physfield[1+i], 1, physfield[0], 1,
1852  velocity[i], 1);
1853  }
1854  }
1855 
1856  /**
1857  * @brief Compute the temperature \f$ T = p/\rho R \f$.
1858  *
1859  * @param physfield Input physical field.
1860  * @param pressure The pressure field corresponding to physfield.
1861  * @param temperature The resulting temperature \f$ T \f$.
1862  */
1863  void CompressibleFlowSystem::GetTemperature(
1864  const Array<OneD, const Array<OneD, NekDouble> > &physfield,
1865  Array<OneD, NekDouble > &pressure,
1866  Array<OneD, NekDouble > &temperature)
1867  {
1868  const int nq = physfield[0].num_elements();
1869 
1870  Vmath::Vdiv(nq, pressure, 1, physfield[0], 1, temperature, 1);
1871  Vmath::Smul(nq, 1.0/m_gasConstant, temperature, 1, temperature, 1);
1872  }
1873 
1874  /**
1875  * @brief Compute the sound speed \f$ c = sqrt(\gamma p/\rho) \f$.
1876  *
1877  * @param physfield Input physical field.
1878  * @param pressure The pressure field corresponding to physfield.
1879  * @param soundspeed The resulting sound speed \f$ c \f$.
1880  */
1881  void CompressibleFlowSystem::GetSoundSpeed(
1882  const Array<OneD, Array<OneD, NekDouble> > &physfield,
1883  Array<OneD, NekDouble > &pressure,
1884  Array<OneD, NekDouble > &soundspeed)
1885  {
1886  const int nq = m_fields[0]->GetTotPoints();
1887  Vmath::Vdiv (nq, pressure, 1, physfield[0], 1, soundspeed, 1);
1888  Vmath::Smul (nq, m_gamma, soundspeed, 1, soundspeed, 1);
1889  Vmath::Vsqrt(nq, soundspeed, 1, soundspeed, 1);
1890  }
1891 
1892  /**
1893  * @brief Compute the mach number \f$ M = \| \mathbf{v} \|^2 / c \f$.
1894  *
1895  * @param physfield Input physical field.
1896  * @param soundfield The speed of sound corresponding to physfield.
1897  * @param mach The resulting mach number \f$ M \f$.
1898  */
1899  void CompressibleFlowSystem::GetMach(
1900  Array<OneD, Array<OneD, NekDouble> > &physfield,
1901  Array<OneD, NekDouble > &soundspeed,
1902  Array<OneD, NekDouble > &mach)
1903  {
1904  const int nq = m_fields[0]->GetTotPoints();
1905 
1906  Vmath::Vmul(nq, physfield[1], 1, physfield[1], 1, mach, 1);
1907 
1908  for (int i = 1; i < m_spacedim; ++i)
1909  {
1910  Vmath::Vvtvp(nq, physfield[1+i], 1, physfield[1+i], 1,
1911  mach, 1, mach, 1);
1912  }
1913 
1914  Vmath::Vdiv(nq, mach, 1, physfield[0], 1, mach, 1);
1915  Vmath::Vdiv(nq, mach, 1, physfield[0], 1, mach, 1);
1916  Vmath::Vsqrt(nq, mach, 1, mach, 1);
1917 
1918  Vmath::Vdiv(nq, mach, 1, soundspeed, 1, mach, 1);
1919  }
1920 
1921  /**
1922  * @brief Compute the dynamic viscosity using the Sutherland's law
1923  * \f$ \mu = \mu_star * (T / T_star)^3/2 * (T_star + 110) / (T + 110) \f$,
1924  * where: \mu_star = 1.7894 * 10^-5 Kg / (m * s)
1925  * T_star = 288.15 K
1926  *
1927  * @param physfield Input physical field.
1928  * @param mu The resulting dynamic viscosity.
1929  */
1930  void CompressibleFlowSystem::GetDynamicViscosity(
1931  const Array<OneD, const NekDouble> &temperature,
1932  Array<OneD, NekDouble> &mu)
1933  {
1934  const int nPts = temperature.num_elements();
1935  NekDouble mu_star = m_mu;
1936  NekDouble T_star = m_pInf / (m_rhoInf * m_gasConstant);
1937  NekDouble ratio;
1938 
1939  for (int i = 0; i < nPts; ++i)
1940  {
1941  ratio = temperature[i] / T_star;
1942  mu[i] = mu_star * ratio * sqrt(ratio) *
1943  (T_star + 110.0) / (temperature[i] + 110.0);
1944  }
1945  }
1946 
1947  /**
1948  * @brief Calcualte entropy.
1949  */
1950  void CompressibleFlowSystem::GetEntropy(
1951  const Array<OneD, const Array<OneD, NekDouble> > &physfield,
1952  const Array<OneD, const NekDouble> &pressure,
1953  const Array<OneD, const NekDouble> &temperature,
1954  Array<OneD, NekDouble> &entropy)
1955  {
1956  NekDouble entropy_sum = 0.0;
1957  const int npts = m_fields[0]->GetTotPoints();
1958  const NekDouble temp_inf = m_pInf/(m_rhoInf*m_gasConstant);;
1959  Array<OneD, NekDouble> L2entropy(npts, 0.0);
1960 
1961  for (int i = 0; i < npts; ++i)
1962  {
1963  entropy[i] = m_gamma/(m_gamma-1.0)*m_gasConstant*log(temperature[i]/temp_inf) -
1964  m_gasConstant*log(pressure[i]/m_pInf);
1965  }
1966 
1967  Vmath::Vmul(npts,entropy,1,entropy,1,L2entropy,1);
1968 
1969  entropy_sum = Vmath::Vsum(npts, L2entropy, 1);
1970 
1971  entropy_sum = sqrt(entropy_sum);
1972 
1973  std::ofstream m_file( "L2entropy.txt", std::ios_base::app);
1974 
1975  m_file << setprecision(16) << scientific << entropy_sum << endl;
1976  //m_file << Vmath::Vmax(entropy.num_elements(),entropy,1) << endl;
1977 
1978  m_file.close();
1979  }
1980 
1981  /**
1982  * @brief Calculate the maximum timestep subject to CFL restrictions.
1983  */
1984  NekDouble CompressibleFlowSystem::v_GetTimeStep(
1985  const Array<OneD, const Array<OneD, NekDouble> > &inarray)
1986  {
1987  int n;
1988  int nElements = m_fields[0]->GetExpSize();
1989  const Array<OneD, int> ExpOrder = GetNumExpModesPerExp();
1990 
1991  Array<OneD, NekDouble> tstep (nElements, 0.0);
1992  Array<OneD, NekDouble> stdVelocity(nElements);
1993 
1994  // Get standard velocity to compute the time-step limit
1995  GetStdVelocity(inarray, stdVelocity);
1996 
1997  // Factors to compute the time-step limit
1998  NekDouble minLength = 0.0;
1999  NekDouble alpha = MaxTimeStepEstimator();
2000  NekDouble cLambda = 0.2; // Spencer book-317
2001 
2002  // Loop over elements to compute the time-step limit for each element
2003  for(n = 0; n < nElements; ++n)
2004  {
2005  int npoints = m_fields[0]->GetExp(n)->GetTotPoints();
2006  Array<OneD, NekDouble> one2D(npoints, 1.0);
2007  NekDouble Area = m_fields[0]->GetExp(n)->Integral(one2D);
2008 
2009  minLength = sqrt(Area);
2010  if (m_fields[0]->GetExp(n)->as<LocalRegions::TriExp>())
2011  {
2012  minLength *= 2.0;
2013  }
2014 
2015  tstep[n] = m_cflSafetyFactor * alpha * minLength
2016  / (stdVelocity[n] * cLambda
2017  * (ExpOrder[n] - 1) * (ExpOrder[n] - 1));
2018  }
2019 
2020  // Get the minimum time-step limit and return the time-step
2021  NekDouble TimeStep = Vmath::Vmin(nElements, tstep, 1);
2022  m_comm->AllReduce(TimeStep, LibUtilities::ReduceMin);
2023  return TimeStep;
2024  }
2025 
2026  /**
2027  * @brief Compute the advection velocity in the standard space
2028  * for each element of the expansion.
2029  *
2030  * @param inarray Momentum field.
2031  * @param stdV Standard velocity field.
2032  */
2033  void CompressibleFlowSystem::GetStdVelocity(
2034  const Array<OneD, const Array<OneD, NekDouble> > &inarray,
2035  Array<OneD, NekDouble> &stdV)
2036  {
2037  int nTotQuadPoints = GetTotPoints();
2038  int n_element = m_fields[0]->GetExpSize();
2039  int nBCEdgePts = 0;
2040 
2041  // Getting the velocity vector on the 2D normal space
2042  Array<OneD, Array<OneD, NekDouble> > velocity (m_spacedim);
2043  Array<OneD, Array<OneD, NekDouble> > stdVelocity(m_spacedim);
2044  Array<OneD, NekDouble> pressure (nTotQuadPoints);
2045  Array<OneD, NekDouble> soundspeed (nTotQuadPoints);
2046  LibUtilities::PointsKeyVector ptsKeys;
2047 
2048  // Zero output array
2049  Vmath::Zero(stdV.num_elements(), stdV, 1);
2050 
2051  for (int i = 0; i < m_spacedim; ++i)
2052  {
2053  velocity [i] = Array<OneD, NekDouble>(nTotQuadPoints);
2054  stdVelocity[i] = Array<OneD, NekDouble>(nTotQuadPoints, 0.0);
2055  }
2056 
2057  GetVelocityVector(inarray, velocity);
2058  GetPressure (inarray, velocity, pressure);
2059  GetSoundSpeed (inarray, pressure, soundspeed);
2060 
2061  for(int el = 0; el < n_element; ++el)
2062  {
2063  ptsKeys = m_fields[0]->GetExp(el)->GetPointsKeys();
2064 
2065  // Possible bug: not multiply by jacobian??
2066  const SpatialDomains::GeomFactorsSharedPtr metricInfo =
2067  m_fields[0]->GetExp(el)->GetGeom()->GetMetricInfo();
2068  const Array<TwoD, const NekDouble> &gmat =
2069  m_fields[0]->GetExp(el)->GetGeom()->GetMetricInfo()
2070  ->GetDerivFactors(ptsKeys);
2071 
2072  int nq = m_fields[0]->GetExp(el)->GetTotPoints();
2073 
2074  if(metricInfo->GetGtype() == SpatialDomains::eDeformed)
2075  {
2076  // d xi/ dx = gmat = 1/J * d x/d xi
2077  for (int i = 0; i < m_spacedim; ++i)
2078  {
2079  Vmath::Vmul(nq, gmat[i], 1, velocity[0], 1,
2080  stdVelocity[i], 1);
2081  for (int j = 1; j < m_spacedim; ++j)
2082  {
2083  Vmath::Vvtvp(nq, gmat[m_spacedim*j+i], 1, velocity[j],
2084  1, stdVelocity[i], 1, stdVelocity[i], 1);
2085  }
2086  }
2087  }
2088  else
2089  {
2090  for (int i = 0; i < m_spacedim; ++i)
2091  {
2092  Vmath::Smul(nq, gmat[i][0], velocity[0], 1,
2093  stdVelocity[i], 1);
2094  for (int j = 1; j < m_spacedim; ++j)
2095  {
2096  Vmath::Svtvp(nq, gmat[m_spacedim*j+i][0], velocity[j],
2097  1, stdVelocity[i], 1, stdVelocity[i], 1);
2098  }
2099  }
2100  }
2101 
2102  for (int i = 0; i < nq; ++i)
2103  {
2104  NekDouble pntVelocity = 0.0;
2105  for (int j = 0; j < m_spacedim; ++j)
2106  {
2107  pntVelocity += stdVelocity[j][i]*stdVelocity[j][i];
2108  }
2109  pntVelocity = sqrt(pntVelocity) + soundspeed[nBCEdgePts];
2110  if (pntVelocity > stdV[el])
2111  {
2112  stdV[el] = pntVelocity;
2113  }
2114  nBCEdgePts++;
2115  }
2116  }
2117  }
2118 
2119  /**
2120  * @brief Set the denominator to compute the time step when a cfl
2121  * control is employed. This function is no longer used but is still
2122  * here for being utilised in the future.
2123  *
2124  * @param n Order of expansion element by element.
2125  */
2126  NekDouble CompressibleFlowSystem::GetStabilityLimit(int n)
2127  {
2128  ASSERTL0(n <= 20, "Illegal modes dimension for CFL calculation "
2129  "(P has to be less then 20)");
2130 
2131  NekDouble CFLDG[21] = { 2.0000, 6.0000, 11.8424, 19.1569,
2132  27.8419, 37.8247, 49.0518, 61.4815,
2133  75.0797, 89.8181, 105.6700, 122.6200,
2134  140.6400, 159.7300, 179.8500, 201.0100,
2135  223.1800, 246.3600, 270.5300, 295.6900,
2136  321.8300}; //CFLDG 1D [0-20]
2137  NekDouble CFL = 0.0;
2138 
2139  if (m_projectionType == MultiRegions::eDiscontinuous)
2140  {
2141  CFL = CFLDG[n];
2142  }
2143  else
2144  {
2145  ASSERTL0(false, "Continuous Galerkin stability coefficients "
2146  "not introduced yet.");
2147  }
2148 
2149  return CFL;
2150  }
2151 
2152  /**
2153  * @brief Compute the vector of denominators to compute the time step
2154  * when a cfl control is employed. This function is no longer used but
2155  * is still here for being utilised in the future.
2156  *
2157  * @param ExpOrder Order of expansion element by element.
2158  */
2159  Array<OneD, NekDouble> CompressibleFlowSystem::GetStabilityLimitVector(
2160  const Array<OneD,int> &ExpOrder)
2161  {
2162  int i;
2163  Array<OneD,NekDouble> returnval(m_fields[0]->GetExpSize(), 0.0);
2164  for (i =0; i<m_fields[0]->GetExpSize(); i++)
2165  {
2166  returnval[i] = GetStabilityLimit(ExpOrder[i]);
2167  }
2168  return returnval;
2169  }
2170 
2171  void CompressibleFlowSystem::GetSensor(
2172  const Array<OneD, const Array<OneD, NekDouble> > &physarray,
2173  Array<OneD, NekDouble> &Sensor,
2174  Array<OneD, NekDouble> &SensorKappa)
2175  {
2176  int e, NumModesElement, nQuadPointsElement;
2177  int nTotQuadPoints = GetTotPoints();
2178  int nElements = m_fields[0]->GetExpSize();
2179 
2180  // Find solution (SolP) at p = P;
2181  // The input array (physarray) is the solution at p = P;
2182 
2183  Array<OneD,int> ExpOrderElement = GetNumExpModesPerExp();
2184 
2185  Array<OneD, NekDouble> SolP(nTotQuadPoints,0.0);
2186  Array<OneD, NekDouble> SolPmOne(nTotQuadPoints,0.0);
2187  Array<OneD, NekDouble> SolNorm(nTotQuadPoints,0.0);
2188 
2189  Vmath::Vcopy(nTotQuadPoints,physarray[0],1,SolP,1);
2190 
2191  int CoeffsCount = 0;
2192 
2193  for (e = 0; e < nElements; e++)
2194  {
2195  NumModesElement = ExpOrderElement[e];
2196 
2197  int nQuadPointsElement = m_fields[0]->GetExp(e)->GetTotPoints();
2198  int nCoeffsElement = m_fields[0]->GetExp(e)->GetNcoeffs();
2199  int numCutOff = NumModesElement - 1;
2200 
2201  // Set-up of the Orthogonal basis for a Quadrilateral element which is
2202  // needed to obtain thesolution at P = p - 1;
2203 
2204  Array<OneD, NekDouble> SolPElementPhys(nQuadPointsElement,0.0);
2205  Array<OneD, NekDouble> SolPElementCoeffs(nCoeffsElement,0.0);
2206 
2207  Array<OneD, NekDouble> SolPmOneElementPhys(nQuadPointsElement,0.0);
2208  Array<OneD, NekDouble> SolPmOneElementCoeffs(nCoeffsElement,0.0);
2209 
2210  // create vector the save the solution points per element at P = p;
2211 
2212  for (int i = 0; i < nQuadPointsElement; i++)
2213  {
2214  SolPElementPhys[i] = SolP[CoeffsCount+i];
2215  }
2216 
2217  m_fields[0]->GetExp(e)->FwdTrans(SolPElementPhys,
2218  SolPElementCoeffs);
2219 
2220  // ReduceOrderCoeffs reduces the polynomial order of the solution that
2221  // is represented by the coeffs given as an inarray. This is done by
2222  // projecting the higher order solution onto the orthogonal basis and
2223  // padding the higher order coefficients with zeros.
2224 
2225  m_fields[0]->GetExp(e)->ReduceOrderCoeffs(numCutOff,
2226  SolPElementCoeffs,
2227  SolPmOneElementCoeffs);
2228 
2229  m_fields[0]->GetExp(e)->BwdTrans(SolPmOneElementCoeffs,
2230  SolPmOneElementPhys);
2231 
2232  for (int i = 0; i < nQuadPointsElement; i++)
2233  {
2234  SolPmOne[CoeffsCount+i] = SolPmOneElementPhys[i];
2235  }
2236 
2237  NekDouble SolPmeanNumerator = 0.0;
2238  NekDouble SolPmeanDenumerator = 0.0;
2239 
2240  // Determining the norm of the numerator of the Sensor
2241 
2242  Vmath::Vsub(nQuadPointsElement,
2243  SolPElementPhys, 1,
2244  SolPmOneElementPhys, 1,
2245  SolNorm, 1);
2246 
2247  Vmath::Vmul(nQuadPointsElement,
2248  SolNorm, 1,
2249  SolNorm, 1,
2250  SolNorm, 1);
2251 
2252  for (int i = 0; i < nQuadPointsElement; i++)
2253  {
2254  SolPmeanNumerator += SolNorm[i];
2255  SolPmeanDenumerator += SolPElementPhys[i];
2256  }
2257 
2258  for (int i = 0; i < nQuadPointsElement; ++i)
2259  {
2260  Sensor[CoeffsCount+i] = sqrt(SolPmeanNumerator/nQuadPointsElement)
2261  /sqrt(SolPmeanDenumerator/nQuadPointsElement);
2262 
2263  Sensor[CoeffsCount+i] = log10(Sensor[CoeffsCount+i]);
2264  }
2265  CoeffsCount += nQuadPointsElement;
2266  }
2267 
2268  CoeffsCount = 0.0;
2269 
2270  for (e = 0; e < nElements; e++)
2271  {
2272  NumModesElement = ExpOrderElement[e];
2273  NekDouble ThetaS = m_mu0;
2274  NekDouble Phi0 = m_Skappa;
2275  NekDouble DeltaPhi = m_Kappa;
2276  nQuadPointsElement = m_fields[0]->GetExp(e)->GetTotPoints();
2277 
2278  for (int i = 0; i < nQuadPointsElement; i++)
2279  {
2280  if (Sensor[CoeffsCount+i] <= (Phi0 - DeltaPhi))
2281  {
2282  SensorKappa[CoeffsCount+i] = 0;
2283  }
2284  else if(Sensor[CoeffsCount+i] >= (Phi0 + DeltaPhi))
2285  {
2286  SensorKappa[CoeffsCount+i] = ThetaS;
2287  }
2288  else if(abs(Sensor[CoeffsCount+i]-Phi0) < DeltaPhi)
2289  {
2290  SensorKappa[CoeffsCount+i] = ThetaS/2*(1+sin(M_PI*
2291  (Sensor[CoeffsCount+i]-Phi0)
2292  /(2*DeltaPhi)));
2293  }
2294  }
2295 
2296  CoeffsCount += nQuadPointsElement;
2297  }
2298 
2299  }
2300 
2301  void CompressibleFlowSystem::GetForcingTerm(
2302  const Array<OneD, const Array<OneD, NekDouble> > &inarray,
2303  Array<OneD, Array<OneD, NekDouble> > outarrayForcing)
2304  {
2305  const int nPts = m_fields[0]->GetTotPoints();
2306  const int nvariables = m_fields.num_elements();
2307  const int nElements = m_fields[0]->GetExpSize();
2308 
2309  Array<OneD, NekDouble> Sensor(nPts, 0.0);
2310  Array<OneD, NekDouble> SensorKappa(nPts, 0.0);
2311  Array <OneD, NekDouble > Lambda(nPts, 0.0);
2312  Array <OneD, NekDouble > Tau(nPts, 1.0);
2313  Array <OneD, NekDouble > soundspeed(nPts, 0.0);
2314  Array <OneD, NekDouble > pressure(nPts, 0.0);
2315  Array <OneD, NekDouble > temperature(nPts, 0.0);
2316  Array <OneD, NekDouble > absVelocity(nPts, 0.0);
2317  Array <OneD, NekDouble > hel(nPts, 0.0);
2318  Array <OneD, NekDouble > h_minmin(m_spacedim, 0.0);
2319 
2320  Array<OneD,int> pOrderElmt = GetNumExpModesPerExp();
2321  Array<OneD, NekDouble> pOrder (nPts, 0.0);
2322  // Thermodynamic related quantities
2323  GetPressure(inarray, pressure);
2324  GetTemperature(inarray, pressure, temperature);
2325  GetSoundSpeed(inarray, pressure, soundspeed);
2326  GetAbsoluteVelocity(inarray, absVelocity);
2327  GetSensor(inarray, Sensor, SensorKappa);
2328 
2329  // Determine the maximum wavespeed
2330  Vmath::Vadd(nPts, absVelocity, 1, soundspeed, 1, Lambda, 1);
2331 
2332  NekDouble LambdaMax = Vmath::Vmax(nPts, Lambda, 1);
2333 
2334  int PointCount = 0;
2335 
2336  for (int e = 0; e < nElements; e++)
2337  {
2338  int nQuadPointsElement = m_fields[0]->GetExp(e)->GetTotPoints();
2339 
2340  for (int n = 0; n < nQuadPointsElement; n++)
2341  {
2342  pOrder[n + PointCount] = pOrderElmt[e];
2343 
2344  Tau[n + PointCount] = 1.0/(m_C1*pOrder[n + PointCount]*LambdaMax); // order 1.0e-06
2345 
2346  outarrayForcing[nvariables-1][n + PointCount] = 1/Tau[n + PointCount]*(m_hFactor*LambdaMax/pOrder[n + PointCount]*SensorKappa[n + PointCount]-inarray[nvariables-1][n + PointCount]);
2347  }
2348  PointCount += nQuadPointsElement;
2349  }
2350  }
2351 
2352  void CompressibleFlowSystem::GetElementDimensions(
2353  Array<OneD, Array<OneD, NekDouble> > &outarray,
2354  Array<OneD, NekDouble > &hmin)
2355  {
2356  // So far, this function is only implemented for quads
2357  const int nElements = m_fields[0]->GetExpSize();
2358 
2359  SpatialDomains::HexGeomSharedPtr ElHexGeom;
2360 
2361  NekDouble hx = 0.0;
2362  NekDouble hy = 0.0;
2363 
2364  for (int e = 0; e < nElements; e++)
2365  {
2366  NekDouble nedges = m_fields[0]->GetExp(e)->GetNedges();
2367  Array <OneD, NekDouble> L1(nedges, 0.0);
2368 
2369  for (int j = 0; j < nedges; ++j)
2370  {
2371 
2372  NekDouble x0 = 0.0;
2373  NekDouble y0 = 0.0;
2374  NekDouble z0 = 0.0;
2375 
2376  NekDouble x1 = 0.0;
2377  NekDouble y1 = 0.0;
2378  NekDouble z1 = 0.0;
2379 
2380  if (boost::dynamic_pointer_cast<SpatialDomains::QuadGeom>(m_fields[0]->GetExp(e)->GetGeom()))
2381  {
2382  SpatialDomains::QuadGeomSharedPtr ElQuadGeom = boost::dynamic_pointer_cast<SpatialDomains::QuadGeom>(m_fields[0]->GetExp(e)->GetGeom());
2383 
2384  ElQuadGeom->GetEdge(j)->GetVertex(0)->GetCoords(x0,y0,z0);
2385  ElQuadGeom->GetEdge(j)->GetVertex(1)->GetCoords(x1,y1,z1);
2386 
2387  L1[j] = sqrt(pow((x0-x1),2)+pow((y0-y1),2)+pow((z0-z1),2));
2388  }
2389  else
2390  {
2391  ASSERTL0(false, "GetElementDimensions() is only implemented for quadrilateral elements");
2392  }
2393  }
2394  // determine the minimum length in x and y direction
2395  // still have to find a better estimate when dealing with unstructured meshes
2396  if(boost::dynamic_pointer_cast<SpatialDomains::QuadGeom>(m_fields[0]->GetExp(e)->GetGeom()))
2397  {
2398  hx = min(L1[0], L1[2]);
2399  hy = min(L1[1], L1[3]);
2400 
2401  outarray[0][e] = hx;
2402  outarray[1][e] = hy;
2403  }
2404  }
2405 
2406  hmin[0] = Vmath::Vmin(outarray[0].num_elements(), outarray[0], 1);
2407  hmin[1] = Vmath::Vmin(outarray[1].num_elements(), outarray[1], 1);
2408  }
2409 
2410  void CompressibleFlowSystem::GetAbsoluteVelocity(
2411  const Array<OneD, const Array<OneD, NekDouble> > &inarray,
2412  Array<OneD, NekDouble> &Vtot)
2413  {
2414  int nTotQuadPoints = GetTotPoints();
2415 
2416  // Getting the velocity vector on the 2D normal space
2417  Array<OneD, Array<OneD, NekDouble> > velocity (m_spacedim);
2418 
2419  Vmath::Zero(Vtot.num_elements(), Vtot, 1);
2420 
2421  for (int i = 0; i < m_spacedim; ++i)
2422  {
2423  velocity[i] = Array<OneD, NekDouble>(nTotQuadPoints);
2424  }
2425 
2426  GetVelocityVector(inarray, velocity);
2427 
2428  for (int i = 0; i < m_spacedim; ++i)
2429  {
2430  Vmath::Vvtvp(nTotQuadPoints,
2431  velocity[i], 1,
2432  velocity[i], 1,
2433  Vtot, 1,
2434  Vtot, 1);
2435  }
2436 
2437  Vmath::Vsqrt(Vtot.num_elements(),Vtot,1,Vtot,1);
2438  }
2439 
2440  void CompressibleFlowSystem::GetSmoothArtificialViscosity(
2441  const Array<OneD, Array<OneD, NekDouble> > &physfield,
2442  Array<OneD, NekDouble > &eps_bar)
2443  {
2444  int nvariables = physfield.num_elements();
2445  int nPts = m_fields[0]->GetTotPoints();
2446 
2447  Array<OneD, NekDouble > pressure (nPts, 0.0);
2448  Array<OneD, NekDouble > temperature (nPts, 0.0);
2449  Array <OneD, NekDouble > sensor (nPts, 0.0);
2450  Array <OneD, NekDouble > SensorKappa (nPts, 0.0);
2451  Array <OneD, NekDouble > absVelocity (nPts, 0.0);
2452  Array <OneD, NekDouble > soundspeed (nPts, 0.0);
2453  Array <OneD, NekDouble > Lambda (nPts, 0.0);
2454  Array <OneD, NekDouble > mu_var (nPts, 0.0);
2455  Array <OneD, NekDouble > h_minmin (m_spacedim, 0.0);
2456  Vmath::Zero(nPts, eps_bar, 1);
2457 
2458  // Thermodynamic related quantities
2459  GetPressure(physfield, pressure);
2460  GetTemperature(physfield, pressure, temperature);
2461  GetSoundSpeed(physfield, pressure, soundspeed);
2462  GetAbsoluteVelocity(physfield, absVelocity);
2463  GetSensor(physfield, sensor, SensorKappa);
2464 
2465  // Determine the maximum wavespeed
2466  Vmath::Vadd(nPts, absVelocity, 1, soundspeed, 1, Lambda, 1);
2467 
2468  // Determine hbar = hx_i/h
2469  Array<OneD,int> pOrderElmt = GetNumExpModesPerExp();
2470 
2471  NekDouble ThetaH = m_FacH;
2472  NekDouble ThetaL = m_FacL;
2473 
2474  NekDouble Phi0 = (ThetaH+ThetaL)/2;
2475  NekDouble DeltaPhi = ThetaH-Phi0;
2476 
2477  Vmath::Zero(eps_bar.num_elements(), eps_bar, 1);
2478 
2479  /*Vmath::Smul(eps_bar.num_elements(),
2480  m_eps_max,
2481  &physfield[nvariables-1][0], 1,
2482  &eps_bar[0], 1);*/
2483 
2484  for (int e = 0; e < eps_bar.num_elements(); e++)
2485  {
2486  if (physfield[nvariables-1][e] <= (Phi0 - DeltaPhi))
2487  {
2488  eps_bar[e] = 0;
2489  }
2490  else if(physfield[nvariables-1][e] >= (Phi0 + DeltaPhi))
2491  {
2492  eps_bar[e] = m_mu0;
2493  }
2494  else if(abs(physfield[nvariables-1][e]-Phi0) < DeltaPhi)
2495  {
2496  eps_bar[e] = m_mu0/2*(1+sin(M_PI*
2497  (physfield[nvariables-1][e]-Phi0)/(2*DeltaPhi)));
2498  }
2499  }
2500 
2501  }
2502 
2503  void CompressibleFlowSystem::GetArtificialDynamicViscosity(
2504  const Array<OneD, Array<OneD, NekDouble> > &physfield,
2505  Array<OneD, NekDouble > &mu_var)
2506  {
2507  const int nElements = m_fields[0]->GetExpSize();
2508 
2509  int PointCount = 0;
2510  int nTotQuadPoints = GetTotPoints();
2511 
2512  Array<OneD, NekDouble> S_e (nTotQuadPoints, 0.0);
2513  Array<OneD, NekDouble> se (nTotQuadPoints, 0.0);
2514  Array<OneD, NekDouble> Sensor (nTotQuadPoints, 0.0);
2515  Array<OneD, NekDouble> SensorKappa(nTotQuadPoints, 0.0);
2516  Array<OneD, NekDouble> absVelocity(nTotQuadPoints, 0.0);
2517  Array<OneD, NekDouble> soundspeed (nTotQuadPoints, 0.0);
2518  Array<OneD, NekDouble> pressure (nTotQuadPoints, 0.0);
2519 
2520  GetAbsoluteVelocity(physfield, absVelocity);
2521  GetPressure (physfield, pressure);
2522  GetSoundSpeed (physfield, pressure, soundspeed);
2523  GetSensor (physfield, Sensor, SensorKappa);
2524 
2525  Array<OneD, int> pOrderElmt = GetNumExpModesPerExp();
2526  Array<OneD, NekDouble> Lambda(nTotQuadPoints, 1.0);
2527  Vmath::Vadd(nTotQuadPoints, absVelocity, 1, soundspeed, 1, Lambda, 1);
2528 
2529  for (int e = 0; e < nElements; e++)
2530  {
2531  // Threshold value specified in C. Biottos thesis. Based on a 1D
2532  // shock tube problem S_k = log10(1/p^4). See G.E. Barter and
2533  // D.L. Darmofal. Shock Capturing with PDE-based artificial
2534  // diffusion for DGFEM: Part 1 Formulation, Journal of Computational
2535  // Physics 229 (2010) 1810-1827 for further reference
2536 
2537  // Adjustable depending on the coarsness of the mesh. Might want to
2538  // move this variable into the session file
2539 
2540  int nQuadPointsElement = m_fields[0]->GetExp(e)->GetTotPoints();
2541  Array <OneD, NekDouble> one2D(nQuadPointsElement, 1.0);
2542 
2543  for (int n = 0; n < nQuadPointsElement; n++)
2544  {
2545  NekDouble mu_0 = m_mu0;
2546 
2547  if (Sensor[n+PointCount] < (m_Skappa-m_Kappa))
2548  {
2549  mu_var[n+PointCount] = 0;
2550  }
2551  else if(Sensor[n+PointCount] >= (m_Skappa-m_Kappa)
2552  && Sensor[n+PointCount] <= (m_Skappa+m_Kappa))
2553  {
2554  mu_var[n+PointCount] = mu_0*(0.5*(1+sin(
2555  M_PI*(Sensor[n+PointCount]-m_Skappa-m_Kappa)/(2*m_Kappa))));
2556  }
2557  else if(Sensor[n+PointCount] > (m_Skappa+m_Kappa))
2558  {
2559  mu_var[n+PointCount] = mu_0;
2560  }
2561  }
2562 
2563  PointCount += nQuadPointsElement;
2564  }
2565  }
2566 
2567  void CompressibleFlowSystem::SetVarPOrderElmt(
2568  const Array<OneD, const Array<OneD, NekDouble> > &physfield,
2569  Array<OneD, NekDouble > &PolyOrder)
2570  {
2571  int e;
2572  NekDouble s_ds, s_sm, s_fl;
2573 
2574  int nElements = m_fields[0]->GetExpSize();
2575  int npts = m_fields[0]->GetTotPoints();
2576 
2577  Array<OneD, NekDouble > Sensor (npts, 0.0);
2578  Array<OneD, NekDouble > SensorKappa (npts, 0.0);
2579  Array<OneD, NekDouble > se (npts,0.0);
2580 
2581  GetSensor(physfield, Sensor, SensorKappa);
2582 
2583  int nQuadPointsElement = 0;
2584  int npCount = 0;
2585  int MinOrder = 2;
2586  int MaxOrder = 12;
2587  int MinOrderShock = 4;
2588 
2589  std::ofstream m_file( "VariablePComposites.txt", std::ios_base::app);
2590  for (int e = 0; e < nElements; e++)
2591  {
2592  m_file << "<C ID=\"" << e+1 << "\"> Q[" << e << "] </C>"<< endl;
2593  }
2594  m_file.close();
2595 
2596  std::ofstream m_file2( "VariablePExpansions.txt", std::ios_base::app);
2597 
2598  for (e = 0; e < nElements; e++)
2599  {
2600  nQuadPointsElement = m_fields[0]->GetExp(e)->GetTotPoints();
2601 
2602  // Define thresholds
2603  // Ideally, these threshold values could be given in the Session File
2604 
2605  s_ds = -5.0;
2606  //s_ds = s_0*log10(PolyOrder[e]);
2607  s_sm = -6;
2608  s_fl = -7;
2609 
2610 
2611  for (int i = 0; i < nQuadPointsElement; i++)
2612  {
2613  se[npCount + i] = (Sensor[npCount + i]);
2614 
2615  if (se[npCount + i] > s_ds)
2616  {
2617  if (PolyOrder[npCount + i] > MinOrderShock)
2618  {
2619  PolyOrder[npCount + i] = PolyOrder[npCount + i] - 1;
2620  }
2621  else if(PolyOrder[e] < MinOrderShock)
2622  {
2623  PolyOrder[npCount + i] = PolyOrder[npCount + i] + 1;
2624  }
2625 
2626  }
2627  else if(se[npCount + i] > s_sm && se[npCount + i] < s_ds)
2628  {
2629  if (PolyOrder[npCount + i] < MaxOrder)
2630  {
2631  PolyOrder[npCount + i] = PolyOrder[npCount + i] + 2;
2632  }
2633  }
2634  else if(se[npCount + i] > s_fl && se[npCount + i] < s_sm)
2635  {
2636  PolyOrder[npCount + i] = PolyOrder[npCount + i] + 1;
2637  }
2638  else if(se[npCount + i] < s_fl)
2639  {
2640  if (PolyOrder[npCount + i] > MinOrder)
2641  {
2642  PolyOrder[npCount + i] = PolyOrder[npCount + i];
2643  }
2644  }
2645  }
2646  m_file2 << "<E COMPOSITE= \"C[" << e+1
2647  << "]\" NUMMODES=\"" << PolyOrder[npCount + 1]
2648  << "\" TYPE=\"MODIFIED\" FIELDS=\"rho,rhou,rhov,rhow,E\" />"
2649  << endl;
2650  npCount += nQuadPointsElement;
2651  }
2652 
2653  m_file2.close();
2654  }
2655 
2656  void CompressibleFlowSystem::v_ExtraFldOutput(
2657  std::vector<Array<OneD, NekDouble> > &fieldcoeffs,
2658  std::vector<std::string> &variables)
2659  {
2660  const int nPhys = m_fields[0]->GetNpoints();
2661  const int nCoeffs = m_fields[0]->GetNcoeffs();
2662  Array<OneD, Array<OneD, NekDouble> > tmp(m_fields.num_elements());
2663 
2664  for (int i = 0; i < m_fields.num_elements(); ++i)
2665  {
2666  tmp[i] = m_fields[i]->GetPhys();
2667  }
2668 
2669  Array<OneD, NekDouble> pressure(nPhys), soundspeed(nPhys), mach(nPhys), sensor(nPhys), SensorKappa(nPhys), smooth(nPhys);
2670 
2671  GetPressure (tmp, pressure);
2672  GetSoundSpeed(tmp, pressure, soundspeed);
2673  GetMach (tmp, soundspeed, mach);
2674  GetSensor (tmp, sensor, SensorKappa);
2675  GetSmoothArtificialViscosity (tmp, smooth);
2676 
2677  Array<OneD, NekDouble> pFwd(nCoeffs), sFwd(nCoeffs), mFwd(nCoeffs), sensFwd(nCoeffs), smoothFwd(nCoeffs);
2678 
2679  m_fields[0]->FwdTrans(pressure, pFwd);
2680  m_fields[0]->FwdTrans(soundspeed, sFwd);
2681  m_fields[0]->FwdTrans(mach, mFwd);
2682  m_fields[0]->FwdTrans(sensor, sensFwd);
2683  m_fields[0]->FwdTrans(smooth, smoothFwd);
2684 
2685  variables.push_back ("p");
2686  variables.push_back ("a");
2687  variables.push_back ("Mach");
2688  variables.push_back ("Sensor");
2689  variables.push_back ("SmoothVisc");
2690  fieldcoeffs.push_back(pFwd);
2691  fieldcoeffs.push_back(sFwd);
2692  fieldcoeffs.push_back(mFwd);
2693  fieldcoeffs.push_back(sensFwd);
2694  fieldcoeffs.push_back(smoothFwd);
2695  }
2696 }
2697 
Generated on Wed Mar 4 2015 16:33:13 for Nektar++ by   doxygen 1.8.1.2