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,
112 int nDim = fields[0]->GetCoordim(0);
113 int nScalars = inarray.num_elements();
114 int nPts = fields[0]->GetTotPoints();
115 int nCoeffs = fields[0]->GetNcoeffs();
116 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]);
225 int nTracePts = fields[0]->GetTrace()->GetTotPoints();
226 int nScalars = inarray.num_elements();
227 int nDim = fields[0]->GetCoordim(0);
233 for(i = 0; i < nDim; ++i)
244 for (i = 0; i < nScalars; ++i)
249 fields[i]->GetFwdBwdTracePhys(inarray[i], Fwd[i], Bwd[i]);
250 fields[0]->GetTrace()->Upwind(Vn, Fwd[i], Bwd[i], numflux[i]);
254 if (fields[0]->GetBndCondExpansions().num_elements())
262 for (i = 0; i < nScalars; ++i)
265 numflux[i], 1, numericalFluxO1[j][i], 1);
283 int nBndEdgePts, nBndEdges, nBndRegions;
285 int nTracePts = fields[0]->GetTrace()->GetTotPoints();
286 int nScalars = inarray.num_elements();
296 for (i = 0; i < nScalars; ++i)
301 fields[i]->ExtractTracePhys(inarray[i], uplus[i]);
305 for (i = 0; i < nScalars-1; ++i)
310 nBndRegions = fields[i+1]->
311 GetBndCondExpansions().num_elements();
312 for (j = 0; j < nBndRegions; ++j)
314 nBndEdges = fields[i+1]->
315 GetBndCondExpansions()[j]->GetExpSize();
316 for (e = 0; e < nBndEdges; ++e)
318 nBndEdgePts = fields[i+1]->
319 GetBndCondExpansions()[j]->GetExp(e)->GetTotPoints();
322 GetBndCondExpansions()[j]->GetPhys_Offset(e);
324 id2 = fields[0]->GetTrace()->
325 GetPhys_Offset(fields[0]->GetTraceMap()->
326 GetBndCondTraceToGlobalTraceMap(cnt++));
329 if (boost::iequals(fields[i]->GetBndConditions()[j]->
330 GetUserDefined(),
"WallViscous") ||
331 boost::iequals(fields[i]->GetBndConditions()[j]->
332 GetUserDefined(),
"WallAdiabatic"))
335 &scalarVariables[i][id2], 1);
340 else if (fields[i]->GetBndConditions()[j]->
341 GetBoundaryConditionType() ==
345 &(fields[i+1]->GetBndCondExpansions()[j]->
346 UpdatePhys())[id1], 1,
347 &(fields[0]->GetBndCondExpansions()[j]->
348 UpdatePhys())[id1], 1,
349 &scalarVariables[i][id2], 1);
353 if (fields[i]->GetBndConditions()[j]->
354 GetBoundaryConditionType() ==
358 &scalarVariables[i][id2], 1,
359 &penaltyfluxO1[i][id2], 1);
363 else if ((fields[i]->GetBndConditions()[j])->
364 GetBoundaryConditionType() ==
369 &penaltyfluxO1[i][id2], 1);
374 &scalarVariables[i][id2], 1,
375 &scalarVariables[i][id2], 1,
392 nBndRegions = fields[nScalars]->
393 GetBndCondExpansions().num_elements();
394 for (j = 0; j < nBndRegions; ++j)
396 nBndEdges = fields[nScalars]->
397 GetBndCondExpansions()[j]->GetExpSize();
398 for (e = 0; e < nBndEdges; ++e)
400 nBndEdgePts = fields[nScalars]->
401 GetBndCondExpansions()[j]->GetExp(e)->GetTotPoints();
403 id1 = fields[nScalars]->
404 GetBndCondExpansions()[j]->GetPhys_Offset(e);
406 id2 = fields[0]->GetTrace()->
407 GetPhys_Offset(fields[0]->GetTraceMap()->
408 GetBndCondTraceToGlobalTraceMap(cnt++));
411 if (boost::iequals(fields[i]->GetBndConditions()[j]->
412 GetUserDefined(),
"WallViscous"))
416 &scalarVariables[nScalars-1][id2], 1);
420 else if (fields[i]->GetBndConditions()[j]->
421 GetBoundaryConditionType() ==
427 GetBndCondExpansions()[j]->
430 GetBndCondExpansions()[j]->
432 &scalarVariables[nScalars-1][id2], 1);
436 &scalarVariables[nScalars-1][id2], 1,
438 &scalarVariables[nScalars-1][id2], 1);
442 &scalarVariables[nScalars-1][id2], 1,
443 &scalarVariables[nScalars-1][id2], 1);
447 if (fields[nScalars]->GetBndConditions()[j]->
448 GetBoundaryConditionType() ==
451 fields[nScalars]->GetBndConditions()[j]
452 ->GetUserDefined(),
"WallAdiabatic"))
455 &scalarVariables[nScalars-1][id2], 1,
456 &penaltyfluxO1[nScalars-1][id2], 1);
461 else if (((fields[nScalars]->GetBndConditions()[j])->
462 GetBoundaryConditionType() ==
464 boost::iequals(fields[nScalars]->GetBndConditions()[j]->
465 GetUserDefined(),
"WallAdiabatic"))
468 &uplus[nScalars-1][id2], 1,
469 &penaltyfluxO1[nScalars-1][id2], 1);
487 int nTracePts = fields[0]->GetTrace()->GetTotPoints();
488 int nVariables = fields.num_elements();
489 int nDim = fields[0]->GetCoordim(0);
500 for(i = 0; i < nDim; ++i)
509 for (i = 1; i < nVariables; ++i)
512 for (j = 0; j < nDim; ++j)
515 fields[i]->GetFwdBwdTracePhys(qfield[j][i], qFwd, qBwd);
518 fields[i]->GetTrace()->Upwind(Vn, qBwd, qFwd, qfluxtemp);
525 if (fields[0]->GetBndCondExpansions().num_elements())
550 int nBndEdges, nBndEdgePts;
554 int nTracePts = fields[0]->GetTrace()->GetTotPoints();
555 int nBndRegions = fields[var]->GetBndCondExpansions().num_elements();
561 fields[var]->ExtractTracePhys(qfield, qtemp);
564 for (i = 0; i < nBndRegions; ++i)
567 nBndEdges = fields[var]->
568 GetBndCondExpansions()[i]->GetExpSize();
571 for (e = 0; e < nBndEdges; ++e)
573 nBndEdgePts = fields[var]->
574 GetBndCondExpansions()[i]->GetExp(e)->GetTotPoints();
576 id2 = fields[0]->GetTrace()->
577 GetPhys_Offset(fields[0]->GetTraceMap()->
578 GetBndCondTraceToGlobalTraceMap(cnt++));
582 if(fields[var]->GetBndConditions()[i]->
584 && !boost::iequals(fields[var]->GetBndConditions()[i]->
585 GetUserDefined(),
"WallAdiabatic"))
590 &penaltyflux[id2], 1);
594 else if((fields[var]->GetBndConditions()[i])->
598 "Neumann bcs not implemented for LDGNS");
609 else if(boost::iequals(fields[var]->GetBndConditions()[i]->
610 GetUserDefined(),
"WallAdiabatic"))
622 &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_WeakPenaltyO2(const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const int var, const int dir, const Array< OneD, const NekDouble > &qfield, Array< OneD, NekDouble > &penaltyflux)
Imposes appropriate bcs for the 2nd order derivatives.
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_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.
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)
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.
virtual void v_WeakPenaltyO1(const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &penaltyfluxO1)
Imposes appropriate bcs for the 1st order derivatives.
std::string m_ViscosityType
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)
Calculate weak DG Diffusion in the LDG form for the Navier-Stokes (NS) equations: ...
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.