# -*- coding: utf-8 -*-
import numpy as np
from pyfr.solvers.baseadvecdiff import BaseAdvectionDiffusionElements
from pyfr.solvers.euler.elements import BaseFluidElements
[docs]class NavierStokesElements(BaseFluidElements, BaseAdvectionDiffusionElements):
# Use the density field for shock sensing
shockvar = 'rho'
[docs] @staticmethod
def grad_con_to_pri(cons, grad_cons, cfg):
rho, *rhouvw = cons[:-1]
grad_rho, *grad_rhouvw, grad_E = grad_cons
# Divide momentum components by ρ
uvw = [rhov / rho for rhov in rhouvw]
# Velocity gradients: ∇u⃗ = 1/ρ·[∇(ρu⃗) - u⃗ ⊗ ∇ρ]
grad_uvw = [(grad_rhov - v*grad_rho) / rho
for grad_rhov, v in zip(grad_rhouvw, uvw)]
# Pressure gradient: ∇p = (γ - 1)·[∇E - 1/2*(u⃗·∇(ρu⃗) - ρu⃗·∇u⃗)]
gamma = cfg.getfloat('constants', 'gamma')
grad_p = grad_E - 0.5*(np.einsum('ijk,iljk->ljk', uvw, grad_rhouvw) +
np.einsum('ijk,iljk->ljk', rhouvw, grad_uvw))
grad_p *= (gamma - 1)
return [grad_rho] + grad_uvw + [grad_p]
[docs] def set_backend(self, *args, **kwargs):
super().set_backend(*args, **kwargs)
# Can elide interior flux calculations at p = 0
if self.basis.order == 0:
return
# Register our flux kernels
self._be.pointwise.register('pyfr.solvers.navstokes.kernels.tflux')
self._be.pointwise.register('pyfr.solvers.navstokes.kernels.tfluxlin')
# Handle shock capturing and Sutherland's law
shock_capturing = self.cfg.get('solver', 'shock-capturing')
visc_corr = self.cfg.get('solver', 'viscosity-correction', 'none')
if visc_corr not in {'sutherland', 'none'}:
raise ValueError('Invalid viscosity-correction option')
# Template parameters for the flux kernels
tplargs = {
'ndims': self.ndims,
'nvars': self.nvars,
'nverts': len(self.basis.linspts),
'c': self.cfg.items_as('constants', float),
'jac_exprs': self.basis.jac_exprs,
'shock_capturing': shock_capturing,
'visc_corr': visc_corr
}
# Helpers
c, l = 'curved', 'linear'
r, s = self._mesh_regions, self._slice_mat
av = self.artvisc
if c in r and 'flux' not in self.antialias:
self.kernels['tdisf_curved'] = lambda uin: self._be.kernel(
'tflux', tplargs=tplargs, dims=[self.nupts, r[c]],
u=s(self.scal_upts[uin], c), f=s(self._vect_upts, c),
artvisc=s(av, c), smats=self.curved_smat_at('upts')
)
elif c in r:
self.kernels['tdisf_curved'] = lambda: self._be.kernel(
'tflux', tplargs=tplargs, dims=[self.nqpts, r[c]],
u=s(self._scal_qpts, c), f=s(self._vect_qpts, c),
artvisc=s(av, c), smats=self.curved_smat_at('qpts')
)
if l in r and 'flux' not in self.antialias:
self.kernels['tdisf_linear'] = lambda uin: self._be.kernel(
'tfluxlin', tplargs=tplargs, dims=[self.nupts, r[l]],
u=s(self.scal_upts[uin], l), f=s(self._vect_upts, l),
artvisc=s(av, l), verts=self.ploc_at('linspts', l),
upts=self.upts
)
elif l in r:
self.kernels['tdisf_linear'] = lambda: self._be.kernel(
'tfluxlin', tplargs=tplargs, dims=[self.nqpts, r[l]],
u=s(self._scal_qpts, l), f=s(self._vect_qpts, l),
artvisc=s(av, l), verts=self.ploc_at('linspts', l),
upts=self.qpts
)