"""flux variable placeholder and related operators (Flux class definition)"""
from typing import Any, Tuple
import os
if os.environ.get('__USE_CUDA', '0') == '1':
import cupy as np
from cupy import ndarray
else:
import numpy as np
from numpy import ndarray
from g3d.defs import G3DWorld
from g3d.defs import G3DScheme
from g3d.defs import FMPTypes
from g3d.defs import ShockDampingMode
from g3d.defs import BCTypes
from g3d.defs import WALLTypes
from g3d.defs import IOFLOWTypes
from g3d.defs import EQNMode
from g3d.defs import DIMMode
from g3d.gas import Gas
from g3d.flow.placeholders import ConsField
from g3d.flow.domain import BoundaryConditionPlaceHolder
from g3d.flow.domain import InterfacePlaceHolder
from g3d.flow.domain import Domain
from g3d.flow.geo import Mesh
from g3d.flow.io import G3DIO
from g3d.flow.solver.models import LogLaw
from g3d.flow.solver.rhs import MOC
from g3d.flow.solver.rhs.fluxmodifiers import FluxScale
from g3d.flow.solver.rhs.fluxmodifiers import Q1DSource
from g3d.flow.solver.rhs.fluxmodifiers import Q3DSource
from g3d.flow.solver.rhs.fluxmodifiers import AxiMod
from g3d.flow.solver.rhs.fluxmodifiers import IRSM
from g3d.flow.domain.boco import BoCo
from g3d.flow.placeholders import FMP
[docs]
class Flux(ConsField):
def __init__(self, *args: Any, **kwargs: Any) -> None:
super().__init__(*args, **kwargs)
g3d_world = G3DWorld()
self.__gas = Gas()
self.__domain = Domain()
self.__mesh = Mesh()
self.__boundary_faces = BoundaryConditionPlaceHolder()
self.__connectivity_faces = InterfacePlaceHolder()
self.__moc = MOC()
self.__log_law = LogLaw()
self.__dt = np.zeros(0)
self.__voli = np.zeros(0)
self.__vol = np.zeros(0)
self.__flux_modifiers = []
# pre-scale modifiers
if g3d_world.options.dim_mode is DIMMode.Q1D:
q1d_source_fm = Q1DSource()
self.__flux_modifiers.append(q1d_source_fm)
if g3d_world.options.dim_mode is DIMMode.Q3D:
q3d_source_fm = Q3DSource()
self.__flux_modifiers.append(q3d_source_fm)
if g3d_world.options.dim_mode is DIMMode.G2D_AXISYMM:
axisymm_fm = AxiMod()
self.__flux_modifiers.append(axisymm_fm)
# flux scale
scale_fm = FluxScale()
self.__flux_modifiers.append(scale_fm)
# post-scale modifiers
if g3d_world.options.use_irsm and g3d_world.num.time.cfl > g3d_world.num.irsm.cfl_crit:
irsm_fm = IRSM()
self.__flux_modifiers.append(irsm_fm)
self.__define_local_vectorized_functions()
[docs]
def setup(self, **kwargs: Any) -> None:
"""
:param kwargs:
"""
super().setup(**kwargs)
self.__allocate_flux_placeholders()
def __allocate_flux_placeholders(self, n_faces: int = None) -> None:
g3d_world = G3DWorld()
if n_faces is None:
n_faces = g3d_world.sizes.n_faces
# generic placeholders
self.__tmp_1 = np.zeros(n_faces)
self.__tmp_2 = np.zeros(n_faces)
self.__tmp_3 = np.zeros(n_faces)
self.__tmp_4 = np.zeros(n_faces)
# allocate arrays used for calculation of convective fluxes
self.__r_l = np.zeros(n_faces)
self.__u_l = np.zeros(n_faces)
self.__v_l = np.zeros(n_faces)
self.__w_l = np.zeros(n_faces)
self.__p_l = np.zeros(n_faces)
self.__r_r = np.zeros(n_faces)
self.__u_r = np.zeros(n_faces)
self.__v_r = np.zeros(n_faces)
self.__w_r = np.zeros(n_faces)
self.__p_r = np.zeros(n_faces)
self.__r_a = np.zeros(n_faces)
self.__u_a = np.zeros(n_faces)
self.__v_a = np.zeros(n_faces)
self.__w_a = np.zeros(n_faces)
self.__c_a = np.zeros(n_faces)
self.__ev1 = np.zeros(n_faces)
self.__ev4 = np.zeros(n_faces)
self.__ev5 = np.zeros(n_faces)
self.__sxl = np.zeros(n_faces)
self.__syl = np.zeros(n_faces)
self.__szl = np.zeros(n_faces)
self.__sh = np.zeros(n_faces)
self.__shl = np.zeros(n_faces)
self.__sxh = np.zeros(n_faces)
self.__syh = np.zeros(n_faces)
self.__sx = np.zeros(n_faces)
self.__sy = np.zeros(n_faces)
self.__c2 = np.zeros(n_faces)
self.__r1 = np.zeros(n_faces)
self.__r2 = np.zeros(n_faces)
self.__r3 = np.zeros(n_faces)
self.__r4 = np.zeros(n_faces)
self.__r5 = np.zeros(n_faces)
self.__u_m = np.zeros(n_faces)
self.__v_m = np.zeros(n_faces)
self.__w_m = np.zeros(n_faces)
self.__p_m = np.zeros(n_faces)
self.__ro_m = np.zeros(n_faces)
self.__ru_m = np.zeros(n_faces)
self.__rv_m = np.zeros(n_faces)
self.__rw_m = np.zeros(n_faces)
self.__et_m = np.zeros(n_faces)
self.__c_filt = np.zeros(n_faces)
self.__s_vel = np.zeros(n_faces)
self.__f_ro = np.zeros(n_faces)
self.__f_ru = np.zeros(n_faces)
self.__f_rv = np.zeros(n_faces)
self.__f_rw = np.zeros(n_faces)
self.__f_et = np.zeros(n_faces)
self.__tvd_div = np.zeros(n_faces)
self.__tvd_alfa = np.zeros(n_faces)
self.__tvd_beta = np.zeros(n_faces)
self.__allocate_optional_placeholders(n_faces)
def __allocate_optional_placeholders(self, n_faces: int = None) -> None:
g3d_world = G3DWorld()
if n_faces is None:
n_faces = g3d_world.sizes.n_faces
# allocate face-gradient placeholders (if needed)
if (g3d_world.options.eqn_mode != EQNMode.EULER or
g3d_world.num.scheme.shock_damping_mode is ShockDampingMode.DUCROSS_FACE_GRADIENTS or
g3d_world.options.dim_mode is DIMMode.G2D_AXISYMM):
self.__q_diff = np.zeros((n_faces, 3))
self.__dudx = np.zeros(n_faces)
self.__dudy = np.zeros(n_faces)
self.__dudz = np.zeros(n_faces)
self.__dvdx = np.zeros(n_faces)
self.__dvdy = np.zeros(n_faces)
self.__dvdz = np.zeros(n_faces)
self.__dwdx = np.zeros(n_faces)
self.__dwdy = np.zeros(n_faces)
self.__dwdz = np.zeros(n_faces)
else:
self.__q_diff = np.zeros((0, 3))
self.__dudx = np.zeros(0)
self.__dudy = np.zeros(0)
self.__dudz = np.zeros(0)
self.__dvdx = np.zeros(0)
self.__dvdy = np.zeros(0)
self.__dvdz = np.zeros(0)
self.__dwdx = np.zeros(0)
self.__dwdy = np.zeros(0)
self.__dwdz = np.zeros(0)
# allocate arrays used for calculation of diffusive fluxes (if needed)
self.__dtdx = []
self.__dtdy = []
self.__dtdz = []
self.__tau_xx = []
self.__tau_xy = []
self.__tau_xz = []
self.__tau_yy = []
self.__tau_yz = []
self.__tau_zz = []
if g3d_world.options.eqn_mode != EQNMode.EULER:
self.__dtdx = np.zeros(n_faces)
self.__dtdy = np.zeros(n_faces)
self.__dtdz = np.zeros(n_faces)
self.__tau_xx = np.zeros(n_faces)
self.__tau_xy = np.zeros(n_faces)
self.__tau_xz = np.zeros(n_faces)
self.__tau_yy = np.zeros(n_faces)
self.__tau_yz = np.zeros(n_faces)
self.__tau_zz = np.zeros(n_faces)
[docs]
def change_eqn_mode(self) -> None:
"""
"""
self.__allocate_optional_placeholders()
[docs]
def associate(self, dom: Domain, dt: ndarray) -> None:
"""
:param dom:
:param dt:
"""
self.__domain = dom
self.__mesh = dom.mesh
self.__gas = dom.gas
self.__boundary_faces = dom.boundary_faces
self.__connectivity_faces = dom.connectivity_faces
self.__dt = dt
self.__voli = self.__mesh.get_inverse_volume()
self.__vol = self.__mesh.get_volume()
self.__moc.associate(self.__gas)
self.__log_law.associate(self.__gas)
for fm in self.__flux_modifiers:
fm.associate(self, dt, dom)
self.__define_local_vectorized_functions()
def __define_local_vectorized_functions(self) -> None:
# self.__v_evaluate_log_law = np.vectorize(self.__log_law.evaluate)
# self.__v_evaluate_moc = np.vectorize(self.__moc.evaluate)
# self.__v_gas_get_temperature = np.vectorize(self.__gas.get_temperature)
# self.__v_gas_get_density = np.vectorize(self.__gas.get_density)
# self.__v_gas_get_pressure = np.vectorize(self.__gas.get_pressure)
# self.__v_gas_get_speed_of_sound = np.vectorize(self.__gas.get_speed_of_sound)
# self.__v_gas_get_energy_from_primitive = np.vectorize(self.__gas.get_energy_from_primitive)
pass
[docs]
def evaluate(self) -> None:
"""
"""
# nullify flux
self.nullify()
# internal face fluxes
self.__internal_face_fluxes()
# connectivity face fluxes
self.__connectivity_face_fluxes()
# boundary face fluxes
self.__boundary_face_fluxes()
# evaluate flux modifiers
for fm in self.__flux_modifiers:
fm.evaluate()
def __internal_face_fluxes(self) -> None:
# internal face fluxes
g3d_world = G3DWorld()
internal = self.__domain.internal_faces
for tmp in internal.get_FMPs():
self.__conv_flux(tmp)
if g3d_world.options.eqn_mode != EQNMode.EULER:
self.__diff_flux(tmp)
def __connectivity_face_fluxes(self) -> None:
# connectivity face fluxes
g3d_world = G3DWorld()
for interface in self.__connectivity_faces.get_interface_array():
A, B = interface.get_face_defs()
for intf in A:
tmp = intf.FMP
self.__conv_flux(tmp)
if g3d_world.options.eqn_mode != EQNMode.EULER:
self.__diff_flux(tmp)
def __boundary_face_fluxes(self) -> None:
# boundary face fluxes
g3d_world = G3DWorld()
for bc in self.__boundary_faces.get_boundary_array():
for boco in bc.get_face_defs():
tmp = boco.FMP
if boco.face_type is BCTypes.STREAM_TUBE_WALLS:
# add streamtube wall pressure force
self.__add_streamtube_pressure_forces(tmp)
else:
self.__update_boundary_ghost_cells(boco)
self.__conv_flux(tmp)
if g3d_world.options.eqn_mode != EQNMode.EULER:
self.__diff_flux(tmp, int(tmp.size / 2)) # diffusive flux for face one step in from boundary
if boco.face_type is BCTypes.WALLS:
if boco.it is WALLTypes.VISCOUS_RESOLVED:
self.__wall_diff_flux_resolved(boco)
elif boco.it is WALLTypes.VISCOUS_WALL_FUNCTION:
self.__wall_diff_flux_wall_function(boco)
elif boco.face_type is BCTypes.SYMM:
self.__symm_diff_flux(boco)
def __update_boundary_ghost_cells(self, boco: BoCo) -> None:
fmp = boco.FMP
face_type = boco.face_type
it = boco.it
idx1, idx2, idx3 = np.int32([0, boco.get_face_count(), fmp.size])
ro, ru, rv, rw, et = self.__domain.get_cons_arrays_with_ghosts()
p, u, v, w, t, c = self.__domain.get_auxv_arrays_with_ghosts()
kp, kc = fmp.get_pointers()
s, sl = fmp.get_conv_metrics()
if (face_type is BCTypes.WALLS or face_type is BCTypes.SYMM or
(face_type is BCTypes.IOFLOW and it is IOFLOWTypes.UNDEFINED)):
self.__s_vel[idx1:idx2] = (u[kc[idx1:idx2, 2]] * s[idx2:idx3, 0] +
v[kc[idx1:idx2, 2]] * s[idx2:idx3, 1] +
w[kc[idx1:idx2, 2]] * s[idx2:idx3, 2]) / np.power(sl[idx2:idx3], 2)
p[kc[idx1:idx2, 3]] = p[kc[idx1:idx2, 2]]
u[kc[idx1:idx2, 3]] = u[kc[idx1:idx2, 2]] - 2.0 * self.__s_vel[idx1:idx2] * s[idx2:idx3, 0]
v[kc[idx1:idx2, 3]] = v[kc[idx1:idx2, 2]] - 2.0 * self.__s_vel[idx1:idx2] * s[idx2:idx3, 1]
w[kc[idx1:idx2, 3]] = w[kc[idx1:idx2, 2]] - 2.0 * self.__s_vel[idx1:idx2] * s[idx2:idx3, 2]
t[kc[idx1:idx2, 3]] = t[kc[idx1:idx2, 2]]
c[kc[idx1:idx2, 3]] = c[kc[idx1:idx2, 2]]
ro[kc[idx1:idx2, 3]] = ro[kc[idx1:idx2, 2]]
ru[kc[idx1:idx2, 3]] = (
ru[kc[idx1:idx2, 2]] - 2.0 * ro[kc[idx1:idx2, 2]] * self.__s_vel[idx1:idx2] * s[idx2:idx3, 0])
rv[kc[idx1:idx2, 3]] = (
rv[kc[idx1:idx2, 2]] - 2.0 * ro[kc[idx1:idx2, 2]] * self.__s_vel[idx1:idx2] * s[idx2:idx3, 1])
rw[kc[idx1:idx2, 3]] = (
rw[kc[idx1:idx2, 2]] - 2.0 * ro[kc[idx1:idx2, 2]] * self.__s_vel[idx1:idx2] * s[idx2:idx3, 2])
et[kc[idx1:idx2, 3]] = et[kc[idx1:idx2, 2]]
self.__s_vel[idx1:idx2] = (u[kc[idx1:idx2, 1]] * s[idx2:idx3, 0] +
v[kc[idx1:idx2, 1]] * s[idx2:idx3, 1] +
w[kc[idx1:idx2, 1]] * s[idx2:idx3, 2]) / np.power(sl[idx2:idx3], 2)
p[kc[idx2:idx3, 3]] = p[kc[idx1:idx2, 1]]
u[kc[idx2:idx3, 3]] = u[kc[idx1:idx2, 1]] - 2.0 * self.__s_vel[idx1:idx2] * s[idx2:idx3, 0]
v[kc[idx2:idx3, 3]] = v[kc[idx1:idx2, 1]] - 2.0 * self.__s_vel[idx1:idx2] * s[idx2:idx3, 1]
w[kc[idx2:idx3, 3]] = w[kc[idx1:idx2, 1]] - 2.0 * self.__s_vel[idx1:idx2] * s[idx2:idx3, 2]
t[kc[idx2:idx3, 3]] = t[kc[idx1:idx2, 1]]
c[kc[idx2:idx3, 3]] = c[kc[idx1:idx2, 1]]
ro[kc[idx2:idx3, 3]] = ro[kc[idx1:idx2, 1]]
ru[kc[idx2:idx3, 3]] = ru[kc[idx1:idx2, 1]] - 2.0 * ro[kc[idx1:idx2, 1]] * self.__s_vel[idx1:idx2] * s[
idx2:idx3, 0]
rv[kc[idx2:idx3, 3]] = rv[kc[idx1:idx2, 1]] - 2.0 * ro[kc[idx1:idx2, 1]] * self.__s_vel[idx1:idx2] * s[
idx2:idx3, 1]
rw[kc[idx2:idx3, 3]] = rw[kc[idx1:idx2, 1]] - 2.0 * ro[kc[idx1:idx2, 1]] * self.__s_vel[idx1:idx2] * s[
idx2:idx3, 2]
et[kc[idx2:idx3, 3]] = et[kc[idx1:idx2, 1]]
return
elif face_type is BCTypes.IOFLOW:
if it is IOFLOWTypes.ABSORBING:
# MOC in/out flow
self.__moc_update_metrics(s, sl, idx1, idx2)
self.__r_r[idx1:idx2] = boco.ro
self.__u_r[idx1:idx2] = boco.u
self.__v_r[idx1:idx2] = boco.v
self.__w_r[idx1:idx2] = boco.w
self.__p_r[idx1:idx2] = boco.p
self.__r_a[idx1:idx2] = 0.5 * (ro[kc[idx1:idx2, 2]] + boco.ro)
self.__u_a[idx1:idx2] = 0.5 * (u[kc[idx1:idx2, 2]] + boco.u)
self.__v_a[idx1:idx2] = 0.5 * (v[kc[idx1:idx2, 2]] + boco.v)
self.__w_a[idx1:idx2] = 0.5 * (w[kc[idx1:idx2, 2]] + boco.w)
self.__c_a[idx1:idx2] = 0.5 * (c[kc[idx1:idx2, 2]] + boco.c)
self.__r_l[idx1:idx2] = 2. * ro[kc[idx1:idx2, 2]] - ro[kc[idx1:idx2, 1]]
self.__u_l[idx1:idx2] = 2. * u[kc[idx1:idx2, 2]] - u[kc[idx1:idx2, 1]]
self.__v_l[idx1:idx2] = 2. * v[kc[idx1:idx2, 2]] - v[kc[idx1:idx2, 1]]
self.__w_l[idx1:idx2] = 2. * w[kc[idx1:idx2, 2]] - w[kc[idx1:idx2, 1]]
self.__p_l[idx1:idx2] = 2. * p[kc[idx1:idx2, 2]] - p[kc[idx1:idx2, 1]]
self.__moc_evaluate(s, sl, idx1, idx2)
ro[kc[idx1:idx2, 3]] = self.__ro_m[idx1:idx2]
ru[kc[idx1:idx2, 3]] = self.__ru_m[idx1:idx2]
rv[kc[idx1:idx2, 3]] = self.__rv_m[idx1:idx2]
rw[kc[idx1:idx2, 3]] = self.__rw_m[idx1:idx2]
et[kc[idx1:idx2, 3]] = self.__et_m[idx1:idx2]
p[kc[idx1:idx2, 3]] = self.__p_m[idx1:idx2]
u[kc[idx1:idx2, 3]] = self.__u_m[idx1:idx2]
v[kc[idx1:idx2, 3]] = self.__v_m[idx1:idx2]
w[kc[idx1:idx2, 3]] = self.__w_m[idx1:idx2]
t[kc[idx1:idx2, 3]] = self.__gas.get_temperature(p[kc[idx1:idx2, 3]], ro[kc[idx1:idx2, 3]])
c[kc[idx1:idx2, 3]] = self.__gas.get_speed_of_sound(t[kc[idx1:idx2, 3]])
self.__r_l[idx1:idx2] = 3. * ro[kc[idx1:idx2, 2]] - 2. * ro[kc[idx1:idx2, 1]]
self.__u_l[idx1:idx2] = 3. * u[kc[idx1:idx2, 2]] - 2. * u[kc[idx1:idx2, 1]]
self.__v_l[idx1:idx2] = 3. * v[kc[idx1:idx2, 2]] - 2. * v[kc[idx1:idx2, 1]]
self.__w_l[idx1:idx2] = 3. * w[kc[idx1:idx2, 2]] - 2. * w[kc[idx1:idx2, 1]]
self.__p_l[idx1:idx2] = 3. * p[kc[idx1:idx2, 2]] - 2. * p[kc[idx1:idx2, 1]]
self.__moc_evaluate(s, sl, idx1, idx2)
ro[kc[idx2:idx3, 3]] = self.__ro_m[idx1:idx2]
ru[kc[idx2:idx3, 3]] = self.__ru_m[idx1:idx2]
rv[kc[idx2:idx3, 3]] = self.__rv_m[idx1:idx2]
rw[kc[idx2:idx3, 3]] = self.__rw_m[idx1:idx2]
et[kc[idx2:idx3, 3]] = self.__et_m[idx1:idx2]
p[kc[idx2:idx3, 3]] = self.__p_m[idx1:idx2]
u[kc[idx2:idx3, 3]] = self.__u_m[idx1:idx2]
v[kc[idx2:idx3, 3]] = self.__v_m[idx1:idx2]
w[kc[idx2:idx3, 3]] = self.__w_m[idx1:idx2]
t[kc[idx2:idx3, 3]] = self.__gas.get_temperature(p[kc[idx2:idx3, 3]], ro[kc[idx2:idx3, 3]])
c[kc[idx2:idx3, 3]] = self.__gas.get_speed_of_sound(t[kc[idx2:idx3, 3]])
return
elif it is IOFLOWTypes.SUBSONIC_INFLOW:
# subsonic inflow
vne = fmp.get_tmp_storage(0)
vns = fmp.get_tmp_storage(1)
vns[idx1:idx2] = (boco.u * s[idx2:idx3, 0] +
boco.v * s[idx2:idx3, 1] +
boco.w * s[idx2:idx3, 2]) / sl[idx2:idx3]
vne[idx1:idx2] = (u[kc[idx1:idx2, 2]] * s[idx2:idx3, 0] +
v[kc[idx1:idx2, 2]] * s[idx2:idx3, 1] +
w[kc[idx1:idx2, 2]] * s[idx2:idx3, 2]) / sl[idx2:idx3]
u[kc[idx1:idx2, 3]] = u[kc[idx2:idx3, 3]] = (
(vne[idx1:idx2] - vns[idx1:idx2]) * s[idx2:idx3, 0] / sl[idx2:idx3] + boco.u)
v[kc[idx1:idx2, 3]] = v[kc[idx2:idx3, 3]] = (
(vne[idx1:idx2] - vns[idx1:idx2]) * s[idx2:idx3, 1] / sl[idx2:idx3] + boco.v)
w[kc[idx1:idx2, 3]] = w[kc[idx2:idx3, 3]] = (
(vne[idx1:idx2] - vns[idx1:idx2]) * s[idx2:idx3, 2] / sl[idx2:idx3] + boco.w)
t[kc[idx1:idx2, 3]] = t[kc[idx2:idx3, 3]] = boco.t0 - 0.5 * (
u[kc[idx1:idx2, 3]] * u[kc[idx1:idx2, 3]] +
v[kc[idx1:idx2, 3]] * v[kc[idx1:idx2, 3]] +
w[kc[idx1:idx2, 3]] * w[kc[idx1:idx2, 3]]) / self.__gas.Cp
p[kc[idx1:idx2, 3]] = p[kc[idx2:idx3, 3]] = (
boco.p0 * np.power(t[kc[idx1:idx2, 3]] / boco.t0, self.__gas.g_over_gm1))
c[kc[idx1:idx2, 3]] = c[kc[idx2:idx3, 3]] = (
self.__gas.get_speed_of_sound(t[kc[idx1:idx2, 3]]))
ro[kc[idx1:idx2, 3]] = ro[kc[idx2:idx3, 3]] = (
self.__gas.get_density(p[kc[idx1:idx2, 3]], t[kc[idx1:idx2, 3]]))
ru[kc[idx1:idx2, 3]] = ru[kc[idx2:idx3, 3]] = ro[kc[idx1:idx2, 3]] * u[kc[idx1:idx2, 3]]
rv[kc[idx1:idx2, 3]] = rv[kc[idx2:idx3, 3]] = ro[kc[idx1:idx2, 3]] * v[kc[idx1:idx2, 3]]
rw[kc[idx1:idx2, 3]] = rw[kc[idx2:idx3, 3]] = ro[kc[idx1:idx2, 3]] * w[kc[idx1:idx2, 3]]
et[kc[idx1:idx2, 3]] = et[kc[idx2:idx3, 3]] = self.__gas.get_energy_from_primitive(
ro[kc[idx1:idx2, 3]],
p[kc[idx1:idx2, 3]],
u[kc[idx1:idx2, 3]],
v[kc[idx1:idx2, 3]],
w[kc[idx1:idx2, 3]])
return
elif it is IOFLOWTypes.SUBSONIC_OUTFLOW:
# subsonic outflow
p[kc[idx1:idx2, 3]] = p[kc[idx2:idx3, 3]] = boco.p
u[kc[idx1:idx2, 3]] = u[kc[idx2:idx3, 3]] = u[kc[idx1:idx2, 2]]
v[kc[idx1:idx2, 3]] = v[kc[idx2:idx3, 3]] = v[kc[idx1:idx2, 2]]
w[kc[idx1:idx2, 3]] = w[kc[idx2:idx3, 3]] = w[kc[idx1:idx2, 2]]
t[kc[idx1:idx2, 3]] = t[kc[idx2:idx3, 3]] = t[kc[idx1:idx2, 2]]
c[kc[idx1:idx2, 3]] = c[kc[idx2:idx3, 3]] = c[kc[idx1:idx2, 2]]
ro[kc[idx1:idx2, 3]] = ro[kc[idx2:idx3, 3]] = ro[kc[idx1:idx2, 2]]
ru[kc[idx1:idx2, 3]] = ru[kc[idx2:idx3, 3]] = ru[kc[idx1:idx2, 2]]
rv[kc[idx1:idx2, 3]] = rv[kc[idx2:idx3, 3]] = rv[kc[idx1:idx2, 2]]
rw[kc[idx1:idx2, 3]] = rw[kc[idx2:idx3, 3]] = rw[kc[idx1:idx2, 2]]
et[kc[idx1:idx2, 3]] = et[kc[idx2:idx3, 3]] = self.__gas.get_energy_from_primitive(
ro[kc[idx1:idx2, 3]],
p[kc[idx1:idx2, 3]],
u[kc[idx1:idx2, 3]],
v[kc[idx1:idx2, 3]],
w[kc[idx1:idx2, 3]])
return
elif it is IOFLOWTypes.SUPERSONIC_INFLOW:
# supersonic inflow
p[kc[idx1:idx2, 3]] = p[kc[idx2:idx3, 3]] = boco.p
u[kc[idx1:idx2, 3]] = u[kc[idx2:idx3, 3]] = boco.u
v[kc[idx1:idx2, 3]] = v[kc[idx2:idx3, 3]] = boco.v
w[kc[idx1:idx2, 3]] = w[kc[idx2:idx3, 3]] = boco.w
t[kc[idx1:idx2, 3]] = t[kc[idx2:idx3, 3]] = boco.t
c[kc[idx1:idx2, 3]] = c[kc[idx2:idx3, 3]] = boco.c
ro[kc[idx1:idx2, 3]] = ro[kc[idx2:idx3, 3]] = boco.ro
ru[kc[idx1:idx2, 3]] = ru[kc[idx2:idx3, 3]] = boco.ru
rv[kc[idx1:idx2, 3]] = rv[kc[idx2:idx3, 3]] = boco.rv
rw[kc[idx1:idx2, 3]] = rw[kc[idx2:idx3, 3]] = boco.rw
et[kc[idx1:idx2, 3]] = et[kc[idx2:idx3, 3]] = boco.et
return
elif it is IOFLOWTypes.SUPERSONIC_OUTFLOW:
# supersonic outflow
p[kc[idx1:idx2, 3]] = p[kc[idx2:idx3, 3]] = p[kc[idx1:idx2, 2]]
u[kc[idx1:idx2, 3]] = u[kc[idx2:idx3, 3]] = u[kc[idx1:idx2, 2]]
v[kc[idx1:idx2, 3]] = v[kc[idx2:idx3, 3]] = v[kc[idx1:idx2, 2]]
w[kc[idx1:idx2, 3]] = w[kc[idx2:idx3, 3]] = w[kc[idx1:idx2, 2]]
t[kc[idx1:idx2, 3]] = t[kc[idx2:idx3, 3]] = t[kc[idx1:idx2, 2]]
c[kc[idx1:idx2, 3]] = c[kc[idx2:idx3, 3]] = c[kc[idx1:idx2, 2]]
ro[kc[idx1:idx2, 3]] = ro[kc[idx2:idx3, 3]] = ro[kc[idx1:idx2, 2]]
ru[kc[idx1:idx2, 3]] = ru[kc[idx2:idx3, 3]] = ru[kc[idx1:idx2, 2]]
rv[kc[idx1:idx2, 3]] = rv[kc[idx2:idx3, 3]] = rv[kc[idx1:idx2, 2]]
rw[kc[idx1:idx2, 3]] = rw[kc[idx2:idx3, 3]] = rw[kc[idx1:idx2, 2]]
et[kc[idx1:idx2, 3]] = et[kc[idx2:idx3, 3]] = et[kc[idx1:idx2, 2]]
return
elif it is IOFLOWTypes.MIXED_INFLOW_OUTFLOW:
# mixed inflow/outflow
vne = fmp.get_tmp_storage(0)
vns = fmp.get_tmp_storage(1)
vne[idx1:idx2] = (u[kc[idx1:idx2, 2]] * s[idx2:idx3, 0] +
v[kc[idx1:idx2, 2]] * s[idx2:idx3, 1] +
w[kc[idx1:idx2, 2]] * s[idx2:idx3, 2]) / sl[idx2:idx3]
n_faces = idx2 - idx1
face_indices = np.where(
np.logical_and(vne[idx1:idx2] <= 0., np.abs(vne[idx1:idx2]) <= c[kc[idx1:idx2, 2]]))[0] + idx1
if len(face_indices):
n_faces -= len(face_indices)
# subsonic inflow
face_indices_shift = face_indices + idx2
vns[face_indices] = (boco.u * s[face_indices_shift, 0] +
boco.v * s[face_indices_shift, 1] +
boco.w * s[face_indices_shift, 2]) / sl[face_indices_shift]
u[kc[face_indices, 3]] = u[kc[face_indices_shift, 3]] = (
(vne[face_indices] - vns[face_indices]) * s[face_indices_shift, 0] /
sl[face_indices_shift] + boco.u)
v[kc[face_indices, 3]] = v[kc[face_indices_shift, 3]] = (
(vne[face_indices] - vns[face_indices]) * s[face_indices_shift, 1] /
sl[face_indices_shift] + boco.v)
w[kc[face_indices, 3]] = w[kc[face_indices_shift, 3]] = (
(vne[face_indices] - vns[face_indices]) * s[face_indices_shift, 2] /
sl[face_indices_shift] + boco.w)
t[kc[face_indices, 3]] = t[kc[face_indices_shift, 3]] = boco.t0 - 0.5 * (
u[kc[face_indices, 3]] * u[kc[face_indices, 3]] +
v[kc[face_indices, 3]] * v[kc[face_indices, 3]] +
w[kc[face_indices, 3]] * w[kc[face_indices, 3]]) / self.__gas.Cp
p[kc[face_indices, 3]] = p[kc[face_indices_shift, 3]] = (
boco.p0 * np.power(t[kc[face_indices, 3]] / boco.t0, self.__gas.g_over_gm1))
c[kc[face_indices, 3]] = c[kc[face_indices_shift, 3]] = self.__gas.get_speed_of_sound(
t[kc[face_indices, 3]])
ro[kc[face_indices, 3]] = ro[kc[face_indices_shift, 3]] = (
self.__gas.get_density(p[kc[face_indices, 3]], t[kc[face_indices, 3]]))
ru[kc[face_indices, 3]] = ru[kc[face_indices_shift, 3]] = (
ro[kc[face_indices, 3]] * u[kc[face_indices, 3]])
rv[kc[face_indices, 3]] = rv[kc[face_indices_shift, 3]] = (
ro[kc[face_indices, 3]] * v[kc[face_indices, 3]])
rw[kc[face_indices, 3]] = rw[kc[face_indices_shift, 3]] = (
ro[kc[face_indices, 3]] * w[kc[face_indices, 3]])
et[kc[face_indices, 3]] = et[kc[face_indices_shift, 3]] = self.__gas.get_energy_from_primitive(
ro[kc[face_indices, 3]],
p[kc[face_indices, 3]],
u[kc[face_indices, 3]],
v[kc[face_indices, 3]],
w[kc[face_indices, 3]])
if n_faces == 0:
return
if n_faces > 0:
face_indices = np.where(np.logical_and(vne[idx1:idx2] <= 0.,
np.abs(vne[idx1:idx2]) > c[kc[idx1:idx2, 2]]))[0] + idx1
if len(face_indices):
n_faces -= len(face_indices)
face_indices_shift = face_indices + idx2
# supersonic inflow
p[kc[face_indices, 3]] = p[kc[face_indices_shift, 3]] = boco.p
u[kc[face_indices, 3]] = u[kc[face_indices_shift, 3]] = boco.u
v[kc[face_indices, 3]] = v[kc[face_indices_shift, 3]] = boco.v
w[kc[face_indices, 3]] = w[kc[face_indices_shift, 3]] = boco.w
t[kc[face_indices, 3]] = t[kc[face_indices_shift, 3]] = boco.t
c[kc[face_indices, 3]] = c[kc[face_indices_shift, 3]] = boco.c
ro[kc[face_indices, 3]] = ro[kc[face_indices_shift, 3]] = boco.ro
ru[kc[face_indices, 3]] = ru[kc[face_indices_shift, 3]] = boco.ru
rv[kc[face_indices, 3]] = rv[kc[face_indices_shift, 3]] = boco.rv
rw[kc[face_indices, 3]] = rw[kc[face_indices_shift, 3]] = boco.rw
et[kc[face_indices, 3]] = et[kc[face_indices_shift, 3]] = boco.et
if n_faces == 0:
return
if n_faces > 0:
face_indices = (
np.where(np.logical_and(
vne[idx1:idx2] > 0., np.abs(vne[idx1:idx2]) <= c[kc[idx1:idx2, 2]]))[0] + idx1)
if len(face_indices):
n_faces -= len(face_indices)
face_indices_shift = face_indices + idx2
# subsonic outflow
p[kc[face_indices, 3]] = p[kc[face_indices_shift, 3]] = boco.p
u[kc[face_indices, 3]] = u[kc[face_indices_shift, 3]] = u[kc[face_indices, 2]]
v[kc[face_indices, 3]] = v[kc[face_indices_shift, 3]] = v[kc[face_indices, 2]]
w[kc[face_indices, 3]] = w[kc[face_indices_shift, 3]] = w[kc[face_indices, 2]]
t[kc[face_indices, 3]] = t[kc[face_indices_shift, 3]] = t[kc[face_indices, 2]]
c[kc[face_indices, 3]] = c[kc[face_indices_shift, 3]] = c[kc[face_indices, 2]]
ro[kc[face_indices, 3]] = ro[kc[face_indices_shift, 3]] = ro[kc[face_indices, 2]]
ru[kc[face_indices, 3]] = ru[kc[face_indices_shift, 3]] = ru[kc[face_indices, 2]]
rv[kc[face_indices, 3]] = rv[kc[face_indices_shift, 3]] = rv[kc[face_indices, 2]]
rw[kc[face_indices, 3]] = rw[kc[face_indices_shift, 3]] = rw[kc[face_indices, 2]]
et[kc[face_indices, 3]] = et[kc[face_indices_shift, 3]] = self.__gas.get_energy_from_primitive(
ro[kc[face_indices, 3]],
p[kc[face_indices, 3]],
u[kc[face_indices, 3]],
v[kc[face_indices, 3]],
w[kc[face_indices, 3]])
if n_faces == 0:
return
if n_faces > 0:
face_indices = (
np.where(np.logical_and(
vne[idx1:idx2] > 0., np.abs(vne[idx1:idx2]) > c[kc[idx1:idx2, 2]]))[0] + idx1)
if len(face_indices):
n_faces -= len(face_indices)
face_indices_shift = face_indices + idx2
# supersonic outflow
p[kc[face_indices, 3]] = p[kc[face_indices_shift, 3]] = p[kc[face_indices, 2]]
u[kc[face_indices, 3]] = u[kc[face_indices_shift, 3]] = u[kc[face_indices, 2]]
v[kc[face_indices, 3]] = v[kc[face_indices_shift, 3]] = v[kc[face_indices, 2]]
w[kc[face_indices, 3]] = w[kc[face_indices_shift, 3]] = w[kc[face_indices, 2]]
t[kc[face_indices, 3]] = t[kc[face_indices_shift, 3]] = t[kc[face_indices, 2]]
c[kc[face_indices, 3]] = c[kc[face_indices_shift, 3]] = c[kc[face_indices, 2]]
ro[kc[face_indices, 3]] = ro[kc[face_indices_shift, 3]] = ro[kc[face_indices, 2]]
ru[kc[face_indices, 3]] = ru[kc[face_indices_shift, 3]] = ru[kc[face_indices, 2]]
rv[kc[face_indices, 3]] = rv[kc[face_indices_shift, 3]] = rv[kc[face_indices, 2]]
rw[kc[face_indices, 3]] = rw[kc[face_indices_shift, 3]] = rw[kc[face_indices, 2]]
et[kc[face_indices, 3]] = et[kc[face_indices_shift, 3]] = et[kc[face_indices, 2]]
if n_faces == 0:
return
return
elif it is IOFLOWTypes.MIXED_SUBSONIC_SUPERSONIC_OUTFLOW:
# mixed subsonic/supersonic outflow
vne = fmp.get_tmp_storage(0)
# vns = fmp.get_tmp_storage(1)
vne[idx1:idx2] = (u[kc[idx1:idx2, 2]] * s[idx2:idx3, 0] +
v[kc[idx1:idx2, 2]] * s[idx2:idx3, 1] +
w[kc[idx1:idx2, 2]] * s[idx2:idx3, 2]) / sl[idx2:idx3]
n_faces = idx2 - idx1
face_indices = np.where(np.abs(vne[idx1:idx2]) <= c[kc[idx1:idx2, 2]])[0] + idx1
if len(face_indices):
n_faces -= len(face_indices)
face_indices_shift = face_indices + idx2
# subsonic outflow
p[kc[face_indices, 3]] = p[kc[face_indices_shift, 3]] = boco.p
u[kc[face_indices, 3]] = u[kc[face_indices_shift, 3]] = u[kc[face_indices, 2]]
v[kc[face_indices, 3]] = v[kc[face_indices_shift, 3]] = v[kc[face_indices, 2]]
w[kc[face_indices, 3]] = w[kc[face_indices_shift, 3]] = w[kc[face_indices, 2]]
t[kc[face_indices, 3]] = t[kc[face_indices_shift, 3]] = t[kc[face_indices, 2]]
c[kc[face_indices, 3]] = c[kc[face_indices_shift, 3]] = c[kc[face_indices, 2]]
ro[kc[face_indices, 3]] = ro[kc[face_indices_shift, 3]] = ro[kc[face_indices, 2]]
ru[kc[face_indices, 3]] = ru[kc[face_indices_shift, 3]] = ru[kc[face_indices, 2]]
rv[kc[face_indices, 3]] = rv[kc[face_indices_shift, 3]] = rv[kc[face_indices, 2]]
rw[kc[face_indices, 3]] = rw[kc[face_indices_shift, 3]] = rw[kc[face_indices, 2]]
et[kc[face_indices, 3]] = et[kc[face_indices_shift, 3]] = (self.__gas.get_energy_from_primitive(
ro[kc[face_indices, 3]],
p[kc[face_indices, 3]],
u[kc[face_indices, 3]],
v[kc[face_indices, 3]],
w[kc[face_indices, 3]]))
if n_faces == 0:
return
if n_faces > 0:
face_indices = np.where(np.abs(vne[idx1:idx2]) > c[kc[idx1:idx2, 2]])[0] + idx1
if len(face_indices):
n_faces -= len(face_indices)
face_indices_shift = face_indices + idx2
# supersonic outflow
p[kc[face_indices, 3]] = p[kc[face_indices_shift, 3]] = p[kc[face_indices, 2]]
u[kc[face_indices, 3]] = u[kc[face_indices_shift, 3]] = u[kc[face_indices, 2]]
v[kc[face_indices, 3]] = v[kc[face_indices_shift, 3]] = v[kc[face_indices, 2]]
w[kc[face_indices, 3]] = w[kc[face_indices_shift, 3]] = w[kc[face_indices, 2]]
t[kc[face_indices, 3]] = t[kc[face_indices_shift, 3]] = t[kc[face_indices, 2]]
c[kc[face_indices, 3]] = c[kc[face_indices_shift, 3]] = c[kc[face_indices, 2]]
ro[kc[face_indices, 3]] = ro[kc[face_indices_shift, 3]] = ro[kc[face_indices, 2]]
ru[kc[face_indices, 3]] = ru[kc[face_indices_shift, 3]] = ru[kc[face_indices, 2]]
rv[kc[face_indices, 3]] = rv[kc[face_indices_shift, 3]] = rv[kc[face_indices, 2]]
rw[kc[face_indices, 3]] = rw[kc[face_indices_shift, 3]] = rw[kc[face_indices, 2]]
et[kc[face_indices, 3]] = et[kc[face_indices_shift, 3]] = et[kc[face_indices, 2]]
if n_faces == 0:
return
return
elif face_type is BCTypes.SING:
p[kc[idx1:idx2, 3]] = p[kc[idx1:idx2, 2]]
u[kc[idx1:idx2, 3]] = u[kc[idx1:idx2, 2]]
v[kc[idx1:idx2, 3]] = v[kc[idx1:idx2, 2]]
w[kc[idx1:idx2, 3]] = w[kc[idx1:idx2, 2]]
t[kc[idx1:idx2, 3]] = t[kc[idx1:idx2, 2]]
c[kc[idx1:idx2, 3]] = c[kc[idx1:idx2, 2]]
ro[kc[idx1:idx2, 3]] = ro[kc[idx1:idx2, 2]]
ru[kc[idx1:idx2, 3]] = ru[kc[idx1:idx2, 2]]
rv[kc[idx1:idx2, 3]] = rv[kc[idx1:idx2, 2]]
rw[kc[idx1:idx2, 3]] = rw[kc[idx1:idx2, 2]]
et[kc[idx1:idx2, 3]] = et[kc[idx1:idx2, 2]]
return
def __calculate_face_gradients(self,
Q: ndarray, fmp: FMP,
idx1: int = 0, idx2: int = None) -> Tuple[ndarray, ndarray, ndarray]:
kp, kc = fmp.get_pointers()
sx, sy, sz, si, voli = fmp.get_diff_metrics()
if idx2 is None:
idx2 = fmp.size
self.__q_diff[idx1:idx2, 0] = Q[kc[idx1:idx2, 2]] - Q[kc[idx1:idx2, 1]]
self.__q_diff[idx1:idx2, 1] = fmp.fctj[idx1:idx2] * (
Q[kc[idx1:idx2, 6]] + Q[kc[idx1:idx2, 7]] - Q[kc[idx1:idx2, 4]] - Q[kc[idx1:idx2, 5]])
self.__q_diff[idx1:idx2, 2] = fmp.fctk[idx1:idx2] * (
Q[kc[idx1:idx2, 10]] + Q[kc[idx1:idx2, 11]] - Q[kc[idx1:idx2, 8]] - Q[kc[idx1:idx2, 9]])
return (self.__q_diff[idx1:idx2, 0] * sx[idx1:idx2, 0] +
self.__q_diff[idx1:idx2, 1] * sx[idx1:idx2, 1] +
self.__q_diff[idx1:idx2, 2] * sx[idx1:idx2, 2]) * voli[idx1:idx2], \
(self.__q_diff[idx1:idx2, 0] * sy[idx1:idx2, 0] +
self.__q_diff[idx1:idx2, 1] * sy[idx1:idx2, 1] +
self.__q_diff[idx1:idx2, 2] * sy[idx1:idx2, 2]) * voli[idx1:idx2], \
(self.__q_diff[idx1:idx2, 0] * sz[idx1:idx2, 0] +
self.__q_diff[idx1:idx2, 1] * sz[idx1:idx2, 1] +
self.__q_diff[idx1:idx2, 2] * sz[idx1:idx2, 2]) * voli[idx1:idx2]
def __moc_update_metrics(self, s: ndarray, sl: ndarray, idx1: int, idx2: int) -> None:
self.__sxl[idx1:idx2] = s[idx1:idx2, 0] / sl[idx1:idx2]
self.__syl[idx1:idx2] = s[idx1:idx2, 1] / sl[idx1:idx2]
self.__szl[idx1:idx2] = s[idx1:idx2, 2] / sl[idx1:idx2]
self.__sh[:] = 0.
self.__sh[idx1:idx2] = np.sqrt(np.power(s[idx1:idx2, 0], 2) + np.power(s[idx1:idx2, 1], 2))
self.__shl[idx1:idx2] = self.__sh[idx1:idx2] / sl[idx1:idx2]
self.__sxh[idx1:idx2] = 1.
self.__syh[idx1:idx2] = 0.
self.__sx[idx1:idx2] = s[idx1:idx2, 0]
self.__sy[idx1:idx2] = s[idx1:idx2, 1]
mask = self.__sh != 0.
self.__sxh[mask] = self.__sx[mask] / self.__sh[mask]
self.__syh[mask] = self.__sy[mask] / self.__sh[mask]
def __moc_calculate_characteristic_variables(self, s: ndarray, sl: ndarray, idx1: int, idx2: int) -> None:
# compute eigenvalues (convection + two pressure waves)
self.__ev1[idx1:idx2] = (
self.__u_a[idx1:idx2] * s[idx1:idx2, 0] +
self.__v_a[idx1:idx2] * s[idx1:idx2, 1] +
self.__w_a[idx1:idx2] * s[idx1:idx2, 2])
self.__ev4[idx1:idx2] = self.__ev1[idx1:idx2] + (self.__c_a[idx1:idx2] * sl[idx1:idx2])
self.__ev5[idx1:idx2] = self.__ev1[idx1:idx2] - (self.__c_a[idx1:idx2] * sl[idx1:idx2])
# compute characteristic variables
self.__c2[idx1:idx2] = np.power(self.__c_a[idx1:idx2], 2)
self.__tmp_1[idx1:idx2] = 0.5 * self.__r_a[idx1:idx2] / self.__c_a[idx1:idx2] * (
self.__u_r[idx1:idx2] * self.__sxl[idx1:idx2] +
self.__v_r[idx1:idx2] * self.__syl[idx1:idx2] +
self.__w_r[idx1:idx2] * self.__szl[idx1:idx2])
self.__tmp_2[idx1:idx2] = 0.5 * self.__r_a[idx1:idx2] / self.__c_a[idx1:idx2] * (
self.__u_l[idx1:idx2] * self.__sxl[idx1:idx2] +
self.__v_l[idx1:idx2] * self.__syl[idx1:idx2] +
self.__w_l[idx1:idx2] * self.__szl[idx1:idx2])
self.__r1[idx1:idx2] = ((self.__r_r[idx1:idx2] - self.__p_r[idx1:idx2] / self.__c2[idx1:idx2]) *
(self.__ev1[idx1:idx2] < 0) +
(self.__r_l[idx1:idx2] - self.__p_l[idx1:idx2] / self.__c2[idx1:idx2]) *
(self.__ev1[idx1:idx2] >= 0))
self.__r2[idx1:idx2] = ((self.__sxh[idx1:idx2] * self.__v_r[idx1:idx2] -
self.__syh[idx1:idx2] * self.__u_r[idx1:idx2]) *
(self.__ev1[idx1:idx2] < 0) +
(self.__sxh[idx1:idx2] * self.__v_l[idx1:idx2] -
self.__syh[idx1:idx2] * self.__u_l[idx1:idx2]) *
(self.__ev1[idx1:idx2] >= 0))
self.__r3[idx1:idx2] = ((self.__shl[idx1:idx2] * self.__w_r[idx1:idx2] - self.__szl[idx1:idx2] *
(self.__sxh[idx1:idx2] * self.__u_r[idx1:idx2] + self.__syh[idx1:idx2] *
self.__v_r[idx1:idx2])) * (self.__ev1[idx1:idx2] < 0) +
(self.__shl[idx1:idx2] * self.__w_l[idx1:idx2] - self.__szl[idx1:idx2] *
(self.__sxh[idx1:idx2] * self.__u_l[idx1:idx2] + self.__syh[idx1:idx2] *
self.__v_l[idx1:idx2])) * (self.__ev1[idx1:idx2] >= 0))
self.__r4[idx1:idx2] = ((0.5 * self.__p_r[idx1:idx2] / self.__c2[idx1:idx2] + self.__tmp_1[idx1:idx2]) *
(self.__ev4[idx1:idx2] < 0) +
(0.5 * self.__p_l[idx1:idx2] / self.__c2[idx1:idx2] + self.__tmp_2[idx1:idx2]) *
(self.__ev4[idx1:idx2] >= 0))
self.__r5[idx1:idx2] = ((0.5 * self.__p_r[idx1:idx2] / self.__c2[idx1:idx2] - self.__tmp_1[idx1:idx2]) *
(self.__ev5[idx1:idx2] < 0) +
(0.5 * self.__p_l[idx1:idx2] / self.__c2[idx1:idx2] - self.__tmp_2[idx1:idx2]) *
(self.__ev5[idx1:idx2] >= 0))
def __moc_convert_back_to_primitive(self, idx1: int, idx2: int) -> None:
# transform back to primitive variables
self.__tmp_1[idx1:idx2] = self.__r4[idx1:idx2] + self.__r5[idx1:idx2]
self.__tmp_2[idx1:idx2] = (self.__c_a[idx1:idx2] /
self.__r_a[idx1:idx2] * (self.__r4[idx1:idx2] - self.__r5[idx1:idx2]))
self.__ro_m[idx1:idx2] = self.__r1[idx1:idx2] + self.__tmp_1[idx1:idx2]
self.__u_m[idx1:idx2] = (
self.__sxl[idx1:idx2] * self.__tmp_2[idx1:idx2] -
self.__sxh[idx1:idx2] * self.__szl[idx1:idx2] * self.__r3[idx1:idx2] -
self.__syh[idx1:idx2] * self.__r2[idx1:idx2])
self.__v_m[idx1:idx2] = (
self.__syl[idx1:idx2] * self.__tmp_2[idx1:idx2] -
self.__syh[idx1:idx2] * self.__szl[idx1:idx2] * self.__r3[idx1:idx2] +
self.__sxh[idx1:idx2] * self.__r2[idx1:idx2])
self.__w_m[idx1:idx2] = (
self.__szl[idx1:idx2] * self.__tmp_2[idx1:idx2] +
self.__shl[idx1:idx2] * self.__r3[idx1:idx2])
self.__p_m[idx1:idx2] = self.__c2[idx1:idx2] * self.__tmp_1[idx1:idx2]
self.__ru_m[idx1:idx2] = self.__ro_m[idx1:idx2] * self.__u_m[idx1:idx2]
self.__rv_m[idx1:idx2] = self.__ro_m[idx1:idx2] * self.__v_m[idx1:idx2]
self.__rw_m[idx1:idx2] = self.__ro_m[idx1:idx2] * self.__w_m[idx1:idx2]
self.__et_m[idx1:idx2] = (
self.__p_m[idx1:idx2] * self.__gas.one_over_gm1 +
0.5 * (self.__ru_m[idx1:idx2] * self.__u_m[idx1:idx2] +
self.__rv_m[idx1:idx2] * self.__v_m[idx1:idx2] +
self.__rw_m[idx1:idx2] * self.__w_m[idx1:idx2]))
def __moc_evaluate(self, s: ndarray, sl: ndarray, idx1: int, idx2: int) -> None:
self.__moc_calculate_characteristic_variables(s, sl, idx1, idx2)
self.__moc_convert_back_to_primitive(idx1, idx2)
def __tvd_ro_loop(self, ro: ndarray, kc: ndarray, idx1: int, idx2: int) -> None:
for idx in range(idx1, idx2):
if self.__ev1[idx] < 0:
div = ro[kc[idx, 1]] - ro[kc[idx, 3]]
if np.abs(div) < 1.e-10:
if div < 0.:
div = -1.e-10
else:
div = 1.e-10
alfa = (ro[kc[idx, 2]] - ro[kc[idx, 3]]) / div
beta = np.max([(self.__r_r[idx] - ro[kc[idx, 3]]) / div, alfa])
if 0. <= alfa < 1.0:
if alfa < 0.1:
beta = np.min([beta, 10. * alfa])
else:
beta = np.min([beta, 1.])
else:
beta = np.min([beta, alfa])
self.__r_r[idx] = ro[kc[idx, 3]] + beta * div
self.__r1[idx] = self.__r_r[idx] - self.__p_r[idx] / self.__c2[idx]
else:
div = ro[kc[idx, 2]] - ro[kc[idx, 0]]
if np.abs(div) < 1.e-10:
if div < 0.:
div = -1.e-10
else:
div = 1.e-10
alfa = (ro[kc[idx, 1]] - ro[kc[idx, 0]]) / div
beta = np.max([(self.__r_l[idx] - ro[kc[idx, 0]]) / div, alfa])
if 0. <= alfa < 1.0:
if alfa < 0.1:
beta = np.min([beta, 10. * alfa])
else:
beta = np.min([beta, 1.])
else:
beta = np.min([beta, alfa])
self.__r_l[idx] = ro[kc[idx, 0]] + beta * div
self.__r1[idx] = self.__r_l[idx] - self.__p_l[idx] / self.__c2[idx]
def __tvd_ro(self, ro: ndarray, kc: ndarray, idx1: int, idx2: int) -> None:
self.__tvd_div[idx1:idx2] = ((ro[kc[idx1:idx2, 1]] - ro[kc[idx1:idx2, 3]]) * (self.__ev1[idx1:idx2] < 0) +
(ro[kc[idx1:idx2, 2]] - ro[kc[idx1:idx2, 0]]) * (self.__ev1[idx1:idx2] >= 0))
self.__tvd_div[idx1:idx2] = np.where((self.__tvd_div[idx1:idx2] > -1.e-10) & (self.__tvd_div[idx1:idx2] < 0.),
-1.e-10, self.__tvd_div[idx1:idx2])
self.__tvd_div[idx1:idx2] = np.where((self.__tvd_div[idx1:idx2] < 1.e-10) & (self.__tvd_div[idx1:idx2] >= 0.),
1.e-10, self.__tvd_div[idx1:idx2])
self.__tvd_alfa[idx1:idx2] = ((ro[kc[idx1:idx2, 2]] - ro[kc[idx1:idx2, 3]]) /
self.__tvd_div[idx1:idx2] * (self.__ev1[idx1:idx2] < 0) +
(ro[kc[idx1:idx2, 1]] - ro[kc[idx1:idx2, 0]]) /
self.__tvd_div[idx1:idx2] * (self.__ev1[idx1:idx2] >= 0))
self.__tvd_beta[idx1:idx2] = ((self.__r_r[idx1:idx2] - ro[kc[idx1:idx2, 3]]) /
self.__tvd_div[idx1:idx2] * (self.__ev1[idx1:idx2] < 0) +
(self.__r_l[idx1:idx2] - ro[kc[idx1:idx2, 0]]) /
self.__tvd_div[idx1:idx2] * (self.__ev1[idx1:idx2] >= 0))
self.__tvd_beta[idx1:idx2] = np.maximum(self.__tvd_beta[idx1:idx2], self.__tvd_alfa[idx1:idx2])
self.__tvd_alfa[idx1:idx2] = np.where((self.__tvd_alfa[idx1:idx2] >= 0.) & (self.__tvd_alfa[idx1:idx2] < 0.1),
10. * self.__tvd_alfa[idx1:idx2], self.__tvd_alfa[idx1:idx2])
self.__tvd_alfa[idx1:idx2] = np.where((self.__tvd_alfa[idx1:idx2] >= 0.1) & (self.__tvd_alfa[idx1:idx2] < 1.0),
1.0, self.__tvd_alfa[idx1:idx2])
self.__tvd_beta[idx1:idx2] = np.minimum(self.__tvd_beta[idx1:idx2], self.__tvd_alfa[idx1:idx2])
self.__r_r[idx1:idx2] = ((ro[kc[idx1:idx2, 3]] + self.__tvd_beta[idx1:idx2] * self.__tvd_div[idx1:idx2]) *
(self.__ev1[idx1:idx2] < 0) + self.__r_r[idx1:idx2] * (self.__ev1[idx1:idx2] >= 0))
self.__r_l[idx1:idx2] = (self.__r_l[idx1:idx2] * (self.__ev1[idx1:idx2] < 0) +
(ro[kc[idx1:idx2, 0]] + self.__tvd_beta[idx1:idx2] * self.__tvd_div[idx1:idx2]) *
(self.__ev1[idx1:idx2] >= 0))
self.__r1[idx1:idx2] = ((self.__r_r[idx1:idx2] - self.__p_r[idx1:idx2] / self.__c2[idx1:idx2]) *
(self.__ev1[idx1:idx2] < 0) +
(self.__r_l[idx1:idx2] - self.__p_l[idx1:idx2] / self.__c2[idx1:idx2]) *
(self.__ev1[idx1:idx2] >= 0))
def __conv_flux(self, fmp: FMP, n_faces: int = 0, StartIndex: int = 0) -> None:
# calculate face convective flux contribution
g3d_scheme = G3DScheme()
if n_faces == 0:
n_faces = fmp.size - StartIndex
n_faces, idx1, idx2 = np.int32([n_faces, StartIndex, StartIndex + n_faces])
ro, ru, rv, rw, et = self.__domain.get_cons_arrays_with_ghosts()
p, u, v, w, t, c = self.__domain.get_auxv_arrays_with_ghosts()
kp, kc = fmp.get_pointers()
s, sl = fmp.get_conv_metrics()
coefs = g3d_scheme.get_coefs()
# compute left- and right-biased primitive variables
self.__r_l[idx1:idx2] = (coefs[0] * ro[kc[idx1:idx2, 0]] + coefs[1] * ro[kc[idx1:idx2, 1]] +
coefs[2] * ro[kc[idx1:idx2, 2]] + coefs[3] * ro[kc[idx1:idx2, 3]])
self.__u_l[idx1:idx2] = (coefs[0] * u[kc[idx1:idx2, 0]] + coefs[1] * u[kc[idx1:idx2, 1]] +
coefs[2] * u[kc[idx1:idx2, 2]] + coefs[3] * u[kc[idx1:idx2, 3]])
self.__v_l[idx1:idx2] = (coefs[0] * v[kc[idx1:idx2, 0]] + coefs[1] * v[kc[idx1:idx2, 1]] +
coefs[2] * v[kc[idx1:idx2, 2]] + coefs[3] * v[kc[idx1:idx2, 3]])
self.__w_l[idx1:idx2] = (coefs[0] * w[kc[idx1:idx2, 0]] + coefs[1] * w[kc[idx1:idx2, 1]] +
coefs[2] * w[kc[idx1:idx2, 2]] + coefs[3] * w[kc[idx1:idx2, 3]])
self.__p_l[idx1:idx2] = (coefs[0] * p[kc[idx1:idx2, 0]] + coefs[1] * p[kc[idx1:idx2, 1]] +
coefs[2] * p[kc[idx1:idx2, 2]] + coefs[3] * p[kc[idx1:idx2, 3]])
self.__r_r[idx1:idx2] = (coefs[3] * ro[kc[idx1:idx2, 0]] + coefs[2] * ro[kc[idx1:idx2, 1]] +
coefs[1] * ro[kc[idx1:idx2, 2]] + coefs[0] * ro[kc[idx1:idx2, 3]])
self.__u_r[idx1:idx2] = (coefs[3] * u[kc[idx1:idx2, 0]] + coefs[2] * u[kc[idx1:idx2, 1]] +
coefs[1] * u[kc[idx1:idx2, 2]] + coefs[0] * u[kc[idx1:idx2, 3]])
self.__v_r[idx1:idx2] = (coefs[3] * v[kc[idx1:idx2, 0]] + coefs[2] * v[kc[idx1:idx2, 1]] +
coefs[1] * v[kc[idx1:idx2, 2]] + coefs[0] * v[kc[idx1:idx2, 3]])
self.__w_r[idx1:idx2] = (coefs[3] * w[kc[idx1:idx2, 0]] + coefs[2] * w[kc[idx1:idx2, 1]] +
coefs[1] * w[kc[idx1:idx2, 2]] + coefs[0] * w[kc[idx1:idx2, 3]])
self.__p_r[idx1:idx2] = (coefs[3] * p[kc[idx1:idx2, 0]] + coefs[2] * p[kc[idx1:idx2, 1]] +
coefs[1] * p[kc[idx1:idx2, 2]] + coefs[0] * p[kc[idx1:idx2, 3]])
# compute averaged quantities for local linearization
self.__r_a[idx1:idx2] = 0.5 * (ro[kc[idx1:idx2, 1]] + ro[kc[idx1:idx2, 2]])
self.__u_a[idx1:idx2] = 0.5 * (u[kc[idx1:idx2, 1]] + u[kc[idx1:idx2, 2]])
self.__v_a[idx1:idx2] = 0.5 * (v[kc[idx1:idx2, 1]] + v[kc[idx1:idx2, 2]])
self.__w_a[idx1:idx2] = 0.5 * (w[kc[idx1:idx2, 1]] + w[kc[idx1:idx2, 2]])
self.__c_a[idx1:idx2] = 0.5 * (c[kc[idx1:idx2, 1]] + c[kc[idx1:idx2, 2]])
# evaluate characteristic variables using MOC
self.__moc_update_metrics(s, sl, idx1, idx2)
self.__moc_calculate_characteristic_variables(s, sl, idx1, idx2)
if g3d_scheme.use_density_tvd_limiter:
self.__tvd_ro(ro, kc, idx1, idx2)
self.__moc_convert_back_to_primitive(idx1, idx2)
# self.__moc_evaluate(s, sl, idx1, idx2)
# add extra artificial damping
self.__c_filt[idx1:idx2] = 0.
if g3d_scheme.shock_damp_coef > 0.:
self.__tmp_1[idx1:idx2] = (
np.abs(p[kc[idx1:idx2, 0]] - 2. * p[kc[idx1:idx2, 1]] + p[kc[idx1:idx2, 2]]) /
np.abs(p[kc[idx1:idx2, 0]] + 2. * p[kc[idx1:idx2, 1]] + p[kc[idx1:idx2, 2]]))
self.__tmp_2[idx1:idx2] = (
np.abs(p[kc[idx1:idx2, 1]] - 2. * p[kc[idx1:idx2, 2]] + p[kc[idx1:idx2, 3]]) /
np.abs(p[kc[idx1:idx2, 1]] + 2. * p[kc[idx1:idx2, 2]] + p[kc[idx1:idx2, 3]]))
if fmp.type != FMPTypes.BOUNDARY and g3d_scheme.shock_damping_mode.value > ShockDampingMode.JAMESON.value:
if g3d_scheme.shock_damping_mode is ShockDampingMode.DUCROSS_FACE_GRADIENTS:
# evaluate Ducross sensor weight factors using face-based velocity gradients
self.__dudx[idx1:idx2], self.__dudy[idx1:idx2], self.__dudz[idx1:idx2] = (
self.__calculate_face_gradients(u, fmp, idx1, idx2))
self.__dvdx[idx1:idx2], self.__dvdy[idx1:idx2], self.__dvdz[idx1:idx2] = (
self.__calculate_face_gradients(v, fmp, idx1, idx2))
self.__dwdx[idx1:idx2], self.__dwdy[idx1:idx2], self.__dwdz[idx1:idx2] = (
self.__calculate_face_gradients(w, fmp, idx1, idx2))
# divergence
self.__tmp_3[idx1:idx2] = (self.__dudx[idx1:idx2] +
self.__dvdy[idx1:idx2] +
self.__dwdz[idx1:idx2]) ** 2
# cross product
self.__tmp_4[idx1:idx2] = (((self.__dwdy[idx1:idx2] - self.__dvdz[idx1:idx2]) ** 2 +
(self.__dudz[idx1:idx2] - self.__dwdx[idx1:idx2]) ** 2 +
(self.__dvdx[idx1:idx2] - self.__dudy[idx1:idx2]) ** 2) ** 0.5 +
1.e-20 * np.ones(len(self.__dudx[idx1:idx2])))
self.__tmp_1[idx1:idx2] = (self.__tmp_1[idx1:idx2] * self.__tmp_3[idx1:idx2] /
(self.__tmp_3[idx1:idx2] + self.__tmp_4[idx1:idx2]))
self.__tmp_2[idx1:idx2] = (self.__tmp_2[idx1:idx2] * self.__tmp_3[idx1:idx2] /
(self.__tmp_3[idx1:idx2] + self.__tmp_4[idx1:idx2]))
elif g3d_scheme.shock_damping_mode is ShockDampingMode.DUCROSS_CELL_GRADIENTS:
# evaluate Ducross sensor weight factors using cell-based velocity gradients
dudx, dudy, dudz, dvdx, dvdy, dvdz, dwdx, dwdy, dwdz = self.__domain.get_vgrad_arrays()
# divergence 1
self.__tmp_3[idx1:idx2] = (dudx[kc[idx1:idx2, 1]] +
dvdy[kc[idx1:idx2, 1]] +
dwdz[kc[idx1:idx2, 1]]) ** 2
# cross product 1
self.__tmp_4[idx1:idx2] = (((dwdy[kc[idx1:idx2, 1]] - dvdz[kc[idx1:idx2, 1]]) ** 2 +
(dudz[kc[idx1:idx2, 1]] - dwdx[kc[idx1:idx2, 1]]) ** 2 +
(dvdx[kc[idx1:idx2, 1]] - dudy[kc[idx1:idx2, 1]]) ** 2) ** 0.5 +
1.e-20 * np.ones(len(self.__c_filt[idx1:idx2])))
self.__tmp_1[idx1:idx2] = (self.__tmp_1[idx1:idx2] * self.__tmp_3[idx1:idx2] /
(self.__tmp_3[idx1:idx2] + self.__tmp_4[idx1:idx2]))
# divergence 2
self.__tmp_3[idx1:idx2] = (dudx[kc[idx1:idx2, 2]] +
dvdy[kc[idx1:idx2, 2]] +
dwdz[kc[idx1:idx2, 2]]) ** 2
# cross produce 2
self.__tmp_4[idx1:idx2] = (((dwdy[kc[idx1:idx2, 2]] - dvdz[kc[idx1:idx2, 2]]) ** 2 +
(dudz[kc[idx1:idx2, 2]] - dwdx[kc[idx1:idx2, 2]]) ** 2 +
(dvdx[kc[idx1:idx2, 2]] - dudy[kc[idx1:idx2, 2]]) ** 2) ** 0.5 +
1.e-20 * np.ones(len(self.__c_filt[idx1:idx2])))
self.__tmp_2[idx1:idx2] = (self.__tmp_2[idx1:idx2] * self.__tmp_3[idx1:idx2] /
(self.__tmp_3[idx1:idx2] + self.__tmp_4[idx1:idx2]))
self.__c_filt[idx1:idx2] = (g3d_scheme.shock_damp_coef *
np.maximum(self.__tmp_1[idx1:idx2], self.__tmp_2[idx1:idx2]) * sl[idx1:idx2] *
(np.abs(self.__u_a[idx1:idx2]) + self.__c_a[idx1:idx2]))
if g3d_scheme.entropy_damp_coef > 0.:
self.__tmp_1[idx1:idx2] = (
np.abs(ro[kc[idx1:idx2, 0]] - 2. * ro[kc[idx1:idx2, 1]] + ro[kc[idx1:idx2, 2]]) /
np.abs(ro[kc[idx1:idx2, 0]] + 2. * ro[kc[idx1:idx2, 1]] + ro[kc[idx1:idx2, 2]]))
self.__tmp_2[idx1:idx2] = (
np.abs(ro[kc[idx1:idx2, 1]] - 2. * ro[kc[idx1:idx2, 2]] + ro[kc[idx1:idx2, 3]]) /
np.abs(ro[kc[idx1:idx2, 1]] + 2. * ro[kc[idx1:idx2, 2]] + ro[kc[idx1:idx2, 3]]))
self.__c_filt[idx1:idx2] = np.maximum(self.__c_filt[idx1:idx2], g3d_scheme.entropy_damp_coef * np.maximum(
self.__tmp_1[idx1:idx2], self.__tmp_2[idx1:idx2]) * sl[idx1:idx2] * np.abs(self.__u_a[idx1:idx2]))
# compute mass, momentum and energy fluxes
self.__tmp_1[idx1:idx2] = (self.__ru_m[idx1:idx2] * s[idx1:idx2, 0] +
self.__rv_m[idx1:idx2] * s[idx1:idx2, 1] +
self.__rw_m[idx1:idx2] * s[idx1:idx2, 2])
self.__f_ro[idx1:idx2] = (self.__tmp_1[idx1:idx2] -
self.__c_filt[idx1:idx2] * (ro[kc[idx1:idx2, 2]] - ro[kc[idx1:idx2, 1]]))
self.__f_ru[idx1:idx2] = (self.__tmp_1[idx1:idx2] * self.__u_m[idx1:idx2] +
self.__p_m[idx1:idx2] * s[idx1:idx2, 0] -
self.__c_filt[idx1:idx2] * (ru[kc[idx1:idx2, 2]] - ru[kc[idx1:idx2, 1]]))
self.__f_rv[idx1:idx2] = (self.__tmp_1[idx1:idx2] * self.__v_m[idx1:idx2] +
self.__p_m[idx1:idx2] * s[idx1:idx2, 1] -
self.__c_filt[idx1:idx2] * (rv[kc[idx1:idx2, 2]] - rv[kc[idx1:idx2, 1]]))
self.__f_rw[idx1:idx2] = (self.__tmp_1[idx1:idx2] * self.__w_m[idx1:idx2] +
self.__p_m[idx1:idx2] * s[idx1:idx2, 2] -
self.__c_filt[idx1:idx2] * (rw[kc[idx1:idx2, 2]] - rw[kc[idx1:idx2, 1]]))
self.__f_et[idx1:idx2] = ((self.__tmp_1[idx1:idx2] * (self.__et_m[idx1:idx2] + self.__p_m[idx1:idx2]) /
self.__ro_m[idx1:idx2]) - self.__c_filt[idx1:idx2] *
(et[kc[idx1:idx2, 2]] - et[kc[idx1:idx2, 1]]))
# add calculated convective fluxes
self._ro[kc[idx1:idx2, 1]] -= self.__f_ro[idx1:idx2]
self._ru[kc[idx1:idx2, 1]] -= self.__f_ru[idx1:idx2]
self._rv[kc[idx1:idx2, 1]] -= self.__f_rv[idx1:idx2]
self._rw[kc[idx1:idx2, 1]] -= self.__f_rw[idx1:idx2]
self._et[kc[idx1:idx2, 1]] -= self.__f_et[idx1:idx2]
self._ro[kc[idx1:idx2, 2]] += self.__f_ro[idx1:idx2]
self._ru[kc[idx1:idx2, 2]] += self.__f_ru[idx1:idx2]
self._rv[kc[idx1:idx2, 2]] += self.__f_rv[idx1:idx2]
self._rw[kc[idx1:idx2, 2]] += self.__f_rw[idx1:idx2]
self._et[kc[idx1:idx2, 2]] += self.__f_et[idx1:idx2]
def __calculate_shear_stress_components(self, idx1: int, idx2: int) -> None:
visc = self.__gas.viscosity
self.__tau_xx[idx1:idx2] = visc * (4. * self.__dudx[idx1:idx2] - 2. *
(self.__dvdy[idx1:idx2] + self.__dwdz[idx1:idx2])) / 3.
self.__tau_yy[idx1:idx2] = visc * (4. * self.__dvdy[idx1:idx2] - 2. *
(self.__dudx[idx1:idx2] + self.__dwdz[idx1:idx2])) / 3.
self.__tau_zz[idx1:idx2] = visc * (4. * self.__dwdz[idx1:idx2] - 2. *
(self.__dudx[idx1:idx2] + self.__dvdy[idx1:idx2])) / 3.
self.__tau_xy[idx1:idx2] = visc * (self.__dudy[idx1:idx2] + self.__dvdx[idx1:idx2])
self.__tau_xz[idx1:idx2] = visc * (self.__dudz[idx1:idx2] + self.__dwdx[idx1:idx2])
self.__tau_yz[idx1:idx2] = visc * (self.__dvdz[idx1:idx2] + self.__dwdy[idx1:idx2])
def __diff_flux(self, fmp: FMP, n_faces: int = 0, StartIndex: int = 0) -> None:
# calculate internal face diffusive flux contribution
if n_faces == 0:
n_faces = fmp.size - StartIndex
n_faces, idx1, idx2 = np.int32([n_faces, StartIndex, StartIndex + n_faces])
kp, kc = fmp.get_pointers()
p, u, v, w, t, c = self.__domain.get_auxv_arrays_with_ghosts()
sx, sy, sz, si, voli = fmp.get_diff_metrics()
# compute velocity gradients at face center
self.__dudx[idx1:idx2], self.__dudy[idx1:idx2], self.__dudz[idx1:idx2] = (
self.__calculate_face_gradients(u, fmp, idx1, idx2))
self.__dvdx[idx1:idx2], self.__dvdy[idx1:idx2], self.__dvdz[idx1:idx2] = (
self.__calculate_face_gradients(v, fmp, idx1, idx2))
self.__dwdx[idx1:idx2], self.__dwdy[idx1:idx2], self.__dwdz[idx1:idx2] = (
self.__calculate_face_gradients(w, fmp, idx1, idx2))
self.__dtdx[idx1:idx2], self.__dtdy[idx1:idx2], self.__dtdz[idx1:idx2] = (
self.__calculate_face_gradients(t, fmp, idx1, idx2))
# calculate face velocities
self.__u_a[idx1:idx2] = 0.5 * (u[kc[idx1:idx2, 2]] + u[kc[idx1:idx2, 1]])
self.__v_a[idx1:idx2] = 0.5 * (v[kc[idx1:idx2, 2]] + v[kc[idx1:idx2, 1]])
self.__w_a[idx1:idx2] = 0.5 * (w[kc[idx1:idx2, 2]] + w[kc[idx1:idx2, 1]])
# calculate viscous stress tensor components
self.__calculate_shear_stress_components(idx1, idx2)
# calculate diffusive flux contribution
self.__tmp_1[idx1:idx2] = self.__gas.Cp * self.__gas.viscosity / self.__gas.Pr * (
self.__dtdx[idx1:idx2] * si[idx1:idx2, 0] +
self.__dtdy[idx1:idx2] * si[idx1:idx2, 1] +
self.__dtdz[idx1:idx2] * si[idx1:idx2, 2])
self.__f_ru[idx1:idx2] = -(
self.__tau_xx[idx1:idx2] * si[idx1:idx2, 0] +
self.__tau_xy[idx1:idx2] * si[idx1:idx2, 1] +
self.__tau_xz[idx1:idx2] * si[idx1:idx2, 2])
self.__f_rv[idx1:idx2] = -(
self.__tau_xy[idx1:idx2] * si[idx1:idx2, 0] +
self.__tau_yy[idx1:idx2] * si[idx1:idx2, 1] +
self.__tau_yz[idx1:idx2] * si[idx1:idx2, 2])
self.__f_rw[idx1:idx2] = -(
self.__tau_xz[idx1:idx2] * si[idx1:idx2, 0] +
self.__tau_yz[idx1:idx2] * si[idx1:idx2, 1] +
self.__tau_zz[idx1:idx2] * si[idx1:idx2, 2])
self.__f_et[idx1:idx2] = (
self.__f_ru[idx1:idx2] * self.__u_a[idx1:idx2] +
self.__f_rv[idx1:idx2] * self.__v_a[idx1:idx2] +
self.__f_rw[idx1:idx2] * self.__w_a[idx1:idx2] + self.__tmp_1[idx1:idx2])
# add viscous fluxes
self._ru[kc[idx1:idx2, 1]] -= self.__f_ru[idx1:idx2]
self._rv[kc[idx1:idx2, 1]] -= self.__f_rv[idx1:idx2]
self._rw[kc[idx1:idx2, 1]] -= self.__f_rw[idx1:idx2]
self._et[kc[idx1:idx2, 1]] -= self.__f_et[idx1:idx2]
self._ru[kc[idx1:idx2, 2]] += self.__f_ru[idx1:idx2]
self._rv[kc[idx1:idx2, 2]] += self.__f_rv[idx1:idx2]
self._rw[kc[idx1:idx2, 2]] += self.__f_rw[idx1:idx2]
self._et[kc[idx1:idx2, 2]] += self.__f_et[idx1:idx2]
def __wall_diff_flux_resolved(self, boco: BoCo) -> None:
# calculate wall diffusive flux contribution (fully resolved boundary layer)
fmp = boco.FMP
wall_velocity = boco.get_velocity()
idx1, idx2 = np.int32([boco.get_face_count(), fmp.size])
u, v, w = self.__domain.get_auxv_velocity_arrays_with_ghosts()
sx, sy, sz, si, voli = fmp.get_diff_metrics()
kp, kc = fmp.get_pointers()
# compute velocity gradients at face center
self.__q_diff[idx1:idx2, 0] = 2. * (wall_velocity[0] - u[kc[idx1:idx2, 1]])
self.__q_diff[idx1:idx2, 1] = 2. * (wall_velocity[1] - v[kc[idx1:idx2, 1]])
self.__q_diff[idx1:idx2, 2] = 2. * (wall_velocity[2] - w[kc[idx1:idx2, 1]])
self.__dudx[idx1:idx2] = self.__q_diff[idx1:idx2, 0] * sx[idx1:idx2, 0] * voli[idx1:idx2]
self.__dudy[idx1:idx2] = self.__q_diff[idx1:idx2, 0] * sy[idx1:idx2, 0] * voli[idx1:idx2]
self.__dudz[idx1:idx2] = self.__q_diff[idx1:idx2, 0] * sz[idx1:idx2, 0] * voli[idx1:idx2]
self.__dvdx[idx1:idx2] = self.__q_diff[idx1:idx2, 1] * sx[idx1:idx2, 0] * voli[idx1:idx2]
self.__dvdy[idx1:idx2] = self.__q_diff[idx1:idx2, 1] * sy[idx1:idx2, 0] * voli[idx1:idx2]
self.__dvdz[idx1:idx2] = self.__q_diff[idx1:idx2, 1] * sz[idx1:idx2, 0] * voli[idx1:idx2]
self.__dwdx[idx1:idx2] = self.__q_diff[idx1:idx2, 2] * sx[idx1:idx2, 0] * voli[idx1:idx2]
self.__dwdy[idx1:idx2] = self.__q_diff[idx1:idx2, 2] * sy[idx1:idx2, 0] * voli[idx1:idx2]
self.__dwdz[idx1:idx2] = self.__q_diff[idx1:idx2, 2] * sz[idx1:idx2, 0] * voli[idx1:idx2]
# calculate viscous stress tensor components
self.__calculate_shear_stress_components(idx1, idx2)
# calculate diffusive flux contribution
self.__f_ru[idx1:idx2] = -(
self.__tau_xx[idx1:idx2] * si[idx1:idx2, 0] +
self.__tau_xy[idx1:idx2] * si[idx1:idx2, 1] +
self.__tau_xz[idx1:idx2] * si[idx1:idx2, 2])
self.__f_rv[idx1:idx2] = -(
self.__tau_xy[idx1:idx2] * si[idx1:idx2, 0] +
self.__tau_yy[idx1:idx2] * si[idx1:idx2, 1] +
self.__tau_yz[idx1:idx2] * si[idx1:idx2, 2])
self.__f_rw[idx1:idx2] = -(
self.__tau_xz[idx1:idx2] * si[idx1:idx2, 0] +
self.__tau_yz[idx1:idx2] * si[idx1:idx2, 1] +
self.__tau_zz[idx1:idx2] * si[idx1:idx2, 2])
self.__f_et[idx1:idx2] = (
wall_velocity[0] * self.__f_ru[idx1:idx2] +
wall_velocity[1] * self.__f_rv[idx1:idx2] +
wall_velocity[2] * self.__f_rw[idx1:idx2])
# add viscous fluxes
self._ru[kc[idx1:idx2, 1]] -= self.__f_ru[idx1:idx2]
self._rv[kc[idx1:idx2, 1]] -= self.__f_rv[idx1:idx2]
self._rw[kc[idx1:idx2, 1]] -= self.__f_rw[idx1:idx2]
self._et[kc[idx1:idx2, 1]] -= self.__f_et[idx1:idx2]
def __wall_diff_flux_wall_function(self, boco: BoCo) -> None:
# calculate wall diffusive flux contribution (wall function approximation)
fmp = boco.FMP
wall_velocity = boco.get_velocity()
ro = self.__domain.get_cons_array_with_ghosts('ro')
u, v, w = self.__domain.get_auxv_velocity_arrays_with_ghosts()
kp, kc = fmp.get_pointers()
s, sl = fmp.get_conv_metrics()
idx1, idx2 = np.int32([boco.get_face_count(), fmp.size])
y_norm = fmp.get_wall_normal()
u_tang = fmp.get_tmp_storage(0)
upr = fmp.get_tmp_storage(1)
vpr = fmp.get_tmp_storage(2)
wpr = fmp.get_tmp_storage(3)
tauw = fmp.get_tmp_storage(4)
self.__s_vel[idx1:idx2] = (u[kc[idx1:idx2, 1]] * s[idx1:idx2, 0] +
v[kc[idx1:idx2, 1]] * s[idx1:idx2, 1] +
w[kc[idx1:idx2, 1]] * s[idx1:idx2, 2]) / np.power(sl[idx1:idx2], 2)
upr[idx1:idx2] = u[kc[idx1:idx2, 1]] - (self.__s_vel[idx1:idx2] * s[idx1:idx2, 0] - wall_velocity[0])
vpr[idx1:idx2] = v[kc[idx1:idx2, 1]] - (self.__s_vel[idx1:idx2] * s[idx1:idx2, 1] - wall_velocity[1])
wpr[idx1:idx2] = w[kc[idx1:idx2, 1]] - (self.__s_vel[idx1:idx2] * s[idx1:idx2, 2] - wall_velocity[2])
# tangential projected velocity
u_tang[idx1:idx2] = np.maximum(np.sqrt(
upr[idx1:idx2] * upr[idx1:idx2] +
vpr[idx1:idx2] * vpr[idx1:idx2] +
wpr[idx1:idx2] * wpr[idx1:idx2]), 1.e-20)
# estimate the wall-shear stress
tauw[idx1:idx2] = self.__log_law.evaluate(u_tang[idx1:idx2], ro[kc[idx1:idx2, 1]], y_norm[idx1:idx2])
# modify the projected tangential wall velocity
upr[idx1:idx2] *= tauw[idx1:idx2] * sl[idx1:idx2] / u_tang[idx1:idx2]
vpr[idx1:idx2] *= tauw[idx1:idx2] * sl[idx1:idx2] / u_tang[idx1:idx2]
wpr[idx1:idx2] *= tauw[idx1:idx2] * sl[idx1:idx2] / u_tang[idx1:idx2]
self._ru[kc[idx1:idx2, 1]] -= upr[idx1:idx2]
self._rv[kc[idx1:idx2, 1]] -= vpr[idx1:idx2]
self._rw[kc[idx1:idx2, 1]] -= wpr[idx1:idx2]
self._et[kc[idx1:idx2, 1]] -= (
upr[idx1:idx2] * wall_velocity[0] +
vpr[idx1:idx2] * wall_velocity[1] +
wpr[idx1:idx2] * wall_velocity[2])
def __symm_diff_flux(self, boco: BoCo) -> None:
# calculate diff metrics for symmetry molecules in solver setup
fmp = boco.FMP
idx1, idx2 = np.int32([fmp.size / 2, fmp.size])
p, u, v, w, t, c = self.__domain.get_auxv_arrays_with_ghosts()
kp, kc = fmp.get_pointers()
s, sl = fmp.get_conv_metrics()
p[kc[idx1:idx2, 2]] = p[kc[idx1:idx2, 1]]
t[kc[idx1:idx2, 2]] = t[kc[idx1:idx2, 1]]
c[kc[idx1:idx2, 2]] = c[kc[idx1:idx2, 1]]
self.__s_vel[idx1:idx2] = (u[kc[idx1:idx2, 1]] * s[idx1:idx2, 0] +
v[kc[idx1:idx2, 1]] * s[idx1:idx2, 1] +
w[kc[idx1:idx2, 1]] * s[idx1:idx2, 2]) / np.power(sl[idx1:idx2], 2)
u[kc[idx1:idx2, 2]] = u[kc[idx1:idx2, 1]] - 2. * self.__s_vel[idx1:idx2] * s[idx1:idx2, 0]
v[kc[idx1:idx2, 2]] = v[kc[idx1:idx2, 1]] - 2. * self.__s_vel[idx1:idx2] * s[idx1:idx2, 1]
w[kc[idx1:idx2, 2]] = w[kc[idx1:idx2, 1]] - 2. * self.__s_vel[idx1:idx2] * s[idx1:idx2, 2]
self.__diff_flux(fmp, int(fmp.size / 2), int(fmp.size / 2))
def __add_streamtube_pressure_forces(self, fmp: FMP, n_faces: int = 0, StartIndex: int = 0) -> None:
# add stream tube surface forces (applicable for Q3D simulations only)
if n_faces == 0:
n_faces = fmp.size - StartIndex
n_faces, idx1, idx2 = np.int32([n_faces, StartIndex, StartIndex + n_faces])
p = self.__domain.get_auxv_array_with_ghosts('p')
kp, kc = fmp.get_pointers()
s, sl = fmp.get_conv_metrics()
self._ru[kc[idx1:idx2, 1]] -= (p[kc[idx1:idx2, 1]] * s[idx1:idx2, 0])
self._rv[kc[idx1:idx2, 1]] -= (p[kc[idx1:idx2, 1]] * s[idx1:idx2, 1])
self._rw[kc[idx1:idx2, 1]] -= (p[kc[idx1:idx2, 1]] * s[idx1:idx2, 2])
[docs]
def flux_check(self) -> None:
"""
"""
NFlux = np.zeros(self.size)
internal = self.__domain.internal_faces
for tmp in internal.get_FMPs():
self.__flux_check(tmp, NFlux)
# connectivity face fluxes
for interface in self.__connectivity_faces.get_interface_array():
A, B = interface.get_face_defs()
for intf in A:
self.__flux_check(intf.FMP, NFlux)
# boundary face fluxes
for bc in self.__boundary_faces.get_boundary_array():
for boco in bc.get_face_defs():
self.__flux_check(boco.FMP, NFlux)
idx1, idx2 = self.get_range()
print(np.min(NFlux[idx1:idx2]), np.argmin(NFlux[idx1:idx2]))
print(np.max(NFlux[idx1:idx2]), np.argmax(NFlux[idx1:idx2]))
io = G3DIO()
io.export_scalar(self.__domain, NFlux, variable_name='NFlux', name='flux_check')
[docs]
def diff_flux_check(self) -> None:
"""
"""
NFlux = np.zeros(self.size)
internal = self.__domain.internal_faces
for tmp in internal.get_FMPs():
self.__flux_check(tmp, NFlux)
# connectivity face fluxes
for interface in self.__connectivity_faces.get_interface_array():
A, B = interface.get_face_defs()
for intf in A:
self.__flux_check(intf.FMP, NFlux)
# boundary face fluxes
for bc in self.__boundary_faces.get_boundary_array():
# note: singularities will give unbalanced output since the boundary face is not included in FMP
for boco in bc.get_face_defs():
tmp = boco.FMP
self.__flux_check(tmp, NFlux, int(tmp.size / 2))
# diffusive flux for face one step in from boundary
if boco.face_type is BCTypes.WALLS:
if boco.it is WALLTypes.VISCOUS_RESOLVED:
self.__flux_check(tmp, NFlux, int(tmp.size / 2), int(tmp.size / 2))
elif boco.it is WALLTypes.VISCOUS_WALL_FUNCTION:
self.__flux_check(tmp, NFlux, int(tmp.size / 2), int(tmp.size / 2))
elif boco.face_type is BCTypes.SYMM:
self.__flux_check(tmp, NFlux, int(tmp.size / 2), int(tmp.size / 2))
elif boco.face_type is BCTypes.IOFLOW:
# no diffusive fluxes for ioflow boundary face. added here to get a balanced output
self.__flux_check(tmp, NFlux, int(tmp.size / 2), int(tmp.size / 2))
idx1, idx2 = self.get_range()
print(np.min(NFlux[idx1:idx2]), np.argmin(NFlux[idx1:idx2]))
print(np.max(NFlux[idx1:idx2]), np.argmax(NFlux[idx1:idx2]))
io = G3DIO()
io.export_scalar(self.__domain, NFlux, variable_name='NFlux', name='diff_flux_check')
@staticmethod
def __flux_check(fmp: FMP, NF: ndarray, n_faces: int = 0, StartIndex: int = 0) -> None:
if n_faces == 0:
n_faces = fmp.size - StartIndex
idx1 = StartIndex
idx2 = StartIndex + n_faces
kp, kc = fmp.get_pointers()
s, sl = fmp.get_conv_metrics()
f_tmp = np.round((s[idx1:idx2, 0] + s[idx1:idx2, 1] + s[idx1:idx2, 2]) / sl[idx1:idx2])
NF[kc[idx1:idx2, 1]] = NF[kc[idx1:idx2, 1]] - f_tmp
NF[kc[idx1:idx2, 2]] = NF[kc[idx1:idx2, 2]] + f_tmp
[docs]
def get_shock_capturing_field(self) -> ndarray:
"""
:return:
"""
F_TMP = np.zeros(self._size)
# internal faces
internal = self.__domain.internal_faces
for tmp in internal.get_FMPs():
self.__evaluate_shock_capturing(tmp, F_TMP)
# connectivity faces
for interface in self.__connectivity_faces.get_interface_array():
A, B = interface.get_face_defs()
for intf in A:
tmp = intf.FMP
self.__evaluate_shock_capturing(tmp, F_TMP)
# boundary faces
for bc in self.__boundary_faces.get_boundary_array():
for boco in bc.get_face_defs():
tmp = boco.FMP
self.__update_boundary_ghost_cells(boco)
self.__evaluate_shock_capturing(tmp, F_TMP)
return F_TMP
def __evaluate_shock_capturing(self, fmp: FMP, F: ndarray, n_faces: int = 0, StartIndex: int = 0) -> None:
g3d_scheme = G3DScheme()
StartIndex = int(StartIndex)
n_faces = int(n_faces)
if n_faces == 0:
n_faces = fmp.size - StartIndex
idx1 = StartIndex
idx2 = int(StartIndex + n_faces)
ro, ru, rv, rw, et = self.__domain.get_cons_arrays_with_ghosts()
p, u, v, w, t, c = self.__domain.get_auxv_arrays_with_ghosts()
kp, kc = fmp.get_pointers()
s, sl = fmp.get_conv_metrics()
# compute averaged quantities (cell face estimates)
self.__r_a[idx1:idx2] = 0.5 * (ro[kc[idx1:idx2, 1]] + ro[kc[idx1:idx2, 2]])
self.__u_a[idx1:idx2] = 0.5 * (u[kc[idx1:idx2, 1]] + u[kc[idx1:idx2, 2]])
self.__v_a[idx1:idx2] = 0.5 * (v[kc[idx1:idx2, 1]] + v[kc[idx1:idx2, 2]])
self.__w_a[idx1:idx2] = 0.5 * (w[kc[idx1:idx2, 1]] + w[kc[idx1:idx2, 2]])
self.__c_a[idx1:idx2] = 0.5 * (c[kc[idx1:idx2, 1]] + c[kc[idx1:idx2, 2]])
# evaluate shock capturing parameters
self.__c_filt[idx1:idx2] = 0.
if g3d_scheme.shock_damp_coef > 0.:
self.__tmp_1[idx1:idx2] = (np.abs(p[kc[idx1:idx2, 0]] - 2. * p[kc[idx1:idx2, 1]] + p[kc[idx1:idx2, 2]]) /
np.abs(p[kc[idx1:idx2, 0]] + 2. * p[kc[idx1:idx2, 1]] + p[kc[idx1:idx2, 2]]))
self.__tmp_2[idx1:idx2] = (np.abs(p[kc[idx1:idx2, 1]] - 2. * p[kc[idx1:idx2, 2]] + p[kc[idx1:idx2, 3]]) /
np.abs(p[kc[idx1:idx2, 1]] + 2. * p[kc[idx1:idx2, 2]] + p[kc[idx1:idx2, 3]]))
if fmp.type != FMPTypes.BOUNDARY and g3d_scheme.shock_damping_mode.value > ShockDampingMode.JAMESON.value:
if g3d_scheme.shock_damping_mode is ShockDampingMode.DUCROSS_FACE_GRADIENTS:
# evaluate Ducross sensor weight factors using face-based velocity gradients
self.__dudx[idx1:idx2], self.__dudy[idx1:idx2], self.__dudz[idx1:idx2] = \
self.__calculate_face_gradients(u, fmp, idx1, idx2)
self.__dvdx[idx1:idx2], self.__dvdy[idx1:idx2], self.__dvdz[idx1:idx2] = \
self.__calculate_face_gradients(v, fmp, idx1, idx2)
self.__dwdx[idx1:idx2], self.__dwdy[idx1:idx2], self.__dwdz[idx1:idx2] = \
self.__calculate_face_gradients(w, fmp, idx1, idx2)
divergence = (self.__dudx[idx1:idx2] + self.__dvdy[idx1:idx2] + self.__dwdz[idx1:idx2]) ** 2
cross_product = (((self.__dwdy[idx1:idx2] - self.__dvdz[idx1:idx2]) ** 2 +
(self.__dudz[idx1:idx2] - self.__dwdx[idx1:idx2]) ** 2 +
(self.__dvdx[idx1:idx2] - self.__dudy[idx1:idx2]) ** 2) ** 0.5 +
1.e-20 * np.ones(len(self.__dudx[idx1:idx2])))
self.__tmp_1[idx1:idx2] = self.__tmp_1[idx1:idx2] * divergence / (divergence + cross_product)
self.__tmp_2[idx1:idx2] = self.__tmp_2[idx1:idx2] * divergence / (divergence + cross_product)
elif g3d_scheme.shock_damping_mode is ShockDampingMode.DUCROSS_CELL_GRADIENTS:
# evaluate Ducross sensor weight factors using cell-based velocity gradients
dudx, dudy, dudz, dvdx, dvdy, dvdz, dwdx, dwdy, dwdz = self.__domain.get_vgrad_arrays()
divergence_a = (dudx[kc[idx1:idx2, 1]] + dvdy[kc[idx1:idx2, 1]] + dwdz[kc[idx1:idx2, 1]]) ** 2
cross_product_a = (((dwdy[kc[idx1:idx2, 1]] - dvdz[kc[idx1:idx2, 1]]) ** 2 +
(dudz[kc[idx1:idx2, 1]] - dwdx[kc[idx1:idx2, 1]]) ** 2 +
(dvdx[kc[idx1:idx2, 1]] - dudy[kc[idx1:idx2, 1]]) ** 2) ** 0.5 +
1.e-20 * np.ones(len(self.__c_filt[idx1:idx2])))
divergence_b = (dudx[kc[idx1:idx2, 2]] + dvdy[kc[idx1:idx2, 2]] + dwdz[kc[idx1:idx2, 2]]) ** 2
cross_product_b = (((dwdy[kc[idx1:idx2, 2]] - dvdz[kc[idx1:idx2, 2]]) ** 2 +
(dudz[kc[idx1:idx2, 2]] - dwdx[kc[idx1:idx2, 2]]) ** 2 +
(dvdx[kc[idx1:idx2, 2]] - dudy[kc[idx1:idx2, 2]]) ** 2) ** 0.5 +
1.e-20 * np.ones(len(self.__c_filt[idx1:idx2])))
self.__tmp_1[idx1:idx2] = self.__tmp_1[idx1:idx2] * divergence_a / (divergence_a + cross_product_a)
self.__tmp_2[idx1:idx2] = self.__tmp_2[idx1:idx2] * divergence_b / (divergence_b + cross_product_b)
self.__c_filt[idx1:idx2] = (g3d_scheme.shock_damp_coef *
np.maximum(self.__tmp_1[idx1:idx2], self.__tmp_2[idx1:idx2]) * sl[idx1:idx2] * (
np.abs(self.__u_a[idx1:idx2]) + self.__c_a[idx1:idx2]))
F[kc[idx1:idx2, 1]] = np.maximum(F[kc[idx1:idx2, 1]], self.__c_filt[idx1:idx2])
F[kc[idx1:idx2, 2]] = np.maximum(F[kc[idx1:idx2, 2]], self.__c_filt[idx1:idx2])