Frank14f/Frank_LBM
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

241104_vortex_lamb

Browse Source
This commit is contained in:
Frank14f
2024-11-04 18:10:36 +08:00
parent acf2b36b8c
commit 5f6337078f
54 changed files with 3292 additions and 883 deletions
Download Patch File Download Diff File Expand all files Collapse all files
CelerisLab/__pycache__/driver.cpython-310.pyc
BIN
View File
Binary file not shown.
CelerisLab/driver.py
+102 -2
View File
@@ -2,6 +2,8 @@
import pycuda.driver as cuda
import numpy as np
import struct
from scipy.special import jv, expi
from typing import List, Tuple, Union, Optional
from . import utils
@@ -13,7 +15,7 @@ SOLID = 0b00000010
GAS = 0b00000100
INTERFACE = 0b00001000
SENSOR = 0b00010000
V_TAYLOR = np.int32(1)
class FlowField:
def __init__(
@@ -94,12 +96,14 @@ class FlowField:
self.flag = np.ones(self.FIELD_SIZE, dtype=np.uint8)
self.indx = np.zeros(self.FIELD_SIZE, dtype=np.int32)
self.delta_curve = np.zeros(0, dtype=self.DATA_TYPE)
self.vortex_config = np.zeros(7, dtype=float)
self.ddf_gpu = cuda.mem_alloc(self.ddf.nbytes)
self.temp_gpu = cuda.mem_alloc(self.ddf.nbytes)
self.flag_gpu = cuda.mem_alloc(self.flag.nbytes)
self.indx_gpu = cuda.mem_alloc(self.indx.nbytes)
self.delta_gpu = cuda.mem_alloc(1)
self.vortex_gpu = cuda.mem_alloc(self.vortex_config.nbytes)
self.objects = {}
self.action = np.zeros(0, dtype=self.DATA_TYPE)
@@ -181,7 +185,7 @@ class FlowField:
self.action_gpu = cuda.mem_alloc(self.action.nbytes)
self.obs = np.zeros(len(self.objects) * self.DIM, dtype=self.DATA_TYPE)
if hasattr(self, "force_gpu"):
if hasattr(self, "obs_gpu"):
self.obs_gpu.free()
self.obs_gpu = cuda.mem_alloc(self.obs.nbytes)
@@ -236,12 +240,108 @@ class FlowField:
self.ptx = cuda.module_from_file(compiler.kernel_path("kernel.ptx"))
self.step = self.ptx.get_function("OneStep")
def add_vortex(self, center: Tuple[float, float, float], radius: float, strength: float, direction: float, type: str):
x_c, y_c, z_c = center
if (
x_c - radius <= 0
or x_c + radius >= self.FIELD_SHAPE[0] - 1
or y_c - radius <= 0
or y_c + radius >= self.FIELD_SHAPE[1] - 1
):
raise ValueError("Vortex is out of bounds.")
if type not in ["lamb", "oseen", "taylor"]:
raise ValueError("Vortex type" + type + " is not supported.")
x = np.linspace(-x_c, self.FIELD_SHAPE[0] - 1 - x_c, self.FIELD_SHAPE[0])
y = np.linspace(-y_c, self.FIELD_SHAPE[1] - 1 - y_c, self.FIELD_SHAPE[1])
X, Y = np.meshgrid(x, y)
r = np.sqrt(X**2 + Y**2)
nu = self.field_config.viscosity
theta = np.arctan2(Y, X)
psi = np.zeros_like(r)
if type == "lamb":
b = 3.831705970207512
n = b / radius
u0 = strength
inside = r <= radius
outside = r > radius
psi[inside] = (2 * u0 / n / jv(0, b) * jv(1, n * r[inside]) - u0 * r[inside]) * np.sin(theta[inside])
psi[outside] = -u0 * radius**2 / r[outside] * np.sin(theta[outside])
u_vor = np.gradient(psi, axis=0)
v_vor = -np.gradient(psi, axis=1)
p_vor = -2 * (np.gradient(v_vor, axis=1) - np.gradient(u_vor, axis=0)) * psi - (u_vor**2 + v_vor**2) / 2
elif type == "oseen":
# 4 nu t = radius^2 / 4
kappa = 2 * np.pi * radius **2 * strength
u_vor = - kappa / (2 * np.pi * r) * (1 - np.exp(-4 * r**2 / radius**2)) * np.sin(theta)
v_vor = kappa / (2 * np.pi * r) * (1 - np.exp(-4 * r**2 / radius**2)) * np.cos(theta)
zeta = 4 * r**2 / radius**2
p_vor = -kappa**2 / 8 / np.pi**2 / r**2 * (-2 * zeta * (expi(-zeta) - expi(-2 * zeta)) + (1 - np.exp(-zeta))**2)
elif type == "taylor":
# 4 nu t = radius^2
M = strength * np.pi * radius**4 / 8 / nu
u_vor = - M * r * 4 * nu / radius**4 * np.exp(-r**2 / radius**2) * np.sin(theta)
v_vor = M * r * 4 * nu / radius**4 * np.exp(-r**2 / radius**2) * np.cos(theta)
p_vor = -4 * M**2 * nu**2 * np.exp(-2 * r**2 / radius**2) / np.pi**2 / radius**6
cuda.memcpy_dtoh(self.ddf, self.ddf_gpu)
ddf_temp = self.ddf.copy().reshape((self.LATTICE, self.FIELD_SHAPE[1], self.FIELD_SHAPE[0])).transpose(2, 1, 0)
u_ddf = ddf_temp[:, :, 1] + ddf_temp[:, :, 5] + ddf_temp[:, :, 8] - ddf_temp[:, :, 3] - ddf_temp[:, :, 6] - ddf_temp[:, :, 7]
v_ddf = ddf_temp[:, :, 2] + ddf_temp[:, :, 5] + ddf_temp[:, :, 6] - ddf_temp[:, :, 4] - ddf_temp[:, :, 7] - ddf_temp[:, :, 8]
p_ddf = np.sum(ddf_temp, axis=2) / 3
for i in range(self.FIELD_SHAPE[0]):
for j in range(self.FIELD_SHAPE[1]):
k = i + j * self.FIELD_SHAPE[0]
if (j == 0 or j == self.FIELD_SHAPE[1] - 1) or (i == 0 or i == self.FIELD_SHAPE[0] - 1):
continue
else:
for e in range(self.LATTICE):
u = u_ddf[i, j] + u_vor[j, i]
v = v_ddf[i, j] + v_vor[j, i]
p = p_ddf[i, j] + p_vor[j, i]
eu = self.E[e][0] * u + self.E[e][1] * v
u2 = u ** 2 + v ** 2
self.ddf[k + e * self.FIELD_SIZE] = self.WW[e] * (3 * p + 3 * eu + 4.5 * eu ** 2 - 1.5 * u2)
cuda.memcpy_htod(self.ddf_gpu, self.ddf)
# def add_vortex_gpu(self, center: Tuple[float, float, float], radius: float, strength: float, direction: float, type: str):
# x_c, y_c, z_c = center
# if (
# x_c - radius <= 0
# or x_c + radius >= self.FIELD_SHAPE[0] - 1
# or y_c - radius <= 0
# or y_c + radius >= self.FIELD_SHAPE[1] - 1
# ):
# raise ValueError("Vortex is out of bounds.")
# if type not in ["lamb", "oseen", "taylor"]:
# raise ValueError("Vortex type" + type + " is not supported.")
# add_vortex = self.ptx.get_function("AddVortex")
# self.vortex_config[0:3] = np.array(center, dtype=float)
# self.vortex_config[3] = radius
# self.vortex_config[4] = strength
# self.vortex_config[5] = direction
# if type == "taylor":
# self.vortex_config[6] =
def run(self, num_steps: int, action_target: np.ndarray):
if (
action_target.size != len(self.objects)
or action_target.dtype != self.DATA_TYPE
):
raise ValueError("action data type or size does not match the objects.")
elif len(self.objects) == 0:
raise ValueError("No objects have been added to the flow field.")
weight = 0.1
stream = cuda.Stream()
CelerisLab/kernels/kernel.cu
+32
View File
@@ -187,4 +187,36 @@ extern "C"
f[k + i * totalCells] = f_share[threadIdx.x + i * NT];
}
}
// __global__ void AddVortex(LBtype *f, int32_t *config)
// {
// __shared__ LBtype f_share[NT * NQ];
// int x, y, k;
// LBtype u, v, u_vor, v_vor;
// Index_lattice(x, y, k);
// int totalCells = NX * NY;
// for (int i = 0; i < NQ; i++)
// {
// f_share[threadIdx.x + i * NT] = f[k + i * totalCells];
// }
// __syncthreads();
// u = f_share[threadIdx.x + 1 * NT] - f_share[threadIdx.x + 3 * NT] + f_share[threadIdx.x + 5 * NT] - f_share[threadIdx.x + 6 * NT] - f_share[threadIdx.x + 7 * NT] + f_share[threadIdx.x + 8 * NT];
// v = f_share[threadIdx.x + 2 * NT] - f_share[threadIdx.x + 4 * NT] + f_share[threadIdx.x + 5 * NT] + f_share[threadIdx.x + 6 * NT] - f_share[threadIdx.x + 7 * NT] - f_share[threadIdx.x + 8 * NT];
// if type & V_TAYLOR
// {
// u_vor = -2 * PI * U0 * sin(2 * PI * x / NX) * sin(2 * PI * y / NY);
// v_vor = 2 * PI * U0 * cos(2 * PI * x / NX) * cos(2 * PI * y / NY);
// }
// else
// {
// u_vor = 0;
// v_vor = 0;
// }
// }
}
CelerisLab/kernels/kernel.ptx
+298 -307
View File
File diff suppressed because it is too large Load Diff
CelerisLab/kernels/macros.h
+9 -6
View File
@@ -9,15 +9,15 @@
// flow parameters
#define LBtype float
#define UX 12
#define UY 20
#define UX 10
#define UY 16
#define UZ 1
#define NX 1536
#define NY 640
#define NX 1280
#define NY 512
#define NZ 1
#define DIM 2
#define NQ 9
#define VIS 0.006
#define VIS 0.004
#define RHO 1.0
#define U0 0.01
@@ -29,6 +29,9 @@
#define INTERFACE 0b00001000
#define SENSOR 0b00010000
// vortex type
#define V_TAYLOR 0b00000001
// variables
#define N_OBJS 7
#define N_OBJS 6
// #define N_SENS 2