"""Domain class definition"""
import os
from typing import Any, Tuple, Union
from g3d.mesh.geo import IJKRegion
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 G3DIRSM
from g3d.defs import G3DSizes
from g3d.defs import BCTypes
from g3d.defs import GradientMode
from g3d.defs import INTFTypes
from g3d.defs import DIMMode
from g3d.defs import EQNMode
from g3d.defs import FaceType
from g3d.defs import clean_path
from g3d.defs import TurbulenceModelMode
from g3d.defs import G3DEntity
from g3d.flow.geo import Mesh
from g3d.flow.geo import Q1DPipe
from g3d.gas import Gas
from g3d.flow.vars import Cons
from g3d.flow.vars import Auxv
from g3d.flow.placeholders import ConsField
from g3d.flow.placeholders import AuxvField
from g3d.flow.placeholders import PostField
from g3d.flow.placeholders import VGradField
from g3d.flow.placeholders import DynamicTurbulenceModelVariables
from g3d.flow.placeholders import FMPArray
from g3d.flow.placeholders import MatrixRows
from g3d.flow.operators import IRSMOperator
from g3d.flow.domain import BoCo
from g3d.flow.domain import BoundaryConditionPlaceHolder
from g3d.flow.domain import InterfacePlaceHolder
from g3d.flow.io import G3DIO
from g3d.flow.domain import TurbulenceModel
from g3d.flow.domain import Smagorinsky
[docs]
class Domain(G3DEntity):
""" """
def __init__(self, *args: Any, **kwargs: Any) -> None:
self._gas = Gas()
"""Gas model"""
self._turb = TurbulenceModel()
"""Turbulence model"""
self._mesh = Mesh()
"""Computational grid: block-structured multi-block mesh"""
self._cons = ConsField()
"""Conservative flow variables container (solver variables: ro, ru, rv, rw, et)"""
self._auxv = AuxvField()
"""Primitive flow variables container (solver auxiliary variables: p, u, v, w, t, c)"""
self._post = PostField()
"""Postprocessing variables container (Po, To, M, s)"""
self._vgrad = VGradField()
"""Velocity gradient container (dudx, dudy, dudz, dvdx, dvdy, dvdz, dwdx, dwdy, dwdz)"""
self._dtmv = DynamicTurbulenceModelVariables()
"""Turbulence model variables container (visct)"""
self._boundary_faces = BoundaryConditionPlaceHolder()
"""Boundary condition container"""
self._connectivity_faces = InterfacePlaceHolder()
"""Block interface connectivity container"""
self._internal_faces = FMPArray()
"""Mesh-block internal face container"""
super().__init__(*args, **kwargs)
def _set_default(self) -> None:
""" """
super()._set_default()
self._tag = 'g3d_domain'
self._bc_mods = []
self.__io = G3DIO()
self.__initsol = Cons()
self.__deallocate_internal_faces()
self.__IRSM_operator = IRSMOperator()
self.__reset_face_counts()
@property
def mesh(self) -> Mesh:
"""Computational mesh instance
Returns
-------
Mesh
mesh object
"""
return self._mesh
@property
def gas(self) -> Gas:
"""Gas model instance
Returns
-------
Gas
gas model
"""
return self._gas
@gas.setter
def gas(self, g: Gas) -> None:
self._gas = g
@property
def turb(self) -> TurbulenceModel:
"""Turbulence model instance
Returns
-------
TurbulenceModel
turbulence model
"""
return self._turb
@property
def cons(self) -> ConsField:
"""Conservative variables placeholder
Returns
-------
ConsField
conservative variables object
"""
return self._cons
@property
def auxv(self) -> AuxvField:
"""Auxiliary variables placeholder
Returns
-------
AuxvField
auxiliary variables object
"""
return self._auxv
@property
def post(self) -> PostField:
"""Postprocessing variables placeholder
Returns
-------
PostField
postprocessing variables object
"""
return self._post
@property
def vgrad(self) -> VGradField:
"""Velocity gradients placeholder
Returns
-------
VGradField
velocity gradients object
"""
return self._vgrad
@property
def dtmv(self) -> DynamicTurbulenceModelVariables:
"""Dynamic turbulence model variables placeholder
Returns
-------
DynamicTurbulenceModelVariables
dynamic turbulence model variables object
"""
return self._dtmv
@property
def bc_mods(self) -> list:
"""Boundary condition modifiers array
Returns
-------
list
array of boundary condition modifiers
"""
return self._bc_mods
@property
def boundary_faces(self) -> BoundaryConditionPlaceHolder:
"""Boundary conditions container
Returns
-------
BoundaryConditionPlaceHolder
boundary conditions placeholder object
"""
return self._boundary_faces
@property
def connectivity_faces(self) -> InterfacePlaceHolder:
"""Grid connectivity container (block interfaces)
Returns
-------
InterfacePlaceHolder
interface placeholder object
"""
return self._connectivity_faces
@property
def internal_faces(self) -> FMPArray:
"""Internal grid faces container
Returns
-------
FMPArray
flux molecule array object
"""
return self._internal_faces
[docs]
def get_max_FMP_size(self) -> int:
"""Gets the maximum flux molecule array size
Returns
-------
int
size
"""
return max([self._internal_faces.get_max_FMP_size(),
self._boundary_faces.get_max_FMP_size(),
self._connectivity_faces.get_max_FMP_size()])
[docs]
def cell_get_block(self, kc) -> int:
"""Get block number associated with specified cell index
Parameters
----------
kc : int
cell index
Returns
-------
int
block number
"""
return self._mesh.cell_get_block(kc)
[docs]
def node_get_block(self, kp) -> int:
"""Get block number associated with specified node index
Parameters
----------
kp : int
node index
"""
return self._mesh.node_get_block(kp)
[docs]
def cell_get_ijk(self, kc) -> Tuple[int, int, int]:
"""Gets cell indices in computational space (i, j, k)
Parameters
----------
kc : int
cell index
"""
return self._mesh.cell_get_ijk(kc)
[docs]
def node_get_ijk(self, kp) -> Tuple[int, int, int]:
"""Gets node indices in computational space (i, j, k)
Parameters
----------
kp : int
node index
"""
return self._mesh.node_get_ijk(kp)
[docs]
def get_block_sizes(self, im) -> Tuple[int, int, int]:
"""Gets block node counts in three directions (ni, nj, nk)
Parameters
----------
im : int
block number
Returns
-------
tuple
block node counts (ni, nj, nk)
"""
return self._mesh.blocks[im].get_block_sizes()
[docs]
def get_node(self, im, i, j, k) -> int:
"""Gets node index based on provided data
Parameters
----------
im : int
block number
i : int
node index in i-direction
j : int
node index in j-direction
k : int
node index in k-direction
Returns
-------
int
node index lp
"""
return self._mesh.blocks[im].get_node(i, j, k)
[docs]
def get_cell(self, im, i, j, k) -> int:
"""Gets cell index based on provided data
Parameters
----------
im : int
block number
i : int
cell index in i-direction (lowest node i-index in cell)
j : int
cell index in j-direction (lowest node j-index in cell)
k : int
cell index in k-direction (lowest node k-index in cell)
Returns
-------
int
cell index lc
"""
return self._mesh.blocks[im].get_cell(i, j, k)
[docs]
def get_cell_centers(self) -> Tuple[ndarray, ...]:
"""Gets cell center coordinates, (x, y, z)-vectors, for all cell in the computational domain
Returns
-------
tuple
cell center coordinates (tuple of arrays)
"""
return self._mesh.get_cell_centers()
[docs]
def get_IRSM_operator(self) -> IRSMOperator:
"""Gets IRMS operator object
Returns
-------
IRSMOperator
IRSM operator object
"""
return self.__IRSM_operator
[docs]
def get_gas_properties(self) -> Tuple[float, float]:
"""Gets gas model properties (ratio of specific heats (gamma) & gas constant (R))
Returns
-------
tuple
ratio of specific heats (gamma) & gas constant (R)
"""
return self._gas.gamma, self._gas.R
[docs]
def set_gas_properties(self, gamma: float, R: float) -> None:
"""Sets gas model properties (ratio of specific heats (gamma) & gas constant (R))
Parameters
----------
gamma : float
ratio of specific heats
R : float
gas constant
"""
self._gas.gamma = gamma
self._gas.R = R
def __initialize_flow_field(self) -> None:
self._cons.set_from_Cons(self.__initsol)
self.auxvr()
def __region_initsol(self, r: IJKRegion) -> None:
im, i1, i2, j1, j2, k1, k2 = r.get_def()
for k in range(k1 - 1, k2 - 1):
for j in range(j1 - 1, j2 - 1):
for i in range(i1 - 1, i2 - 1):
self._cons.set_by_index_from_Cons(self._mesh.blocks[im - 1].get_node(i, j, k), self.__initsol)
def __block_initsol(self, im: int) -> None:
idx1 = self._mesh.blocks[im - 1].nmc
idx2 = idx1 + self._mesh.blocks[im - 1].n_cells
for ic in range(idx1, idx2):
self._cons.set_by_index_from_Cons(ic, self.__initsol)
[docs]
def initsol_conservative(self, ro: float, ru: float, rv: float, rw: float, et: float) -> None:
"""Initializes flow field using specified flow property values for all cells (conservative variables)
Parameters
----------
ro : float
density [kg/m^3]
ru : float
x-momentum [kg/(m^2 s)]
rv : float
y-momentum [kg/(m^2 s)]
rw : float
z-momentum [kg/(m^2 s)]
et : float
total energy times density [J/m^3]
"""
self.__initsol.set(ro, ru, rv, rw, et)
self.__initialize_flow_field()
[docs]
def initsol_primitive(self, p: float, u: float, v: float, w: float, t: float) -> None:
"""Initializes flow field using specified flow property values for all cells (primitive variables)
Parameters
----------
p : float
pressure [Pa]
u : float
x-velocity [m/s]
v : float
y-velocity [m/s]
w : float
z-velocity [m/s]
t : float
temperature [K]
"""
self.__initsol.set_from_primitive(p, u, v, w, t, self._gas)
self.__initialize_flow_field()
[docs]
def initsol_from_Cons(self, q: Cons) -> None:
"""Initialize flow field variables from specified conservative variables container
Parameters
----------
q : Cons
conservative variables container object
"""
self.__initsol.copy(q)
self.__initialize_flow_field()
[docs]
def initsol_from_Auxv(self, p: Auxv) -> None:
"""Initialize flow field variables from specified primitive variables container
Parameters
----------
p : Auxv
primitive variables container object
"""
q = p.get_Cons(self._gas)
self.__initsol.copy(q)
self.__initialize_flow_field()
[docs]
def initsol_region_from_Cons(self, r: IJKRegion, q: Cons) -> None:
"""Initialize flow field variables in specified region using specified conservative variables object
Parameters
----------
r : IJKRegion
computational space region
q : Cons
conservative variables container object
"""
self.__initsol.copy(q)
self.__region_initsol(r)
[docs]
def initsol_region_from_Auxv(self, r: IJKRegion, p: Auxv) -> None:
"""Initialize flow field variables in specified region using specified primitive variables object
Parameters
----------
r : IJKRegion
computational space region
p : Auxv
primitive variables container object
"""
q = p.get_Cons(self._gas)
self.__initsol.copy(q)
self.__region_initsol(r)
[docs]
def initsol_block_from_Cons(self, im: int, q: Cons) -> None:
"""Initialize flow field variables in specified mesh block using specified conservative variables object
Parameters
----------
im : int
block number
q : Cons
conservative variables container object
"""
self.__initsol.copy(q)
self.__block_initsol(im)
[docs]
def initsol_block_from_Auxv(self, im: int, p: Auxv) -> None:
"""Initialize flow field variables in specified mesh block using specified primitive variables object
Parameters
----------
im : int
block number
p : Auxv
primitive variables container object
"""
q = p.get_Cons(self._gas)
self.__initsol.copy(q)
self.__block_initsol(im)
[docs]
def initsol_linear(self, ql: Cons, qr: Cons, direction: str = 'x') -> None:
"""Initialize flow field variables using linear interpolation between specified left and right values
(conservative variable container objects)
Parameters
----------
ql : Cons
conservative variable container objects (left-end values)
qr : Cons
conservative variable container objects (right-end values)
direction :
interpolation direction, (Default value = 'x')
"""
# TODO maybe change from string input to DIRECTION Enum?
g3d_world = G3DWorld()
ro, ru, rv, rw, et = self.get_cons_arrays()
xc, yc, zc = self.get_cell_centers()
if direction == 'x':
x_min, x_max = self._mesh.get_x_min_x_max()
length = float(x_max) - float(x_min)
for ic in range(g3d_world.sizes.n_cells):
f = (xc[ic] - float(x_min)) / length
ro[ic] = ql.ro + (qr.ro - ql.ro) * f
ru[ic] = ql.ru + (qr.ru - ql.ru) * f
rv[ic] = ql.rv + (qr.rv - ql.rv) * f
rw[ic] = ql.rw + (qr.rw - ql.rw) * f
et[ic] = ql.et + (qr.et - ql.et) * f
elif direction == 'y':
y_min, y_max = self._mesh.get_y_min_y_max()
length = float(y_max) - float(y_min)
for ic in range(g3d_world.sizes.n_cells):
f = (yc[ic] - float(y_min)) / length
ro[ic] = ql.ro + (qr.ro - ql.ro) * f
ru[ic] = ql.ru + (qr.ru - ql.ru) * f
rv[ic] = ql.rv + (qr.rv - ql.rv) * f
rw[ic] = ql.rw + (qr.rw - ql.rw) * f
et[ic] = ql.et + (qr.et - ql.et) * f
elif direction == 'z':
z_min, z_max = self._mesh.get_z_min_z_max()
length = float(z_max) - float(z_min)
for ic in range(g3d_world.sizes.n_cells):
f = (zc[ic] - float(z_min)) / length
ro[ic] = ql.ro + (qr.ro - ql.ro) * f
ru[ic] = ql.ru + (qr.ru - ql.ru) * f
rv[ic] = ql.rv + (qr.rv - ql.rv) * f
rw[ic] = ql.rw + (qr.rw - ql.rw) * f
et[ic] = ql.et + (qr.et - ql.et) * f
else:
print('specified interpolation direction (%s) not supported' % direction)
print('must be one of ''x'', ''y'', or ''z''')
# conservative variables (solver variables)
[docs]
def get_cons_arrays(self) -> Tuple[ndarray, ...]:
"""Gets conservative variable arrays (solver variables)
Returns
-------
tuple
tuple of conservative variable arrays (ro, ru, rv, rw, et)
"""
return self._cons.get_arrays()
[docs]
def get_cons_array(self, var: str = 'ro') -> ndarray:
"""Gets specified conservative variable array (solver variable)
Parameters
----------
var : str
conservative variable (any of 'ro', 'ru', 'rv', 'rw', 'et'), defaults to 'ro'
var: str :
(Default value = 'ro')
Returns
-------
ndarray
variable array
"""
return self._cons.get_array(var)
[docs]
def get_cons_arrays_with_ghosts(self) -> Tuple[ndarray, ...]:
"""Gets conservative variable arrays (solver variables) including ghost cell values
Returns
-------
tuple
tuple of conservative variable arrays (ro, ru, rv, rw, et)
"""
return self._cons.get_arrays_with_ghosts()
[docs]
def get_cons_array_with_ghosts(self, var: str = 'ro') -> ndarray:
"""Gets specified conservative variable array (solver variable) including ghost cell values
Parameters
----------
var : str
conservative variable (any of 'ro', 'ru', 'rv', 'rw', 'et'), (Default value = 'ro')
Returns
-------
ndarray
variable array
"""
return self._cons.get_array_with_ghosts(var)
# auxiliary variables
[docs]
def get_auxv_arrays(self) -> Tuple[ndarray, ...]:
"""Gets primitive variable arrays (auxiliary solver variables)
Returns
-------
tuple
tuple of primitive variable arrays (p, u, v, w, t, c)
"""
return self._auxv.get_arrays()
[docs]
def get_auxv_array(self, var: str = 'p') -> ndarray:
"""Gets specified primitive variable array (auxiliary solver variable)
Parameters
----------
var : str
primitive variable (any of 'p', 'u', 'v', 'w', 't', 'c), (Default value = 'p')
Returns
-------
ndarray
variable array
"""
return self._auxv.get_array(var)
[docs]
def get_auxv_velocity_arrays(self) -> Tuple[ndarray, ...]:
"""Gets velocity vectors auxiliary solver variables placeholder
Returns
-------
tuple
tuple of velocity arrays (u, v, w)
"""
return self._auxv.get_velocity_arrays()
[docs]
def get_auxv_arrays_with_ghosts(self) -> Tuple[ndarray, ...]:
"""Gets primitive variable arrays (auxiliary solver variables) including ghost cell values
Returns
-------
tuple
tuple of primitive variable arrays (p, u, v, w, t, c)
"""
return self._auxv.get_arrays_with_ghosts()
[docs]
def get_auxv_array_with_ghosts(self, var: str = 'p') -> ndarray:
"""Gets specified primitive variable array (auxiliary solver variable) including ghost cell values
Parameters
----------
var : str
primitive variable (any of 'p', 'u', 'v', 'w', 't', 'c'), (Default value = 'p')
Returns
-------
ndarray
variable array
"""
return self._auxv.get_array_with_ghosts(var)
[docs]
def get_auxv_velocity_arrays_with_ghosts(self) -> Tuple[ndarray, ...]:
"""Gets velocity vectors auxiliary solver variables placeholder including ghost cell values
Returns
-------
tuple
tuple of velocity arrays (u, v, w)
"""
return self._auxv.get_velocity_arrays_with_ghosts()
# dynamic turbulence model variables
[docs]
def get_dtmv_arrays(self) -> Tuple[ndarray, ...]:
"""Gets dynamic turbulence model variable arrays
Returns
-------
tuple
tuple of dynamic turbulence model variable arrays (visct)
"""
return self._dtmv.get_arrays()
[docs]
def get_dtmv_array(self, var: str = 'visct') -> ndarray:
"""Gets specified dynamic turbulence model variable array
Parameters
----------
var : str
dynamic turbulence model variable (any of 'visct'), (Default value = 'visct')
Returns
-------
ndarray
variable array
"""
return self._dtmv.get_array(var)
[docs]
def get_dtmv_arrays_with_ghosts(self) -> Tuple[ndarray, ...]:
"""Gets dynamic turbulence model variable arrays including ghost cell values
Returns
-------
tuple
tuple of dynamic turbulence model variable arrays (visct)
"""
return self._dtmv.get_arrays_with_ghosts()
[docs]
def get_dtmv_array_with_ghosts(self, var: str = 'visct') -> ndarray:
"""Gets specified dynamic turbulence model variable array including ghost cell values
Parameters
----------
var : str
dynamic turbulence model variable (any of 'visct'), (Default value = 'visct')
Returns
-------
ndarray
variable array
"""
return self._dtmv.get_array_with_ghosts(var)
[docs]
def get_visct_array(self) -> ndarray:
"""Gets turbulent viscosity array
Returns
-------
ndarray
turbulent viscosity array (visct)
"""
return self._dtmv.get_visct_array()
# postprocessing variables
[docs]
def get_post_arrays(self) -> Tuple[ndarray, ...]:
"""Gets postprocessing variable arrays
Returns
-------
tuple
tuple of postprocessing variable arrays (Po, To, M, s)
"""
return self._post.get_arrays()
[docs]
def get_post_array(self, var: str = 'M') -> ndarray:
"""Gets specified postprocessing variable array
Parameters
----------
var : str
primitive variable (any of 'p0', 't0', 'M', 's'), defaults to 'M'
var: str :
(Default value = 'M')
Returns
-------
ndarray
variable array
"""
return self._post.get_array(var)
[docs]
def get_post_arrays_with_ghosts(self) -> Tuple[ndarray, ...]:
"""Gets postprocessing variable arrays including ghost cell values
Returns
-------
tuple
tuple of postprocessing variable arrays (Po, To, M, s)
"""
return self._post.get_arrays_with_ghosts()
[docs]
def get_post_array_with_ghosts(self, var: str = 'M') -> ndarray:
"""Gets specified postprocessing variable array including ghost cell values
Parameters
----------
var : str
primitive variable (any of 'p0', 't0', 'M', 's'), defaults to 'M'
var: str :
(Default value = 'M')
Returns
-------
ndarray
variable array
"""
return self._post.get_array_with_ghosts(var)
# velocity gradients
[docs]
def get_vgrad_arrays(self) -> Tuple[ndarray, ...]:
"""Gets velocity gradient arrays
Returns
-------
tuple
tuple of velocity gradient arrays (dudx, dudy, dudz, dvdx, dvdy, dvdz, dwdx, dwdy, dwdz)
"""
return self._vgrad.get_arrays()
[docs]
def get_vgrad_Array(self, var: str = 'dudx') -> ndarray:
"""Gets specified velocity gradient array
Parameters
----------
var : str
primitive variable (any of 'dudx', 'dudy', 'dudz', 'dvdx', 'dvdy', 'dvdz', 'dwdx', 'dwdy', 'dwdz'),
defaults to 'dudx'
var: str :
(Default value = 'dudx')
Returns
-------
ndarray
variable array
"""
return self._vgrad.get_array(var)
[docs]
def get_vgrad_arrays_with_ghosts(self) -> Tuple[ndarray, ...]:
"""Gets velocity gradient arrays including ghost cell values
Returns
-------
tuple
tuple of velocity gradient arrays (dudx, dudy, dudz, dvdx, dvdy, dvdz, dwdx, dwdy, dwdz)
"""
return self._vgrad.get_arrays_with_ghosts()
[docs]
def get_vgrad_array_with_ghosts(self, var: str = 'dudx') -> ndarray:
"""Gets specified velocity gradient array including ghost cell values
Parameters
----------
var : str
primitive variable (any of 'dudx', 'dudy', 'dudz', 'dvdx', 'dvdy', 'dvdz', 'dwdx', 'dwdy', 'dwdz'),
(Default value = 'dudx')
Returns
-------
ndarray
variable array
"""
return self._vgrad.get_array_with_ghosts(var)
# node coordinates
[docs]
def get_node_coordinates(self) -> Tuple[ndarray, ...]:
"""Gets node coordinate arrays
Returns
-------
tuple
tuple of node coordinate arrays (x, y, z)
"""
return self._mesh.nodes.get_XYZ()
# boundary conditions and interfaces
[docs]
def bc_update_speed_of_sound(self, boco: BoCo) -> None:
"""Updates boundary condition speed of sound
Parameters
----------
boco : BoCo
boundary condition container
"""
boco.bc_update_speed_of_sound(self._gas)
[docs]
def bc_update_primitive(self, boco: BoCo) -> None:
"""Updates boundary condition primitive variables (p, u, v, w, t, c)
Parameters
----------
boco : BoCo
boundary condition container
"""
boco.bc_update_primitive(self._gas)
[docs]
def bc_update_conservative(self, boco: BoCo) -> None:
"""Updates boundary condition conservative variables (ro, ru, rv, rw, et)
Parameters
----------
boco : BoCo
boundary condition container
"""
boco.bc_update_conservative(self._gas)
[docs]
def auxvr(self) -> None:
"""Updates flow domain primitive variable fields (p, u, v, w, t, c)"""
g3d_world = G3DWorld()
self._auxv.update(self._cons, self._gas)
if g3d_world.options.use_velocity_gradients:
self.__calculate_velocity_gradients()
[docs]
def print_flow_field_extreme_values(self) -> None:
"""Prints extreme values (min/max) of density (ro), velocity (u, v, w), pressure (p), and temperature (t)"""
self._cons.print_extreme_values('ro')
self._auxv.print_extreme_values(['u', 'v', 'w', 'p', 't'])
[docs]
def print_velocity_gradients_extreme_values(self) -> None:
"""
prints velocity gradient extreme values (min/max) of (dudx, dudy, dudz, dvdx, dvdy, dvdz, dwdx, dwdy, dwdz)
"""
self._vgrad.print_extreme_values()
@staticmethod
def __reset_face_counts() -> None:
g3d_sizes = G3DSizes()
g3d_sizes.n_internal_faces = 0
g3d_sizes.n_internal_faces_i = 0
g3d_sizes.n_internal_faces_j = 0
g3d_sizes.n_internal_faces_k = 0
g3d_sizes.n_block_faces = 0
g3d_sizes.n_connectivity_faces = 0
g3d_sizes.n_boundary_faces = 0
def __allocate_flux_molecule_placeholders(self) -> None:
g3d_sizes = G3DSizes()
# -- update face counts
self.__reset_face_counts()
for b in self._mesh.blocks:
ni, nj, nk = b.get_block_sizes()
g3d_sizes.add_internal_faces_i(max((ni - 4) * (nj - 1) * (nk - 1), 0))
g3d_sizes.add_internal_faces_j(max((nj - 4) * (ni - 1) * (nk - 1), 0))
g3d_sizes.add_internal_faces_k(max((nk - 4) * (ni - 1) * (nj - 1), 0))
if ni > 2:
g3d_sizes.add_block_faces(2 * (nj - 1) * (nk - 1))
if nj > 2:
g3d_sizes.add_block_faces(2 * (ni - 1) * (nk - 1))
if nk > 2:
g3d_sizes.add_block_faces(2 * (ni - 1) * (nj - 1))
g3d_sizes.n_internal_faces = (g3d_sizes.n_internal_faces_i +
g3d_sizes.n_internal_faces_j +
g3d_sizes.n_internal_faces_k)
# -- allocate flux molecule placeholders
if g3d_sizes.n_internal_faces_i > 0:
self._internal_faces.add_FMP(size=g3d_sizes.n_internal_faces_i, type=FaceType.I_FACE)
if g3d_sizes.n_internal_faces_j > 0:
self._internal_faces.add_FMP(size=g3d_sizes.n_internal_faces_j, type=FaceType.J_FACE)
if g3d_sizes.n_internal_faces_k > 0:
self._internal_faces.add_FMP(size=g3d_sizes.n_internal_faces_k, type=FaceType.K_FACE)
g3d_sizes.n_boundary_faces = 0
for bc in self._boundary_faces.get_boundary_array():
bc.allocate_FMPs()
g3d_sizes.add_boundary_faces(bc.get_FMP_size())
g3d_sizes.n_connectivity_faces = 0
for intf in self._connectivity_faces.get_interface_array():
intf.allocate_FMPs()
g3d_sizes.add_connectivity_faces(intf.get_FMP_size())
self._connectivity_faces.print()
self._boundary_faces.print()
self.__FMP_size_report()
def __deallocate_internal_faces(self) -> None:
self._internal_faces.deallocate_FMPs()
def __FMP_size_report(self) -> None:
g3d_sizes = G3DSizes()
g3d_world = G3DWorld()
g3d_world.display.header('face counts')
print(' internal faces: %5d' % g3d_sizes.n_internal_faces)
print(' block interfaces: %5d' % g3d_sizes.n_connectivity_faces)
print(' boundary faces: %5d' % g3d_sizes.n_boundary_faces)
face_check = (g3d_sizes.n_block_faces * 2 -
g3d_sizes.n_connectivity_faces / 3 -
g3d_sizes.n_connectivity_faces -
g3d_sizes.n_boundary_faces -
self._boundary_faces.get_face_count(BCTypes.SING))
if face_check != 0:
print(' face check: %d' % face_check)
if face_check > 0:
print(' there are undefined faces')
if face_check < 0:
print(' there are over-defined faces')
print()
[docs]
def export_mesh(self, *args: Any, **kwargs: Any) -> None:
"""Exports mesh to file
Parameters
----------
*args : Any :
arguments passed on to mesh exporter (BinaryFormat), optional
**kwargs : Any :
keyword-value pairs passed on to mesh exporter, optional
valid key values: 'format' binary_format (defaults to BinaryFormat.NATIVE)
"""
self._mesh.export(*args, **kwargs)
[docs]
def export_flow(self, *args: Any, **kwargs: Any) -> None:
"""Exports flow field data to file
Parameters
----------
*args : Any :
arguments passed on to flow data exporter (filename (str)), optional
**kwargs : Any :
keyword-value pairs passed on to flow data exporter, 'filename' -> name (str), optional
"""
name = ''
if len(args) > 0:
for arg in args:
if type(arg) is str:
name = arg
for key, value in kwargs.items():
if key == 'name' or key == 'filename':
name = value
self.auxvr()
self._post.update(self._cons, self._auxv, self._gas)
print()
self.print_flow_field_extreme_values()
print()
if name == '':
g3d_world = G3DWorld()
base_name = clean_path(g3d_world.io.wd + '/' + g3d_world.io.jobname)
sol_file = clean_path(base_name + '.' + str(g3d_world.status.curiter) + '/' + g3d_world.io.jobname)
self.__io.export_flow(self, self._cons, self._auxv, self._post, name=sol_file, xmlname=base_name)
else:
self.__io.export_flow(self, self._cons, self._auxv, self._post, **kwargs)
[docs]
def export_velocity_gradients(self, *args: Any, **kwargs: Any) -> None:
"""Exports velocity gradient data to file
Parameters
----------
*args : Any :
arguments passed on to flow data exporter (filename (str)), optional
**kwargs : Any :
keyword-value pairs passed on to flow data exporter, 'filename' -> name (str), optional
"""
name = ''
if len(args) > 0:
for arg in args:
if type(arg) is str:
name = arg
for key, value in kwargs.items():
if key == 'name' or key == 'filename':
name = value
g3d_world = G3DWorld()
tmp = g3d_world.options.use_velocity_gradients
g3d_world.options.use_velocity_gradients = True
if not tmp:
self._vgrad.setup(n_cells=g3d_world.sizes.n_cells, n_ghosts=g3d_world.sizes.n_ghost_cells)
self._vgrad.allocate_flux_molecule_placeholders(self.get_max_FMP_size())
self.auxvr()
print()
self.print_flow_field_extreme_values()
print()
self.print_velocity_gradients_extreme_values()
print()
if name == '':
base_name = clean_path(g3d_world.io.wd + '/' + g3d_world.io.jobname)
sol_file = clean_path(base_name + '.' + str(g3d_world.status.curiter) + '/' + g3d_world.io.jobname)
self.__io.export_velocity_gradients(self, self._vgrad, name=sol_file, xmlname=base_name)
else:
self.__io.export_velocity_gradients(self, self._vgrad, **kwargs)
g3d_world.options.use_velocity_gradients = tmp
if not tmp:
self._vgrad.setup(size=0)
[docs]
def update_metrics(self) -> None:
"""Updates metric data - cell volumes and face areas"""
self._mesh.update_metrics()
[docs]
def generate_mesh_from_pipe_data(self, q1d_pipe: Q1DPipe) -> None:
"""Generates computational mesh from provided Q1D domain object
Parameters
----------
q1d_pipe : Q1DPipe
quasi-one-dimensional domain object
"""
ni = q1d_pipe.n_cells + 1
nj = 2
nk = 2
self._mesh.add_Block(ni, nj, nk)
self._mesh.allocate_arrays()
self._mesh.get_mesh_coordinates_from_pipe_data(q1d_pipe)
# allocate flow field data placeholders
self.__initial_flow_field_allocation()
self._boundary_faces = BoundaryConditionPlaceHolder()
self._boundary_faces.add_Q1D_boundary_definitions(ni, q1d_pipe.left_end_BC, q1d_pipe.right_end_BC)
[docs]
def import_mesh(self, *args: Any, **kwargs: Any) -> None:
"""Imports mesh from file
Parameters
----------
*args : Any :
arguments (BinaryFormat, file base name (str)), optional
**kwargs : Any :
keyword value pairs, optional
valid key values: 'format' binary_format (defaults to BinaryFormat.NATIVE),
'msh_file' -> filename, 'topo_file' -> filename
"""
g3d_world = G3DWorld()
self._mesh.import_mesh(*args, **kwargs)
# allocate flow field data placeholders
self.__initial_flow_field_allocation()
if g3d_world.options.dim_mode == DIMMode.Q1D:
ni, nj, nk = self._mesh.blocks[0].get_block_sizes()
self._boundary_faces = BoundaryConditionPlaceHolder()
self._boundary_faces.add_Q1D_boundary_definitions(ni)
[docs]
def import_boundary_definitions(self, **kwargs: Any) -> None:
"""Imports boundary definitions from file
Parameters
----------
**kwargs : Any :
keyword-value pairs ('filename' -> name (str)), optional
"""
g3d_world = G3DWorld()
filename = ''
for key, value in kwargs.items():
if key == 'filename':
filename = value
if len(filename) == 0:
filename = clean_path(g3d_world.io.wd + '/' + g3d_world.io.jobname + '.bc')
with open(filename, 'r') as f:
line = f.readline()
while line:
line = line.strip()
for bc in BCTypes:
if line.find(bc.name) == 0:
self._boundary_faces.read_from_file(f, bc)
for intf in INTFTypes:
if line.find(intf.name) == 0:
self._connectivity_faces.read_from_file(f, intf)
line = f.readline()
for bc in self._boundary_faces.get_boundary_array():
for a in bc.get_face_defs():
self._mesh.check_face_def(a)
for intf in self._connectivity_faces.get_interface_array():
A, B = intf.get_face_defs()
for a in A:
self._mesh.check_face_def(a)
for b in B:
self._mesh.check_face_def(b, FaceType.B_SIDE)
[docs]
def update_BC(self, boco: BoCo, **kwargs: Any) -> None:
"""Updates boundary face boundary conditions from provided boundary condition container
Parameters
----------
boco : BoCo
boundary condition container
**kwargs : Any
keyword-vale pair ('verbose' -> value (bool)), optional
"""
verbose = True
for key, value in kwargs.items():
if key == 'verbose':
verbose = value
bc_found = False
for bc in self._boundary_faces.get_boundary_array():
for a in bc.get_face_defs():
if a.tag == boco.tag:
bc_found = True
a.bc_update(boco)
a.update_status()
if verbose:
if bc_found:
print(' boundary patches tagged %s updated' % boco.tag)
else:
print(' there are no surface patches tagged %s' % boco.tag)
[docs]
def BC_mod(self) -> None:
"""Modifies boundary condition data going through specified boundary condition modifiers"""
for mod in self._bc_mods:
mod.evaluate(self._gas)
self.update_BC(mod, verbose=False)
[docs]
def update_BC_family(self, *args: Any, **kwargs: Any) -> None:
"""Updates boundary condition family according to specifications
Parameters
----------
*args : Any
arguments (tag (str), bc_type (BCTypes)), optional
**kwargs : Any
keyword-value pairs ('tag' -> tag (str), 'bc_type' -> bc_type (BCTypes)), optional
"""
tag = ''
bc_type = None
if len(args) > 0:
for arg in args:
if type(arg) is str:
tag = arg
elif type(arg) is BCTypes:
bc_type = arg
for key, value in kwargs.items():
if key == 'tag':
tag = value
elif key == 'bc_type':
bc_type = value
if type(bc_type) is BCTypes:
bc_found = False
for bc in self._boundary_faces.get_boundary_array():
for a in bc.get_face_defs():
if a.tag == tag:
bc_found = True
a.bc_update_family(bc_type)
bc.bc_type = bc_type
if bc_found:
print(' boundary patches tagged %s updated' % tag)
else:
print(' there are no surface patches tagged %s' % tag)
[docs]
def update_BC_type(self, *args: Any, **kwargs: Any) -> None:
"""Updates boundary condition type according to specifications
Parameters
----------
*args : Any
arguments (tag (str), bc_it (int)), optional
**kwargs : Any
keyword-value pairs ('tag' -> tag (str), 'bc_it' -> bc_it (int)), optional
"""
tag = ''
bc_it = None
if len(args) > 0:
for arg in args:
if type(arg) is str:
tag = arg
elif type(arg) is BCTypes:
bc_it = arg
for key, value in kwargs.items():
if key == 'tag':
tag = value
elif key == 'bc_it':
bc_it = value
bc_found = False
for bc in self._boundary_faces.get_boundary_array():
for a in bc.get_face_defs():
if a.tag == tag:
bc_found = True
a.bc_update_type(bc_it)
if bc_found:
print(' boundary patches tagged %s updated' % tag)
else:
print(' there are no surface patches tagged %s' % tag)
def __populate_internal_flux_molecule_placeholder(self) -> None:
g3d_sizes = G3DSizes()
# fctj = 0.
# fctk = 0.
# kc = [0] * 12
# kp = [0] * 12
kc = np.zeros(12, dtype=np.int64)
kp = np.zeros(12, dtype=np.int64)
nb = self._mesh.n_blocks
fmp_idx = 0
if g3d_sizes.n_internal_faces_i > 0:
# I-faces
fmp = self._internal_faces.FMPs[fmp_idx]
face_count = 0
for im in range(nb):
ni, nj, nk = self._mesh.blocks[im].get_block_sizes()
for k in range(nk - 1):
fctk = 0.25
mkm = -1
mkp = 1
if k == 0:
mkm = 0
fctk = 0.5
if k == (nk - 2):
mkp = 0
fctk = 0.5
if k == 0 and k == (nk - 2):
fctk = 0.
for j in range(nj - 1):
fctj = 0.25
mjm = -1
mjp = 1
if j == 0:
mjm = 0
fctj = 0.5
if j == (nj - 2):
mjp = 0
fctj = 0.5
if j == 0 and j == (nj - 2):
fctj = 0.
for i in range(2, (ni - 2)):
kc0 = self._mesh.blocks[im].get_cell(i, j, k)
kp0 = self._mesh.blocks[im].get_node(i, j, k)
kc[0] = kc0 - 2
kc[1] = kc0 - 1
kc[2] = kc0
kc[3] = kc0 + 1
kc[4] = kc0 - 1 + mjm * (ni - 1)
kc[5] = kc0 + mjm * (ni - 1)
kc[6] = kc0 - 1 + mjp * (ni - 1)
kc[7] = kc0 + mjp * (ni - 1)
kc[8] = kc0 - 1 + mkm * (ni - 1) * (nj - 1)
kc[9] = kc0 + mkm * (ni - 1) * (nj - 1)
kc[10] = kc0 - 1 + mkp * (ni - 1) * (nj - 1)
kc[11] = kc0 + mkp * (ni - 1) * (nj - 1)
kp[0] = kp0 - 1
kp[1] = kp0 - 1 + ni
kp[2] = kp0 - 1 + ni * nj
kp[3] = kp0 - 1 + ni + ni * nj
kp[4] = kp0
kp[5] = kp0 + ni
kp[6] = kp0 + ni * nj
kp[7] = kp0 + ni + ni * nj
kp[8] = kp0 + 1
kp[9] = kp0 + 1 + ni
kp[10] = kp0 + 1 + ni * nj
kp[11] = kp0 + 1 + ni + ni * nj
# self._internal_faces.set_FMP_data(fmp_idx, face_count, kp, kc, fctj, fctk)
fmp.set_data(face_count, kp, kc, fctj, fctk)
face_count += 1
fmp_idx += 1
if g3d_sizes.n_internal_faces_j > 0:
# J-faces
fmp = self._internal_faces.FMPs[fmp_idx]
face_count = 0
for im in range(nb):
ni, nj, nk = self._mesh.blocks[im].get_block_sizes()
for k in range(nk - 1):
fctj = 0.25
mkm = -1
mkp = 1
if k == 0:
mkm = 0
fctj = 0.5
if k == (nk - 2):
mkp = 0
fctj = 0.5
if k == 0 and k == (nk - 2):
fctj = 0.
for j in range(2, (nj - 2)):
for i in range(ni - 1):
fctk = 0.25
mim = -1
mip = 1
if i == 0:
mim = 0
fctk = 0.5
if i == (ni - 2):
mip = 0
fctk = 0.5
if i == 0 and i == (ni - 2):
fctk = 0.
kc0 = self._mesh.blocks[im].get_cell(i, j, k)
kp0 = self._mesh.blocks[im].get_node(i, j, k)
kc[0] = kc0 - 2 * (ni - 1)
kc[1] = kc0 - 1 * (ni - 1)
kc[2] = kc0
kc[3] = kc0 + 1 * (ni - 1)
kc[4] = kc0 - 1 * (ni - 1) + mkm * (ni - 1) * (nj - 1)
kc[5] = kc0 + mkm * (ni - 1) * (nj - 1)
kc[6] = kc0 - 1 * (ni - 1) + mkp * (ni - 1) * (nj - 1)
kc[7] = kc0 + mkp * (ni - 1) * (nj - 1)
kc[8] = kc0 - 1 * (ni - 1) + mim
kc[9] = kc0 + mim
kc[10] = kc0 - 1 * (ni - 1) + mip
kc[11] = kc0 + mip
kp[0] = kp0 - ni
kp[1] = kp0 - ni + ni * nj
kp[2] = kp0 - ni + 1
kp[3] = kp0 - ni + ni * nj + 1
kp[4] = kp0
kp[5] = kp0 + ni * nj
kp[6] = kp0 + 1
kp[7] = kp0 + ni * nj + 1
kp[8] = kp0 + ni
kp[9] = kp0 + ni + ni * nj
kp[10] = kp0 + ni + 1
kp[11] = kp0 + ni + ni * nj + 1
# self._internal_faces.set_FMP_data(fmp_idx, face_count, kp, kc, fctj, fctk)
fmp.set_data(face_count, kp, kc, fctj, fctk)
face_count += 1
fmp_idx += 1
if g3d_sizes.n_internal_faces_k > 0:
# K - faces
fmp = self._internal_faces.FMPs[fmp_idx]
face_count = 0
for im in range(nb):
ni, nj, nk = self._mesh.blocks[im].get_block_sizes()
for k in range(2, (nk - 2)):
for j in range(nj - 1):
fctk = 0.25
mjm = -1
mjp = 1
if j == 0:
mjm = 0
fctk = 0.5
if j == (nj - 2):
mjp = 0
fctk = 0.5
if j == 0 and j == (nj - 2):
fctk = 0.
for i in range(ni - 1):
fctj = 0.25
mim = -1
mip = 1
if i == 0:
mim = 0
fctj = 0.5
if i == (ni - 2):
mip = 0
fctj = 0.5
if i == 0 and i == (ni - 2):
fctj = 0.
kc0 = self._mesh.blocks[im].get_cell(i, j, k)
kp0 = self._mesh.blocks[im].get_node(i, j, k)
kc[0] = kc0 - 2 * (ni - 1) * (nj - 1)
kc[1] = kc0 - 1 * (ni - 1) * (nj - 1)
kc[2] = kc0
kc[3] = kc0 + 1 * (ni - 1) * (nj - 1)
kc[4] = kc0 - 1 * (ni - 1) * (nj - 1) + mim
kc[5] = kc0 + mim
kc[6] = kc0 - 1 * (ni - 1) * (nj - 1) + mip
kc[7] = kc0 + mip
kc[8] = kc0 - 1 * (ni - 1) * (nj - 1) + mjm * (ni - 1)
kc[9] = kc0 + mjm * (ni - 1)
kc[10] = kc0 - 1 * (ni - 1) * (nj - 1) + mjp * (ni - 1)
kc[11] = kc0 + mjp * (ni - 1)
kp[0] = kp0 - ni * nj
kp[1] = kp0 - ni * nj + 1
kp[2] = kp0 - ni * nj + ni
kp[3] = kp0 - ni * nj + ni + 1
kp[4] = kp0
kp[5] = kp0 + 1
kp[6] = kp0 + ni
kp[7] = kp0 + ni + 1
kp[8] = kp0 + ni * nj
kp[9] = kp0 + ni * nj + 1
kp[10] = kp0 + ni * nj + ni
kp[11] = kp0 + ni * nj + ni + 1
# self._internal_faces.set_FMP_data(fmp_idx, face_count, kp, kc, fctj, fctk)
fmp.set_data(face_count, kp, kc, fctj, fctk)
face_count += 1
def __populate_boundary_flux_molecule_placeholder(self) -> None:
g3d_sizes = G3DSizes()
# kc = [0] * 12
# kp = [0] * 12
kc = np.zeros(12, dtype=np.int64)
kp = np.zeros(12, dtype=np.int64)
for bc in self._boundary_faces.get_boundary_array():
TYPE = bc.bc_type
for boco in bc.get_face_defs():
face_count = 0
boco_incr = int(boco.get_face_count())
im, i0, j0, k0, idir, jdir, npi, npj = boco.get_def()
l0c, l0p, ldic, ldjc, ldkc, ldip, ldjp, ldkp = self._mesh.extract_face_pointers(boco)
if TYPE == BCTypes.STREAM_TUBE_WALLS:
for i in range(12):
kc[i] = 0
kp[i] = 0
for jp in range(npj - 1):
for ip in range(npi - 1):
lp = l0p + ldip * ip + ldjp * jp
lc = l0c + ldic * ip + ldjc * jp
kc[1] = lc
kp[4] = lp
kp[5] = lp + ldip
kp[6] = lp + ldjp
kp[7] = lp + ldip + ldjp
boco.FMP.set_pointers(face_count, kp, kc)
face_count += 1
else:
for jp in range(npj - 1):
for ip in range(npi - 1):
lp = l0p + ldip * ip + ldjp * jp
lc = l0c + ldic * ip + ldjc * jp
mim, mip, mjm, mjp = self._mesh.check_bounds(lc, ldic, ldjc)
boco.FMP.update_factors(face_count, mim, mip, mjm, mjp)
lgc = ip + (npi - 1) * jp + g3d_sizes.n_cells
lgp = ip + npi * jp + g3d_sizes.n_nodes
kc[0] = lc - 2 * ldkc
kc[1] = lc - 1 * ldkc
kc[2] = lc
kc[3] = lgc
kc[4] = lc - ldkc + ldic * mim
kc[5] = lc + ldic * mim
kc[6] = lc - ldkc + ldic * mip
kc[7] = lc + ldic * mip
kc[8] = lc - ldkc + ldjc * mjm
kc[9] = lc + ldjc * mjm
kc[10] = lc - ldkc + ldjc * mjp
kc[11] = lc + ldjc * mjp
kp[0] = lp - 2 * ldkp
kp[1] = lp - 2 * ldkp + ldip
kp[2] = lp - 2 * ldkp + ldjp
kp[3] = lp - 2 * ldkp + ldip + ldjp
kp[4] = lp - ldkp
kp[5] = lp - ldkp + ldip
kp[6] = lp - ldkp + ldjp
kp[7] = lp - ldkp + ldip + ldjp
kp[8] = lp
kp[9] = lp + ldip
kp[10] = lp + ldjp
kp[11] = lp + ldip + ldjp
boco.FMP.set_pointers(face_count, kp, kc)
if TYPE != BCTypes.SING:
kc[0] = kc[1]
kc[1] = kc[2]
kc[2] = kc[3]
kc[3] = lgc + (npi - 1) * (npj - 1)
kc[4] = kc[5]
kc[5] = lgc + mim
kc[6] = kc[7]
kc[7] = lgc + mip
kc[8] = kc[9]
kc[9] = lgc + (npi - 1) * mjm
kc[10] = kc[11]
kc[11] = lgc + (npi - 1) * mjp
kp[0] = kp[4]
kp[1] = kp[5]
kp[2] = kp[6]
kp[3] = kp[7]
kp[4] = kp[8]
kp[5] = kp[9]
kp[6] = kp[10]
kp[7] = kp[11]
kp[8] = lgp
kp[9] = lgp + 1
kp[10] = lgp + npi
kp[11] = lgp + npi + 1
fctj, fctk = boco.FMP.get_factors(face_count)
boco.FMP.set_data(face_count + boco_incr, kp, kc, fctj, fctk)
face_count += 1
def __populate_connectivity_flux_molecule_placeholder(self) -> None:
# kc = [0] * 12
# kp = [0] * 12
kc = np.zeros(12, dtype=np.int64)
kp = np.zeros(12, dtype=np.int64)
for intf in self._connectivity_faces.get_interface_array():
A, B = intf.get_face_defs()
for i in range(len(A)):
face_count = 0
a = A[i]
b = B[i]
ima, i0a, j0a, k0a, idira, jdira, npia, npja = a.get_def()
# imb, i0b, j0b, k0b, idirb, jdirb, npib, npjb = b.get_def()
l0ac, l0ap, ldiac, ldjac, ldkac, ldiap, ldjap, ldkap = (
self._mesh.extract_face_pointers(a))
l0bc, l0bp, ldibc, ldjbc, ldkbc, ldibp, ldjbp, ldkbp = (
self._mesh.extract_face_pointers(b, FaceType.B_SIDE))
# one step in on the A-side of the interface
for jp in range(npja - 1):
for ip in range(npia - 1):
lap = l0ap + ldiap * ip + ldjap * jp
lac = l0ac + ldiac * ip + ldjac * jp
# lbp = l0bp + ldibp * ip + ldjbp * jp
lbc = l0bc + ldibc * ip + ldjbc * jp
mim, mip, mjm, mjp = self._mesh.check_bounds(lac, ldiac, ldjac)
a.FMP.update_factors(face_count, mim, mip, mjm, mjp)
kc[0] = lac - 2 * ldkac
kc[1] = lac - 1 * ldkac
kc[2] = lac
kc[3] = lbc
kc[4] = lac + ldiac * mim - ldkac
kc[5] = lac + ldiac * mim
kc[6] = lac + ldiac * mip - ldkac
kc[7] = lac + ldiac * mip
kc[8] = lac + ldjac * mjm - ldkac
kc[9] = lac + ldjac * mjm
kc[10] = lac + ldjac * mjp - ldkac
kc[11] = lac + ldjac * mjp
kp[0] = lap - 2 * ldkap
kp[1] = lap - 2 * ldkap + ldiap
kp[2] = lap - 2 * ldkap + ldjap
kp[3] = lap - 2 * ldkap + ldiap + ldjap
kp[4] = lap - ldkap
kp[5] = lap - ldkap + ldiap
kp[6] = lap - ldkap + ldjap
kp[7] = lap - ldkap + ldiap + ldjap
kp[8] = lap
kp[9] = lap + ldiap
kp[10] = lap + ldjap
kp[11] = lap + ldiap + ldjap
a.FMP.set_pointers(face_count, kp, kc)
face_count += 1
# block interface
for jp in range(npja - 1):
for ip in range(npia - 1):
lap = l0ap + ldiap * ip + ldjap * jp
lac = l0ac + ldiac * ip + ldjac * jp
lbp = l0bp + ldibp * ip + ldjbp * jp
lbc = l0bc + ldibc * ip + ldjbc * jp
mim, mip, mjm, mjp = self._mesh.check_bounds_block_interface(lac, ldiac, ldjac, lbc, ldibc,
ldjbc)
a.FMP.update_block_interface_factors(face_count, mim, mip, mjm, mjp)
kc[0] = lac - ldkac
kc[1] = lac
kc[2] = lbc
kc[3] = lbc + ldkbc
kc[4] = lac + ldiac * mim
kc[5] = lbc + ldibc * mim
kc[6] = lac + ldiac * mip
kc[7] = lbc + ldibc * mip
kc[8] = lac + ldjac * mjm
kc[9] = lbc + ldjbc * mjm
kc[10] = lac + ldjac * mjp
kc[11] = lbc + ldjbc * mjp
kp[0] = lap - ldkap
kp[1] = lap - ldkap + ldiap
kp[2] = lap - ldkap + ldjap
kp[3] = lap - ldkap + ldiap + ldjap
kp[4] = lap
kp[5] = lap + ldiap
kp[6] = lap + ldjap
kp[7] = lap + ldiap + ldjap
kp[8] = lbp + ldkbp
kp[9] = lbp + ldkbp + ldibp
kp[10] = lbp + ldkbp + ldjbp
kp[11] = lbp + ldkbp + ldibp + ldjbp
a.FMP.set_pointers(face_count, kp, kc)
face_count += 1
# one step in on the B-side of the interface
for jp in range(npja - 1):
for ip in range(npia - 1):
lap = l0ap + ldiap * ip + ldjap * jp
lac = l0ac + ldiac * ip + ldjac * jp
lbp = l0bp + ldibp * ip + ldjbp * jp
lbc = l0bc + ldibc * ip + ldjbc * jp
mim, mip, mjm, mjp = self._mesh.check_bounds(lbc, ldibc, ldjbc)
a.FMP.update_factors(face_count, mim, mip, mjm, mjp)
kc[0] = lac
kc[1] = lbc
kc[2] = lbc + ldkbc
kc[3] = lbc + 2 * ldkbc
kc[4] = lbc + ldibc * mim
kc[5] = lbc + ldkbc + ldibc * mim
kc[6] = lbc + ldibc * mip
kc[7] = lbc + ldkbc + ldibc * mip
kc[8] = lbc + ldjbc * mjm
kc[9] = lbc + ldkbc + ldjbc * mjm
kc[10] = lbc + ldjbc * mjp
kc[11] = lbc + ldkbc + ldjbc * mjp
kp[0] = lap
kp[1] = lap + ldiap
kp[2] = lap + ldjap
kp[3] = lap + ldiap + ldjap
kp[4] = lbp + ldkbp
kp[5] = lbp + ldkbp + ldibp
kp[6] = lbp + ldkbp + ldjbp
kp[7] = lbp + ldkbp + ldibp + ldjbp
kp[8] = lbp + 2 * ldkbp
kp[9] = lbp + 2 * ldkbp + ldibp
kp[10] = lbp + 2 * ldkbp + ldjbp
kp[11] = lbp + 2 * ldkbp + ldibp + ldjbp
a.FMP.set_pointers(face_count, kp, kc)
face_count += 1
def __update_flux_molecule_metrics(self) -> None:
X, Y, Z = self._mesh.nodes.get_XYZ()
self._internal_faces.update_FMP_metrics(X, Y, Z)
self._internal_faces.update_FMP_diff_metrics(X, Y, Z)
self._boundary_faces.update_FMP_metrics(X, Y, Z)
self._connectivity_faces.update_FMP_metrics(X, Y, Z)
[docs]
def get_volume(self, idx: int = None) -> Union[float, ndarray]:
"""Gets volume for specific cell if cell index is specified otherwise all cell volumes are returned as an array
Parameters
----------
idx : int
cell index, optional (Default value = None)
Returns
-------
float or ndarray
cell volume or array of cell volumes
"""
return self._mesh.get_volume(idx)
[docs]
def get_inverse_volume(self, idx: int = None) -> Union[float, ndarray]:
"""Gets inverse volume for specific cell if cell index is specified otherwise all cell inverse volumes are
returned as an array
Parameters
----------
idx : int
cell index, optional (Default value = None)
Returns
-------
float or ndarray
inverse cell volume or array of inverse cell volumes
"""
return self._mesh.get_inverse_volume(idx)
[docs]
def get_area(self, idx: int = None) -> Union[float, ndarray]:
"""Gets average area for specific cell if cell index is specified otherwise all average cell areas are
returned as an array
Parameters
----------
idx : int
cell index, optional (Default value = None)
Returns
-------
float or ndarray
average cell area or array of average cell areas
"""
return self._mesh.get_area(idx)
[docs]
def get_streamtube_normal(self, idx: int = None) -> Union[float, ndarray]:
"""Gets stream-tube normal for specific cell if cell index is specified otherwise all stream-tube normals are
returned as an array
Parameters
----------
idx : int
cell index, optional (Default value = None)
Returns
-------
float or ndarray
stream-tube normal or array of stream-tube normals
"""
return self._mesh.get_streamtube_normal(idx)
[docs]
def check_FMP_metrics(self) -> Tuple[ndarray, ...]:
"""Checks flux molecule cell metrics for all flux molecules (internal faces, mesh block connectivity faces,
and boundary faces)
Returns
-------
tuple
tuple of arrays with face check data (X, Y, Z)
"""
g3d_sizes = G3DSizes()
X = np.zeros(g3d_sizes.n_cells + g3d_sizes.n_ghost_cells)
Y = np.zeros(g3d_sizes.n_cells + g3d_sizes.n_ghost_cells)
Z = np.zeros(g3d_sizes.n_cells + g3d_sizes.n_ghost_cells)
self._internal_faces.check_FMP_metrics(X, Y, Z, 'I-faces')
self._internal_faces.check_FMP_metrics(X, Y, Z, 'J-faces')
self._internal_faces.check_FMP_metrics(X, Y, Z, 'K-faces')
self._boundary_faces.check_FMP_metrics(X, Y, Z)
self._connectivity_faces.check_FMP_metrics(X, Y, Z)
print(' SX (min) [max]: (%10.5E) [%10.5E]' % (min(X), max(X)))
print(' SY (min) [max]: (%10.5E) [%10.5E]' % (min(Y), max(Y)))
print(' SZ (min) [max]: (%10.5E) [%10.5E]' % (min(Z), max(Z)))
return X, Y, Z
[docs]
def check_cell_metrics(self, ic: int) -> None:
"""Checks cell metrics for cell with specified cell index
Parameters
----------
ic : int
cell index
"""
check_sum = [0., 0., 0.]
print(' internal faces:')
self._internal_faces.check_FMP_cell_metrics(ic, check_sum)
self._boundary_faces.check_FMP_cell_metrics(ic, check_sum)
self._connectivity_faces.check_FMP_cell_metrics(ic, check_sum)
print(' check sums for cell %d' % ic)
print(' ', check_sum)
def __populate_flux_molecule_arrays(self) -> None:
# -- update internal face flux molecule pointer arrays
self.__populate_internal_flux_molecule_placeholder()
# -- update boundary face flux molecule pointer arrays
self.__populate_boundary_flux_molecule_placeholder()
# -- update connectivity face flux molecule pointer arrays
self.__populate_connectivity_flux_molecule_placeholder()
def __initial_flow_field_allocation(self) -> None:
g3d_world = G3DWorld()
n_cells = g3d_world.sizes.n_cells
self._cons.setup(n_cells=n_cells)
self._auxv.setup(n_cells=n_cells)
self._post.setup(n_cells=n_cells)
if g3d_world.options.eqn_mode.value > EQNMode.DNS.value:
self._dtmv.setup(n_cells=n_cells)
if g3d_world.options.use_velocity_gradients:
self._vgrad.setup(n_cells=n_cells)
def __allocate_placeholders(self) -> None:
g3d_world = G3DWorld()
if g3d_world.display.verbose:
self._boundary_faces.ghost_report()
n_cells = g3d_world.sizes.n_cells
g3d_world.sizes.n_ghost_cells = n_ghost_cells = self._boundary_faces.get_max_ghost_cells()
g3d_world.sizes.n_ghost_nodes = self._boundary_faces.get_max_ghost_nodes()
self._mesh.redefine_ghosts()
self._cons.setup(n_cells=n_cells, n_ghosts=n_ghost_cells, backup_data=True)
self._auxv.setup(n_cells=n_cells, n_ghosts=n_ghost_cells, backup_data=True)
self._post.setup(size=n_cells, backup_data=True)
if g3d_world.options.eqn_mode.value > EQNMode.DNS.value:
self._dtmv.setup(n_cells=n_cells, n_ghosts=n_ghost_cells, backup_data=True)
if g3d_world.options.use_velocity_gradients:
self._vgrad.setup(n_cells=n_cells, n_ghosts=n_ghost_cells, backup_data=True)
self._vgrad.allocate_flux_molecule_placeholders(self.get_max_FMP_size())
[docs]
def setup_and_allocate(self) -> None:
"""Sets up and allocates all flow domain placeholders and data containers
updates flux molecule metrics for all flux molecule types (internal faces,
block connectivity faces, and boundary faces)
allocates turbulence model if applicable
"""
g3d_world = G3DWorld()
# -- setup models
if g3d_world.models.turbulence_model is TurbulenceModelMode.ERLEBACHER_SMAGORINSKY:
self._turb = Smagorinsky()
# -- allocate flux molecule placeholders (internal + boundary/connectivity faces)
self.__allocate_flux_molecule_placeholders()
# -- populate flux molecule placeholders
self.__populate_flux_molecule_arrays()
# -- re-allocate coordinate arrays and flow field placeholders
self.__allocate_placeholders()
# -- update face metrics
self.__update_flux_molecule_metrics()
# self.check_FMP_metrics()
if g3d_world.display.verbose:
print('max boundary FMP size: %d' % self._boundary_faces.get_max_FMP_size())
print('max connectivity FMP size: %d' % self._connectivity_faces.get_max_FMP_size())
def __calculate_velocity_gradients(self) -> None:
g3d_world = G3DWorld()
u, v, w = self._auxv.get_velocity_arrays_with_ghosts()
voli = self._mesh.get_inverse_volume()
# nullify gradient field
self._vgrad.nullify()
# add axisymm correction if applicable
if g3d_world.options.dim_mode is DIMMode.G2D_AXISYMM:
area = self._mesh.get_area()
dudy = self.get_vgrad_Array('dudy')
dvdy = self.get_vgrad_Array('dvdy')
dwdy = self.get_vgrad_Array('dwdy')
u, v, w = self.get_auxv_velocity_arrays()
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')
ic = self._mesh.blocks[im].get_cell(i, j, k)
dudy[ic] = -area[ic] * u[ic]
dvdy[ic] = -area[ic] * v[ic]
dwdy[ic] = -area[ic] * w[ic]
# internal face velocity gradient contributions
for tmp in self._internal_faces.get_FMPs():
self._vgrad.add_gradients(u, v, w, tmp)
# connectivity face velocity gradient contributions
for interface in self._connectivity_faces.get_interface_array():
A, B = interface.get_face_defs()
for intf in A:
self._vgrad.add_gradients(u, v, w, intf.FMP)
# boundary face velocity gradient contributions
for bc in self._boundary_faces.get_boundary_array():
for boco in bc.get_face_defs():
TYPE = boco.face_type
# IT = boco.it
tmp = boco.FMP
self._vgrad.add_gradients(u, v, w, tmp, 0, tmp.size / 2)
if TYPE != BCTypes.SING:
if TYPE is BCTypes.SYMM:
self._vgrad.add_gradients(u, v, w, tmp, tmp.size / 2, tmp.size / 2,
GradientMode.BOUNDARY_MIRROR)
else:
self._vgrad.add_gradients(u, v, w, tmp, tmp.size / 2, tmp.size / 2,
GradientMode.BOUNDARY_EXTRAPOLATE)
# scale with cell volume:
self._vgrad.vector_scale(voli)
[docs]
def status_check(self) -> None:
"""Checks if all domain boundaries are defined"""
all_bcs_defined = True
for bc in self._boundary_faces.get_boundary_array():
if bc.bc_type is BCTypes.IOFLOW:
for a in bc.get_face_defs():
if not a.get_status():
all_bcs_defined = False
print(' ioflow boundary patch tagged %s is undefined (inviscid wall will be used)' % a.GetTag())
if all_bcs_defined:
print(' all boundary patches defined')
[docs]
def generate_irsm_matrix(self) -> None:
"""Generates matrix for Implicit Residual Smoothing (IRSM) operator"""
g3dirsm = G3DIRSM()
g3dirsm.update_lambda()
M = MatrixRows()
for im in range(self._mesh.n_blocks):
ni, nj, nk = self._mesh.blocks[im].get_block_sizes()
mkm = 0
mkp = (ni - 1) * (nj - 1)
for k in range(nk - 1):
if k == nk - 2:
mkp = 0
mjm = 0
mjp = ni - 1
for j in range(nj - 1):
if j == nj - 2:
mjp = 0
mim = 0
mip = 1
for i in range(ni - 1):
if i == ni - 2:
mip = 0
lc = self.get_cell(im, i, j, k)
M.add_entry(lc, lc, 6.0)
M.add_entry(lc, lc + mip, -1.0)
M.add_entry(lc, lc - mim, -1.0)
M.add_entry(lc, lc + mjp, -1.0)
M.add_entry(lc, lc - mjm, -1.0)
M.add_entry(lc, lc + mkp, -1.0)
M.add_entry(lc, lc - mkm, -1.0)
mim = 1
mjm = ni - 1
mkm = (ni - 1) * (nj - 1)
for interface in self._connectivity_faces.get_interface_array():
A, B = interface.get_face_defs()
for intf in A:
fmp = intf.get_FMP
kp, kc = fmp.get_pointers()
idx1, idx2 = np.int32([fmp.size / 3, 2 * fmp.size / 3])
for i in range(idx1, idx2):
kc1 = kc[i, 1]
kc2 = kc[i, 2]
M.add_entry(kc1, kc1, 1.0)
M.add_entry(kc1, kc2, -1.0)
M.add_entry(kc2, kc2, 1.0)
M.add_entry(kc2, kc1, -1.0)
M.scale(g3dirsm.irsm_lambda)
M.add_to_diagonals(1.0)
rows, cols, coefs = M.extract_vectors()
self.__IRSM_operator.populate_sparse_matrix(rows, cols, coefs)
self.__IRSM_operator.print()