-
Notifications
You must be signed in to change notification settings - Fork 25
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Merge pull request #63 from hsalehipour/dev
Added external force and fixed dtype issues
- Loading branch information
Showing
46 changed files
with
893 additions
and
343 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1 @@ | ||
| {"y": [0, 0.000301, 0.0012, 0.00271, 0.00482, 0.00752, 0.0108, 0.0147, 0.0192, 0.0243, 0.03, 0.0362, 0.0431, 0.0505, 0.0585, 0.067, 0.0761, 0.0858, 0.096, 0.107, 0.118, 0.13, 0.142, 0.155, 0.169, 0.182, 0.197, 0.212, 0.227, 0.243, 0.259, 0.276, 0.293, 0.31, 0.328, 0.347, 0.366, 0.385, 0.404, 0.424, 0.444, 0.465, 0.486, 0.507, 0.529, 0.55, 0.572, 0.595, 0.617, 0.64, 0.663, 0.686, 0.71, 0.733, 0.757, 0.781, 0.805, 0.829, 0.853, 0.878, 0.902, 0.926, 0.951, 0.975, 1], "y+": [0, 0.053648, 0.21456, 0.48263, 0.85771, 1.3396, 1.9279, 2.6224, 3.4226, 4.328, 5.3381, 6.4523, 7.67, 8.9902, 10.412, 11.936, 13.559, 15.281, 17.102, 19.019, 21.033, 23.141, 25.342, 27.635, 30.019, 32.492, 35.053, 37.701, 40.432, 43.247, 46.143, 49.118, 52.171, 55.3, 58.503, 61.778, 65.123, 68.536, 72.016, 75.559, 79.164, 82.828, 86.55, 90.327, 94.157, 98.037, 101.97, 105.94, 109.96, 114.02, 118.12, 122.25, 126.42, 130.62, 134.84, 139.1, 143.37, 147.67, 151.99, 156.32, 160.66, 165.02, 169.38, 173.75, 178.12], "Umean": [0, 0.053639, 0.21443, 0.48197, 0.85555, 1.3339, 1.9148, 2.5939, 3.3632, 4.2095, 5.1133, 6.0493, 6.9892, 7.9052, 8.7741, 9.579, 10.311, 10.967, 11.55, 12.066, 12.52, 12.921, 13.276, 13.59, 13.87, 14.121, 14.349, 14.557, 14.75, 14.931, 15.101, 15.264, 15.419, 15.569, 15.714, 15.855, 15.993, 16.128, 16.26, 16.389, 16.515, 16.637, 16.756, 16.872, 16.985, 17.094, 17.2, 17.302, 17.4, 17.494, 17.585, 17.672, 17.756, 17.835, 17.911, 17.981, 18.045, 18.103, 18.154, 18.198, 18.235, 18.264, 18.285, 18.297, 18.301], "dUmean/dy": [178, 178, 178, 178, 177, 176, 175, 173, 169, 163, 155, 144, 131, 116, 101, 87.1, 73.9, 62.2, 52.2, 43.8, 36.9, 31.1, 26.4, 22.6, 19.4, 16.9, 14.9, 13.3, 12, 10.9, 10.1, 9.38, 8.79, 8.29, 7.86, 7.49, 7.19, 6.91, 6.63, 6.35, 6.07, 5.81, 5.58, 5.36, 5.14, 4.92, 4.68, 4.45, 4.23, 4.04, 3.85, 3.66, 3.48, 3.28, 3.06, 2.81, 2.54, 2.25, 1.96, 1.67, 1.35, 1.02, 0.673, 0.33, 0], "Wmean": [0, 7.07e-05, 0.000283, 0.000636, 0.00113, 0.00176, 0.00252, 0.00339, 0.00435, 0.00538, 0.00643, 0.00751, 0.00864, 0.00986, 0.0112, 0.0126, 0.0141, 0.0156, 0.017, 0.0181, 0.0186, 0.0184, 0.0176, 0.0163, 0.0149, 0.0135, 0.0124, 0.0116, 0.0107, 0.00966, 0.00843, 0.00695, 0.00519, 0.00329, 0.00145, -0.000284, -0.00177, -0.00292, -0.00377, -0.00445, -0.00497, -0.0054, -0.00594, -0.00681, -0.0082, -0.00996, -0.0119, -0.0139, -0.0163, -0.0191, -0.0225, -0.0263, -0.0306, -0.0354, -0.0405, -0.0455, -0.05, -0.0539, -0.0577, -0.0615, -0.0653, -0.0685, -0.071, -0.0724, -0.0729], "dWmean/dy": [0.235, 0.235, 0.235, 0.234, 0.234, 0.232, 0.228, 0.22, 0.208, 0.194, 0.179, 0.168, 0.164, 0.164, 0.166, 0.167, 0.162, 0.148, 0.121, 0.076, 0.0159, -0.0439, -0.087, -0.107, -0.106, -0.0871, -0.0643, -0.0546, -0.061, -0.0707, -0.0818, -0.0958, -0.108, -0.106, -0.0989, -0.0881, -0.0697, -0.0506, -0.0379, -0.0303, -0.0221, -0.0216, -0.0314, -0.0522, -0.0756, -0.0841, -0.0884, -0.0974, -0.114, -0.136, -0.154, -0.172, -0.196, -0.214, -0.215, -0.199, -0.174, -0.156, -0.155, -0.159, -0.147, -0.118, -0.0788, -0.0387, 0], "Pmean": [6.217e-13, -7.3193e-10, -1.5832e-07, -3.7598e-06, -3.3837e-05, -0.00017683, -0.00065008, -0.001865, -0.0044488, -0.0092047, -0.017023, -0.028777, -0.045228, -0.066952, -0.094281, -0.12724, -0.16551, -0.20842, -0.25498, -0.30396, -0.35398, -0.40362, -0.45163, -0.49698, -0.5388, -0.57639, -0.60919, -0.63686, -0.6593, -0.67652, -0.68867, -0.69613, -0.69928, -0.69854, -0.69444, -0.68744, -0.67802, -0.66675, -0.65429, -0.64131, -0.62817, -0.61487, -0.60122, -0.58703, -0.57221, -0.55678, -0.5409, -0.52493, -0.50917, -0.49371, -0.47867, -0.46421, -0.4505, -0.43759, -0.4255, -0.41436, -0.40444, -0.39595, -0.389, -0.3836, -0.37966, -0.37702, -0.37542, -0.3746, -0.37436]} |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,207 @@ | ||
| import xlb | ||
| import time | ||
| from xlb.compute_backend import ComputeBackend | ||
| from xlb.precision_policy import PrecisionPolicy | ||
| from xlb.helper import create_nse_fields, initialize_eq | ||
| from xlb.operator.stepper import IncompressibleNavierStokesStepper | ||
| from xlb.operator.boundary_condition import RegularizedBC | ||
| from xlb.operator.macroscopic import Macroscopic | ||
| from xlb.operator.boundary_masker import IndicesBoundaryMasker | ||
| from xlb.utils import save_fields_vtk, save_image | ||
| import warp as wp | ||
| import numpy as np | ||
| import jax.numpy as jnp | ||
| import matplotlib.pyplot as plt | ||
| import json | ||
|
|
||
|
|
||
| # helper functions for this benchmark example | ||
| def vonKarman_loglaw_wall(yplus): | ||
| vonKarmanConst = 0.41 | ||
| cplus = 5.5 | ||
| uplus = np.log(yplus) / vonKarmanConst + cplus | ||
| return uplus | ||
|
|
||
|
|
||
| def get_dns_data(): | ||
| """ | ||
| Reference: DNS of Turbulent Channel Flow up to Re_tau=590, 1999, | ||
| Physics of Fluids, vol 11, 943-945. | ||
| https://turbulence.oden.utexas.edu/data/MKM/chan180/profiles/chan180.means | ||
| """ | ||
| file_name = "examples/cfd/data/turbulent_channel_dns_data.json" | ||
| with open(file_name, "r") as file: | ||
| return json.load(file) | ||
|
|
||
|
|
||
| class TurbulentChannel3D: | ||
| def __init__(self, channel_half_width, Re_tau, u_tau, grid_shape, velocity_set, backend, precision_policy): | ||
| # initialize backend | ||
| xlb.init( | ||
| velocity_set=velocity_set, | ||
| default_backend=backend, | ||
| default_precision_policy=precision_policy, | ||
| ) | ||
|
|
||
| self.channel_half_width = channel_half_width | ||
| self.Re_tau = Re_tau | ||
| self.u_tau = u_tau | ||
| self.visc = u_tau * channel_half_width / Re_tau | ||
| self.omega = 1.0 / (3.0 * self.visc + 0.5) | ||
| # DeltaPlus = Re_tau / channel_half_width | ||
| # DeltaPlus = u_tau / nu * Delta where u_tau / nu = Re_tau / channel_half_width | ||
|
|
||
| self.grid_shape = grid_shape | ||
| self.velocity_set = velocity_set | ||
| self.backend = backend | ||
| self.precision_policy = precision_policy | ||
| self.grid, self.f_0, self.f_1, self.missing_mask, self.bc_mask = create_nse_fields(grid_shape) | ||
| self.stepper = None | ||
| self.boundary_conditions = [] | ||
|
|
||
| # Setup the simulation BC, its initial conditions, and the stepper | ||
| self._setup() | ||
|
|
||
| def get_force(self): | ||
| # define the external force | ||
| shape = (self.velocity_set.d,) | ||
| force = np.zeros(shape) | ||
| force[0] = self.Re_tau**2 * self.visc**2 / self.channel_half_width**3 | ||
| return force | ||
|
|
||
| def _setup(self): | ||
| self.setup_boundary_conditions() | ||
| self.setup_boundary_masker() | ||
| self.initialize_fields() | ||
| self.setup_stepper() | ||
|
|
||
| def define_boundary_indices(self): | ||
| # top and bottom sides of the channel are no-slip and the other directions are periodic | ||
| boundingBoxIndices = self.grid.bounding_box_indices(remove_edges=True) | ||
| walls = [boundingBoxIndices["bottom"][i] + boundingBoxIndices["top"][i] for i in range(self.velocity_set.d)] | ||
| return walls | ||
|
|
||
| def setup_boundary_conditions(self): | ||
| walls = self.define_boundary_indices() | ||
| bc_walls = RegularizedBC("velocity", (0.0, 0.0, 0.0), indices=walls) | ||
| self.boundary_conditions = [bc_walls] | ||
|
|
||
| def setup_boundary_masker(self): | ||
| indices_boundary_masker = IndicesBoundaryMasker( | ||
| velocity_set=self.velocity_set, | ||
| precision_policy=self.precision_policy, | ||
| compute_backend=self.backend, | ||
| ) | ||
| self.bc_mask, self.missing_mask = indices_boundary_masker(self.boundary_conditions, self.bc_mask, self.missing_mask) | ||
|
|
||
| def initialize_fields(self): | ||
| shape = (self.velocity_set.d,) + (self.grid_shape) | ||
| np.random.seed(0) | ||
| u_init = np.random.random(shape) | ||
| if self.backend == ComputeBackend.JAX: | ||
| u_init = jnp.full(shape=shape, fill_value=1e-2 * u_init) | ||
| else: | ||
| u_init = wp.array(1e-2 * u_init, dtype=self.precision_policy.compute_precision.wp_dtype) | ||
| self.f_0 = initialize_eq(self.f_0, self.grid, self.velocity_set, self.backend, u=u_init) | ||
|
|
||
| def setup_stepper(self): | ||
| force = self.get_force() | ||
| self.stepper = IncompressibleNavierStokesStepper( | ||
| self.omega, boundary_conditions=self.boundary_conditions, collision_type="KBC", forcing_scheme="exact_difference", force_vector=force | ||
| ) | ||
|
|
||
| def run(self, num_steps, print_interval, post_process_interval=100): | ||
| start_time = time.time() | ||
| for i in range(num_steps): | ||
| self.f_1 = self.stepper(self.f_0, self.f_1, self.bc_mask, self.missing_mask, i) | ||
| self.f_0, self.f_1 = self.f_1, self.f_0 | ||
|
|
||
| if (i + 1) % print_interval == 0: | ||
| elapsed_time = time.time() - start_time | ||
| print(f"Iteration: {i + 1}/{num_steps} | Time elapsed: {elapsed_time:.2f}s") | ||
|
|
||
| if i % post_process_interval == 0 or i == num_steps - 1: | ||
| self.post_process(i) | ||
|
|
||
| def post_process(self, i): | ||
| # Write the results. We'll use JAX backend for the post-processing | ||
| if not isinstance(self.f_0, jnp.ndarray): | ||
| f_0 = wp.to_jax(self.f_0) | ||
| else: | ||
| f_0 = self.f_0 | ||
|
|
||
| macro = Macroscopic( | ||
| compute_backend=ComputeBackend.JAX, | ||
| precision_policy=self.precision_policy, | ||
| velocity_set=xlb.velocity_set.D3Q27(precision_policy=self.precision_policy, backend=ComputeBackend.JAX), | ||
| ) | ||
|
|
||
| rho, u = macro(f_0) | ||
|
|
||
| # compute velocity magnitude | ||
| u_magnitude = (u[0] ** 2 + u[1] ** 2 + u[2] ** 2) ** 0.5 | ||
| fields = {"rho": rho[0], "u_x": u[0], "u_y": u[1], "u_z": u[2], "u_magnitude": u_magnitude} | ||
| save_fields_vtk(fields, timestep=i) | ||
| save_image(fields["u_magnitude"][:, grid_size_y // 2, :], timestep=i) | ||
|
|
||
| # Save monitor plot | ||
| self.plot_uplus(u, i) | ||
| return | ||
|
|
||
| def plot_uplus(self, u, timestep): | ||
| # Compute moving average of drag coefficient, 100, 1000, 10000 | ||
| # mean streamwise velocity in wall units u^+(z) | ||
| # Wall distance in wall units to be used inside output_data | ||
| zz = np.arange(self.grid_shape[-1]) | ||
| zz = np.minimum(zz, zz.max() - zz) | ||
| yplus = zz * self.u_tau / self.visc | ||
| uplus = np.mean(u[0], axis=(0, 1)) / self.u_tau | ||
| uplus_loglaw = vonKarman_loglaw_wall(yplus) | ||
| dns_dic = get_dns_data() | ||
| plt.clf() | ||
| plt.semilogx(yplus, uplus, "r.", yplus, uplus_loglaw, "k:", dns_dic["y+"], dns_dic["Umean"], "b-") | ||
| ax = plt.gca() | ||
| ax.set_xlim([0.1, 300]) | ||
| ax.set_ylim([0, 20]) | ||
| fname = "uplus_" + str(timestep // 10000).zfill(5) + ".png" | ||
| plt.savefig(fname, format="png") | ||
|
|
||
|
|
||
| if __name__ == "__main__": | ||
| # Problem Configuration | ||
| # h: channel half-width | ||
| channel_half_width = 50 | ||
|
|
||
| # Define channel geometry based on h | ||
| grid_size_x = 6 * channel_half_width | ||
| grid_size_y = 3 * channel_half_width | ||
| grid_size_z = 2 * channel_half_width | ||
|
|
||
| # Grid parameters | ||
| grid_shape = (grid_size_x, grid_size_y, grid_size_z) | ||
|
|
||
| # Define flow regime | ||
| # Set up Reynolds number and deduce relaxation time (omega) | ||
| Re_tau = 180 | ||
| u_tau = 0.001 | ||
|
|
||
| # Runtime & backend configurations | ||
| backend = ComputeBackend.WARP | ||
| precision_policy = PrecisionPolicy.FP64FP64 | ||
| velocity_set = xlb.velocity_set.D3Q27(precision_policy=precision_policy, backend=backend) | ||
| num_steps = 10000000 | ||
| print_interval = 100000 | ||
|
|
||
| # Print simulation info | ||
| print("\n" + "=" * 50 + "\n") | ||
| print("Simulation Configuration:") | ||
| print(f"Grid size: {grid_size_x} x {grid_size_y} x {grid_size_z}") | ||
| print(f"Backend: {backend}") | ||
| print(f"Velocity set: {velocity_set}") | ||
| print(f"Precision policy: {precision_policy}") | ||
| print(f"Reynolds number: {Re_tau}") | ||
| print(f"Max iterations: {num_steps}") | ||
| print("\n" + "=" * 50 + "\n") | ||
|
|
||
| simulation = TurbulentChannel3D(channel_half_width, Re_tau, u_tau, grid_shape, velocity_set, backend, precision_policy) | ||
| simulation.run(num_steps, print_interval, post_process_interval=100000) |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Oops, something went wrong.