Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
DiffusionLDGNS.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File: DiffusionLDGNS.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: LDGNS diffusion class.
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
37 #include <iostream>
38 #include <iomanip>
39 
41 
42 namespace Nektar
43 {
44  namespace SolverUtils
45  {
47  RegisterCreatorFunction("LDGNS", DiffusionLDGNS::create);
48 
50  {
51  }
52 
56  {
57  m_session = pSession;
58  m_session->LoadParameter ("Gamma", m_gamma, 1.4);
59  m_session->LoadParameter ("GasConstant", m_gasConstant, 287.058);
60  m_session->LoadParameter ("Twall", m_Twall, 300.15);
61  m_session->LoadSolverInfo("ViscosityType", m_ViscosityType,
62  "Constant");
63  m_session->LoadParameter ("mu", m_mu, 1.78e-05);
64  m_session->LoadParameter ("thermalConductivity",
65  m_thermalConductivity, 0.0257);
66  m_session->LoadParameter ("rhoInf", m_rhoInf, 1.225);
67  m_session->LoadParameter ("pInf", m_pInf, 101325);
68 
69  // Setting up the normals
70  int i;
71  int nDim = pFields[0]->GetCoordim(0);
72  int nTracePts = pFields[0]->GetTrace()->GetTotPoints();
73 
74  m_spaceDim = nDim;
75  if (pSession->DefinesSolverInfo("HOMOGENEOUS"))
76  {
77  m_spaceDim = 3;
78  }
79 
80  m_diffDim = m_spaceDim - nDim;
81 
84  for(i = 0; i < m_spaceDim; ++i)
85  {
86  m_traceVel[i] = Array<OneD, NekDouble> (nTracePts, 0.0);
87  m_traceNormals[i] = Array<OneD, NekDouble> (nTracePts);
88  }
89  pFields[0]->GetTrace()->GetNormals(m_traceNormals);
90  }
91 
92  /**
93  * @brief Calculate weak DG Diffusion in the LDG form for the
94  * Navier-Stokes (NS) equations:
95  *
96  * \f$ \langle\psi, \hat{u}\cdot n\rangle
97  * - \langle\nabla\psi \cdot u\rangle
98  * \langle\phi, \hat{q}\cdot n\rangle -
99  * (\nabla \phi \cdot q) \rangle \f$
100  *
101  * The equations that need a diffusion operator are those related
102  * with the velocities and with the energy.
103  *
104  */
106  const int nConvectiveFields,
108  const Array<OneD, Array<OneD, NekDouble> > &inarray,
109  Array<OneD, Array<OneD, NekDouble> > &outarray,
110  const Array<OneD, Array<OneD, NekDouble> > &pFwd,
111  const Array<OneD, Array<OneD, NekDouble> > &pBwd)
112  {
113  int i, j;
114  int nDim = fields[0]->GetCoordim(0);
115  int nScalars = inarray.num_elements();
116  int nPts = fields[0]->GetTotPoints();
117  int nCoeffs = fields[0]->GetNcoeffs();
118  int nTracePts = fields[0]->GetTrace()->GetTotPoints();
119 
120  Array<OneD, NekDouble> tmp1(nCoeffs);
121  Array<OneD, Array<OneD, NekDouble> > tmp2(nConvectiveFields);
122 
124  numericalFluxO1(m_spaceDim);
126  derivativesO1(m_spaceDim);
127 
128  for (j = 0; j < m_spaceDim; ++j)
129  {
130  numericalFluxO1[j] = Array<OneD, Array<OneD, NekDouble> >(
131  nScalars);
132  derivativesO1[j] = Array<OneD, Array<OneD, NekDouble> >(
133  nScalars);
134 
135  for (i = 0; i < nScalars; ++i)
136  {
137  numericalFluxO1[j][i] = Array<OneD, NekDouble>(
138  nTracePts, 0.0);
139  derivativesO1[j][i] = Array<OneD, NekDouble>(nPts, 0.0);
140  }
141  }
142 
143  // Compute the numerical fluxes for the first order derivatives
144  v_NumericalFluxO1(fields, inarray, numericalFluxO1, pFwd, pBwd);
145 
146  for (j = 0; j < nDim; ++j)
147  {
148  for (i = 0; i < nScalars; ++i)
149  {
150  fields[i]->IProductWRTDerivBase (j, inarray[i], tmp1);
151  Vmath::Neg (nCoeffs, tmp1, 1);
152  fields[i]->AddTraceIntegral (numericalFluxO1[j][i],
153  tmp1);
154  fields[i]->SetPhysState (false);
155  fields[i]->MultiplyByElmtInvMass(tmp1, tmp1);
156  fields[i]->BwdTrans (tmp1, derivativesO1[j][i]);
157  }
158  }
159 
160  // For 3D Homogeneous 1D only take derivatives in 3rd direction
161  if (m_diffDim == 1)
162  {
163  for (i = 0; i < nScalars; ++i)
164  {
165  derivativesO1[2][i] = m_homoDerivs[i];
166  }
167  }
168 
169  // Initialisation viscous tensor
171  (m_spaceDim);
172  Array<OneD, Array<OneD, NekDouble> > viscousFlux(nConvectiveFields);
173 
174  for (j = 0; j < m_spaceDim; ++j)
175  {
177  nScalars+1);
178  for (i = 0; i < nScalars+1; ++i)
179  {
180  m_viscTensor[j][i] = Array<OneD, NekDouble>(nPts, 0.0);
181  }
182  }
183 
184  for (i = 0; i < nConvectiveFields; ++i)
185  {
186  viscousFlux[i] = Array<OneD, NekDouble>(nTracePts, 0.0);
187  }
188 
189  m_fluxVectorNS(inarray, derivativesO1, m_viscTensor);
190 
191  // Compute u from q_{\eta} and q_{\xi}
192  // Obtain numerical fluxes
193  v_NumericalFluxO2(fields, inarray, m_viscTensor, viscousFlux);
194 
195  for (i = 0; i < nConvectiveFields; ++i)
196  {
197  tmp2[i] = Array<OneD, NekDouble>(nCoeffs, 0.0);
198 
199  for (j = 0; j < nDim; ++j)
200  {
201  fields[i]->IProductWRTDerivBase(j, m_viscTensor[j][i], tmp1);
202  Vmath::Vadd(nCoeffs, tmp1, 1, tmp2[i], 1, tmp2[i], 1);
203  }
204 
205  // Evaulate <\phi, \hat{F}\cdot n> - outarray[i]
206  Vmath::Neg (nCoeffs, tmp2[i], 1);
207  fields[i]->AddTraceIntegral (viscousFlux[i], tmp2[i]);
208  fields[i]->SetPhysState (false);
209  fields[i]->MultiplyByElmtInvMass(tmp2[i], tmp2[i]);
210  fields[i]->BwdTrans (tmp2[i], outarray[i]);
211  }
212  }
213 
214  /**
215  * @brief Builds the numerical flux for the 1st order derivatives
216  *
217  */
220  const Array<OneD, Array<OneD, NekDouble> > &inarray,
222  &numericalFluxO1,
223  const Array<OneD, Array<OneD, NekDouble> > &pFwd,
224  const Array<OneD, Array<OneD, NekDouble> > &pBwd)
225  {
226  int i, j;
227  int nTracePts = fields[0]->GetTrace()->GetTotPoints();
228  int nScalars = inarray.num_elements();
229  int nDim = fields[0]->GetCoordim(0);
230 
231  Array<OneD, NekDouble > Vn (nTracePts, 0.0);
232 
233  // Get the normal velocity Vn
234  for(i = 0; i < nDim; ++i)
235  {
236  Vmath::Svtvp(nTracePts, 1.0, m_traceNormals[i], 1,
237  Vn, 1, Vn, 1);
238  }
239 
240  // Store forwards/backwards space along trace space
243  Array<OneD, Array<OneD, NekDouble> > numflux(nScalars);
244 
245  for (i = 0; i < nScalars; ++i)
246  {
247  if (pFwd == NullNekDoubleArrayofArray ||
249  {
250  Fwd = Array<OneD, NekDouble>(nTracePts);
251  Bwd = Array<OneD, NekDouble>(nTracePts);
252  fields[i]->GetFwdBwdTracePhys(inarray[i], Fwd, Bwd);
253  }
254  else
255  {
256  Fwd = pFwd[i];
257  Bwd = pBwd[i];
258  }
259  numflux[i] = Array<OneD, NekDouble>(nTracePts);
260  fields[0]->GetTrace()->Upwind(Vn, Fwd, Bwd, numflux[i]);
261  }
262 
263  // Extract internal values of the scalar variables for Neumann bcs
264  Array< OneD, Array<OneD, NekDouble > > uplus(nScalars);
265 
266  for (i = 0; i < nScalars; ++i)
267  {
268  uplus[i] = Array<OneD, NekDouble>(nTracePts, 0.0);
269  fields[i]->ExtractTracePhys(inarray[i], uplus[i]);
270  }
271 
272  // Modify the values in case of boundary interfaces
273  if (fields[0]->GetBndCondExpansions().num_elements())
274  {
275  v_WeakPenaltyO1(fields, inarray, uplus, numflux);
276  }
277 
278  // Splitting the numerical flux into the dimensions
279  for (j = 0; j < m_spaceDim; ++j)
280  {
281  for (i = 0; i < nScalars; ++i)
282  {
283  Vmath::Vmul(nTracePts, m_traceNormals[j], 1,
284  numflux[i], 1, numericalFluxO1[j][i], 1);
285  }
286  }
287  }
288 
289  /**
290  * @brief Imposes appropriate bcs for the 1st order derivatives
291  *
292  */
295  const Array<OneD, Array<OneD, NekDouble> > &inarray,
296  const Array<OneD, Array<OneD, NekDouble> > &uplus,
297  Array<OneD, Array<OneD, NekDouble> > &penaltyfluxO1)
298  {
299  int cnt;
300  int i, j, e;
301  int id1, id2;
302 
303  int nBndEdgePts, nBndEdges, nBndRegions;
304 
305  int nTracePts = fields[0]->GetTrace()->GetTotPoints();
306  int nScalars = inarray.num_elements();
307 
308  Array<OneD, NekDouble> tmp1(nTracePts, 0.0);
309  Array<OneD, NekDouble> tmp2(nTracePts, 0.0);
310  Array<OneD, NekDouble> Tw(nTracePts, m_Twall);
311 
312  Array< OneD, Array<OneD, NekDouble > > scalarVariables(nScalars);
313 
314  // Extract internal values of the scalar variables for Neumann bcs
315  for (i = 0; i < nScalars; ++i)
316  {
317  scalarVariables[i] = Array<OneD, NekDouble>(nTracePts, 0.0);
318  }
319 
320  // Compute boundary conditions for velocities
321  for (i = 0; i < nScalars-1; ++i)
322  {
323  // Note that cnt has to loop on nBndRegions and nBndEdges
324  // and has to be reset to zero per each equation
325  cnt = 0;
326  nBndRegions = fields[i+1]->
327  GetBndCondExpansions().num_elements();
328  for (j = 0; j < nBndRegions; ++j)
329  {
330  nBndEdges = fields[i+1]->
331  GetBndCondExpansions()[j]->GetExpSize();
332  for (e = 0; e < nBndEdges; ++e)
333  {
334  nBndEdgePts = fields[i+1]->
335  GetBndCondExpansions()[j]->GetExp(e)->GetTotPoints();
336 
337  id1 = fields[i+1]->
338  GetBndCondExpansions()[j]->GetPhys_Offset(e);
339 
340  id2 = fields[0]->GetTrace()->
341  GetPhys_Offset(fields[0]->GetTraceMap()->
342  GetBndCondTraceToGlobalTraceMap(cnt++));
343 
344  // Reinforcing bcs for velocity in case of Wall bcs
345  if (boost::iequals(fields[i]->GetBndConditions()[j]->
346  GetUserDefined(),"WallViscous") ||
347  boost::iequals(fields[i]->GetBndConditions()[j]->
348  GetUserDefined(),"WallAdiabatic"))
349  {
350  Vmath::Zero(nBndEdgePts,
351  &scalarVariables[i][id2], 1);
352 
353  }
354 
355  // Imposing velocity bcs if not Wall
356  else if (fields[i]->GetBndConditions()[j]->
357  GetBoundaryConditionType() ==
359  {
360  Vmath::Vdiv(nBndEdgePts,
361  &(fields[i+1]->GetBndCondExpansions()[j]->
362  UpdatePhys())[id1], 1,
363  &(fields[0]->GetBndCondExpansions()[j]->
364  UpdatePhys())[id1], 1,
365  &scalarVariables[i][id2], 1);
366  }
367 
368  // For Dirichlet boundary condition: uflux = u_bcs
369  if (fields[i]->GetBndConditions()[j]->
370  GetBoundaryConditionType() ==
372  {
373  Vmath::Vcopy(nBndEdgePts,
374  &scalarVariables[i][id2], 1,
375  &penaltyfluxO1[i][id2], 1);
376  }
377 
378  // For Neumann boundary condition: uflux = u_+
379  else if ((fields[i]->GetBndConditions()[j])->
380  GetBoundaryConditionType() ==
382  {
383  Vmath::Vcopy(nBndEdgePts,
384  &uplus[i][id2], 1,
385  &penaltyfluxO1[i][id2], 1);
386  }
387 
388  // Building kinetic energy to be used for T bcs
389  Vmath::Vmul(nBndEdgePts,
390  &scalarVariables[i][id2], 1,
391  &scalarVariables[i][id2], 1,
392  &tmp1[id2], 1);
393 
394  Vmath::Smul(nBndEdgePts, 0.5,
395  &tmp1[id2], 1,
396  &tmp1[id2], 1);
397 
398  Vmath::Vadd(nBndEdgePts,
399  &tmp2[id2], 1,
400  &tmp1[id2], 1,
401  &tmp2[id2], 1);
402  }
403  }
404  }
405 
406  // Compute boundary conditions for temperature
407  cnt = 0;
408  nBndRegions = fields[nScalars]->
409  GetBndCondExpansions().num_elements();
410  for (j = 0; j < nBndRegions; ++j)
411  {
412  nBndEdges = fields[nScalars]->
413  GetBndCondExpansions()[j]->GetExpSize();
414  for (e = 0; e < nBndEdges; ++e)
415  {
416  nBndEdgePts = fields[nScalars]->
417  GetBndCondExpansions()[j]->GetExp(e)->GetTotPoints();
418 
419  id1 = fields[nScalars]->
420  GetBndCondExpansions()[j]->GetPhys_Offset(e);
421 
422  id2 = fields[0]->GetTrace()->
423  GetPhys_Offset(fields[0]->GetTraceMap()->
424  GetBndCondTraceToGlobalTraceMap(cnt++));
425 
426  // Imposing Temperature Twall at the wall
427  if (boost::iequals(fields[i]->GetBndConditions()[j]->
428  GetUserDefined(),"WallViscous"))
429  {
430  Vmath::Vcopy(nBndEdgePts,
431  &Tw[0], 1,
432  &scalarVariables[nScalars-1][id2], 1);
433  }
434  // Imposing Temperature through condition on the Energy
435  // for no wall boundaries (e.g. farfield)
436  else if (fields[i]->GetBndConditions()[j]->
437  GetBoundaryConditionType() ==
439  {
440  // Divide E by rho
441  Vmath::Vdiv(nBndEdgePts,
442  &(fields[nScalars]->
443  GetBndCondExpansions()[j]->
444  GetPhys())[id1], 1,
445  &(fields[0]->
446  GetBndCondExpansions()[j]->
447  GetPhys())[id1], 1,
448  &scalarVariables[nScalars-1][id2], 1);
449 
450  // Subtract kinetic energy to E/rho
451  Vmath::Vsub(nBndEdgePts,
452  &scalarVariables[nScalars-1][id2], 1,
453  &tmp2[id2], 1,
454  &scalarVariables[nScalars-1][id2], 1);
455 
456  // Multiply by constant factor (gamma-1)/R
457  Vmath::Smul(nBndEdgePts, (m_gamma - 1)/m_gasConstant,
458  &scalarVariables[nScalars-1][id2], 1,
459  &scalarVariables[nScalars-1][id2], 1);
460  }
461 
462  // For Dirichlet boundary condition: uflux = u_bcs
463  if (fields[nScalars]->GetBndConditions()[j]->
464  GetBoundaryConditionType() ==
466  !boost::iequals(
467  fields[nScalars]->GetBndConditions()[j]
468  ->GetUserDefined(), "WallAdiabatic"))
469  {
470  Vmath::Vcopy(nBndEdgePts,
471  &scalarVariables[nScalars-1][id2], 1,
472  &penaltyfluxO1[nScalars-1][id2], 1);
473 
474  }
475 
476  // For Neumann boundary condition: uflux = u_+
477  else if (((fields[nScalars]->GetBndConditions()[j])->
478  GetBoundaryConditionType() ==
480  boost::iequals(fields[nScalars]->GetBndConditions()[j]->
481  GetUserDefined(), "WallAdiabatic"))
482  {
483  Vmath::Vcopy(nBndEdgePts,
484  &uplus[nScalars-1][id2], 1,
485  &penaltyfluxO1[nScalars-1][id2], 1);
486 
487  }
488  }
489  }
490  }
491 
492  /**
493  * @brief Build the numerical flux for the 2nd order derivatives
494  *
495  */
498  const Array<OneD, Array<OneD, NekDouble> > &ufield,
501  {
502  int i, j;
503  int nTracePts = fields[0]->GetTrace()->GetTotPoints();
504  int nVariables = fields.num_elements();
505  int nDim = fields[0]->GetCoordim(0);
506 
507  Array<OneD, NekDouble > Vn (nTracePts, 0.0);
508 
509  Array<OneD, NekDouble > qFwd (nTracePts);
510  Array<OneD, NekDouble > qBwd (nTracePts);
511  Array<OneD, NekDouble > qfluxtemp(nTracePts, 0.0);
512 
513  // Get the normal velocity Vn
514  for(i = 0; i < nDim; ++i)
515  {
516  Vmath::Svtvp(nTracePts, 1.0, m_traceNormals[i], 1,
517  Vn, 1, Vn, 1);
518  }
519 
520  Array<OneD, NekDouble > qtemp(nTracePts);
521 
522  // Evaulate Riemann flux
523  // qflux = \hat{q} \cdot u = q \cdot n
524  // Notice: i = 1 (first row of the viscous tensor is zero)
525  for (i = 1; i < nVariables; ++i)
526  {
527  qflux[i] = Array<OneD, NekDouble> (nTracePts, 0.0);
528  for (j = 0; j < nDim; ++j)
529  {
530  // Compute qFwd and qBwd value of qfield in position 'ji'
531  fields[i]->GetFwdBwdTracePhys(qfield[j][i], qFwd, qBwd);
532 
533  // Get Riemann flux of qflux --> LDG implies upwind
534  fields[i]->GetTrace()->Upwind(Vn, qBwd, qFwd, qfluxtemp);
535 
536  // Multiply the Riemann flux by the trace normals
537  Vmath::Vmul(nTracePts, m_traceNormals[j], 1, qfluxtemp, 1,
538  qfluxtemp, 1);
539 
540  // Extract the physical values of the solution at the boundaries
541  fields[i]->ExtractTracePhys(qfield[j][i], qtemp);
542 
543  // Impose weak boundary condition with flux
544  if (fields[0]->GetBndCondExpansions().num_elements())
545  {
546  v_WeakPenaltyO2(fields, i, j, qfield[j][i], qtemp, qfluxtemp);
547  }
548 
549  // Store the final flux into qflux
550  Vmath::Vadd(nTracePts, qfluxtemp, 1, qflux[i], 1,
551  qflux[i], 1);
552  }
553  }
554  }
555 
556 
557  /**
558  * @brief Imposes appropriate bcs for the 2nd order derivatives
559  *
560  */
563  const int var,
564  const int dir,
565  const Array<OneD, const NekDouble> &qfield,
566  const Array<OneD, const NekDouble> &qtemp,
567  Array<OneD, NekDouble> &penaltyflux)
568  {
569  int cnt = 0;
570  int nBndEdges, nBndEdgePts;
571  int i, e;
572  int id2;
573 
574  int nBndRegions = fields[var]->GetBndCondExpansions().num_elements();
575 
576  // Loop on the boundary regions to apply appropriate bcs
577  for (i = 0; i < nBndRegions; ++i)
578  {
579  // Number of boundary regions related to region 'i'
580  nBndEdges = fields[var]->
581  GetBndCondExpansions()[i]->GetExpSize();
582 
583  // Weakly impose bcs by modifying flux values
584  for (e = 0; e < nBndEdges; ++e)
585  {
586  nBndEdgePts = fields[var]->
587  GetBndCondExpansions()[i]->GetExp(e)->GetTotPoints();
588 
589  id2 = fields[0]->GetTrace()->
590  GetPhys_Offset(fields[0]->GetTraceMap()->
591  GetBndCondTraceToGlobalTraceMap(cnt++));
592 
593  // In case of Dirichlet bcs:
594  // uflux = gD
595  if(fields[var]->GetBndConditions()[i]->
596  GetBoundaryConditionType() == SpatialDomains::eDirichlet
597  && !boost::iequals(fields[var]->GetBndConditions()[i]->
598  GetUserDefined(), "WallAdiabatic"))
599  {
600  Vmath::Vmul(nBndEdgePts,
601  &m_traceNormals[dir][id2], 1,
602  &qtemp[id2], 1,
603  &penaltyflux[id2], 1);
604  }
605  // 3.4) In case of Neumann bcs:
606  // uflux = u+
607  else if((fields[var]->GetBndConditions()[i])->
608  GetBoundaryConditionType() == SpatialDomains::eNeumann)
609  {
610  ASSERTL0(false,
611  "Neumann bcs not implemented for LDGNS");
612 
613  /*
614  Vmath::Vmul(nBndEdgePts,
615  &m_traceNormals[dir][id2], 1,
616  &(fields[var]->
617  GetBndCondExpansions()[i]->
618  UpdatePhys())[id1], 1,
619  &penaltyflux[id2], 1);
620  */
621  }
622  else if(boost::iequals(fields[var]->GetBndConditions()[i]->
623  GetUserDefined(), "WallAdiabatic"))
624  {
625  if ((var == m_spaceDim + 1))
626  {
627  Vmath::Zero(nBndEdgePts, &penaltyflux[id2], 1);
628  }
629  else
630  {
631 
632  Vmath::Vmul(nBndEdgePts,
633  &m_traceNormals[dir][id2], 1,
634  &qtemp[id2], 1,
635  &penaltyflux[id2], 1);
636 
637  }
638  }
639  }
640  }
641  }
642  }//end of namespace SolverUtils
643 }//end of namespace Nektar
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Array< OneD, Array< OneD, NekDouble > > m_homoDerivs
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
static DiffusionSharedPtr create(std::string diffType)
virtual void v_Diffuse(const int nConvective, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const Array< OneD, Array< OneD, NekDouble > > &pFwd=NullNekDoubleArrayofArray, const Array< OneD, Array< OneD, NekDouble > > &pBwd=NullNekDoubleArrayofArray)
Calculate weak DG Diffusion in the LDG form for the Navier-Stokes (NS) equations: ...
DiffusionFactory & GetDiffusionFactory()
Definition: Diffusion.cpp:42
void Svtvp(int n, const T alpha, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
svtvp (scalar times vector plus vector): z = alpha*x + y
Definition: Vmath.cpp:485
virtual void v_WeakPenaltyO2(const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const int var, const int dir, const Array< OneD, const NekDouble > &qfield, const Array< OneD, const NekDouble > &qtemp, Array< OneD, NekDouble > &penaltyflux)
Imposes appropriate bcs for the 2nd order derivatives.
virtual void v_InitObject(LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields)
void Vdiv(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x/y.
Definition: Vmath.cpp:241
DiffusionFluxVecCBNS m_fluxVectorNS
Definition: Diffusion.h:130
LibUtilities::SessionReaderSharedPtr m_session
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:51
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_viscTensor
Array< OneD, Array< OneD, NekDouble > > m_traceVel
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
virtual void v_WeakPenaltyO1(const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &inarray, const Array< OneD, Array< OneD, NekDouble > > &uplus, Array< OneD, Array< OneD, NekDouble > > &penaltyfluxO1)
Imposes appropriate bcs for the 1st order derivatives.
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:396
virtual void v_NumericalFluxO2(const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &ufield, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &qfield, Array< OneD, Array< OneD, NekDouble > > &qflux)
Build the numerical flux for the 2nd order derivatives.
virtual void v_NumericalFluxO1(const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &numericalFluxO1, const Array< OneD, Array< OneD, NekDouble > > &pFwd=NullNekDoubleArrayofArray, const Array< OneD, Array< OneD, NekDouble > > &pBwd=NullNekDoubleArrayofArray)
Builds the numerical flux for the 1st order derivatives.
void Vsub(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Subtract vector z = x-y.
Definition: Vmath.cpp:343
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