Source code for g3d.flow.geo.q1dpipe
"""Q1D Pipe placeholder (Q1DPipe class definition)"""
from typing import Any, Union, 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 Q1DPipeType
from g3d.defs import BCTypes
from g3d.defs import G3DEntity
from g3d.defs import FunctionByPass
[docs]
class Q1DPipe(G3DEntity):
"""Definition of class Q1DPipe: quasi-one-dimensional pipe definition container"""
def __init__(self, *args: Any, **kwargs: Any) -> None:
self._pipe_type = Q1DPipeType.CONSTANT_AREA_PIPE
self._left_end_BC = BCTypes.IOFLOW
self._right_end_BC = BCTypes.IOFLOW
self._x1 = 0.
self._xt = 0.33
self._xl = 0.5
self._xi = 0.5
self._x2 = 1.0
self.__L = self._x2 - self._x1
self._r1 = np.sqrt(0.30 / np.pi)
self._rt = np.sqrt(0.25 / np.pi)
self._r2 = np.sqrt(0.50 / np.pi)
self._NC1 = 0
self._NC2 = 0
self._NC3 = 0
self._n_cells = self._NC1 + self._NC2 + self._NC3
self.__n_nodes = self._n_cells + 1
self._delta = np.zeros(self._n_cells)
self._X = np.zeros(self.__n_nodes)
self._R = np.zeros(self.__n_nodes)
self._h1 = -1
self._h2 = -1
self._ht = -1
self._hl = -1
self._hi = -1
self.reset()
super().__init__(*args, **kwargs)
def _set_default(self) -> None:
super()._set_default()
self._tag = 'q1d_pipe'
def _set_from_args(self, *args: Any, **kwargs: Any) -> ndarray:
bps = self._get_function_evaluate_status(**kwargs)
if bps is FunctionByPass.SELF:
args_stat = np.ones(len(args))
else:
args_stat = super()._set_from_args(*args)
if bps in [FunctionByPass.EVALUATE, FunctionByPass.SELF]:
if np.sum(args_stat) > 0:
for i in range(len(args)):
if args_stat[i] == 1:
arg = args[i]
if type(arg) is Q1DPipeType:
args_stat[i] = 0
self._pipe_type = arg
return args_stat
[docs]
def set_from_key_value(self, key: str, value: Any) -> bool:
"""Sets data from key-value pair input
Parameters
----------
key: str :
keyword
value: Any :
value corresponding to keyword
Returns
-------
bool
status - True if key is found in list of valid keywords
"""
if super().set_from_key_value(key, value):
return True
elif key == 'pipe_type':
self._pipe_type = value
return True
elif key == 'left':
self._left_end_BC = value
return True
elif key == 'right':
self._right_end_BC = value
return True
elif key == 'x1':
self._x1 = value
self.__L = self._x2 - self._x1
return True
elif key == 'x2':
self._x2 = value
self.__L = self._x2 - self._x1
return True
elif key == 'xt':
self._xt = value
return True
elif key == 'xl':
self._xl = value
return True
elif key == 'xi':
self._xi = value
return True
elif key == 'r1':
self._r1 = value
return True
elif key == 'rt':
self._rt = value
return True
elif key == 'r2':
self._r2 = value
return True
elif key == 'NC1':
self._NC1 = value
self.update_cell_count()
return True
elif key == 'NC2':
self._NC2 = value
self.update_cell_count()
return True
elif key == 'NC3':
self._NC3 = value
self.update_cell_count()
return True
elif key == 'h1':
self._h1 = value
return True
elif key == 'ht':
self._ht = value
return True
elif key == 'h2':
self._h2 = value
return True
elif key == 'hl':
self._hl = value
return True
elif key == 'hi':
self._hi = value
return True
elif key == 'n_cells':
self.set_cell_count(value)
return True
return False
[docs]
def reset(self) -> None:
"""Resets pipe definition to default condition"""
self._left_end_BC = BCTypes.IOFLOW
self._right_end_BC = BCTypes.IOFLOW
self.reset_cell_counts()
self.update_cell_count()
self._x1 = 0.
self._xt = 0.33
self._xl = 0.5
self._xi = 0.5
self._x2 = 1.0
self.update_pipe_length()
self._r1 = np.sqrt(0.30 / np.pi)
self._rt = np.sqrt(0.25 / np.pi)
self._r2 = np.sqrt(0.50 / np.pi)
self.reset_spacing()
self.reset_placeholders()
[docs]
def reset_cell_counts(self) -> None:
"""Sets all cell counts to zero"""
self._NC1 = 0
self._NC2 = 0
self._NC3 = 0
self.update_cell_count()
[docs]
def update_cell_count(self) -> None:
"""Update cell counts from pipe definition"""
self._n_cells = self._NC1 + self._NC2 + self._NC3
self.__n_nodes = self._n_cells + 1
[docs]
def set_cell_count(self, N) -> None:
"""Set total cell count
Parameters
----------
N : int
number of cells
"""
self._n_cells = N
self.__n_nodes = N + 1
[docs]
def reset_spacing(self) -> None:
"""Resets all spacing parameters (inlet, outlet, throat, interface, local refinement)"""
self._h1 = -1
self._h2 = -1
self._ht = -1
self._hl = -1
self._hi = -1
[docs]
def reset_placeholders(self) -> None:
"""Reset array placeholders (x, R, delta)"""
self._delta = np.zeros(self._n_cells)
self._X = np.zeros(self.__n_nodes)
self._R = np.zeros(self.__n_nodes)
[docs]
def get_arrays(self) -> Tuple[ndarray, ndarray]:
"""Get arrays (X, R)
Returns
-------
tuple
axial coordinate (X) and radius (R) arrays
"""
return self._X, self._R
@property
def pipe_type(self) -> Any:
"""Pipe type ('CONSTANT_AREA_PIPE', 'TWO_REGION_PIPE', 'CONVERGENT_OR_DIVERGENT_PIPE',
'CONVERGENT_DIVERGENT_PIPE')
Returns
-------
Q1DPipeType
defined pipe type
"""
return self._pipe_type
@pipe_type.setter
def pipe_type(self, val: Any) -> None:
if type(val) is Q1DPipeType:
self._pipe_type = val
@property
def left_end_BC(self) -> Any:
"""Left-end boundary condition
Returns
-------
BCTypes
defined left-end boundary condition
"""
return self._left_end_BC
@left_end_BC.setter
def left_end_BC(self, val: Any) -> None:
if type(val) is BCTypes:
self._left_end_BC = val
@property
def right_end_BC(self) -> Any:
"""Right-end boundary condition
Returns
-------
BCTypes
defined right-end boundary condition
"""
return self._right_end_BC
@right_end_BC.setter
def right_end_BC(self, val: Any) -> None:
if type(val) is BCTypes:
self._right_end_BC = val
# AXIAL COORDINATES
@property
def x1(self) -> float:
"""Domain inlet coordinate (x1)
Returns
-------
float
defined inlet coordinate value
"""
return self._x1
@x1.setter
def x1(self, val: float) -> None:
self._x1 = val
@property
def x2(self) -> float:
"""Domain outlet coordinate (x2)
Returns
-------
float
defined outlet coordinate value
"""
return self._x2
@x2.setter
def x2(self, val: float) -> None:
self._x2 = val
@property
def xt(self) -> float:
"""Throat coordinate (xt)
Returns
-------
float
defined throat coordinate value
"""
return self._xt
@xt.setter
def xt(self, val: float) -> None:
self._xt = val
@property
def xl(self) -> float:
"""Local refinement coordinate (xl)
Returns
-------
float
defined local refinement coordinate value
"""
return self._xl
@xl.setter
def xl(self, val: float) -> None:
self._xl = val
@property
def xi(self) -> float:
"""Interface coordinate (xi)
Returns
-------
float
defined interface coordinate value
"""
return self._xi
@xi.setter
def xi(self, val: float) -> None:
self._xi = val
# PIPE RADIUS
@property
def r1(self) -> float:
"""Domain inlet radius (r1)
Returns
-------
float
defined inlet radius value
"""
return self._r1
@r1.setter
def r1(self, val: float) -> None:
self._r1 = val
@property
def r2(self) -> float:
"""Domain outlet radius (r2)
Returns
-------
float
defined outlet radius value
"""
return self._r2
@r2.setter
def r2(self, val: float) -> None:
self._r2 = val
@property
def rt(self) -> float:
"""Throat radius (rt)
Returns
-------
float
defined throat radius value
"""
return self._rt
@rt.setter
def rt(self, val: float) -> None:
self._rt = val
# CELL WIDTH
@property
def h1(self) -> Union[float, int]:
"""Domain inlet spacing (h1)
Returns
-------
float
defined inlet spacing value
"""
return self._h1
@h1.setter
def h1(self, val: Union[float, int]) -> None:
self._h1 = val
@property
def h2(self) -> Union[float, int]:
"""Domain outlet spacing (h2)
Returns
-------
float
defined outlet spacing value
"""
return self._h2
@h2.setter
def h2(self, val: Union[float, int]) -> None:
self._h2 = val
@property
def ht(self) -> Union[float, int]:
"""Throat spacing (ht)
Returns
-------
float
defined throat spacing value
"""
return self._ht
@ht.setter
def ht(self, val: Union[float, int]) -> None:
self._ht = val
@property
def hl(self) -> Union[float, int]:
"""Local refinement spacing (hl)
Returns
-------
float
defined local refinement spacing value
"""
return self._hl
@hl.setter
def hl(self, val: Union[float, int]) -> None:
self._hl = val
@property
def hi(self) -> Union[float, int]:
"""Interface spacing (hl)
Returns
-------
float
defined interface spacing value
"""
return self._hi
@hi.setter
def hi(self, val: Union[float, int]) -> None:
self._hi = val
# REGION CELL COUNTS
@property
def NC1(self) -> int:
"""Number of cells in region 1 - from inlet to throat or from inlet to interface
Returns
-------
int
cell count (region 1)
"""
return self._NC1
@NC1.setter
def NC1(self, N: int) -> None:
self._NC1 = N
self.update_cell_count()
@property
def NC2(self) -> int:
"""Number of cells in region 2 - from interface/throat to local refinement or outlet
Returns
-------
int
cell count (region 2)
"""
return self._NC2
@NC2.setter
def NC2(self, N: int) -> None:
self._NC2 = N
self.update_cell_count()
@property
def NC3(self) -> int:
"""Number of cells in region 3 - from local refinement to outlet
Returns
-------
int
cell count (region 3)
"""
return self._NC3
@NC3.setter
def NC3(self, N: int) -> None:
self._NC3 = N
self.update_cell_count()
@property
def n_cells(self) -> int:
"""Total number of cells
Returns
-------
int
total cell count
"""
return self._n_cells
@n_cells.setter
def n_cells(self, N: int) -> None:
self.set_cell_count(N)
# SETTINGS COPY FUNCTIONALITY
[docs]
def copy(self, p: 'Q1DPipe') -> None:
"""Copies data from prototype Q1DPipe object
Parameters
----------
p: Q1DPipe :
prototype quasi-one-dimensional pipe (Q1DPipe) object
"""
super().copy(p)
self._pipe_type = p.pipe_type
self._NC1 = p.NC1
self._NC2 = p.NC2
self._NC3 = p.NC3
self._x1 = p.x1
self._xt = p.xt
self._xl = p.xl
self._xi = p.xi
self._x2 = p.x2
self._r1 = p.r1
self._rt = p.rt
self._r2 = p.r2
self._h1 = p.h1
self._h2 = p.h2
self._ht = p.ht
self._hl = p.hl
self._hi = p.hi
self._delta = []
# CELL DISTRIBUTION FUNCTIONALITY
[docs]
def update_throat_coordinate(self) -> None:
"""Calculate throat coordinate unless specified (length / 3)"""
self.update_pipe_length()
self._xt = self.__L / 3. + self._x1
[docs]
def update_throat_coordinate_uniform(self) -> None:
"""Calculate throat coordinate based on equidistant cells"""
self.update_pipe_length()
dx = self.__L / self._n_cells
self._NC1 = np.floor(self._n_cells / 3)
self._xt = self._x1 + self._NC1 * dx
[docs]
def update_pipe_length(self) -> None:
"""Calculate total domain length"""
self.__L = self._x2 - self._x1
def __one_sided_stretching(self,
x1: float,
x2: float,
h: Union[float, int],
NC: int,
idx0: int = 0,
sign: Union[float, int] = 1,
alpha: float = 0.5) -> None: # CURV_1
float_n_cells = np.float64(NC)
x_prev = x1
for i in range(1, NC + 1):
xsi = np.float64(i) / float_n_cells
GF = (np.exp(alpha * xsi) - 1.0 - alpha * xsi) / (np.exp(alpha) - 1.0 - alpha)
G10 = 1.0 - GF
G11 = xsi - GF
G2 = GF
dx = G10 * x1 + G11 * h * float_n_cells + G2 * x2 - x_prev
x_prev += dx
self._delta[np.int32(sign * i - 0.5 * (sign + 1) + idx0)] = dx
def __two_sided_stretching(self,
x1: float,
x2: float,
h1: Union[float, int],
h2: Union[float, int],
NC: int, idx0: int = 0,
alpha: float = 0.5) -> None: # CURV_2
float_n_cells = np.float64(NC)
x_prev = x1
L = x2 - x1
DT1 = h1 / L * float_n_cells
DT2 = h2 / L * float_n_cells
for i in range(1, NC + 1):
xsi = np.float64(i) / float_n_cells
if alpha != 0.:
if xsi < 0.5:
GF = 0.5 * (np.exp(2.0 * alpha * xsi) - 1.0 - 2.0 * alpha * xsi) / \
(np.exp(alpha) - 1.0 - alpha)
else:
GF = 1.0 - 0.5 * (np.exp(2.0 * alpha * (1.0 - xsi)) - 1.0 - 2.0 * alpha * (1.0 - xsi)) / \
(np.exp(alpha) - 1.0 - alpha)
# G10 = 1.0 - GF
G20 = GF
G11 = xsi * (1.0 - GF)
G21 = (xsi - 1.0) * GF
T = G11 * DT1 + G20 + G21 * DT2
else:
# G10 = 1.0 - 3.0 * xsi**2.0 + 2.0 * xsi**3.0
G11 = xsi - 2.0 * xsi ** 2.0 + xsi ** 3.0
G20 = 3.0 * xsi ** 2.0 - 2.0 * xsi ** 3.0
G21 = -xsi ** 2.0 + xsi ** 3.0
T = G11 * DT1 + G20 + G21 * DT2
dx = (1.0 - T) * x1 + T * x2 - x_prev
x_prev += dx
self._delta[i - 1 + idx0] = dx
[docs]
def extract_coordinates(self) -> None:
"""Calculate node coordinates based on domain definition"""
self.update_pipe_length()
dx = self.__L / self._n_cells
self.reset_placeholders()
# pipe_type == 0 : constant-area pipe
# pipe_type == 1 : constant-area pipe with two flow regions (shock tube)
# pipe_type == 2 : convergent or divergent pipe (varying cross-section area)
# pipe_type == 3 : convergent-divergent nozzle (automatic setup)
# pipe_type == 4 : convergent-divergent nozzle (fixed throat location, specified spacing)
# pipe_type == 5 : convergent-divergent nozzle (fixed throat location and a
# specified estimated shock location, specified spacing)
# Q1DPipeType = Enum('Q1DPipeType', [
# 'CONSTANT_AREA_PIPE',
# 'TWO_REGION_PIPE',
# 'CONVERGENT_OR_DIVERGENT_PIPE',
# 'CONVERGENT_DIVERGENT_PIPE'])
use_cell_spacing_control = (self._h1 != -1 or self._h2 != -1)
if self._pipe_type == Q1DPipeType.CONVERGENT_DIVERGENT_PIPE and self._ht != -1:
use_cell_spacing_control = True
if self._pipe_type == Q1DPipeType.TWO_REGION_PIPE and self._hi != -1:
use_cell_spacing_control = True
if self._hl != -1:
use_cell_spacing_control = True
if (self._pipe_type == Q1DPipeType.CONSTANT_AREA_PIPE or
self._pipe_type == Q1DPipeType.CONVERGENT_OR_DIVERGENT_PIPE or
self._pipe_type == Q1DPipeType.CONVERGENT_DIVERGENT_PIPE):
if not use_cell_spacing_control:
self._delta[:] = dx
if self._pipe_type == Q1DPipeType.CONVERGENT_DIVERGENT_PIPE:
self._NC1 = np.int32(np.floor(self._n_cells * (self._xt - self._x1) / self.__L))
self._NC2 = self._n_cells - self._NC1
self.__two_sided_stretching(self._x1, self._xt, dx, dx, self._NC1)
self.__two_sided_stretching(self._xt, self._x2, dx, dx, self._NC2, self._NC1)
else: # spacing control activated
control_count = 0
if self._h2 != -1:
control_count += 1
if self._hl != -1:
control_count += 1
if self._pipe_type == Q1DPipeType.CONVERGENT_DIVERGENT_PIPE:
if self._ht != -1:
control_count += 1
if self._xl <= self._xt and self._hl != -1: # local mesh refinement in divergent section only
self._hl = -1
control_count -= 1
# distribute cells
if control_count == 1:
if self._ht != -1:
self.__one_sided_stretching(self._x1, self._xt, self._ht, self._NC1, self._NC1, -1)
self.__one_sided_stretching(self._xt, self._x2, self._ht, self._NC2, self._NC1)
else:
for i in range(self._NC1):
self._delta[i] = dx
if self._hl != -1:
self.__two_sided_stretching(self._xt, self._xl, dx, self._hl, self._NC2, self._NC1)
self.__one_sided_stretching(self._xl, self._x2, self._hl, self._NC3,
self._NC1 + self._NC2)
if self._h2 != -1:
self.__two_sided_stretching(self._xt, self._x2, dx, self._h2, self._NC2, self._NC1)
if control_count == 2:
if self._ht != -1 and self._h2 != -1:
self.__one_sided_stretching(self._x1, self._xt, self._ht, self._NC1, self._NC1, -1)
self.__two_sided_stretching(self._xt, self._x2, self._ht, self._h2, self._NC2, self._NC1)
if self._hl != -1:
if self._ht != -1:
self.__one_sided_stretching(self._x1, self._xt, self._ht, self._NC1, self._NC1, -1)
self.__two_sided_stretching(self._xt, self._xl, self._ht, self._hl, self._NC2,
self._NC1)
self.__one_sided_stretching(self._xl, self._x2, self._hl, self._NC3,
self._NC1 + self._NC2)
if self._h2 != -1:
for i in range(self._NC1):
self._delta[i] = dx
self.__two_sided_stretching(self._xt, self._xl, dx, self._hl, self._NC2, self._NC1)
self.__two_sided_stretching(self._xl, self._x2, self._hl,
self._h2, self._NC3, self._NC1 + self._NC2)
if control_count == 3:
self.__one_sided_stretching(self._x1, self._xt, self._ht, self._NC1, self._NC1, -1)
self.__two_sided_stretching(self._xt, self._xl, self._ht, self._hl, self._NC2, self._NC1)
self.__two_sided_stretching(self._xl, self._x2, self._hl, self._h2, self._NC3,
self._NC1 + self._NC2)
else: # pipe_type == 0 || pipe_type == 2
if self._h1 != -1:
control_count += 1
NC = self._n_cells
if self._hl != -1:
if self._xl >= self._x2 or self._xl <= self._x1:
self._hl = -1
control_count -= 1
if control_count == 1:
if self._h1 != -1:
self.__one_sided_stretching(self._x1, self._x2, self._h1, NC)
if self._h2 != -1:
self.__one_sided_stretching(self._x1, self._x2, self._h2, NC, NC, -1)
if self._hl != -1:
self.__one_sided_stretching(self._x1, self._xl, self._hl, self._NC1, self._NC1, -1)
self.__one_sided_stretching(self._xl, self._x2, self._hl, self._NC2, self._NC1)
if control_count == 2:
if self._h1 != -1 and self._h2 != -1:
self.__two_sided_stretching(self._x1, self._x2, self._h1, self._h2, NC)
if self._h1 != -1 and self._hl != -1:
self.__two_sided_stretching(self._x1, self._xl, self._h1, self._hl, self._NC1)
self.__one_sided_stretching(self._xl, self._x2, self._hl, self._NC2, self._NC1)
if self._h2 != -1 and self._hl != -1:
self.__one_sided_stretching(self._x1, self._xl, self._hl, self._NC1, self._NC1, -1)
self.__two_sided_stretching(self._xl, self._x2, self._hl, self._h2, self._NC2, self._NC1)
if control_count == 3:
self.__two_sided_stretching(self._x1, self._xl, self._h1, self._hl, self._NC1)
self.__two_sided_stretching(self._xl, self._x2, self._hl, self._h2, self._NC2, self._NC1)
elif self._pipe_type == Q1DPipeType.TWO_REGION_PIPE:
if self._hi != -1:
self.__one_sided_stretching(self._x1, self._xi, self._hi, self._NC1, self._NC1, -1)
self.__one_sided_stretching(self._xi, self._x2, self._hi, self._NC2, self._NC1)
else:
self.__two_sided_stretching(self._x1, self._xi, dx, dx, self._NC1)
self.__two_sided_stretching(self._xi, self._x2, dx, dx, self._NC2, self._NC1)
# get axial coordinates for all nodes
self._X[0] = self._x1
for i in range(1, self.__n_nodes):
self._X[i] = self._X[i - 1] + self._delta[i - 1]
self._X[-1] = self._x2
# calculate radial coordinates for all nodes
r = self._r1
# x = self._x1
for i in range(self.__n_nodes):
x = self._X[i]
if self._pipe_type != Q1DPipeType.CONSTANT_AREA_PIPE and self._pipe_type != Q1DPipeType.TWO_REGION_PIPE:
if self._pipe_type == Q1DPipeType.CONVERGENT_OR_DIVERGENT_PIPE:
r = 0.5 * (self._r1 + self._r2) + 0.5 * (self._r1 - self._r2) * np.cos(
np.pi * x / (self._x2 - self._x1))
else:
if x <= self._xt:
r = 0.5 * (self._r1 + self._rt) + \
0.5 * (self._r1 - self._rt) * np.cos(np.pi * x / (self._xt - self._x1))
else:
r = 0.5 * (self._rt + self._r2) + \
0.5 * (self._r2 - self._rt) * np.cos(np.pi + np.pi / (self._x2 - self._xt) * (x - self._xt))
self._R[i] = r