64 m_session->LoadParameter (
"thermalConductivity",
71 int nDim = pFields[0]->GetCoordim(0);
72 int nTracePts = pFields[0]->GetTrace()->GetTotPoints();
75 if (pSession->DefinesSolverInfo(
"HOMOGENEOUS"))
89 pFields[0]->GetTrace()->GetNormals(m_traceNormals);
106 const int nConvectiveFields,
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();
135 for (i = 0; i < nScalars; ++i)
146 for (j = 0; j < nDim; ++j)
148 for (i = 0; i < nScalars; ++i)
150 fields[i]->IProductWRTDerivBase (j, inarray[i], tmp1);
152 fields[i]->AddTraceIntegral (numericalFluxO1[j][i],
154 fields[i]->SetPhysState (
false);
155 fields[i]->MultiplyByElmtInvMass(tmp1, tmp1);
156 fields[i]->BwdTrans (tmp1, derivativesO1[j][i]);
163 for (i = 0; i < nScalars; ++i)
178 for (i = 0; i < nScalars+1; ++i)
184 for (i = 0; i < nConvectiveFields; ++i)
195 for (i = 0; i < nConvectiveFields; ++i)
199 for (j = 0; j < nDim; ++j)
201 fields[i]->IProductWRTDerivBase(j,
m_viscTensor[j][i], tmp1);
202 Vmath::Vadd(nCoeffs, tmp1, 1, tmp2[i], 1, 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]);
227 int nTracePts = fields[0]->GetTrace()->GetTotPoints();
228 int nScalars = inarray.num_elements();
229 int nDim = fields[0]->GetCoordim(0);
234 for(i = 0; i < nDim; ++i)
245 for (i = 0; i < nScalars; ++i)
252 fields[i]->GetFwdBwdTracePhys(inarray[i], Fwd, Bwd);
260 fields[0]->GetTrace()->Upwind(Vn, Fwd, Bwd, numflux[i]);
266 for (i = 0; i < nScalars; ++i)
269 fields[i]->ExtractTracePhys(inarray[i], uplus[i]);
273 if (fields[0]->GetBndCondExpansions().num_elements())
281 for (i = 0; i < nScalars; ++i)
284 numflux[i], 1, numericalFluxO1[j][i], 1);
303 int nBndEdgePts, nBndEdges, nBndRegions;
305 int nTracePts = fields[0]->GetTrace()->GetTotPoints();
306 int nScalars = inarray.num_elements();
315 for (i = 0; i < nScalars; ++i)
321 for (i = 0; i < nScalars-1; ++i)
326 nBndRegions = fields[i+1]->
327 GetBndCondExpansions().num_elements();
328 for (j = 0; j < nBndRegions; ++j)
330 nBndEdges = fields[i+1]->
331 GetBndCondExpansions()[j]->GetExpSize();
332 for (e = 0; e < nBndEdges; ++e)
334 nBndEdgePts = fields[i+1]->
335 GetBndCondExpansions()[j]->GetExp(e)->GetTotPoints();
338 GetBndCondExpansions()[j]->GetPhys_Offset(e);
340 id2 = fields[0]->GetTrace()->
341 GetPhys_Offset(fields[0]->GetTraceMap()->
342 GetBndCondTraceToGlobalTraceMap(cnt++));
345 if (boost::iequals(fields[i]->GetBndConditions()[j]->
346 GetUserDefined(),
"WallViscous") ||
347 boost::iequals(fields[i]->GetBndConditions()[j]->
348 GetUserDefined(),
"WallAdiabatic"))
351 &scalarVariables[i][id2], 1);
356 else if (fields[i]->GetBndConditions()[j]->
357 GetBoundaryConditionType() ==
361 &(fields[i+1]->GetBndCondExpansions()[j]->
362 UpdatePhys())[id1], 1,
363 &(fields[0]->GetBndCondExpansions()[j]->
364 UpdatePhys())[id1], 1,
365 &scalarVariables[i][id2], 1);
369 if (fields[i]->GetBndConditions()[j]->
370 GetBoundaryConditionType() ==
374 &scalarVariables[i][id2], 1,
375 &penaltyfluxO1[i][id2], 1);
379 else if ((fields[i]->GetBndConditions()[j])->
380 GetBoundaryConditionType() ==
385 &penaltyfluxO1[i][id2], 1);
390 &scalarVariables[i][id2], 1,
391 &scalarVariables[i][id2], 1,
408 nBndRegions = fields[nScalars]->
409 GetBndCondExpansions().num_elements();
410 for (j = 0; j < nBndRegions; ++j)
412 nBndEdges = fields[nScalars]->
413 GetBndCondExpansions()[j]->GetExpSize();
414 for (e = 0; e < nBndEdges; ++e)
416 nBndEdgePts = fields[nScalars]->
417 GetBndCondExpansions()[j]->GetExp(e)->GetTotPoints();
419 id1 = fields[nScalars]->
420 GetBndCondExpansions()[j]->GetPhys_Offset(e);
422 id2 = fields[0]->GetTrace()->
423 GetPhys_Offset(fields[0]->GetTraceMap()->
424 GetBndCondTraceToGlobalTraceMap(cnt++));
427 if (boost::iequals(fields[i]->GetBndConditions()[j]->
428 GetUserDefined(),
"WallViscous"))
432 &scalarVariables[nScalars-1][id2], 1);
436 else if (fields[i]->GetBndConditions()[j]->
437 GetBoundaryConditionType() ==
443 GetBndCondExpansions()[j]->
446 GetBndCondExpansions()[j]->
448 &scalarVariables[nScalars-1][id2], 1);
452 &scalarVariables[nScalars-1][id2], 1,
454 &scalarVariables[nScalars-1][id2], 1);
458 &scalarVariables[nScalars-1][id2], 1,
459 &scalarVariables[nScalars-1][id2], 1);
463 if (fields[nScalars]->GetBndConditions()[j]->
464 GetBoundaryConditionType() ==
467 fields[nScalars]->GetBndConditions()[j]
468 ->GetUserDefined(),
"WallAdiabatic"))
471 &scalarVariables[nScalars-1][id2], 1,
472 &penaltyfluxO1[nScalars-1][id2], 1);
477 else if (((fields[nScalars]->GetBndConditions()[j])->
478 GetBoundaryConditionType() ==
480 boost::iequals(fields[nScalars]->GetBndConditions()[j]->
481 GetUserDefined(),
"WallAdiabatic"))
484 &uplus[nScalars-1][id2], 1,
485 &penaltyfluxO1[nScalars-1][id2], 1);
503 int nTracePts = fields[0]->GetTrace()->GetTotPoints();
504 int nVariables = fields.num_elements();
505 int nDim = fields[0]->GetCoordim(0);
514 for(i = 0; i < nDim; ++i)
525 for (i = 1; i < nVariables; ++i)
528 for (j = 0; j < nDim; ++j)
531 fields[i]->GetFwdBwdTracePhys(qfield[j][i], qFwd, qBwd);
534 fields[i]->GetTrace()->Upwind(Vn, qBwd, qFwd, qfluxtemp);
541 fields[i]->ExtractTracePhys(qfield[j][i], qtemp);
544 if (fields[0]->GetBndCondExpansions().num_elements())
570 int nBndEdges, nBndEdgePts;
574 int nBndRegions = fields[var]->GetBndCondExpansions().num_elements();
577 for (i = 0; i < nBndRegions; ++i)
580 nBndEdges = fields[var]->
581 GetBndCondExpansions()[i]->GetExpSize();
584 for (e = 0; e < nBndEdges; ++e)
586 nBndEdgePts = fields[var]->
587 GetBndCondExpansions()[i]->GetExp(e)->GetTotPoints();
589 id2 = fields[0]->GetTrace()->
590 GetPhys_Offset(fields[0]->GetTraceMap()->
591 GetBndCondTraceToGlobalTraceMap(cnt++));
595 if(fields[var]->GetBndConditions()[i]->
597 && !boost::iequals(fields[var]->GetBndConditions()[i]->
598 GetUserDefined(),
"WallAdiabatic"))
603 &penaltyflux[id2], 1);
607 else if((fields[var]->GetBndConditions()[i])->
611 "Neumann bcs not implemented for LDGNS");
622 else if(boost::iequals(fields[var]->GetBndConditions()[i]->
623 GetUserDefined(),
"WallAdiabatic"))
635 &penaltyflux[id2], 1);
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Array< OneD, Array< OneD, NekDouble > > m_homoDerivs
#define ASSERTL0(condition, msg)
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()
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
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.
DiffusionFluxVecCBNS m_fluxVectorNS
LibUtilities::SessionReaderSharedPtr m_session
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
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.
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.
NekDouble m_thermalConductivity
void Neg(int n, T *x, const int incx)
Negate x = -x.
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.
std::string m_ViscosityType
static Array< OneD, Array< OneD, NekDouble > > NullNekDoubleArrayofArray
void Zero(int n, T *x, const int incx)
Zero vector.
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
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.
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.