from typing import Any
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 EQNMode
from g3d.defs import DIMMode
from g3d.defs import TIMESTEPTypes
from g3d.gas import Gas
from g3d.flow.domain import Domain
from g3d.flow.geo import Mesh
from g3d.flow.placeholders import ConsField, AuxvField, PostField, VGradField
from g3d.flow.solver.rhs import Flux
from g3d.defs import G3DEntity
[docs]
class G3DStep(G3DEntity):
def __init__(self, *args: Any, **kwargs: Any) -> None:
super().__init__(*args, **kwargs)
self._dt_min = 1.e20
self._dt_max = -1.e20
self._cfl_min = 1.e20
self._cfl_max = -1.e20
self._gas = Gas()
self._domain = Domain()
self._mesh = Mesh()
self._cons = ConsField()
self._auxv = AuxvField()
self._post = PostField()
self._vgrad = VGradField()
self._dcons = Flux()
self._dcons_pointers = self._dcons.get_pointers()
self.__sum3 = np.zeros(0)
self._allocate()
def _set_default(self) -> None:
super()._set_default()
self._tag = 'g3d step'
def _allocate(self) -> None:
g3d_world = G3DWorld()
n_cells = g3d_world.sizes.n_cells
n_ghost_cells = g3d_world.sizes.n_ghost_cells
self._dcons.setup(n_cells=n_cells, n_ghosts=n_ghost_cells)
self.__six = np.zeros(n_cells)
self.__siy = np.zeros(n_cells)
self.__siz = np.zeros(n_cells)
self.__sjx = np.zeros(n_cells)
self.__sjy = np.zeros(n_cells)
self.__sjz = np.zeros(n_cells)
self.__skx = np.zeros(n_cells)
self.__sky = np.zeros(n_cells)
self.__skz = np.zeros(n_cells)
self.__sum1 = np.zeros(n_cells)
self.__sum2 = np.zeros(n_cells)
if g3d_world.options.eqn_mode is not EQNMode.EULER:
self.__sum3 = np.zeros(n_cells)
self._s_max = np.zeros(n_cells)
self._tsf = np.zeros(n_cells)
self._dt = np.zeros(n_cells)
self._cfl_array = np.zeros(n_cells)
[docs]
def change_eqn_mode(self) -> None:
"""
"""
g3d_world = G3DWorld()
if (g3d_world.options.eqn_mode is not EQNMode.EULER) and len(self.__sum3) != g3d_world.sizes.n_cells:
self.__sum3 = np.zeros(g3d_world.sizes.n_cells)
self.__get_cell_metrics()
self._dcons.change_eqn_mode()
[docs]
def setup_and_allocate(self) -> None:
"""
"""
self._allocate()
[docs]
def associate(self, dom: Domain) -> None:
"""
:param dom:
"""
self._domain = dom
self._gas = dom.gas
self._mesh = dom.mesh
self._cons = self._domain.cons
self._auxv = self._domain.auxv
self._post = self._domain.post
self._vgrad = self._domain.vgrad
self._dcons.associate(dom, self._dt)
self.__get_cell_metrics()
self._dcons_pointers = self._dcons.get_pointers()
def __get_cell_metrics(self) -> None:
g3d_world = G3DWorld()
x, y, z = self._mesh.nodes.get_XYZ()
for im in range(self._mesh.n_blocks):
ni, nj, nk = self._mesh.blocks[im].get_block_sizes()
i = np.arange(ni - 1)
j = np.arange(nj - 1)
k = np.arange(nk - 1)
i, j, k = np.meshgrid(i, j, k, indexing='ij')
l0 = self._mesh.blocks[im].get_node(i, j, k)
ic = self._mesh.blocks[im].get_cell(i, j, k)
xp_0 = x[l0]
xp_1 = x[l0 + 1]
xp_2 = x[l0 + ni]
xp_3 = x[l0 + ni + 1]
xp_4 = x[l0 + ni * nj]
xp_5 = x[l0 + ni * nj + 1]
xp_6 = x[l0 + ni * nj + ni]
xp_7 = x[l0 + ni * nj + ni + 1]
yp_0 = y[l0]
yp_1 = y[l0 + 1]
yp_2 = y[l0 + ni]
yp_3 = y[l0 + ni + 1]
yp_4 = y[l0 + ni * nj]
yp_5 = y[l0 + ni * nj + 1]
yp_6 = y[l0 + ni * nj + ni]
yp_7 = y[l0 + ni * nj + ni + 1]
zp_0 = z[l0]
zp_1 = z[l0 + 1]
zp_2 = z[l0 + ni]
zp_3 = z[l0 + ni + 1]
zp_4 = z[l0 + ni * nj]
zp_5 = z[l0 + ni * nj + 1]
zp_6 = z[l0 + ni * nj + ni]
zp_7 = z[l0 + ni * nj + ni + 1]
if g3d_world.options.dim_mode is DIMMode.G2D_AXISYMM:
self.__six[ic] = 0.25 * ((yp_0 + yp_2) * (yp_2 - yp_0) + (yp_1 + yp_3) * (yp_3 - yp_1))
self.__siy[ic] = -0.25 * ((yp_0 + yp_2) * (xp_2 - xp_0) + (yp_1 + yp_3) * (xp_3 - xp_1))
self.__sjx[ic] = -0.25 * ((yp_0 + yp_1) * (yp_1 - yp_0) + (yp_2 + yp_3) * (yp_3 - yp_2))
self.__sjy[ic] = 0.25 * ((yp_0 + yp_1) * (xp_1 - xp_0) + (yp_2 + yp_3) * (xp_3 - xp_2))
else:
xi = .25 * ((xp_1 - xp_0) +
(xp_3 - xp_2) +
(xp_5 - xp_4) +
(xp_7 - xp_6))
yi = .25 * ((yp_1 - yp_0) +
(yp_3 - yp_2) +
(yp_5 - yp_4) +
(yp_7 - yp_6))
zi = .25 * ((zp_1 - zp_0) +
(zp_3 - zp_2) +
(zp_5 - zp_4) +
(zp_7 - zp_6))
xj = .25 * ((xp_2 - xp_0) +
(xp_3 - xp_1) +
(xp_6 - xp_4) +
(xp_7 - xp_5))
yj = .25 * ((yp_2 - yp_0) +
(yp_3 - yp_1) +
(yp_6 - yp_4) +
(yp_7 - yp_5))
zj = .25 * ((zp_2 - zp_0) +
(zp_3 - zp_1) +
(zp_6 - zp_4) +
(zp_7 - zp_5))
xk = .25 * ((xp_4 - xp_0) +
(xp_5 - xp_1) +
(xp_6 - xp_2) +
(xp_7 - xp_3))
yk = .25 * ((yp_4 - yp_0) +
(yp_5 - yp_1) +
(yp_6 - yp_2) +
(yp_7 - yp_3))
zk = .25 * ((zp_4 - zp_0) +
(zp_5 - zp_1) +
(zp_6 - zp_2) +
(zp_7 - zp_3))
self.__six[ic] = (yj * zk) - (zj * yk)
self.__siy[ic] = (zj * xk) - (xj * zk)
self.__siz[ic] = (xj * yk) - (yj * xk)
self.__sjx[ic] = (yk * zi) - (zk * yi)
self.__sjy[ic] = (zk * xi) - (xk * zi)
self.__sjz[ic] = (xk * yi) - (yk * xi)
self.__skx[ic] = (yi * zj) - (zi * yj)
self.__sky[ic] = (zi * xj) - (xi * zj)
self.__skz[ic] = (xi * yj) - (yi * xj)
# f4 = si dot si
# f5 = sj dot sj
# f6 = sk dot sk
# f7 = si dot sj
# f8 = si dot sk
# f9 = sj dot sk
f4 = (self.__six[ic] * self.__six[ic] +
self.__siy[ic] * self.__siy[ic] +
self.__siz[ic] * self.__siz[ic])
f5 = (self.__sjx[ic] * self.__sjx[ic] +
self.__sjy[ic] * self.__sjy[ic] +
self.__sjz[ic] * self.__sjz[ic])
f6 = (self.__skx[ic] * self.__skx[ic] +
self.__sky[ic] * self.__sky[ic] +
self.__skz[ic] * self.__skz[ic])
f7 = 2.0 * np.abs(self.__six[ic] * self.__sjx[ic] +
self.__siy[ic] * self.__sjy[ic] +
self.__siz[ic] * self.__sjz[ic])
f8 = 2.0 * np.abs(self.__six[ic] * self.__skx[ic] +
self.__siy[ic] * self.__sky[ic] +
self.__siz[ic] * self.__skz[ic])
f9 = 2.0 * np.abs(self.__sjx[ic] * self.__skx[ic] +
self.__sjy[ic] * self.__sky[ic] +
self.__sjz[ic] * self.__skz[ic])
# sum2 = f1 + f2 + f3 + f4 + f5 + f6 + f7 + f8 + f9
# sum3 = 3 * (f4 + f5 + f6) + sum2
self.__sum2[ic] = f4 + f5 + f6 + f7 + f8 + f9
if g3d_world.options.eqn_mode is not EQNMode.EULER:
self.__sum3[ic] = 3.0 * (f4 + f5 + f6) + self.__sum2[ic]
[docs]
def initialize(self) -> None:
"""
"""
g3d_world = G3DWorld()
if g3d_world.num.time.time_mode is TIMESTEPTypes.TIME_ACCURATE_CONSTANT_TIME_STEP:
self._dt[0:g3d_world.sizes.n_cells] = g3d_world.num.time.dtime
def _calculate_s_max(self) -> None:
g3d_world = G3DWorld()
p, u, v, w, t, c = self._domain.get_auxv_arrays()
ro, ru, rv, rw, et = self._domain.get_cons_arrays()
visct = self._domain.get_visct_array()
voli = self._mesh.get_inverse_volume()
idx1 = 0
idx2 = g3d_world.sizes.n_cells
# sum1 = (vel dot si) + (vel dot sj) + (vel dot sk)
self.__sum1[idx1:idx2] = (np.abs(u[idx1:idx2] * self.__six[idx1:idx2] +
v[idx1:idx2] * self.__siy[idx1:idx2] +
w[idx1:idx2] * self.__siz[idx1:idx2]) +
np.abs(u[idx1:idx2] * self.__sjx[idx1:idx2] +
v[idx1:idx2] * self.__sjy[idx1:idx2] +
w[idx1:idx2] * self.__sjz[idx1:idx2]) +
np.abs(u[idx1:idx2] * self.__skx[idx1:idx2] +
v[idx1:idx2] * self.__sky[idx1:idx2] +
w[idx1:idx2] * self.__skz[idx1:idx2]))
self._s_max[idx1:idx2] = (self.__sum1[idx1:idx2] +
c[idx1:idx2] * np.power(self.__sum2[idx1:idx2], 0.5)) * voli[idx1:idx2]
if g3d_world.options.eqn_mode is not EQNMode.EULER:
if g3d_world.options.eqn_mode.value > EQNMode.DNS.value:
self._s_max[idx1:idx2] = np.maximum(self._s_max[idx1:idx2],
2. * (self._gas.viscosity + visct[idx1:idx2]) *
self.__sum3[idx1:idx2] * np.power(voli[idx1:idx2], 2) /
ro[idx1:idx2])
else:
self._s_max[idx1:idx2] = np.maximum(self._s_max[idx1:idx2],
2. * self._gas.viscosity *
self.__sum3[idx1:idx2] * np.power(voli[idx1:idx2], 2) /
ro[idx1:idx2])
def _ltsp(self) -> None:
g3d_world = G3DWorld()
idx1 = 0
idx2 = g3d_world.sizes.n_cells
self._calculate_s_max()
self._tsf[idx1:idx2] = 1. / self._s_max[idx1:idx2]
if g3d_world.num.time.time_mode is TIMESTEPTypes.LOCAL_TIME_STEPPING:
# local time stepping
self._dt[idx1:idx2] = g3d_world.num.time.cfl * self._tsf[idx1:idx2]
if g3d_world.display.verbose:
self._dt_min = float(np.min(self._dt[idx1:idx2]))
self._dt_max = float(np.max(self._dt[idx1:idx2]))
print(' delta time (min) [max]: (%10.5E) [%10.5E]' % (self._dt_min, self._dt_max))
elif g3d_world.num.time.time_mode is TIMESTEPTypes.TIME_ACCURATE_CFL:
# time accurate (CFL-based)
tsf_min = float(np.min(self._tsf[idx1:idx2]))
dt_tmp = g3d_world.num.time.cfl * tsf_min
g3d_world.status.time = g3d_world.status.time + dt_tmp
self._dt[idx1:idx2] = dt_tmp
if g3d_world.display.verbose:
self._cfl_array[idx1:idx2] = self._dt[idx1:idx2] / self._tsf[idx1:idx2]
self._cfl_min = float(np.min(self._cfl_array[idx1:idx2]))
self._cfl_max = float(np.max(self._cfl_array[idx1:idx2]))
print(' CFL (min) [max]: (%10.5E) [%10.5E]' % (self._cfl_min, self._cfl_max))
elif g3d_world.num.time.time_mode is TIMESTEPTypes.TIME_ACCURATE_CONSTANT_TIME_STEP:
# time accurate (specified time step)
g3d_world.status.time = g3d_world.status.time + g3d_world.num.time.dtime
if g3d_world.display.verbose:
self._cfl_array[idx1:idx2] = self._dt[idx1:idx2] / self._tsf[idx1:idx2]
self._cfl_min = float(np.min(self._cfl_array[idx1:idx2]))
self._cfl_max = float(np.max(self._cfl_array[idx1:idx2]))
print(' CFL (min) [max]: (%10.5E) [%10.5E]' % (self._cfl_min, self._cfl_max))
[docs]
def fwd(self) -> None:
"""
"""
# base step class does not have an associated forward method
pass
[docs]
def flux_check(self) -> None:
"""
"""
self._dcons.flux_check()
[docs]
def diff_flux_check(self) -> None:
"""
"""
self._dcons.diff_flux_check()
[docs]
def flux_test(self) -> None: # test function
"""
"""
self._domain.status_check()
self._dcons.evaluate()
self._dcons.print_extreme_values()
[docs]
def ltsp_test(self) -> None: # test function
"""
"""
self._domain.status_check()
self._ltsp()
[docs]
def get_shock_capturing_field(self) -> ndarray:
"""
:return:
"""
return self._dcons.get_shock_capturing_field()