Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
NavierStokesCFE.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File NavierStokesCFE.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: Navier Stokes equations in conservative variables
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
37 
38 using namespace std;
39 
40 namespace Nektar
41 {
42  string NavierStokesCFE::className =
44  "NavierStokesCFE", NavierStokesCFE::create,
45  "NavierStokes equations in conservative variables.");
46 
47  NavierStokesCFE::NavierStokesCFE(
49  : CompressibleFlowSystem(pSession)
50  {
51  }
52 
54  {
55 
56  }
57 
58  /**
59  * @brief Initialization object for CompressibleFlowSystem class.
60  */
62  {
64 
65  string diffName;
66  m_session->LoadSolverInfo("DiffusionType", diffName, "LDGNS");
67 
69  .CreateInstance(diffName, diffName);
70 
72  {
73  m_diffusion->SetFluxVectorNS(
75  this);
76  }
77  else
78  {
79  m_diffusion->SetFluxVectorNS(&NavierStokesCFE::
80  GetViscousFluxVector, this);
81  }
82 
83  // Concluding initialisation of diffusion operator
84  m_diffusion->InitObject (m_session, m_fields);
85  }
86 
88  const Array<OneD, const Array<OneD, NekDouble> > &inarray,
89  Array<OneD, Array<OneD, NekDouble> > &outarray,
90  const Array<OneD, Array<OneD, NekDouble> > &pFwd,
91  const Array<OneD, Array<OneD, NekDouble> > &pBwd)
92  {
93  int i;
94  int nvariables = inarray.num_elements();
95  int npoints = GetNpoints();
96  int nTracePts = GetTraceTotPoints();
97 
98  Array<OneD, Array<OneD, NekDouble> > outarrayDiff(nvariables);
99 
100  Array<OneD, Array<OneD, NekDouble> > inarrayDiff(nvariables-1);
101  Array<OneD, Array<OneD, NekDouble> > inFwd(nvariables-1);
102  Array<OneD, Array<OneD, NekDouble> > inBwd(nvariables-1);
103 
104  for (i = 0; i < nvariables; ++i)
105  {
106  outarrayDiff[i] = Array<OneD, NekDouble>(npoints);
107  }
108 
109  for (i = 0; i < nvariables-1; ++i)
110  {
111  inarrayDiff[i] = Array<OneD, NekDouble>(npoints);
112  inFwd[i] = Array<OneD, NekDouble>(nTracePts);
113  inBwd[i] = Array<OneD, NekDouble>(nTracePts);
114  }
115 
116  // Extract pressure
117  // (use inarrayDiff[0] as a temporary storage for the pressure)
118  m_varConv->GetPressure(inarray, inarrayDiff[0]);
119 
120  // Extract temperature
121  m_varConv->GetTemperature(inarray, inarrayDiff[0],
122  inarrayDiff[nvariables-2]);
123 
124  // Extract velocities
125  m_varConv->GetVelocityVector(inarray, inarrayDiff);
126 
127  // Repeat calculation for trace space
128  if (pFwd == NullNekDoubleArrayofArray ||
130  {
133  }
134  else
135  {
136  m_varConv->GetPressure(pFwd, inFwd[0]);
137  m_varConv->GetPressure(pBwd, inBwd[0]);
138 
139  m_varConv->GetTemperature(pFwd, inFwd[0],
140  inFwd[nvariables-2]);
141  m_varConv->GetTemperature(pBwd, inBwd[0],
142  inBwd[nvariables-2]);
143 
144  m_varConv->GetVelocityVector(pFwd, inFwd);
145  m_varConv->GetVelocityVector(pBwd, inBwd);
146  }
147 
148  // Diffusion term in physical rhs form
149  m_diffusion->Diffuse(nvariables, m_fields, inarrayDiff, outarrayDiff,
150  inFwd, inBwd);
151 
152  for (i = 0; i < nvariables; ++i)
153  {
154  Vmath::Vadd(npoints,
155  outarrayDiff[i], 1,
156  outarray[i], 1,
157  outarray[i], 1);
158  }
159  }
160 
161  /**
162  * @brief Return the flux vector for the LDG diffusion problem.
163  * \todo Complete the viscous flux vector
164  */
166  const Array<OneD, Array<OneD, NekDouble> > &physfield,
167  Array<OneD, Array<OneD, Array<OneD, NekDouble> > > &derivativesO1,
168  Array<OneD, Array<OneD, Array<OneD, NekDouble> > > &viscousTensor)
169  {
170  int i, j;
171  int nVariables = m_fields.num_elements();
172  int nPts = physfield[0].num_elements();
173 
174  // Stokes hypothesis
175  const NekDouble lambda = -2.0/3.0;
176 
177  // Auxiliary variables
178  Array<OneD, NekDouble > mu (nPts, 0.0);
179  Array<OneD, NekDouble > thermalConductivity(nPts, 0.0);
180  Array<OneD, NekDouble > divVel (nPts, 0.0);
181 
182  // Variable viscosity through the Sutherland's law
183  if (m_ViscosityType == "Variable")
184  {
185  m_varConv->GetDynamicViscosity(physfield[nVariables-2], mu);
186  NekDouble tRa = m_Cp / m_Prandtl;
187  Vmath::Smul(nPts, tRa, mu, 1, thermalConductivity, 1);
188  }
189  else
190  {
191  Vmath::Fill(nPts, m_mu, mu, 1);
193  thermalConductivity, 1);
194  }
195 
196  // Velocity divergence
197  for (j = 0; j < m_spacedim; ++j)
198  {
199  Vmath::Vadd(nPts, divVel, 1, derivativesO1[j][j], 1,
200  divVel, 1);
201  }
202 
203  // Velocity divergence scaled by lambda * mu
204  Vmath::Smul(nPts, lambda, divVel, 1, divVel, 1);
205  Vmath::Vmul(nPts, mu, 1, divVel, 1, divVel, 1);
206 
207  // Viscous flux vector for the rho equation = 0
208  for (i = 0; i < m_spacedim; ++i)
209  {
210  Vmath::Zero(nPts, viscousTensor[i][0], 1);
211  }
212 
213  // Viscous stress tensor (for the momentum equations)
214  for (i = 0; i < m_spacedim; ++i)
215  {
216  for (j = i; j < m_spacedim; ++j)
217  {
218  Vmath::Vadd(nPts, derivativesO1[i][j], 1,
219  derivativesO1[j][i], 1,
220  viscousTensor[i][j+1], 1);
221 
222  Vmath::Vmul(nPts, mu, 1,
223  viscousTensor[i][j+1], 1,
224  viscousTensor[i][j+1], 1);
225 
226  if (i == j)
227  {
228  // Add divergence term to diagonal
229  Vmath::Vadd(nPts, viscousTensor[i][j+1], 1,
230  divVel, 1,
231  viscousTensor[i][j+1], 1);
232  }
233  else
234  {
235  // Copy to make symmetric
236  Vmath::Vcopy(nPts, viscousTensor[i][j+1], 1,
237  viscousTensor[j][i+1], 1);
238  }
239  }
240  }
241 
242  // Terms for the energy equation
243  for (i = 0; i < m_spacedim; ++i)
244  {
245  Vmath::Zero(nPts, viscousTensor[i][m_spacedim+1], 1);
246  // u_j * tau_ij
247  for (j = 0; j < m_spacedim; ++j)
248  {
249  Vmath::Vvtvp(nPts, physfield[j], 1,
250  viscousTensor[i][j+1], 1,
251  viscousTensor[i][m_spacedim+1], 1,
252  viscousTensor[i][m_spacedim+1], 1);
253  }
254  // Add k*T_i
255  Vmath::Vvtvp(nPts, thermalConductivity, 1,
256  derivativesO1[i][m_spacedim], 1,
257  viscousTensor[i][m_spacedim+1], 1,
258  viscousTensor[i][m_spacedim+1], 1);
259  }
260  }
261 
262  /**
263  * @brief Return the flux vector for the LDG diffusion problem.
264  * \todo Complete the viscous flux vector
265  */
267  const Array<OneD, Array<OneD, NekDouble> > &physfield,
268  Array<OneD, Array<OneD, Array<OneD, NekDouble> > > &derivativesO1,
269  Array<OneD, Array<OneD, Array<OneD, NekDouble> > > &viscousTensor)
270  {
271  int i, j;
272  int nVariables = m_fields.num_elements();
273  // Factor to rescale 1d points in dealiasing.
274  NekDouble OneDptscale = 2;
275  // Get number of points to dealias a cubic non-linearity
276  int nPts = m_fields[0]->Get1DScaledTotPoints(OneDptscale);
277  int nPts_orig = physfield[0].num_elements();
278 
279  // Stokes hypothesis
280  const NekDouble lambda = -2.0/3.0;
281 
282  // Auxiliary variables
283  Array<OneD, NekDouble > mu (nPts, 0.0);
284  Array<OneD, NekDouble > thermalConductivity(nPts, 0.0);
285  Array<OneD, NekDouble > divVel (nPts, 0.0);
286 
287  // Variable viscosity through the Sutherland's law
288  if (m_ViscosityType == "Variable")
289  {
290  m_varConv->GetDynamicViscosity(physfield[nVariables-2], mu);
291  NekDouble tRa = m_Cp / m_Prandtl;
292  Vmath::Smul(nPts, tRa, mu, 1, thermalConductivity, 1);
293  }
294  else
295  {
296  Vmath::Fill(nPts, m_mu, mu, 1);
298  thermalConductivity, 1);
299  }
300 
301  // Interpolate inputs and initialise interpolated output
304  deriv_interp(m_spacedim);
306  out_interp(m_spacedim);
307  for (i = 0; i < m_spacedim; ++i)
308  {
309  // Interpolate velocity
310  vel_interp[i] = Array<OneD, NekDouble> (nPts);
311  m_fields[0]->PhysInterp1DScaled(
312  OneDptscale, physfield[i], vel_interp[i]);
313 
314  // Interpolate derivatives
315  deriv_interp[i] = Array<OneD,Array<OneD,NekDouble> > (m_spacedim+1);
316  for (j = 0; j < m_spacedim+1; ++j)
317  {
318  deriv_interp[i][j] = Array<OneD, NekDouble> (nPts);
319  m_fields[0]->PhysInterp1DScaled(
320  OneDptscale, derivativesO1[i][j], deriv_interp[i][j]);
321  }
322 
323  // Output (start from j=1 since flux is zero for rho)
324  out_interp[i] = Array<OneD,Array<OneD,NekDouble> > (m_spacedim+2);
325  for (j = 1; j < m_spacedim+2; ++j)
326  {
327  out_interp[i][j] = Array<OneD, NekDouble> (nPts);
328  }
329  }
330 
331  // Velocity divergence
332  for (j = 0; j < m_spacedim; ++j)
333  {
334  Vmath::Vadd(nPts, divVel, 1, deriv_interp[j][j], 1,
335  divVel, 1);
336  }
337 
338  // Velocity divergence scaled by lambda * mu
339  Vmath::Smul(nPts, lambda, divVel, 1, divVel, 1);
340  Vmath::Vmul(nPts, mu, 1, divVel, 1, divVel, 1);
341 
342  // Viscous flux vector for the rho equation = 0 (no need to dealias)
343  for (i = 0; i < m_spacedim; ++i)
344  {
345  Vmath::Zero(nPts_orig, viscousTensor[i][0], 1);
346  }
347 
348  // Viscous stress tensor (for the momentum equations)
349  for (i = 0; i < m_spacedim; ++i)
350  {
351  for (j = i; j < m_spacedim; ++j)
352  {
353  Vmath::Vadd(nPts, deriv_interp[i][j], 1,
354  deriv_interp[j][i], 1,
355  out_interp[i][j+1], 1);
356 
357  Vmath::Vmul(nPts, mu, 1,
358  out_interp[i][j+1], 1,
359  out_interp[i][j+1], 1);
360 
361  if (i == j)
362  {
363  // Add divergence term to diagonal
364  Vmath::Vadd(nPts, out_interp[i][j+1], 1,
365  divVel, 1,
366  out_interp[i][j+1], 1);
367  }
368  else
369  {
370  // Make symmetric
371  out_interp[j][i+1] = out_interp[i][j+1];
372  }
373  }
374  }
375 
376  // Terms for the energy equation
377  for (i = 0; i < m_spacedim; ++i)
378  {
379  Vmath::Zero(nPts, out_interp[i][m_spacedim+1], 1);
380  // u_j * tau_ij
381  for (j = 0; j < m_spacedim; ++j)
382  {
383  Vmath::Vvtvp(nPts, vel_interp[j], 1,
384  out_interp[i][j+1], 1,
385  out_interp[i][m_spacedim+1], 1,
386  out_interp[i][m_spacedim+1], 1);
387  }
388  // Add k*T_i
389  Vmath::Vvtvp(nPts, thermalConductivity, 1,
390  deriv_interp[i][m_spacedim], 1,
391  out_interp[i][m_spacedim+1], 1,
392  out_interp[i][m_spacedim+1], 1);
393  }
394 
395  // Project to original space
396  for (i = 0; i < m_spacedim; ++i)
397  {
398  for (j = 1; j < m_spacedim+2; ++j)
399  {
400  m_fields[0]->PhysGalerkinProjection1DScaled(
401  OneDptscale,
402  out_interp[i][j],
403  viscousTensor[i][j]);
404  }
405  }
406  }
407 }
virtual void v_InitObject()
Initialization object for CompressibleFlowSystem class.
virtual void v_InitObject()
Initialization object for CompressibleFlowSystem class.
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
VariableConverterSharedPtr m_varConv
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition: Vmath.cpp:46
DiffusionFactory & GetDiffusionFactory()
Definition: Diffusion.cpp:42
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:442
STL namespace.
bool m_specHP_dealiasing
Flag to determine if dealisising is usde for the Spectral/hp element discretisation.
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.
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:51
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:213
int m_spacedim
Spatial dimension (>= expansion dim).
double NekDouble
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.
EquationSystemFactory & GetEquationSystemFactory()
SolverUtils::DiffusionSharedPtr m_diffusion
SOLVER_UTILS_EXPORT int GetTraceTotPoints()
SOLVER_UTILS_EXPORT int GetNpoints()
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
virtual void v_DoDiffusion(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const Array< OneD, Array< OneD, NekDouble > > &pFwd, const Array< OneD, Array< OneD, NekDouble > > &pBwd)
static Array< OneD, Array< OneD, NekDouble > > NullNekDoubleArrayofArray
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:373
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1061
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:299
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:183
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215