From bd8580c242f4ad100f7dbf4df0d2c3f39ba91e23 Mon Sep 17 00:00:00 2001 From: Hesam Salehipour Date: Tue, 22 Aug 2023 10:30:52 -0400 Subject: [PATCH 1/6] modified the channel3d for validation --- examples/channel3d.py | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/examples/channel3d.py b/examples/channel3d.py index d60ba36..097ef19 100644 --- a/examples/channel3d.py +++ b/examples/channel3d.py @@ -47,7 +47,7 @@ class turbulentChannel(KBCSimForced): def set_boundary_conditions(self): # top and bottom sides of the channel are no-slip and the other directions are periodic wall = np.concatenate((self.boundingBoxIndices['bottom'], self.boundingBoxIndices['top'])) - self.BCs.append(BounceBack(tuple(wall.T), self.grid_info, self.precision_policy)) + self.BCs.append(Regularized(tuple(wall.T), self.grid_info, self.precision_policy, 'velocity', np.zeros((wall.shape[0], 3)))) return def initialize_macroscopic_fields(self): @@ -89,8 +89,11 @@ def output_data(self, **kwargs): dns_dic = get_dns_data() plt.clf() plt.semilogx(yplus, uplus,'r.', yplus, uplus_loglaw, 'k:', dns_dic['y+'], dns_dic['Umean'], 'b-') - fname = "uplus_" + str(timestep).zfill(4) + '.png' - plt.savefig(fname) + ax = plt.gca() + ax.set_xlim([0.1, 300]) + ax.set_ylim([0, 20]) + fname = "uplus_" + str(timestep//10000).zfill(5) + '.pdf' + plt.savefig(fname, format='pdf') fields = {"rho": rho, "u_x": u[..., 0], "u_y": u[..., 1], "u_z": u[..., 2]} save_fields_vtk(timestep, fields) @@ -109,7 +112,7 @@ def output_data(self, **kwargs): # define flow regime Re_tau = 180 - u_tau = 0.004 + u_tau = 0.001 DeltaPlus = Re_tau/h # DeltaPlus = u_tau / nu * Delta where u_tau / nu = Re_tau/h visc = u_tau * h / Re_tau omega = 1.0 / (3.0 * visc + 0.5) @@ -124,4 +127,4 @@ def output_data(self, **kwargs): os.system("rm -rf ./*.vtk && rm -rf ./*.png") sim = turbulentChannel(lattice, omega, nx, ny, nz, precision=precision, optimize=False) - sim.run(4000000, print_iter=20000, io_iter=20000) + sim.run(10000000, print_iter=100000, io_iter=500000) From be3b9f11add91364453b6b410c9c4181c331d69e Mon Sep 17 00:00:00 2001 From: Hesam Salehipour Date: Tue, 22 Aug 2023 10:51:00 -0400 Subject: [PATCH 2/6] modified channel3d exaple for validation --- .gitignore | 1 + examples/CFD/channel3d.py | 85 ++++++++++++++------------------------- 2 files changed, 32 insertions(+), 54 deletions(-) diff --git a/.gitignore b/.gitignore index 85e8687..9997d33 100644 --- a/.gitignore +++ b/.gitignore @@ -10,6 +10,7 @@ *zraw *mhd *png +*pdf # Exclusions !jax-lbm.png !airfoil.png diff --git a/examples/CFD/channel3d.py b/examples/CFD/channel3d.py index 521baa0..097ef19 100644 --- a/examples/CFD/channel3d.py +++ b/examples/CFD/channel3d.py @@ -1,25 +1,9 @@ -""" -This script performs a 3D simulation of turbulent channel flow using the lattice Boltzmann method (LBM). -Turbulent channel flow, also known as plane Couette flow, is a fundamental case in the study of wall-bounded turbulent flows. - -In this example you'll be introduced to the following concepts: - -1. Lattice: A D3Q27 lattice is used, which is a three-dimensional lattice model with 27 discrete velocity directions. This type of lattice allows for a more precise representation of fluid flow in three dimensions. - -2. Initial Conditions: The initial conditions for the flow are randomly generated, and the populations are initialized to be the solution of an advection-diffusion equation. - -3. Boundary Conditions: Bounce back boundary conditions are applied at the top and bottom walls, simulating a no-slip condition typical for wall-bounded flows. - -4. External Force: An external force is applied to drive the flow. - -""" - from src.boundary_conditions import * from jax.config import config from src.utils import * import numpy as np from src.lattice import LatticeD3Q27 -from src.models import KBCSim, AdvectionDiffusionBGK +from src.models import KBCSimForced, AdvectionDiffusionBGK import jax.numpy as jnp import os import matplotlib.pyplot as plt @@ -27,11 +11,14 @@ # Use 8 CPU devices # os.environ["XLA_FLAGS"] = '--xla_force_host_platform_device_count=8' import jax - +jax.config.update("jax_array", True) # disable JIt compilation - +# jax.config.update('jax_disable_jit', True) jax.config.update('jax_enable_x64', True) + +precision = "f64/f64" + def vonKarman_loglaw_wall(yplus): vonKarmanConst = 0.41 cplus = 5.5 @@ -55,39 +42,34 @@ def get_dns_data(): } return dns_dic -class turbulentChannel(KBCSim): - def __init__(self, **kwargs): - super().__init__(**kwargs) +class turbulentChannel(KBCSimForced): def set_boundary_conditions(self): # top and bottom sides of the channel are no-slip and the other directions are periodic wall = np.concatenate((self.boundingBoxIndices['bottom'], self.boundingBoxIndices['top'])) - self.BCs.append(BounceBack(tuple(wall.T), self.gridInfo, self.precisionPolicy)) + self.BCs.append(Regularized(tuple(wall.T), self.grid_info, self.precision_policy, 'velocity', np.zeros((wall.shape[0], 3)))) return def initialize_macroscopic_fields(self): - rho = self.precisionPolicy.cast_to_output(1.0) + rho = self.precision_policy.cast_to_output(1.0) u = self.distributed_array_init((self.nx, self.ny, self.nz, self.dim), - self.precisionPolicy.compute_dtype, initVal=1e-2 * np.random.random((self.nx, self.ny, self.nz, self.dim))) - u = self.precisionPolicy.cast_to_output(u) + self.precision_policy.compute_dtype, initVal=1e-2 * np.random.random((self.nx, self.ny, self.nz, self.dim))) + u = self.precision_policy.cast_to_output(u) return rho, u def initialize_populations(self, rho, u): omegaADE = 1.0 - lattice = LatticeD3Q27(precision) - - kwargs = {'lattice': lattice, 'nx': self.nx, 'ny': self.ny, 'nz': self.nz, 'precision': precision, 'omega': omegaADE, 'vel': u} - ADE = AdvectionDiffusionBGK(**kwargs) + ADE = AdvectionDiffusionBGK(u, lattice, omegaADE, self.nx, self.ny, self.nz, precision=precision) ADE.initialize_macroscopic_fields = self.initialize_macroscopic_fields print("Initializing the distribution functions using the specified macroscopic fields....") - f = ADE.run(50000) + f = ADE.run(50000, print_iter=0, io_iter=0) return f def get_force(self): # define the external force force = np.zeros((self.nx, self.ny, self.nz, 3)) force[..., 0] = Re_tau**2 * visc**2 / h**3 - return self.precisionPolicy.cast_to_output(force) + return self.precision_policy.cast_to_output(force) def output_data(self, **kwargs): rho = np.array(kwargs["rho"]) @@ -107,47 +89,42 @@ def output_data(self, **kwargs): dns_dic = get_dns_data() plt.clf() plt.semilogx(yplus, uplus,'r.', yplus, uplus_loglaw, 'k:', dns_dic['y+'], dns_dic['Umean'], 'b-') - fname = "uplus_" + str(timestep).zfill(4) + '.png' - plt.savefig(fname) - fields = {"rho": rho[..., 0], "u_x": u[..., 0], "u_y": u[..., 1], "u_z": u[..., 2]} + ax = plt.gca() + ax.set_xlim([0.1, 300]) + ax.set_ylim([0, 20]) + fname = "uplus_" + str(timestep//10000).zfill(5) + '.pdf' + plt.savefig(fname, format='pdf') + fields = {"rho": rho, "u_x": u[..., 0], "u_y": u[..., 1], "u_z": u[..., 2]} save_fields_vtk(timestep, fields) if __name__ == "__main__": - precision = "f64/f64" lattice = LatticeD3Q27(precision) + # h: channel half-width - h = 10 - # Define channel geometry based on h + h = 50 + + # define channel geometry based on h nx = 6*h ny = 3*h nz = 2*h - # Define flow regime + # define flow regime Re_tau = 180 - u_tau = 0.004 - # DeltaPlus = Re_tau/h # DeltaPlus = u_tau / nu * Delta where u_tau / nu = Re_tau/h + u_tau = 0.001 + DeltaPlus = Re_tau/h # DeltaPlus = u_tau / nu * Delta where u_tau / nu = Re_tau/h visc = u_tau * h / Re_tau omega = 1.0 / (3.0 * visc + 0.5) - # Wall distance in wall units to be used inside output_data + # wall distance in wall units to be used inside output_data zz = np.arange(nz) zz = np.minimum(zz, zz.max() - zz) yplus = zz * u_tau / visc print("omega = ", omega) + assert omega < 2.0, "omega must be less than 2.0" os.system("rm -rf ./*.vtk && rm -rf ./*.png") - kwargs = { - 'lattice': lattice, - 'omega': omega, - 'nx': nx, - 'ny': ny, - 'nz': nz, - 'precision': precision, - 'io_rate': 20000, - 'print_info_rate': 20000 - } - sim = turbulentChannel(**kwargs) - sim.run(4000000) + sim = turbulentChannel(lattice, omega, nx, ny, nz, precision=precision, optimize=False) + sim.run(10000000, print_iter=100000, io_iter=500000) From 916531a798f6f844fb08d41c2c1d263178981d38 Mon Sep 17 00:00:00 2001 From: Hesam Salehipour Date: Tue, 22 Aug 2023 11:48:41 -0400 Subject: [PATCH 3/6] modified channel3d exaple for validation --- examples/CFD/channel3d.py | 70 +++++++++++++++++++++++++++------------ 1 file changed, 49 insertions(+), 21 deletions(-) diff --git a/examples/CFD/channel3d.py b/examples/CFD/channel3d.py index 097ef19..b356ac9 100644 --- a/examples/CFD/channel3d.py +++ b/examples/CFD/channel3d.py @@ -1,9 +1,25 @@ +""" +This script performs a 3D simulation of turbulent channel flow using the lattice Boltzmann method (LBM). +Turbulent channel flow, also known as plane Couette flow, is a fundamental case in the study of wall-bounded turbulent flows. + +In this example you'll be introduced to the following concepts: + +1. Lattice: A D3Q27 lattice is used, which is a three-dimensional lattice model with 27 discrete velocity directions. This type of lattice allows for a more precise representation of fluid flow in three dimensions. + +2. Initial Conditions: The initial conditions for the flow are randomly generated, and the populations are initialized to be the solution of an advection-diffusion equation. + +3. Boundary Conditions: Bounce back boundary conditions are applied at the top and bottom walls, simulating a no-slip condition typical for wall-bounded flows. + +4. External Force: An external force is applied to drive the flow. + +""" + from src.boundary_conditions import * from jax.config import config from src.utils import * import numpy as np from src.lattice import LatticeD3Q27 -from src.models import KBCSimForced, AdvectionDiffusionBGK +from src.models import KBCSim, AdvectionDiffusionBGK import jax.numpy as jnp import os import matplotlib.pyplot as plt @@ -11,13 +27,10 @@ # Use 8 CPU devices # os.environ["XLA_FLAGS"] = '--xla_force_host_platform_device_count=8' import jax -jax.config.update("jax_array", True) -# disable JIt compilation -# jax.config.update('jax_disable_jit', True) -jax.config.update('jax_enable_x64', True) +# disable JIt compilation -precision = "f64/f64" +jax.config.update('jax_enable_x64', True) def vonKarman_loglaw_wall(yplus): vonKarmanConst = 0.41 @@ -42,34 +55,39 @@ def get_dns_data(): } return dns_dic -class turbulentChannel(KBCSimForced): +class turbulentChannel(KBCSim): + def __init__(self, **kwargs): + super().__init__(**kwargs) def set_boundary_conditions(self): # top and bottom sides of the channel are no-slip and the other directions are periodic wall = np.concatenate((self.boundingBoxIndices['bottom'], self.boundingBoxIndices['top'])) - self.BCs.append(Regularized(tuple(wall.T), self.grid_info, self.precision_policy, 'velocity', np.zeros((wall.shape[0], 3)))) + self.BCs.append(Regularized(tuple(wall.T), self.gridInfo, self.precisionPolicy, 'velocity', np.zeros((wall.shape[0], 3)))) return def initialize_macroscopic_fields(self): - rho = self.precision_policy.cast_to_output(1.0) + rho = self.precisionPolicy.cast_to_output(1.0) u = self.distributed_array_init((self.nx, self.ny, self.nz, self.dim), - self.precision_policy.compute_dtype, initVal=1e-2 * np.random.random((self.nx, self.ny, self.nz, self.dim))) - u = self.precision_policy.cast_to_output(u) + self.precisionPolicy.compute_dtype, initVal=1e-2 * np.random.random((self.nx, self.ny, self.nz, self.dim))) + u = self.precisionPolicy.cast_to_output(u) return rho, u def initialize_populations(self, rho, u): omegaADE = 1.0 - ADE = AdvectionDiffusionBGK(u, lattice, omegaADE, self.nx, self.ny, self.nz, precision=precision) + lattice = LatticeD3Q27(precision) + + kwargs = {'lattice': lattice, 'nx': self.nx, 'ny': self.ny, 'nz': self.nz, 'precision': precision, 'omega': omegaADE, 'vel': u} + ADE = AdvectionDiffusionBGK(**kwargs) ADE.initialize_macroscopic_fields = self.initialize_macroscopic_fields print("Initializing the distribution functions using the specified macroscopic fields....") - f = ADE.run(50000, print_iter=0, io_iter=0) + f = ADE.run(50000) return f def get_force(self): # define the external force force = np.zeros((self.nx, self.ny, self.nz, 3)) force[..., 0] = Re_tau**2 * visc**2 / h**3 - return self.precision_policy.cast_to_output(force) + return self.precisionPolicy.cast_to_output(force) def output_data(self, **kwargs): rho = np.array(kwargs["rho"]) @@ -94,37 +112,47 @@ def output_data(self, **kwargs): ax.set_ylim([0, 20]) fname = "uplus_" + str(timestep//10000).zfill(5) + '.pdf' plt.savefig(fname, format='pdf') - fields = {"rho": rho, "u_x": u[..., 0], "u_y": u[..., 1], "u_z": u[..., 2]} + fields = {"rho": rho[..., 0], "u_x": u[..., 0], "u_y": u[..., 1], "u_z": u[..., 2]} save_fields_vtk(timestep, fields) if __name__ == "__main__": + precision = "f64/f64" lattice = LatticeD3Q27(precision) # h: channel half-width h = 50 - # define channel geometry based on h + # Define channel geometry based on h nx = 6*h ny = 3*h nz = 2*h - # define flow regime + # Define flow regime Re_tau = 180 u_tau = 0.001 DeltaPlus = Re_tau/h # DeltaPlus = u_tau / nu * Delta where u_tau / nu = Re_tau/h visc = u_tau * h / Re_tau omega = 1.0 / (3.0 * visc + 0.5) - # wall distance in wall units to be used inside output_data + # Wall distance in wall units to be used inside output_data zz = np.arange(nz) zz = np.minimum(zz, zz.max() - zz) yplus = zz * u_tau / visc print("omega = ", omega) - assert omega < 2.0, "omega must be less than 2.0" os.system("rm -rf ./*.vtk && rm -rf ./*.png") - sim = turbulentChannel(lattice, omega, nx, ny, nz, precision=precision, optimize=False) - sim.run(10000000, print_iter=100000, io_iter=500000) + kwargs = { + 'lattice': lattice, + 'omega': omega, + 'nx': nx, + 'ny': ny, + 'nz': nz, + 'precision': precision, + 'io_rate': 500000, + 'print_info_rate': 100000 + } + sim = turbulentChannel(**kwargs) + sim.run(10000000) From 17a314b58b283661df69c558716d815e85b945aa Mon Sep 17 00:00:00 2001 From: Hesam Salehipour Date: Tue, 19 Sep 2023 13:52:43 -0400 Subject: [PATCH 4/6] improved the t-g example to showcase the error with various precision/computation pairs --- examples/CFD/taylor_green_vortex.py | 109 ++++++++++++++++------------ 1 file changed, 62 insertions(+), 47 deletions(-) diff --git a/examples/CFD/taylor_green_vortex.py b/examples/CFD/taylor_green_vortex.py index 9434d68..1071025 100644 --- a/examples/CFD/taylor_green_vortex.py +++ b/examples/CFD/taylor_green_vortex.py @@ -6,16 +6,13 @@ from src.boundary_conditions import * -from jax.config import config from src.utils import * import numpy as np from src.lattice import LatticeD2Q9 from src.models import BGKSim, KBCSim, AdvectionDiffusionBGK -import jax.numpy as jnp -from jax.experimental import mesh_utils -from jax.sharding import PositionalSharding import os import matplotlib.pyplot as plt +import json # Use 8 CPU devices # os.environ["XLA_FLAGS"] = '--xla_force_host_platform_device_count=8' @@ -25,8 +22,8 @@ jax.config.update('jax_enable_x64', True) def taylor_green_initial_fields(xx, yy, u0, rho0, nu, time): - ux = -u0 * np.cos(xx) * np.sin(yy) * np.exp(-2 * nu * time) - uy = u0 * np.sin(xx) * np.cos(yy) * np.exp(-2 * nu * time) + ux = u0 * np.sin(xx) * np.cos(yy) * np.exp(-2 * nu * time) + uy = -u0 * np.cos(xx) * np.sin(yy) * np.exp(-2 * nu * time) rho = 1.0 - rho0 * u0 ** 2 / 12. * (np.cos(2. * xx) + np.cos(2. * yy)) * np.exp(-4 * nu * time) return ux, uy, np.expand_dims(rho, axis=-1) @@ -51,7 +48,7 @@ def initialize_populations(self, rho, u): ADE = AdvectionDiffusionBGK(**kwargs) ADE.initialize_macroscopic_fields = self.initialize_macroscopic_fields print("Initializing the distribution functions using the specified macroscopic fields....") - f = ADE.run(20000) + f = ADE.run(int(20000*32/nx)) return f def output_data(self, **kwargs): @@ -64,49 +61,67 @@ def output_data(self, **kwargs): time = timestep * (kx**2 + ky**2)/2. ux_th, uy_th, rho_th = taylor_green_initial_fields(xx, yy, vel_ref, 1, visc, time) vel_err_L2 = np.sqrt(np.sum((u[..., 0]-ux_th)**2 + (u[..., 1]-uy_th)**2) / np.sum(ux_th**2 + uy_th**2)) - print("error= {:07.6f}".format(vel_err_L2)) + rho_err_L2 = np.sqrt(np.sum((rho - rho_th)**2) / np.sum(rho_th**2)) + print("Vel error= {:07.6f}, Pressure error= {:07.6f}".format(vel_err_L2, rho_err_L2)) if timestep == endTime: ErrL2ResList.append(vel_err_L2) + ErrL2ResListRho.append(rho_err_L2) # save_image(timestep, u) if __name__ == "__main__": - precision = "f64/f64" - lattice = LatticeD2Q9(precision) - - resList = [32, 64, 128, 256, 512] - ErrL2ResList = [] - - for nx in resList: - print("Running at nx = ny = {:07.6f}".format(nx)) - ny = nx - twopi = 2.0 * np.pi - coord = np.array([(i, j) for i in range(nx) for j in range(ny)]) - xx, yy = coord[:, 0], coord[:, 1] - kx, ky = twopi / nx, twopi / ny - xx = xx.reshape((nx, ny)) * kx - yy = yy.reshape((nx, ny)) * ky - - Re = 1000.0 - vel_ref = 0.04*32/nx - - visc = vel_ref * nx / Re - omega = 1.0 / (3.0 * visc + 0.5) - print("omega = ", omega) - os.system("rm -rf ./*.vtk && rm -rf ./*.png") - kwargs = { - 'lattice': lattice, - 'omega': omega, - 'nx': nx, - 'ny': ny, - 'nz': 0, - 'precision': precision, - 'io_rate': 500, - 'print_info_rate': 500 - } - sim = TaylorGreenVortex(**kwargs) - endTime = int(20000*nx/32.0) - sim.run(endTime) - plt.loglog(resList, ErrL2ResList, '-o') - plt.loglog(resList, 1e-3*(np.array(resList)/128)**(-2), '--') - plt.savefig('Error.png'); plt.savefig('Error.pdf', format='pdf') + precision_list = ["f32/f32", "f64/f32", "f64/f64"] + resList = [32, 64, 128, 256, 512, 1024] + result_dict = dict.fromkeys(precision_list) + result_dict['resolution_list'] = resList + + for precision in precision_list: + lattice = LatticeD2Q9(precision) + ErrL2ResList = [] + ErrL2ResListRho = [] + result_dict[precision] = dict.fromkeys(['vel_error', 'rho_error']) + for nx in resList: + print("Running at nx = ny = {:07.6f}".format(nx)) + ny = nx + twopi = 2.0 * np.pi + coord = np.array([(i, j) for i in range(nx) for j in range(ny)]) + xx, yy = coord[:, 0], coord[:, 1] + kx, ky = twopi / nx, twopi / ny + xx = xx.reshape((nx, ny)) * kx + yy = yy.reshape((nx, ny)) * ky + + Re = 1600.0 + vel_ref = 0.04*32/nx + + visc = vel_ref * nx / Re + omega = 1.0 / (3.0 * visc + 0.5) + print("omega = ", omega) + os.system("rm -rf ./*.vtk && rm -rf ./*.png") + kwargs = { + 'lattice': lattice, + 'omega': omega, + 'nx': nx, + 'ny': ny, + 'nz': 0, + 'precision': precision, + 'io_rate': 5000, + 'print_info_rate': 1000 + } + sim = TaylorGreenVortex(**kwargs) + tc = 2.0/(2. * visc * (kx**2 + ky**2)) + endTime = int(0.05*tc) + sim.run(endTime) + result_dict[precision]['vel_error'] = ErrL2ResList + result_dict[precision]['rho_error'] = ErrL2ResListRho + + with open('data.json', 'w') as fp: + json.dump(result_dict, fp) + + # plt.loglog(resList, ErrL2ResList, '-o') + # plt.loglog(resList, 1e-3*(np.array(resList)/128)**(-2), '--') + # plt.savefig('ErrorVel.png'); plt.savefig('ErrorVel.pdf', format='pdf') + + # plt.figure() + # plt.loglog(resList, ErrL2ResListRho, '-o') + # plt.loglog(resList, 1e-3*(np.array(resList)/128)**(-2), '--') + # plt.savefig('ErrorRho.png'); plt.savefig('ErrorRho.pdf', format='pdf') From 088e7647692ac66db332324637cf74d5c5495f33 Mon Sep 17 00:00:00 2001 From: Hesam Salehipour Date: Tue, 19 Sep 2023 15:58:34 -0400 Subject: [PATCH 5/6] fixed a bug in moving bounce-back and added velocity to the halfway bounceback --- src/boundary_conditions.py | 18 ++++++++++++------ 1 file changed, 12 insertions(+), 6 deletions(-) diff --git a/src/boundary_conditions.py b/src/boundary_conditions.py index baad633..c33c556 100644 --- a/src/boundary_conditions.py +++ b/src/boundary_conditions.py @@ -445,8 +445,8 @@ def apply(self, fout, fin, time): """ indices, vel = self.update_function(time) c = jnp.array(self.lattice.c, dtype=self.precisionPolicy.compute_dtype) - cu = 3.0 * jnp.dot(vel, c) - return fout.at[indices].set(fin[indices][..., self.lattice.opp_indices] - 6.0 * self.lattice.w * cu) + cu = 6.0 * self.lattice.w * jnp.dot(vel, c) + return fout.at[indices].set(fin[indices][..., self.lattice.opp_indices] - cu) class BounceBackHalfway(BoundaryCondition): @@ -468,12 +468,13 @@ class BounceBackHalfway(BoundaryCondition): isSolid : bool Whether the boundary condition represents a solid boundary. For this class, it is True. """ - def __init__(self, indices, gridInfo, precision_policy): + def __init__(self, indices, gridInfo, precision_policy, vel=None): super().__init__(indices, gridInfo, precision_policy) self.name = "BounceBackHalfway" self.implementationStep = "PostStreaming" self.needsExtraConfiguration = True self.isSolid = True + self.vel = vel def configure(self, boundaryBitmask): """ @@ -523,7 +524,12 @@ def apply(self, fout, fin): nbd = len(self.indices[0]) bindex = np.arange(nbd)[:, None] fbd = fout[self.indices] - fbd = fbd.at[bindex, self.imissing].set(fin[self.indices][bindex, self.iknown]) + if self.vel is not None: + c = jnp.array(self.lattice.c, dtype=self.precisionPolicy.compute_dtype) + cu = 6.0 * self.lattice.w * jnp.dot(self.vel, c) + fbd = fbd.at[bindex, self.imissing].set(fin[self.indices][bindex, self.iknown] - cu[bindex, self.iknown]) + else: + fbd = fbd.at[bindex, self.imissing].set(fin[self.indices][bindex, self.iknown]) return fbd @@ -926,8 +932,8 @@ def configure(self, boundaryBitmask): Parameters ---------- - connectivity_bitmask : np.ndarray - The connectivity bitmask for the lattice. + boundaryBitmask : np.ndarray + The connectivity bitmask for the boundary voxels. """ shiftDir = ~boundaryBitmask[:, self.lattice.opp_indices] idx = np.array(self.indices).T From 21268fbd4229e84f3a645cd15d4017dd5472f97b Mon Sep 17 00:00:00 2001 From: Hesam Salehipour Date: Tue, 19 Sep 2023 18:09:26 -0400 Subject: [PATCH 6/6] added the description for the velocity attribute of the halfway bounceback --- src/boundary_conditions.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/boundary_conditions.py b/src/boundary_conditions.py index c33c556..73cf55d 100644 --- a/src/boundary_conditions.py +++ b/src/boundary_conditions.py @@ -467,6 +467,8 @@ class BounceBackHalfway(BoundaryCondition): Whether the boundary condition needs extra configuration before it can be applied. For this class, it is True. isSolid : bool Whether the boundary condition represents a solid boundary. For this class, it is True. + vel : array-like + The prescribed value of velocity vector for the boundary condition. No-slip BC is assumed if vel=None (default). """ def __init__(self, indices, gridInfo, precision_policy, vel=None): super().__init__(indices, gridInfo, precision_policy)