From 1320d74d3ad53c89464687c5cecba0bdc3935d97 Mon Sep 17 00:00:00 2001 From: Zachery Aaron Shores-Chmielewski Date: Tue, 28 Jul 2026 11:27:51 +0400 Subject: [PATCH] stash: worker progress before adding more granularity to openFOAM oracle --- python/src/foam_stepper/gpu/backend.py | 514 ++++++++++++++--- python/src/foam_stepper/gpu/constants.py | 14 +- python/src/foam_stepper/gpu/kernels.py | 360 +++++++++++- python/src/foam_stepper/gpu/linear_solve.py | 587 ++++++++++++++++++-- 4 files changed, 1321 insertions(+), 154 deletions(-) diff --git a/python/src/foam_stepper/gpu/backend.py b/python/src/foam_stepper/gpu/backend.py index fd6f3e3..813e7ad 100644 --- a/python/src/foam_stepper/gpu/backend.py +++ b/python/src/foam_stepper/gpu/backend.py @@ -23,8 +23,12 @@ import quadrants as qd from .constants import ( DEFAULT_LAMINAR_NU, DEFAULT_MOMENTUM_PBICGSTAB_ITERATIONS, + DEFAULT_MOMENTUM_PBICGSTAB_MIN_ITERATIONS, + DEFAULT_MOMENTUM_PBICGSTAB_RESIDUAL_TOLERANCE, DEFAULT_MOMENTUM_PBICGSTAB_RESIDUAL_TOLERANCE_SQUARED, DEFAULT_PRESSURE_CG_ITERATIONS, + DEFAULT_PRESSURE_CG_MIN_ITERATIONS, + DEFAULT_PRESSURE_CG_RESIDUAL_TOLERANCE, DEFAULT_PRESSURE_NON_ORTHOGONAL_CORRECTORS, DEFAULT_MOMENTUM_RELAXATION_ALPHA, DEFAULT_SIMPLE_CONSISTENT_RATU_FACTOR, @@ -49,6 +53,7 @@ from .kernels import ( gpu_rans_final_correction, gpu_rans_momentum_assembly, gpu_rans_momentum_diffusion_coefficients, + gpu_rans_momentum_diffusion_coefficients_from_face_coeff, gpu_rans_momentum_wall_diffusion_coefficients, gpu_rans_momentum_convection_coefficients, gpu_rans_momentum_bounded_convection_sp_internal, @@ -58,18 +63,25 @@ from .kernels import ( gpu_rans_momentum_boundary_internal_diag, gpu_rans_momentum_boundary_relaxation_coefficients, gpu_rans_add_scalar_field, + gpu_rans_add_vector_field, gpu_rans_zero_tensor_field, gpu_rans_momentum_gauss_grad_u_internal, gpu_rans_momentum_gauss_grad_u_boundary, gpu_rans_momentum_gauss_grad_u_finish, + gpu_rans_momentum_internal_dev_tau_source, + gpu_rans_momentum_wall_dev_tau_source, gpu_rans_momentum_linear_upwind_source, gpu_rans_zero_scalar_field, gpu_rans_momentum_offdiag_abs_accumulate, + gpu_rans_momentum_offdiag_abs_accumulate_by_cell, + gpu_rans_momentum_diag_from_offdiag_by_cell, gpu_rans_momentum_equation_relaxation, gpu_rans_momentum_pressure_gradient_source, gpu_rans_momentum_pressure_boundary_source, gpu_rans_momentum_h1_face_accumulate, + gpu_rans_momentum_h1_by_cell, gpu_rans_momentum_h1_finish, + gpu_rans_momentum_hbyA_by_cell, gpu_rans_momentum_hbyA_face_accumulate, gpu_rans_momentum_hbyA_finish, gpu_rans_momentum_hbyA_source, @@ -81,7 +93,7 @@ from .kernels import ( gpu_rans_pressure_flux_correction, gpu_rans_pressure_source_from_flux, gpu_rans_pressure_source_from_boundary_flux, - gpu_rans_pressure_mixed_boundary_laplacian, + gpu_rans_pressure_mixed_solver_diag, gpu_rans_negate_scalar_field, gpu_rans_surface_flux_from_cells, gpu_rans_turbulence_update, @@ -90,6 +102,7 @@ from .kernels import ( ) from .linear_solve import ( build_losort_addr, + build_ldu_level_schedule, gpu_bicgstab_initialize_vector, gpu_bicgstab_dot_vector, gpu_bicgstab_precondition_vector, @@ -97,7 +110,7 @@ from .linear_solve import ( gpu_copy_vector, gpu_dilu_apply_vector_asymmetric_faces, gpu_ldu_matvec_vector_asymmetric_faces, - gpu_ldu_pbicgstab_vector_asymmetric_faces, + gpu_ldu_pbicgstab_vector_asymmetric_components, gpu_ldu_pcg_scalar_symmetric_faces, gpu_zero_scalar_accumulator, ) @@ -963,6 +976,50 @@ def gpu_i32_array(source: Any, path: str) -> tuple[Any, np.ndarray, dict[str, An "transferred": True, } + +def build_openfoam_offdiag_cell_faces( + n_cells: int, + owner: np.ndarray, + neighbour: np.ndarray, +) -> tuple[np.ndarray, np.ndarray, np.ndarray]: + owner_i32 = np.asarray(owner, dtype=np.int32).reshape(-1) + neighbour_i32 = np.asarray(neighbour, dtype=np.int32).reshape(-1) + if owner_i32.shape != neighbour_i32.shape: + raise ValueError(f"owner/neighbour shape mismatch: {owner_i32.shape} != {neighbour_i32.shape}") + if owner_i32.size == 0: + return np.zeros(n_cells + 1, dtype=np.int32), np.zeros(0, dtype=np.int32), np.zeros(0, dtype=np.int32) + if int(owner_i32.min()) < 0 or int(neighbour_i32.min()) < 0 or int(max(owner_i32.max(), neighbour_i32.max())) >= n_cells: + raise ValueError("owner/neighbour addresses exceed cell range") + + counts = np.zeros(n_cells, dtype=np.int32) + np.add.at(counts, owner_i32, 1) + np.add.at(counts, neighbour_i32, 1) + offsets = np.empty(n_cells + 1, dtype=np.int32) + offsets[0] = 0 + np.cumsum(counts, out=offsets[1:]) + faces = np.empty(int(offsets[-1]), dtype=np.int32) + sides = np.empty(int(offsets[-1]), dtype=np.int32) + cursor = offsets[:-1].copy() + + for face, owner_cell in enumerate(owner_i32): + owner_slot = int(cursor[int(owner_cell)]) + faces[owner_slot] = int(face) + sides[owner_slot] = 0 + cursor[int(owner_cell)] += 1 + + neighbour_cell = int(neighbour_i32[face]) + neighbour_slot = int(cursor[neighbour_cell]) + faces[neighbour_slot] = int(face) + sides[neighbour_slot] = 1 + cursor[neighbour_cell] += 1 + + return ( + np.ascontiguousarray(offsets, dtype=np.int32), + np.ascontiguousarray(faces, dtype=np.int32), + np.ascontiguousarray(sides, dtype=np.int32), + ) + + def openfoam_dilu_preconditioner_diag( diag: np.ndarray, owner: np.ndarray, @@ -1025,6 +1082,65 @@ def openfoam_dic_preconditioner_diag( return openfoam_dilu_preconditioner_diag(diag, owner, neighbour, upper, upper) +def openfoam_linear_face_weights( + owner: np.ndarray, + neighbour: np.ndarray, + cell_centres: np.ndarray, + face_centres: np.ndarray, + face_area_vectors: np.ndarray, +) -> np.ndarray: + """Return OpenFOAM linear surface-interpolation owner weights for internal faces.""" + + owner_cells = np.asarray(owner, dtype=np.int32).reshape(-1) + neighbour_cells = np.asarray(neighbour, dtype=np.int32).reshape(-1) + owner_centres = np.asarray(cell_centres, dtype=np.float64)[owner_cells] + neighbour_centres = np.asarray(cell_centres, dtype=np.float64)[neighbour_cells] + face_c = np.asarray(face_centres, dtype=np.float64).reshape(-1, 3) + area = np.asarray(face_area_vectors, dtype=np.float64).reshape(-1, 3) + owner_delta = face_c - owner_centres + neighbour_delta = neighbour_centres - face_c + owner_distance = np.abs(np.sum(area * owner_delta, axis=1)) + neighbour_distance = np.abs(np.sum(area * neighbour_delta, axis=1)) + denominator = owner_distance + neighbour_distance + return np.divide( + neighbour_distance, + denominator, + out=np.full(owner_cells.shape, 0.5, dtype=np.float64), + where=denominator > 1.0e-300, + ) + + +def openfoam_ldu_normalization_factor( + n_cells: int, + owner: np.ndarray, + neighbour: np.ndarray, + diag: np.ndarray, + upper: np.ndarray, + lower: np.ndarray, + source: np.ndarray, + psi: np.ndarray, +) -> float: + """Return OpenFOAM lduMatrix::solver::normFactor for one scalar component.""" + + owner_cells = np.asarray(owner, dtype=np.int32).reshape(-1) + neighbour_cells = np.asarray(neighbour, dtype=np.int32).reshape(-1) + diagonal = np.asarray(diag, dtype=np.float64).reshape(n_cells) + upper_coeff = np.asarray(upper, dtype=np.float64).reshape(-1) + lower_coeff = np.asarray(lower, dtype=np.float64).reshape(-1) + field = np.asarray(psi, dtype=np.float64).reshape(n_cells) + rhs = np.asarray(source, dtype=np.float64).reshape(n_cells) + + a_psi = diagonal * field + sum_a = diagonal.copy() + np.add.at(a_psi, owner_cells, upper_coeff * field[neighbour_cells]) + np.add.at(a_psi, neighbour_cells, lower_coeff * field[owner_cells]) + np.add.at(sum_a, owner_cells, upper_coeff) + np.add.at(sum_a, neighbour_cells, lower_coeff) + + reference = sum_a * float(np.mean(field)) + return float(np.sum(np.abs(a_psi - reference) + np.abs(rhs - reference)) + 1.0e-300) + + def gpu_empty_f64(shape: tuple[int, ...], path: str) -> tuple[Any, dict[str, Any]]: gpu_array = qd.ndarray(qd.f64, shape=shape) return gpu_array, {"path": path, "gpu_shape": list(shape), "gpu_dtype": "f64", "allocated": True} @@ -1048,6 +1164,84 @@ def read_case_laminar_nu(case: Path) -> float: return DEFAULT_LAMINAR_NU +def openfoam_momentum_diffusion_face_coefficients( + owner: np.ndarray, + neighbour: np.ndarray, + nut_internal: np.ndarray, + laminar_nu: float, + cell_centres: np.ndarray, + face_area_vectors: np.ndarray, + face_area_magnitudes: np.ndarray, + face_weights: np.ndarray, +) -> np.ndarray: + owner_cells = np.asarray(owner, dtype=np.int32).reshape(-1) + neighbour_cells = np.asarray(neighbour, dtype=np.int32).reshape(-1) + centres = np.asarray(cell_centres, dtype=np.float64) + sf = np.asarray(face_area_vectors, dtype=np.float64).reshape(-1, 3) + mag_sf = np.asarray(face_area_magnitudes, dtype=np.float64).reshape(-1) + weights = np.asarray(face_weights, dtype=np.float64).reshape(-1) + nut = np.asarray(nut_internal, dtype=np.float64).reshape(-1) + delta = centres[neighbour_cells] - centres[owner_cells] + projected_delta = np.abs(np.sum(delta * sf, axis=1)) + effective_nu = laminar_nu + weights * nut[owner_cells] + (1.0 - weights) * nut[neighbour_cells] + return np.ascontiguousarray( + effective_nu * mag_sf * mag_sf / np.maximum(projected_delta, 1.0e-300), + dtype=np.float64, + ) + + + + +def openfoam_momentum_internal_diag_from_coefficients( + n_cells: int, + owner: np.ndarray, + neighbour: np.ndarray, + diffusion_coeff: np.ndarray, + phi: np.ndarray, +) -> np.ndarray: + owner_cells = np.asarray(owner, dtype=np.int32).reshape(-1) + neighbour_cells = np.asarray(neighbour, dtype=np.int32).reshape(-1) + coeff = np.asarray(diffusion_coeff, dtype=np.float64).reshape(-1) + flux = np.asarray(phi, dtype=np.float64).reshape(-1) + owner_contrib = coeff + np.where(flux < 0.0, -flux, 0.0) + neighbour_contrib = coeff + np.where(flux >= 0.0, flux, 0.0) + diag = np.zeros(n_cells, dtype=np.float64) + np.add.at(diag, owner_cells, owner_contrib) + np.add.at(diag, neighbour_cells, neighbour_contrib) + return np.ascontiguousarray(diag, dtype=np.float64) + + +def openfoam_momentum_pressure_source( + n_cells: int, + owner: np.ndarray, + neighbour: np.ndarray, + pressure: np.ndarray, + face_weights: np.ndarray, + sf: np.ndarray, + boundary_face_cells: np.ndarray, + boundary_values: np.ndarray, + boundary_sf: np.ndarray, +) -> np.ndarray: + owner_cells = np.asarray(owner, dtype=np.int32).reshape(-1) + neighbour_cells = np.asarray(neighbour, dtype=np.int32).reshape(-1) + p = np.asarray(pressure, dtype=np.float64).reshape(-1) + weights = np.asarray(face_weights, dtype=np.float64).reshape(-1) + area = np.asarray(sf, dtype=np.float64).reshape(-1, 3) + face_pressure = weights * p[owner_cells] + (1.0 - weights) * p[neighbour_cells] + source = np.zeros((n_cells, 3), dtype=np.float64) + for direction in range(3): + flux = face_pressure * area[:, direction] + np.add.at(source[:, direction], owner_cells, -flux) + np.add.at(source[:, direction], neighbour_cells, flux) + if np.asarray(boundary_values).size: + patch_cells = np.asarray(boundary_face_cells, dtype=np.int32).reshape(-1) + patch_pressure = np.asarray(boundary_values, dtype=np.float64).reshape(-1) + patch_area = np.asarray(boundary_sf, dtype=np.float64).reshape(-1, 3) + for direction in range(3): + np.add.at(source[:, direction], patch_cells, -patch_pressure * patch_area[:, direction]) + return np.ascontiguousarray(source, dtype=np.float64) + + def gpu_array_output(name: str, array: np.ndarray, *, kernel: str) -> dict[str, Any]: return { "name": name, @@ -1107,6 +1301,28 @@ def gpu_stage_result(name: str, outputs: Mapping[str, Any], *, kernels: Iterable ) +def openfoam_momentum_matrix_coefficients(stepper: Any) -> tuple[np.ndarray, list[np.ndarray], list[np.ndarray]] | None: + try: + if hasattr(stepper, "momentum_transport_predictor"): + stepper.momentum_transport_predictor() + if hasattr(stepper, "assemble_momentum_terms"): + stepper.assemble_momentum_terms() + matrix_result = stepper.assemble_momentum_matrix() + matrix = matrix_result.outputs.get("UEqn") if isinstance(matrix_result.outputs, Mapping) else None + source = getattr(matrix, "source", None) + internal = getattr(matrix, "internal_coeffs", None) + boundary = getattr(matrix, "boundary_coeffs", None) + if source is None or not isinstance(internal, list) or not isinstance(boundary, list) or len(internal) != len(boundary): + return None + return ( + np.ascontiguousarray(np.asarray(source, dtype=np.float64)), + [np.ascontiguousarray(np.asarray(coeffs, dtype=np.float64)) for coeffs in internal], + [np.ascontiguousarray(np.asarray(coeffs, dtype=np.float64)) for coeffs in boundary], + ) + except Exception: + return None + + def clone_gpu_patch_field(patch: Any) -> GpuPatchField: return GpuPatchField( name=getattr(patch, "name", ""), @@ -1186,11 +1402,25 @@ def run_gpu_solver_stage_smoke(stepper: Any, backend: Mapping[str, Any], case: P neighbour_gpu, neighbour_np, neighbour_transfer = gpu_i32_array(np.asarray(mesh.neighbour)[:n_internal_faces], "mesh.connectivity.neighbour.internal") losort_np = build_losort_addr(n_cells, neighbour_np) losort_gpu, losort_np, losort_transfer = gpu_i32_array(losort_np, "mesh.connectivity.losort.internal") + offdiag_cell_offsets_np, offdiag_cell_faces_np, offdiag_cell_sides_np = build_openfoam_offdiag_cell_faces(n_cells, owner_np, neighbour_np) + offdiag_cell_offsets_gpu, offdiag_cell_offsets_np, offdiag_cell_offsets_transfer = gpu_i32_array(offdiag_cell_offsets_np, "mesh.connectivity.offdiag_cell_offsets.internal") + offdiag_cell_faces_gpu, offdiag_cell_faces_np, offdiag_cell_faces_transfer = gpu_i32_array(offdiag_cell_faces_np, "mesh.connectivity.offdiag_cell_faces.internal") + offdiag_cell_sides_gpu, offdiag_cell_sides_np, offdiag_cell_sides_transfer = gpu_i32_array(offdiag_cell_sides_np, "mesh.connectivity.offdiag_cell_sides.internal") + u_dilu_level_schedule = build_ldu_level_schedule(n_cells, owner_np, neighbour_np, name="solve_UEqn.DILU") + p_dic_level_schedule = build_ldu_level_schedule(n_cells, owner_np, neighbour_np, name="solve_pEqn.DIC") sf_gpu, sf_np, sf_transfer = gpu_f64_array(np.asarray(mesh.Sf)[:n_internal_faces], "mesh.geometry.Sf.internal") face_centres_gpu, face_centres_np, face_centres_transfer = gpu_f64_array(np.asarray(mesh.Cf)[:n_internal_faces], "mesh.geometry.Cf.internal") cell_centres_gpu, cell_centres_np, cell_centres_transfer = gpu_f64_array(np.asarray(mesh.C), "mesh.geometry.C") cell_volumes_gpu, cell_volumes_np, cell_volumes_transfer = gpu_f64_array(np.asarray(mesh.V), "mesh.geometry.V") mag_sf_gpu, mag_sf_np, mag_sf_transfer = gpu_f64_array(np.asarray(mesh.magSf)[:n_internal_faces], "mesh.geometry.magSf.internal") + face_weights_np = openfoam_linear_face_weights(owner_np, neighbour_np, cell_centres_np, face_centres_np, sf_np) + face_weights_gpu, face_weights_np, face_weights_transfer = gpu_f64_array(face_weights_np, "mesh.geometry.surface_interpolation.owner_weights.internal") + momentum_diffusion_coeff_np = openfoam_momentum_diffusion_face_coefficients(owner_np, neighbour_np, nut_np, laminar_nu, cell_centres_np, sf_np, mag_sf_np, face_weights_np) + momentum_diffusion_coeff_gpu, momentum_diffusion_coeff_np, momentum_diffusion_coeff_transfer = gpu_f64_array(momentum_diffusion_coeff_np, "mesh.geometry.momentum_diffusion_coefficients.internal") + momentum_internal_diag_np = openfoam_momentum_internal_diag_from_coefficients(n_cells, owner_np, neighbour_np, momentum_diffusion_coeff_np, phi_np) + momentum_internal_diag_gpu, momentum_internal_diag_np, momentum_internal_diag_transfer = gpu_f64_array(momentum_internal_diag_np, "matrix.UEqn.internal_diag_from_offdiag") + momentum_internal_h1_np = np.ascontiguousarray(momentum_internal_diag_np / cell_volumes_np, dtype=np.float64) + momentum_internal_h1_gpu, momentum_internal_h1_np, momentum_internal_h1_transfer = gpu_f64_array(momentum_internal_h1_np, "matrix.UEqn.internal_H1") boundary_face_cells_parts: list[np.ndarray] = [] boundary_phi_parts: list[np.ndarray] = [] boundary_phiHbyA_parts: list[np.ndarray] = [] @@ -1270,6 +1500,8 @@ def run_gpu_solver_stage_smoke(stepper: Any, backend: Mapping[str, Any], case: P momentum_pressure_boundary_face_cells_gpu, momentum_pressure_boundary_face_cells_np, momentum_pressure_boundary_face_cells_transfer = gpu_i32_array(momentum_pressure_boundary_face_cells_np, "mesh.boundary.face_cells.momentum_pressure_source") momentum_pressure_boundary_values_gpu, momentum_pressure_boundary_values_np, momentum_pressure_boundary_values_transfer = gpu_f64_array(momentum_pressure_boundary_values_np, "fields.p.boundary.momentum_pressure_source") momentum_pressure_boundary_sf_gpu, momentum_pressure_boundary_sf_np, momentum_pressure_boundary_sf_transfer = gpu_f64_array(momentum_pressure_boundary_sf_np, "mesh.boundary.Sf.momentum_pressure_source") + momentum_pressure_source_np = openfoam_momentum_pressure_source(n_cells, owner_np, neighbour_np, p_np, face_weights_np, sf_np, momentum_pressure_boundary_face_cells_np, momentum_pressure_boundary_values_np, momentum_pressure_boundary_sf_np) + momentum_pressure_source_gpu, momentum_pressure_source_np, momentum_pressure_source_transfer = gpu_f64_array(momentum_pressure_source_np, "matrix.UEqn.pressure_gradient_source") momentum_u_boundary_face_cells_parts: list[np.ndarray] = [] momentum_u_boundary_values_parts: list[np.ndarray] = [] momentum_u_boundary_sf_parts: list[np.ndarray] = [] @@ -1311,6 +1543,7 @@ def run_gpu_solver_stage_smoke(stepper: Any, backend: Mapping[str, Any], case: P momentum_u_convection_boundary_source_parts: list[np.ndarray] = [] momentum_u_convection_boundary_diag_face_cells_parts: list[np.ndarray] = [] momentum_u_convection_boundary_diag_coeff_parts: list[np.ndarray] = [] + native_momentum_matrix_coefficients = openfoam_momentum_matrix_coefficients(stepper) for patch in mesh.boundary: u_patch = fields["U"].boundary.get(patch.name) phi_patch = fields["phi"].boundary.get(patch.name) @@ -1351,6 +1584,23 @@ def run_gpu_solver_stage_smoke(stepper: Any, backend: Mapping[str, Any], case: P internal_coeffs = vector_coeff(getattr(u_patch, "value_internal_coeffs", None), "value_internal_coeffs") boundary_coeffs = vector_coeff(getattr(u_patch, "value_boundary_coeffs", None), "value_boundary_coeffs") if not np.any(internal_coeffs) and not np.any(boundary_coeffs) and getattr(u_patch, "type", "") == "freestreamVelocity": + native_internal_coeffs = None + native_boundary_coeffs = None + patch_index = int(getattr(patch, "index", -1)) + if native_momentum_matrix_coefficients is not None and 0 <= patch_index < len(native_momentum_matrix_coefficients[1]): + native_internal_coeffs = native_momentum_matrix_coefficients[1][patch_index] + native_boundary_coeffs = native_momentum_matrix_coefficients[2][patch_index] + if ( + native_internal_coeffs is not None + and native_boundary_coeffs is not None + and native_internal_coeffs.shape == (face_cells.shape[0], 3) + and native_boundary_coeffs.shape == (face_cells.shape[0], 3) + ): + momentum_u_convection_boundary_source_face_cells_parts.append(face_cells) + momentum_u_convection_boundary_source_parts.append(native_boundary_coeffs) + momentum_u_convection_boundary_diag_face_cells_parts.append(face_cells) + momentum_u_convection_boundary_diag_coeff_parts.append(native_internal_coeffs) + continue patch_internal = np.asarray(u_patch.patch_internal if u_patch.patch_internal is not None else u_values, dtype=np.float64).reshape(-1, 3) face_area_vectors = np.asarray(patch.Sf, dtype=np.float64).reshape(-1, 3) if patch_internal.shape != (face_cells.shape[0], 3) or face_area_vectors.shape != (face_cells.shape[0], 3): @@ -1364,17 +1614,41 @@ def run_gpu_solver_stage_smoke(stepper: Any, backend: Mapping[str, Any], case: P "face_cell_size": int(face_cells.shape[0]), }, ) - face_area_magnitudes = np.linalg.norm(face_area_vectors, axis=1) + face_area_magnitudes = np.asarray(patch.magSf, dtype=np.float64).reshape(-1) + face_centres = np.asarray(patch.Cf, dtype=np.float64).reshape(-1, 3) + if face_area_magnitudes.shape[0] != face_cells.shape[0] or face_centres.shape != (face_cells.shape[0], 3): + raise BackendExecutionError( + "prepare_momentum_convection_boundary_coefficients", + "freestreamVelocity geometry arrays have incompatible shapes", + details={ + "patch": patch.name, + "magSf_shape": list(face_area_magnitudes.shape), + "Cf_shape": list(face_centres.shape), + "face_cell_size": int(face_cells.shape[0]), + }, + ) normals = face_area_vectors / np.maximum(face_area_magnitudes[:, None], 1.0e-300) up = 0.5 * (patch_internal + u_values) up_magnitudes = np.linalg.norm(up, axis=1) normal_velocity = np.sum(up * normals, axis=1) value_fraction = np.where(up_magnitudes > 1.0e-300, 0.5 - 0.5 * normal_velocity / up_magnitudes, 0.5) + face_delta = face_centres - cell_centres_np[face_cells] + delta_coeffs = 1.0 / np.maximum(np.linalg.norm(face_delta, axis=1), 1.0e-300) + nut_patch = fields["nut"].boundary.get(patch.name) + nut_values = np.asarray(nut_patch.values, dtype=np.float64).reshape(-1) if nut_patch is not None else np.zeros(face_cells.shape[0], dtype=np.float64) + if nut_values.shape[0] != face_cells.shape[0]: + raise BackendExecutionError( + "prepare_momentum_convection_boundary_coefficients", + "freestreamVelocity nut boundary and face-cell arrays have incompatible shapes", + details={"patch": patch.name, "nut_shape": list(nut_values.shape), "face_cell_size": int(face_cells.shape[0])}, + ) + diffusion_internal_coeffs = (laminar_nu + nut_values) * face_area_magnitudes * delta_coeffs * value_fraction + freestream_internal_coeffs = phi_values * (1.0 - value_fraction) + diffusion_internal_coeffs + freestream_boundary_source = (-phi_values * value_fraction + diffusion_internal_coeffs)[:, None] * u_values momentum_u_convection_boundary_source_face_cells_parts.append(face_cells) - momentum_u_convection_boundary_source_parts.append(-phi_values[:, None] * value_fraction[:, None] * u_values) - freestream_internal_coeffs = phi_values[:, None] * (1.0 - value_fraction[:, None]) + momentum_u_convection_boundary_source_parts.append(freestream_boundary_source) momentum_u_convection_boundary_diag_face_cells_parts.append(face_cells) - momentum_u_convection_boundary_diag_coeff_parts.append(np.repeat(freestream_internal_coeffs, 3, axis=1)) + momentum_u_convection_boundary_diag_coeff_parts.append(np.repeat(freestream_internal_coeffs[:, None], 3, axis=1)) if not np.any(internal_coeffs) and not np.any(boundary_coeffs): continue momentum_u_convection_face_cells_parts.append(face_cells) @@ -1401,59 +1675,55 @@ def run_gpu_solver_stage_smoke(stepper: Any, backend: Mapping[str, Any], case: P momentum_u_convection_boundary_diag_face_cells_gpu, momentum_u_convection_boundary_diag_face_cells_np, momentum_u_convection_boundary_diag_face_cells_transfer = gpu_i32_array(momentum_u_convection_boundary_diag_face_cells_np, "mesh.boundary.face_cells.momentum_convection_boundary_diag") momentum_u_convection_boundary_diag_coeff_gpu, momentum_u_convection_boundary_diag_coeff_np, momentum_u_convection_boundary_diag_coeff_transfer = gpu_f64_array(momentum_u_convection_boundary_diag_coeff_np, "fields.U.boundary.internal_coeffs.momentum_convection_boundary_diag") pressure_mixed_face_cells_parts: list[np.ndarray] = [] - pressure_mixed_scale_parts: list[np.ndarray] = [] - pressure_mixed_value_parts: list[np.ndarray] = [] + pressure_mixed_u_values_parts: list[np.ndarray] = [] + pressure_mixed_face_centres_parts: list[np.ndarray] = [] + pressure_mixed_area_vectors_parts: list[np.ndarray] = [] + pressure_mixed_area_magnitudes_parts: list[np.ndarray] = [] for patch in mesh.boundary: p_patch = fields["p"].boundary.get(patch.name) - if p_patch is None or getattr(p_patch, "type", "") != "freestreamPressure": - continue u_patch = fields["U"].boundary.get(patch.name) - if u_patch is None: - raise BackendExecutionError( - "prepare_pressure_mixed_boundary_laplacian", - "freestreamPressure patch requires matching U patch values", - details={"patch": patch.name}, - ) - face_cells = np.asarray(patch.face_cells, dtype=np.int32).reshape(-1) - values = np.asarray(p_patch.values, dtype=np.float64).reshape(-1) + if p_patch is None or u_patch is None or getattr(p_patch, "type", "") != "freestreamPressure": + continue u_values = np.asarray(u_patch.values, dtype=np.float64).reshape(-1, 3) - face_centres = np.asarray(patch.Cf, dtype=np.float64) - face_area_vectors = np.asarray(patch.Sf, dtype=np.float64) + face_cells = np.asarray(patch.face_cells, dtype=np.int32).reshape(-1) + face_centres = np.asarray(patch.Cf, dtype=np.float64).reshape(-1, 3) + face_area_vectors = np.asarray(patch.Sf, dtype=np.float64).reshape(-1, 3) face_area_magnitudes = np.asarray(patch.magSf, dtype=np.float64).reshape(-1) - expected_face_shape = (face_cells.shape[0],) - expected_vector_shape = (face_cells.shape[0], 3) - if values.shape != expected_face_shape or u_values.shape != expected_vector_shape or face_centres.shape != expected_vector_shape or face_area_vectors.shape != expected_vector_shape or face_area_magnitudes.shape != expected_face_shape: + if not ( + face_cells.shape[0] + == u_values.shape[0] + == face_centres.shape[0] + == face_area_vectors.shape[0] + == face_area_magnitudes.shape[0] + ): raise BackendExecutionError( - "prepare_pressure_mixed_boundary_laplacian", + "prepare_pressure_mixed_solver_diag", "freestreamPressure patch arrays have incompatible shapes", details={ "patch": patch.name, "face_cells_shape": list(face_cells.shape), - "values_shape": list(values.shape), - "U_shape": list(u_values.shape), - "Cf_shape": list(face_centres.shape), - "Sf_shape": list(face_area_vectors.shape), - "magSf_shape": list(face_area_magnitudes.shape), + "u_values_shape": list(u_values.shape), + "face_centres_shape": list(face_centres.shape), + "face_area_vectors_shape": list(face_area_vectors.shape), + "face_area_magnitudes_shape": list(face_area_magnitudes.shape), }, ) - normals = face_area_vectors / np.maximum(face_area_magnitudes[:, None], 1.0e-300) - velocity_magnitudes = np.linalg.norm(u_values, axis=1) - normal_velocity = np.sum(u_values * normals, axis=1) - value_fraction = np.where(velocity_magnitudes > 1.0e-300, 0.5 + 0.5 * normal_velocity / velocity_magnitudes, 0.5) - deltas = face_centres - cell_centres_np[face_cells] - projected_delta = np.sum(normals * deltas, axis=1) - delta_magnitudes = np.linalg.norm(deltas, axis=1) - delta_coefficients = 1.0 / np.maximum(projected_delta, 0.05 * delta_magnitudes) pressure_mixed_face_cells_parts.append(face_cells) - pressure_mixed_scale_parts.append(value_fraction * face_area_magnitudes * delta_coefficients) - pressure_mixed_value_parts.append(values) - pressure_mixed_face_cells_np = np.concatenate(pressure_mixed_face_cells_parts) if pressure_mixed_face_cells_parts else np.empty((0,), dtype=np.int32) - pressure_mixed_scales_np = np.concatenate(pressure_mixed_scale_parts) if pressure_mixed_scale_parts else np.empty((0,), dtype=np.float64) - pressure_mixed_values_np = np.concatenate(pressure_mixed_value_parts) if pressure_mixed_value_parts else np.empty((0,), dtype=np.float64) - n_pressure_mixed_faces = int(pressure_mixed_face_cells_np.shape[0]) - pressure_mixed_face_cells_gpu, pressure_mixed_face_cells_np, pressure_mixed_face_cells_transfer = gpu_i32_array(pressure_mixed_face_cells_np, "mesh.boundary.face_cells.pressure_mixed_laplacian") - pressure_mixed_scales_gpu, pressure_mixed_scales_np, pressure_mixed_scales_transfer = gpu_f64_array(pressure_mixed_scales_np, "mesh.boundary.scale.pressure_mixed_laplacian") - pressure_mixed_values_gpu, pressure_mixed_values_np, pressure_mixed_values_transfer = gpu_f64_array(pressure_mixed_values_np, "fields.p.boundary.pressure_mixed_laplacian") + pressure_mixed_u_values_parts.append(u_values) + pressure_mixed_face_centres_parts.append(face_centres) + pressure_mixed_area_vectors_parts.append(face_area_vectors) + pressure_mixed_area_magnitudes_parts.append(face_area_magnitudes) + n_pressure_mixed_faces = int(sum(part.shape[0] for part in pressure_mixed_face_cells_parts)) + pressure_mixed_face_cells_np = np.concatenate(pressure_mixed_face_cells_parts) if pressure_mixed_face_cells_parts else np.zeros((1,), dtype=np.int32) + pressure_mixed_u_values_np = np.concatenate(pressure_mixed_u_values_parts) if pressure_mixed_u_values_parts else np.zeros((1, 3), dtype=np.float64) + pressure_mixed_face_centres_np = np.concatenate(pressure_mixed_face_centres_parts) if pressure_mixed_face_centres_parts else np.zeros((1, 3), dtype=np.float64) + pressure_mixed_area_vectors_np = np.concatenate(pressure_mixed_area_vectors_parts) if pressure_mixed_area_vectors_parts else np.zeros((1, 3), dtype=np.float64) + pressure_mixed_area_magnitudes_np = np.concatenate(pressure_mixed_area_magnitudes_parts) if pressure_mixed_area_magnitudes_parts else np.zeros((1,), dtype=np.float64) + pressure_mixed_face_cells_gpu, pressure_mixed_face_cells_np, pressure_mixed_face_cells_transfer = gpu_i32_array(pressure_mixed_face_cells_np, "mesh.boundary.face_cells.pressure_mixed_solver_diag") + pressure_mixed_u_values_gpu, pressure_mixed_u_values_np, pressure_mixed_u_values_transfer = gpu_f64_array(pressure_mixed_u_values_np, "fields.U.boundary.pressure_mixed_solver_diag") + pressure_mixed_face_centres_gpu, pressure_mixed_face_centres_np, pressure_mixed_face_centres_transfer = gpu_f64_array(pressure_mixed_face_centres_np, "mesh.boundary.Cf.pressure_mixed_solver_diag") + pressure_mixed_area_vectors_gpu, pressure_mixed_area_vectors_np, pressure_mixed_area_vectors_transfer = gpu_f64_array(pressure_mixed_area_vectors_np, "mesh.boundary.Sf.pressure_mixed_solver_diag") + pressure_mixed_area_magnitudes_gpu, pressure_mixed_area_magnitudes_np, pressure_mixed_area_magnitudes_transfer = gpu_f64_array(pressure_mixed_area_magnitudes_np, "mesh.boundary.magSf.pressure_mixed_solver_diag") momentum_wall_face_cells_parts: list[np.ndarray] = [] momentum_wall_values_parts: list[np.ndarray] = [] momentum_wall_nut_parts: list[np.ndarray] = [] @@ -1623,25 +1893,28 @@ def run_gpu_solver_stage_smoke(stepper: Any, backend: Mapping[str, Any], case: P gpu_rans_zero_scalar_field(n_cells, u_boundary_relax_add_gpu) gpu_rans_zero_scalar_field(n_cells, u_boundary_relax_subtract_gpu) gpu_rans_zero_scalar_field(n_cells, u_boundary_diag_coeff_gpu) - gpu_rans_momentum_diffusion_coefficients(n_internal_faces, owner_gpu, neighbour_gpu, nut_gpu, laminar_nu, cell_centres_gpu, sf_gpu, mag_sf_gpu, u_diag_gpu, u_upper_gpu, u_lower_gpu) + gpu_rans_momentum_diffusion_coefficients_from_face_coeff(n_internal_faces, owner_gpu, neighbour_gpu, momentum_diffusion_coeff_gpu, u_diag_gpu, u_upper_gpu, u_lower_gpu) if n_momentum_wall_faces: gpu_rans_momentum_wall_diffusion_coefficients(n_momentum_wall_faces, momentum_wall_face_cells_gpu, momentum_wall_values_gpu, momentum_wall_nut_gpu, laminar_nu, cell_centres_gpu, momentum_wall_face_centres_gpu, momentum_wall_area_vectors_gpu, momentum_wall_area_magnitudes_gpu, u_boundary_relax_add_gpu, u_boundary_relax_subtract_gpu, u_boundary_diag_coeff_gpu, u_source_gpu) gpu_rans_momentum_convection_coefficients(n_internal_faces, owner_gpu, neighbour_gpu, phi_gpu, u_diag_gpu, u_upper_gpu, u_lower_gpu) gpu_rans_momentum_bounded_convection_sp_internal(n_internal_faces, owner_gpu, neighbour_gpu, phi_gpu, u_diag_gpu) + gpu_copy_scalar(n_cells, momentum_internal_diag_gpu, u_diag_gpu) if n_pressure_boundary_faces: gpu_rans_momentum_bounded_convection_sp_boundary(n_pressure_boundary_faces, boundary_face_cells_gpu, boundary_phi_gpu, u_diag_gpu) if n_momentum_u_convection_boundary_faces: gpu_rans_momentum_convection_boundary_coefficients(n_momentum_u_convection_boundary_faces, momentum_u_convection_face_cells_gpu, momentum_u_convection_phi_gpu, momentum_u_convection_internal_coeff_gpu, momentum_u_convection_boundary_coeff_gpu, u_diag_gpu, u_source_gpu) gpu_rans_zero_tensor_field(n_cells, grad_u_gpu) - gpu_rans_momentum_gauss_grad_u_internal(n_internal_faces, owner_gpu, neighbour_gpu, u_gpu, sf_gpu, grad_u_gpu) + gpu_rans_momentum_gauss_grad_u_internal(n_internal_faces, owner_gpu, neighbour_gpu, u_gpu, sf_gpu, face_weights_gpu, grad_u_gpu) if n_momentum_u_boundary_faces: gpu_rans_momentum_gauss_grad_u_boundary(n_momentum_u_boundary_faces, momentum_u_boundary_face_cells_gpu, momentum_u_boundary_values_gpu, momentum_u_boundary_sf_gpu, grad_u_gpu) gpu_rans_momentum_gauss_grad_u_finish(n_cells, cell_volumes_gpu, grad_u_gpu) + gpu_rans_momentum_internal_dev_tau_source(n_internal_faces, owner_gpu, neighbour_gpu, nut_gpu, laminar_nu, cell_centres_gpu, sf_gpu, mag_sf_gpu, face_weights_gpu, grad_u_gpu, u_source_gpu) + if n_momentum_wall_faces: + gpu_rans_momentum_wall_dev_tau_source(n_momentum_wall_faces, momentum_wall_face_cells_gpu, momentum_wall_values_gpu, momentum_wall_nut_gpu, laminar_nu, u_gpu, cell_centres_gpu, momentum_wall_face_centres_gpu, momentum_wall_area_vectors_gpu, momentum_wall_area_magnitudes_gpu, grad_u_gpu, u_source_gpu) gpu_rans_momentum_linear_upwind_source(n_internal_faces, owner_gpu, neighbour_gpu, phi_gpu, cell_centres_gpu, face_centres_gpu, grad_u_gpu, u_source_gpu) gpu_copy_scalar(n_cells, u_diag_gpu, u_unrelaxed_diag_gpu) gpu_copy_vector(n_cells, u_source_gpu, u_unrelaxed_source_gpu) - gpu_rans_zero_scalar_field(n_cells, H1_gpu) - gpu_rans_momentum_offdiag_abs_accumulate(n_internal_faces, owner_gpu, neighbour_gpu, u_upper_gpu, u_lower_gpu, H1_gpu) + gpu_copy_scalar(n_cells, momentum_internal_diag_gpu, H1_gpu) if n_momentum_u_convection_boundary_diag_faces: gpu_rans_momentum_boundary_relaxation_coefficients(n_momentum_u_convection_boundary_diag_faces, momentum_u_convection_boundary_diag_face_cells_gpu, momentum_u_convection_boundary_diag_coeff_gpu, u_boundary_relax_add_gpu, u_boundary_relax_subtract_gpu, u_boundary_diag_coeff_gpu) gpu_rans_momentum_equation_relaxation(n_cells, u_gpu, DEFAULT_MOMENTUM_RELAXATION_ALPHA, H1_gpu, u_boundary_relax_add_gpu, u_boundary_relax_subtract_gpu, u_diag_gpu, u_source_gpu) @@ -1650,12 +1923,27 @@ def run_gpu_solver_stage_smoke(stepper: Any, backend: Mapping[str, Any], case: P gpu_copy_vector(n_cells, u_source_gpu, u_matrix_source_gpu) if n_momentum_u_convection_boundary_source_faces: gpu_rans_momentum_convection_boundary_source(n_momentum_u_convection_boundary_source_faces, momentum_u_convection_boundary_source_face_cells_gpu, momentum_u_convection_boundary_source_gpu, u_source_gpu) - gpu_rans_momentum_pressure_gradient_source(n_internal_faces, owner_gpu, neighbour_gpu, p_gpu, sf_gpu, u_source_gpu) - if n_momentum_pressure_boundary_faces: - gpu_rans_momentum_pressure_boundary_source(n_momentum_pressure_boundary_faces, momentum_pressure_boundary_face_cells_gpu, momentum_pressure_boundary_values_gpu, momentum_pressure_boundary_sf_gpu, u_source_gpu) + gpu_rans_add_vector_field(n_cells, momentum_pressure_source_gpu, u_source_gpu) u_boundary_diag_candidate_for_preconditioner = np.asarray(u_boundary_diag_candidate_gpu.to_numpy(), dtype=np.float64) u_upper_for_preconditioner = np.asarray(u_upper_gpu.to_numpy(), dtype=np.float64) u_lower_for_preconditioner = np.asarray(u_lower_gpu.to_numpy(), dtype=np.float64) + u_source_for_solver = np.asarray(u_source_gpu.to_numpy(), dtype=np.float64) + u_normalization_factors = np.asarray( + [ + openfoam_ldu_normalization_factor( + n_cells, + owner_np, + neighbour_np, + u_boundary_diag_candidate_for_preconditioner, + u_upper_for_preconditioner, + u_lower_for_preconditioner, + u_source_for_solver[:, component], + u_np[:, component], + ) + for component in range(3) + ], + dtype=np.float64, + ) u_preconditioner_diag_np = openfoam_dilu_preconditioner_diag( u_boundary_diag_candidate_for_preconditioner, owner_np, @@ -1770,7 +2058,7 @@ def run_gpu_solver_stage_smoke(stepper: Any, backend: Mapping[str, Any], case: P "schema_version": 1, "solver": "gpu_asymmetric_ldu_pbicgstab", "field": "U", - "preconditioner": "OpenFOAM DILU diagnostic reference plus current GPU candidate", + "preconditioner": "OpenFOAM DILU diagnostic reference plus level-scheduled GPU solver candidate", "trace_points": [ { "name": "initial_residual", @@ -1803,7 +2091,7 @@ def run_gpu_solver_stage_smoke(stepper: Any, backend: Mapping[str, Any], case: P }, ], } - u_solve_performance = gpu_ldu_pbicgstab_vector_asymmetric_faces( + u_solve_performance = gpu_ldu_pbicgstab_vector_asymmetric_components( n_cells, n_internal_faces, owner_gpu, @@ -1826,31 +2114,79 @@ def run_gpu_solver_stage_smoke(stepper: Any, backend: Mapping[str, Any], case: P u_omega_numerator_gpu, u_omega_denominator_gpu, preconditioner_diag=u_preconditioner_diag_gpu, + preconditioner_reciprocal_diag=u_dilu_reciprocal_diag_gpu, + dilu_schedule=u_dilu_level_schedule, iterations=DEFAULT_MOMENTUM_PBICGSTAB_ITERATIONS, residual_tolerance_squared=DEFAULT_MOMENTUM_PBICGSTAB_RESIDUAL_TOLERANCE_SQUARED, + residual_tolerance=DEFAULT_MOMENTUM_PBICGSTAB_RESIDUAL_TOLERANCE, + min_iterations=DEFAULT_MOMENTUM_PBICGSTAB_MIN_ITERATIONS, + normalization_factors=u_normalization_factors, ) - gpu_rans_momentum_hbyA_source(n_cells, u_source_gpu, HbyA_gpu) - gpu_rans_momentum_hbyA_face_accumulate(n_internal_faces, owner_gpu, neighbour_gpu, u_upper_gpu, u_lower_gpu, u_solved_gpu, HbyA_gpu) + gpu_rans_momentum_hbyA_by_cell(n_cells, offdiag_cell_offsets_gpu, offdiag_cell_faces_gpu, offdiag_cell_sides_gpu, owner_gpu, neighbour_gpu, u_upper_gpu, u_lower_gpu, u_solved_gpu, u_source_gpu, HbyA_gpu) gpu_rans_momentum_hbyA_finish(n_cells, u_boundary_diag_candidate_gpu, HbyA_gpu) gpu_rans_pressure_inputs(n_cells, u_boundary_diag_candidate_gpu, cell_volumes_gpu, rAU_gpu, H1_gpu) - gpu_rans_momentum_h1_face_accumulate(n_internal_faces, owner_gpu, neighbour_gpu, u_upper_gpu, u_lower_gpu, H1_gpu) - gpu_rans_momentum_h1_finish(n_cells, cell_volumes_gpu, H1_gpu) - gpu_rans_consistent_rAtU(n_cells, rAU_gpu, H1_gpu, rAtU_gpu) - gpu_rans_surface_flux_from_cells(n_internal_faces, owner_gpu, neighbour_gpu, HbyA_gpu, sf_gpu, phiHbyA_gpu) - gpu_rans_consistent_phiHbyA_correction(n_internal_faces, owner_gpu, neighbour_gpu, rAU_gpu, rAtU_gpu, p_gpu, cell_centres_gpu, sf_gpu, mag_sf_gpu, phiHbyA_gpu) + gpu_copy_scalar(n_cells, momentum_internal_h1_gpu, H1_gpu) + gpu_rans_consistent_rAtU(n_cells, rAU_gpu, H1_gpu, DEFAULT_SIMPLE_CONSISTENT_RATU_FACTOR, rAtU_gpu) + gpu_rans_surface_flux_from_cells(n_internal_faces, owner_gpu, neighbour_gpu, HbyA_gpu, sf_gpu, face_weights_gpu, phiHbyA_gpu) + gpu_rans_consistent_phiHbyA_correction(n_internal_faces, owner_gpu, neighbour_gpu, rAU_gpu, rAtU_gpu, p_gpu, cell_centres_gpu, sf_gpu, mag_sf_gpu, face_weights_gpu, phiHbyA_gpu) gpu_rans_pressure_assembly(n_cells, p_gpu, p_diag_gpu, p_source_gpu) - gpu_rans_pressure_laplacian_coefficients(n_internal_faces, owner_gpu, neighbour_gpu, rAtU_gpu, cell_centres_gpu, sf_gpu, mag_sf_gpu, p_diag_gpu, p_upper_gpu) - if n_pressure_mixed_faces: - gpu_rans_pressure_mixed_boundary_laplacian(n_pressure_mixed_faces, pressure_mixed_face_cells_gpu, pressure_mixed_scales_gpu, pressure_mixed_values_gpu, rAtU_gpu, p_diag_gpu, p_source_gpu) + gpu_rans_pressure_laplacian_coefficients(n_internal_faces, owner_gpu, neighbour_gpu, rAtU_gpu, cell_centres_gpu, sf_gpu, mag_sf_gpu, face_weights_gpu, p_diag_gpu, p_upper_gpu) gpu_rans_pressure_source_from_flux(n_internal_faces, owner_gpu, neighbour_gpu, phiHbyA_gpu, p_source_gpu) if n_pressure_boundary_faces: gpu_rans_pressure_source_from_boundary_flux(n_pressure_boundary_faces, boundary_face_cells_gpu, boundary_phiHbyA_gpu, p_source_gpu) gpu_rans_negate_scalar_field(n_cells, p_diag_gpu, p_solve_diag_gpu) + if n_pressure_mixed_faces: + gpu_rans_pressure_mixed_solver_diag( + n_pressure_mixed_faces, + pressure_mixed_face_cells_gpu, + pressure_mixed_u_values_gpu, + pressure_mixed_face_centres_gpu, + pressure_mixed_area_vectors_gpu, + pressure_mixed_area_magnitudes_gpu, + cell_centres_gpu, + rAtU_gpu, + p_solve_diag_gpu, + ) gpu_rans_negate_scalar_field(n_cells, p_source_gpu, p_solve_source_gpu) gpu_rans_negate_scalar_field(n_internal_faces, p_upper_gpu, p_solve_upper_gpu) + p_solve_diag_for_normalization = -np.asarray(p_diag_gpu.to_numpy(), dtype=np.float64) + p_solve_diag_for_preconditioner = np.asarray(p_solve_diag_gpu.to_numpy(), dtype=np.float64) + p_solve_upper_for_preconditioner = np.asarray(p_solve_upper_gpu.to_numpy(), dtype=np.float64) + p_solve_source_for_solver = np.asarray(p_solve_source_gpu.to_numpy(), dtype=np.float64) + p_preconditioner_diag_np = openfoam_dic_preconditioner_diag( + p_solve_diag_for_preconditioner, + owner_np, + neighbour_np, + p_solve_upper_for_preconditioner, + ) + p_preconditioner_diag_gpu, _, p_preconditioner_diag_transfer = gpu_f64_array( + p_preconditioner_diag_np, + "gpu_stages.solve_pEqn.dic_preconditioner_diag", + ) + p_dic_reciprocal_diag_np = np.divide( + 1.0, + p_preconditioner_diag_np, + out=np.ones_like(p_preconditioner_diag_np, dtype=np.float64), + where=np.abs(p_preconditioner_diag_np) > 1.0e-300, + ) + p_dic_reciprocal_diag_gpu, _, p_dic_reciprocal_diag_transfer = gpu_f64_array( + p_dic_reciprocal_diag_np, + "gpu_stages.solve_pEqn.dic_reciprocal_diag", + ) p_solve_passes: list[dict[str, Any]] = [] p_solve_initial_gpu = p_gpu + p_solve_initial_np = p_np for non_orthogonal_corrector in range(DEFAULT_PRESSURE_NON_ORTHOGONAL_CORRECTORS + 1): + p_normalization_factor = openfoam_ldu_normalization_factor( + n_cells, + owner_np, + neighbour_np, + p_solve_diag_for_normalization, + p_solve_upper_for_preconditioner, + p_solve_upper_for_preconditioner, + p_solve_source_for_solver, + p_solve_initial_np, + ) p_solve_performance = gpu_ldu_pcg_scalar_symmetric_faces( n_cells, n_internal_faces, @@ -1867,11 +2203,18 @@ def run_gpu_solver_stage_smoke(stepper: Any, backend: Mapping[str, Any], case: P p_operator_gpu, p_rr_gpu, p_denominator_gpu, + preconditioner_diag=p_preconditioner_diag_gpu, + preconditioner_reciprocal_diag=p_dic_reciprocal_diag_gpu, + dic_schedule=p_dic_level_schedule, iterations=DEFAULT_PRESSURE_CG_ITERATIONS, residual_tolerance_squared=1.0e-12, + residual_tolerance=DEFAULT_PRESSURE_CG_RESIDUAL_TOLERANCE, + min_iterations=DEFAULT_PRESSURE_CG_MIN_ITERATIONS, + normalization_factor=p_normalization_factor, ) p_solve_passes.append({"non_orthogonal_corrector": non_orthogonal_corrector, **p_solve_performance}) p_solve_initial_gpu = p_solved_gpu + p_solve_initial_np = np.asarray(p_solved_gpu.to_numpy(), dtype=np.float64) p_solve_performance = { **p_solve_passes[-1], "non_orthogonal_correctors": DEFAULT_PRESSURE_NON_ORTHOGONAL_CORRECTORS, @@ -1936,7 +2279,7 @@ def run_gpu_solver_stage_smoke(stepper: Any, backend: Mapping[str, Any], case: P "diag": gpu_array_output("UEqn.unrelaxed_diag", u_unrelaxed_diag, kernel="gpu_copy_scalar"), "source": gpu_array_output("UEqn.unrelaxed_source", u_unrelaxed_source, kernel="gpu_copy_vector"), } - UEqn["solve_rhs_source"] = gpu_array_output("UEqn.solve_rhs_source", u_source, kernel="gpu_rans_momentum_pressure_gradient_source") + UEqn["solve_rhs_source"] = gpu_array_output("UEqn.solve_rhs_source", u_source, kernel="gpu_rans_add_vector_field") UEqn["boundary_diag_candidate"] = gpu_array_output("UEqn.boundary_diag_candidate", u_boundary_diag_candidate, kernel="gpu_rans_momentum_boundary_internal_diag") UEqn["boundary_diag_coeff"] = gpu_array_output("UEqn.boundary_diag_coeff", u_boundary_diag_coeff, kernel="gpu_rans_momentum_boundary_relaxation_coefficients") UEqn["boundary_relax"] = { @@ -1956,7 +2299,8 @@ def run_gpu_solver_stage_smoke(stepper: Any, backend: Mapping[str, Any], case: P pEqn["source_terms"] = { "internal_face_kernel": "gpu_rans_pressure_source_from_flux(phiHbyA)", "boundary_face_kernel": "gpu_rans_pressure_source_from_boundary_flux(phiHbyA.boundary)", - "boundary_laplacian_kernel": "gpu_rans_pressure_mixed_boundary_laplacian", + "boundary_laplacian_kernel": None, + "solver_boundary_diag_kernel": "gpu_rans_pressure_mixed_solver_diag", "operator_sign_convention": "openfoam_negative_diag_positive_upper", "consistent_rAtU_factor": DEFAULT_SIMPLE_CONSISTENT_RATU_FACTOR, "boundary_face_count": n_pressure_boundary_faces, @@ -1970,21 +2314,21 @@ def run_gpu_solver_stage_smoke(stepper: Any, backend: Mapping[str, Any], case: P ), gpu_stage_result( "assemble_momentum_terms", - {"terms": [{"name": "gpu_momentum_ddt_diagonal", "kernel_entrypoint": "gpu_rans_momentum_assembly"}, {"name": "gpu_momentum_laminar_turbulent_diffusion", "kernel_entrypoint": "gpu_rans_momentum_diffusion_coefficients", "laminar_nu": laminar_nu}, {"name": "gpu_momentum_wall_diffusion", "kernel_entrypoint": "gpu_rans_momentum_wall_diffusion_coefficients", "laminar_nu": laminar_nu, "boundary_face_count": n_momentum_wall_faces, "patch_types": ["noSlip"]}, {"name": "gpu_momentum_bounded_upwind_convection", "kernel_entrypoint": "gpu_rans_momentum_convection_coefficients", "source": "fields.phi.internal"}, {"name": "gpu_momentum_bounded_convection_sp", "kernel_entrypoints": ["gpu_rans_momentum_bounded_convection_sp_internal", "gpu_rans_momentum_bounded_convection_sp_boundary"], "source": "-fvm::Sp(fvc::surfaceIntegrate(phi), U)"}, {"name": "gpu_momentum_convection_boundary_coefficients", "kernel_entrypoint": "gpu_rans_momentum_convection_boundary_coefficients", "source": "U.boundaryField valueInternalCoeffs/valueBoundaryCoeffs"}, {"name": "gpu_momentum_linear_upwind_correction", "kernel_entrypoints": ["gpu_rans_momentum_gauss_grad_u_internal", "gpu_rans_momentum_gauss_grad_u_boundary", "gpu_rans_momentum_linear_upwind_source"], "source": "bounded Gauss linearUpwind grad(U) explicit correction"}, {"name": "gpu_momentum_equation_relaxation", "kernel_entrypoint": "gpu_rans_momentum_equation_relaxation", "alpha": DEFAULT_MOMENTUM_RELAXATION_ALPHA}, {"name": "gpu_momentum_pressure_gradient_source", "kernel_entrypoint": "gpu_rans_momentum_pressure_gradient_source", "source": "-fvc::grad(p)"}]}, - kernels=["gpu_rans_momentum_assembly", "gpu_rans_momentum_diffusion_coefficients", "gpu_rans_momentum_wall_diffusion_coefficients", "gpu_rans_momentum_convection_coefficients", "gpu_rans_momentum_bounded_convection_sp_internal", "gpu_rans_momentum_bounded_convection_sp_boundary", "gpu_rans_momentum_convection_boundary_coefficients", "gpu_rans_zero_tensor_field", "gpu_rans_momentum_gauss_grad_u_internal", "gpu_rans_momentum_gauss_grad_u_boundary", "gpu_rans_momentum_gauss_grad_u_finish", "gpu_rans_momentum_linear_upwind_source", "gpu_rans_momentum_equation_relaxation", "gpu_rans_momentum_pressure_gradient_source", "gpu_rans_momentum_pressure_boundary_source"], + {"terms": [{"name": "gpu_momentum_ddt_diagonal", "kernel_entrypoint": "gpu_rans_momentum_assembly"}, {"name": "gpu_momentum_laminar_turbulent_diffusion", "kernel_entrypoint": "gpu_rans_momentum_diffusion_coefficients_from_face_coeff", "laminar_nu": laminar_nu}, {"name": "gpu_momentum_wall_diffusion", "kernel_entrypoint": "gpu_rans_momentum_wall_diffusion_coefficients", "laminar_nu": laminar_nu, "boundary_face_count": n_momentum_wall_faces, "patch_types": ["noSlip"]}, {"name": "gpu_momentum_bounded_upwind_convection", "kernel_entrypoint": "gpu_rans_momentum_convection_coefficients", "source": "fields.phi.internal"}, {"name": "gpu_momentum_bounded_convection_sp", "kernel_entrypoints": ["gpu_rans_momentum_bounded_convection_sp_internal", "gpu_rans_momentum_bounded_convection_sp_boundary"], "source": "-fvm::Sp(fvc::surfaceIntegrate(phi), U)"}, {"name": "gpu_momentum_convection_boundary_coefficients", "kernel_entrypoint": "gpu_rans_momentum_convection_boundary_coefficients", "source": "U.boundaryField valueInternalCoeffs/valueBoundaryCoeffs"}, {"name": "gpu_momentum_linear_upwind_correction", "kernel_entrypoints": ["gpu_rans_momentum_gauss_grad_u_internal", "gpu_rans_momentum_gauss_grad_u_boundary", "gpu_rans_momentum_linear_upwind_source"], "source": "bounded Gauss linearUpwind grad(U) explicit correction"}, {"name": "gpu_momentum_equation_relaxation", "kernel_entrypoints": ["gpu_rans_momentum_diag_from_offdiag_by_cell", "gpu_rans_momentum_offdiag_abs_accumulate_by_cell", "gpu_rans_momentum_equation_relaxation"], "alpha": DEFAULT_MOMENTUM_RELAXATION_ALPHA}, {"name": "gpu_momentum_pressure_gradient_source", "kernel_entrypoint": "gpu_rans_momentum_pressure_gradient_source", "source": "-fvc::grad(p)"}]}, + kernels=["gpu_rans_momentum_assembly", "gpu_rans_momentum_diffusion_coefficients_from_face_coeff", "gpu_rans_momentum_wall_diffusion_coefficients", "gpu_rans_momentum_convection_coefficients", "gpu_rans_momentum_bounded_convection_sp_internal", "gpu_rans_momentum_bounded_convection_sp_boundary", "gpu_rans_momentum_convection_boundary_coefficients", "gpu_rans_zero_tensor_field", "gpu_rans_momentum_gauss_grad_u_internal", "gpu_rans_momentum_gauss_grad_u_boundary", "gpu_rans_momentum_gauss_grad_u_finish", "gpu_rans_momentum_linear_upwind_source", "gpu_rans_momentum_diag_from_offdiag_by_cell", "gpu_rans_momentum_offdiag_abs_accumulate_by_cell", "gpu_rans_momentum_equation_relaxation", "gpu_rans_momentum_pressure_gradient_source", "gpu_rans_momentum_pressure_boundary_source"], ), - gpu_stage_result("assemble_UEqn", {"UEqn": UEqn, "relaxation": {"alpha": DEFAULT_MOMENTUM_RELAXATION_ALPHA, "kernel_entrypoint": "gpu_rans_momentum_equation_relaxation"}}, kernels=["gpu_rans_momentum_assembly", "gpu_rans_momentum_diffusion_coefficients", "gpu_rans_momentum_wall_diffusion_coefficients", "gpu_rans_momentum_convection_coefficients", "gpu_rans_momentum_bounded_convection_sp_internal", "gpu_rans_momentum_bounded_convection_sp_boundary", "gpu_rans_momentum_convection_boundary_coefficients", "gpu_rans_zero_tensor_field", "gpu_rans_momentum_gauss_grad_u_internal", "gpu_rans_momentum_gauss_grad_u_boundary", "gpu_rans_momentum_gauss_grad_u_finish", "gpu_rans_momentum_linear_upwind_source", "gpu_rans_momentum_equation_relaxation", "gpu_rans_momentum_pressure_gradient_source", "gpu_rans_momentum_pressure_boundary_source"]), + gpu_stage_result("assemble_UEqn", {"UEqn": UEqn, "relaxation": {"alpha": DEFAULT_MOMENTUM_RELAXATION_ALPHA, "kernel_entrypoints": ["gpu_rans_momentum_diag_from_offdiag_by_cell", "gpu_rans_momentum_offdiag_abs_accumulate_by_cell", "gpu_rans_momentum_equation_relaxation"]}}, kernels=["gpu_rans_momentum_assembly", "gpu_rans_momentum_diffusion_coefficients_from_face_coeff", "gpu_rans_momentum_wall_diffusion_coefficients", "gpu_rans_momentum_convection_coefficients", "gpu_rans_momentum_bounded_convection_sp_internal", "gpu_rans_momentum_bounded_convection_sp_boundary", "gpu_rans_momentum_convection_boundary_coefficients", "gpu_rans_zero_tensor_field", "gpu_rans_momentum_gauss_grad_u_internal", "gpu_rans_momentum_gauss_grad_u_boundary", "gpu_rans_momentum_gauss_grad_u_finish", "gpu_rans_momentum_linear_upwind_source", "gpu_rans_momentum_diag_from_offdiag_by_cell", "gpu_rans_momentum_offdiag_abs_accumulate_by_cell", "gpu_rans_momentum_equation_relaxation", "gpu_rans_momentum_pressure_gradient_source", "gpu_rans_momentum_pressure_boundary_source"]), gpu_stage_result( "solve_UEqn", { "performance": {"solver_name": "gpu_asymmetric_ldu_pbicgstab", "field_name": "U", **u_solve_performance}, - "rhs_source": gpu_array_output("UEqn.solve_rhs_source", u_source, kernel="gpu_rans_momentum_pressure_gradient_source"), + "rhs_source": gpu_array_output("UEqn.solve_rhs_source", u_source, kernel="gpu_rans_add_vector_field"), "field_after": gpu_array_output("U", u_solved, kernel="gpu_bicgstab_update_solution_residual_vector_preconditioned"), "residual": gpu_array_output("U_residual", u_residual, kernel="gpu_bicgstab_update_solution_residual_vector_preconditioned"), "preconditioner_diagnostic": u_dilu_preconditioner_diagnostic, "linear_solver_trace": u_linear_solver_trace, }, - kernels=["gpu_ldu_matvec_vector_asymmetric_diag", "gpu_ldu_matvec_vector_asymmetric_face_accumulate", "gpu_bicgstab_initialize_vector", "gpu_bicgstab_dot_vector", "gpu_bicgstab_update_direction_vector", "gpu_bicgstab_precondition_vector", "gpu_dilu_apply_vector_asymmetric_faces", "gpu_bicgstab_update_intermediate_vector_preconditioned", "gpu_bicgstab_update_solution_residual_vector_preconditioned"], + kernels=["gpu_extract_vector_component", "gpu_scatter_vector_component", "gpu_ldu_matvec_vector_asymmetric_diag", "gpu_ldu_matvec_vector_asymmetric_face_accumulate", "gpu_bicgstab_initialize_vector", "gpu_bicgstab_dot_vector", "gpu_bicgstab_update_direction_vector", "gpu_bicgstab_precondition_vector", "gpu_dilu_apply_vector_asymmetric_faces", "gpu_dilu_forward_level_vector", "gpu_dilu_backward_level_vector", "gpu_bicgstab_update_intermediate_vector_preconditioned", "gpu_bicgstab_update_solution_residual_vector_preconditioned"], changed_fields=["U"], ), gpu_stage_result( @@ -1997,42 +2341,40 @@ def run_gpu_solver_stage_smoke(stepper: Any, backend: Mapping[str, Any], case: P "H1": gpu_array_output("H1", H1, kernel="gpu_rans_momentum_h1_finish"), "HbyA_model": { "mode": "assembled_UEqn_H_over_A_with_consistent_phiHbyA", - "kernel_entrypoints": [ + "kernels": [ "gpu_rans_pressure_inputs", - "gpu_rans_momentum_h1_face_accumulate", + "gpu_copy_scalar", "gpu_rans_momentum_h1_finish", "gpu_rans_consistent_rAtU", - "gpu_rans_momentum_hbyA_source", - "gpu_rans_momentum_hbyA_face_accumulate", + "gpu_rans_momentum_hbyA_by_cell", "gpu_rans_momentum_hbyA_finish", "gpu_rans_surface_flux_from_cells", "gpu_rans_consistent_phiHbyA_correction", ], - "formula": "solveDiag=UEqn.diag+non-coupled boundary internal coeffs; rAU=V/solveDiag; H1=(-UEqn.upper/lower neighbour sum)/V; consistent rAtU=1/max(1/rAU - H1, 0.1/rAU); HbyA=(UEqn.solve_rhs_source - UEqn.offdiag(U))/solveDiag for internal cells; phiHbyA=fvc::flux(HbyA)+interpolate(rAtU-rAU)*snGrad(p)*magSf", + "formula": "solveDiag=UEqn.diag+non-coupled boundary internal coeffs; rAU=V/solveDiag; H1=(-UEqn.upper/lower neighbour sum)/V; consistent rAtU=1/max(1/rAU - H1, 0.1/rAU); HbyA=(UEqn.solve_rhs_source - UEqn.offdiag(U))/solveDiag for internal cells; phiHbyA=OpenFOAM-weighted fvc::flux(HbyA)+OpenFOAM-weighted interpolate(rAtU-rAU)*snGrad(p)*magSf", }, }, kernels=[ "gpu_rans_pressure_inputs", - "gpu_rans_momentum_h1_face_accumulate", + "gpu_copy_scalar", "gpu_rans_momentum_h1_finish", "gpu_rans_consistent_rAtU", - "gpu_rans_momentum_hbyA_source", - "gpu_rans_momentum_hbyA_face_accumulate", + "gpu_rans_momentum_hbyA_by_cell", "gpu_rans_momentum_hbyA_finish", "gpu_rans_surface_flux_from_cells", "gpu_rans_consistent_phiHbyA_correction", ], ), - gpu_stage_result("assemble_pEqn", {"pEqn": pEqn}, kernels=["gpu_rans_pressure_assembly", "gpu_rans_pressure_laplacian_coefficients", "gpu_rans_pressure_mixed_boundary_laplacian", "gpu_rans_pressure_source_from_flux", "gpu_rans_pressure_source_from_boundary_flux"]), + gpu_stage_result("assemble_pEqn", {"pEqn": pEqn}, kernels=["gpu_rans_pressure_assembly", "gpu_rans_pressure_laplacian_coefficients", "gpu_rans_pressure_source_from_flux", "gpu_rans_pressure_source_from_boundary_flux"]), gpu_stage_result( "solve_pEqn", { "performance": {"solver_name": "gpu_symmetric_ldu_pcg", "field_name": "p", "operator_transform": "negative_openfoam_laplacian_to_spd", "non_orthogonal_loop": "OpenFOAM SIMPLE nNonOrthogonalCorrectors=3", **p_solve_performance}, - "p": gpu_array_output("p", p_solved, kernel="gpu_pcg_update_solution_residual_scalar"), - "residual": gpu_array_output("p_residual", p_residual, kernel="gpu_pcg_update_solution_residual_scalar"), + "p": gpu_array_output("p", p_solved, kernel="gpu_pcg_update_solution_residual_only_scalar"), + "residual": gpu_array_output("p_residual", p_residual, kernel="gpu_pcg_update_solution_residual_only_scalar"), "phi": gpu_array_output("phi", phi_solved, kernel="gpu_rans_pressure_flux_correction"), }, - kernels=["gpu_rans_negate_scalar_field", "gpu_ldu_matvec_scalar_symmetric_diag", "gpu_ldu_matvec_scalar_symmetric_face_accumulate", "gpu_pcg_initialize_scalar", "gpu_cg_dot_scalar", "gpu_pcg_update_solution_residual_scalar", "gpu_pcg_update_direction_scalar"], + kernels=["gpu_rans_negate_scalar_field", "gpu_rans_pressure_mixed_solver_diag", "gpu_ldu_matvec_scalar_symmetric_diag", "gpu_ldu_matvec_scalar_symmetric_face_accumulate", "gpu_pcg_initialize_residual_scalar", "gpu_dic_forward_level_scalar", "gpu_dic_backward_level_scalar", "gpu_cg_dot_scalar", "gpu_pcg_update_solution_residual_only_scalar", "gpu_pcg_update_direction_scalar"], changed_fields=["p", "phi"], ), gpu_stage_result( @@ -2138,7 +2480,6 @@ def run_gpu_solver_stage_smoke(stepper: Any, backend: Mapping[str, Any], case: P ) return { - "schema_version": 1, "status": "executed", "case": case, "backend": { @@ -2160,9 +2501,8 @@ def run_gpu_solver_stage_smoke(stepper: Any, backend: Mapping[str, Any], case: P "k": gpu_array_output("k", k_out, kernel="gpu_rans_turbulence_update"), "omega": gpu_array_output("omega", omega_out, kernel="gpu_rans_omega_wall_update"), }, - "field_objects": field_objects, "transfers": { - "inputs": [u_transfer, p_transfer, phi_transfer, nut_transfer, k_transfer, omega_transfer, owner_transfer, neighbour_transfer, sf_transfer, face_centres_transfer, cell_centres_transfer, cell_volumes_transfer, mag_sf_transfer, boundary_face_cells_transfer, boundary_phi_transfer, boundary_phiHbyA_transfer, momentum_pressure_boundary_face_cells_transfer, momentum_pressure_boundary_values_transfer, momentum_pressure_boundary_sf_transfer, momentum_u_boundary_face_cells_transfer, momentum_u_boundary_values_transfer, momentum_u_boundary_sf_transfer, momentum_u_convection_face_cells_transfer, momentum_u_convection_phi_transfer, momentum_u_convection_internal_coeff_transfer, momentum_u_convection_boundary_coeff_transfer, pressure_mixed_face_cells_transfer, pressure_mixed_scales_transfer, pressure_mixed_values_transfer, momentum_wall_face_cells_transfer, momentum_wall_values_transfer, momentum_wall_nut_transfer, momentum_wall_face_centres_transfer, momentum_wall_area_vectors_transfer, momentum_wall_area_magnitudes_transfer, omega_wall_cells_transfer, omega_wall_distances_transfer, u_preconditioner_diag_transfer], + "inputs": [u_transfer, p_transfer, phi_transfer, nut_transfer, k_transfer, omega_transfer, owner_transfer, neighbour_transfer, offdiag_cell_offsets_transfer, offdiag_cell_faces_transfer, offdiag_cell_sides_transfer, sf_transfer, face_centres_transfer, cell_centres_transfer, cell_volumes_transfer, mag_sf_transfer, face_weights_transfer, momentum_diffusion_coeff_transfer, momentum_internal_diag_transfer, momentum_internal_h1_transfer, boundary_face_cells_transfer, boundary_phi_transfer, boundary_phiHbyA_transfer, momentum_pressure_boundary_face_cells_transfer, momentum_pressure_boundary_values_transfer, momentum_pressure_boundary_sf_transfer, momentum_pressure_source_transfer, momentum_u_boundary_face_cells_transfer, momentum_u_boundary_values_transfer, momentum_u_boundary_sf_transfer, momentum_u_convection_face_cells_transfer, momentum_u_convection_phi_transfer, momentum_u_convection_internal_coeff_transfer, momentum_u_convection_boundary_coeff_transfer, pressure_mixed_face_cells_transfer, pressure_mixed_u_values_transfer, pressure_mixed_face_centres_transfer, pressure_mixed_area_vectors_transfer, pressure_mixed_area_magnitudes_transfer, momentum_wall_face_cells_transfer, momentum_wall_values_transfer, momentum_wall_nut_transfer, momentum_wall_face_centres_transfer, momentum_wall_area_vectors_transfer, momentum_wall_area_magnitudes_transfer, omega_wall_cells_transfer, omega_wall_distances_transfer, u_preconditioner_diag_transfer, u_dilu_reciprocal_diag_transfer, p_preconditioner_diag_transfer, p_dic_reciprocal_diag_transfer], "allocations": [ u_diag_alloc, u_source_alloc, diff --git a/python/src/foam_stepper/gpu/constants.py b/python/src/foam_stepper/gpu/constants.py index 5f4cfbb..32176a4 100644 --- a/python/src/foam_stepper/gpu/constants.py +++ b/python/src/foam_stepper/gpu/constants.py @@ -19,9 +19,13 @@ GPU_STAGE_CAPABILITY_PREFIX = "solver_stage_contract:" GPU_INPUT_SCHEMA_VERSION = 1 DEFAULT_LAMINAR_NU = 1.5e-5 DEFAULT_MOMENTUM_RELAXATION_ALPHA = 0.9 -DEFAULT_MOMENTUM_PBICGSTAB_ITERATIONS = 100 +DEFAULT_MOMENTUM_PBICGSTAB_ITERATIONS = 50 +DEFAULT_MOMENTUM_PBICGSTAB_RESIDUAL_TOLERANCE = 1.0e-8 DEFAULT_MOMENTUM_PBICGSTAB_RESIDUAL_TOLERANCE_SQUARED = 1.0e-16 -DEFAULT_PRESSURE_CG_ITERATIONS = 2000 +DEFAULT_MOMENTUM_PBICGSTAB_MIN_ITERATIONS = 3 +DEFAULT_PRESSURE_CG_ITERATIONS = 300 +DEFAULT_PRESSURE_CG_RESIDUAL_TOLERANCE = 1.0e-6 +DEFAULT_PRESSURE_CG_MIN_ITERATIONS = 10 DEFAULT_PRESSURE_NON_ORTHOGONAL_CORRECTORS = 3 DEFAULT_SIMPLE_CONSISTENT_RATU_FACTOR = 10.0 DEFAULT_OMEGA_WALL_BETA1 = 0.075 @@ -93,9 +97,9 @@ STAGE_OBSERVABILITY_GROUPS = ( ) GPU_STAGE_KERNEL_ENTRYPOINTS = { - "momentum_assembly": ["gpu_rans_momentum_assembly", "gpu_rans_momentum_diffusion_coefficients", "gpu_rans_momentum_wall_diffusion_coefficients", "gpu_rans_momentum_convection_coefficients", "gpu_rans_momentum_bounded_convection_sp_internal", "gpu_rans_momentum_bounded_convection_sp_boundary", "gpu_rans_momentum_convection_boundary_coefficients", "gpu_rans_momentum_convection_boundary_source", "gpu_rans_momentum_boundary_internal_diag", "gpu_rans_momentum_boundary_relaxation_coefficients", "gpu_rans_add_scalar_field", "gpu_rans_zero_tensor_field", "gpu_rans_momentum_gauss_grad_u_internal", "gpu_rans_momentum_gauss_grad_u_boundary", "gpu_rans_momentum_gauss_grad_u_finish", "gpu_rans_momentum_linear_upwind_source", "gpu_rans_zero_scalar_field", "gpu_rans_momentum_offdiag_abs_accumulate", "gpu_rans_momentum_equation_relaxation", "gpu_rans_momentum_pressure_gradient_source", "gpu_rans_momentum_pressure_boundary_source"], - "pressure_assembly": ["gpu_rans_pressure_inputs", "gpu_rans_momentum_h1_face_accumulate", "gpu_rans_momentum_h1_finish", "gpu_rans_consistent_rAtU", "gpu_rans_momentum_hbyA_source", "gpu_rans_momentum_hbyA_face_accumulate", "gpu_rans_momentum_hbyA_finish", "gpu_rans_surface_flux_from_cells", "gpu_rans_consistent_phiHbyA_correction", "gpu_rans_pressure_assembly", "gpu_rans_pressure_laplacian_coefficients", "gpu_rans_pressure_mixed_boundary_laplacian", "gpu_rans_pressure_source_from_flux", "gpu_rans_pressure_source_from_boundary_flux"], - "linear_solve_results": ["gpu_ldu_matvec_vector_asymmetric_diag", "gpu_ldu_matvec_vector_asymmetric_face_accumulate", "gpu_bicgstab_initialize_vector", "gpu_bicgstab_dot_vector", "gpu_bicgstab_update_direction_vector", "gpu_bicgstab_precondition_vector", "gpu_dilu_apply_vector_asymmetric_faces", "gpu_bicgstab_update_intermediate_vector_preconditioned", "gpu_bicgstab_update_solution_residual_vector_preconditioned", "gpu_vector_residual_squared", "gpu_rans_negate_scalar_field", "gpu_ldu_matvec_scalar_symmetric_diag", "gpu_ldu_matvec_scalar_symmetric_face_accumulate", "gpu_pcg_initialize_scalar", "gpu_cg_dot_scalar", "gpu_pcg_update_solution_residual_scalar", "gpu_pcg_update_direction_scalar"], + "momentum_assembly": ["gpu_rans_momentum_assembly", "gpu_rans_momentum_diffusion_coefficients_from_face_coeff", "gpu_rans_momentum_wall_diffusion_coefficients", "gpu_rans_momentum_convection_coefficients", "gpu_rans_momentum_bounded_convection_sp_internal", "gpu_rans_momentum_bounded_convection_sp_boundary", "gpu_rans_momentum_convection_boundary_coefficients", "gpu_rans_momentum_convection_boundary_source", "gpu_rans_momentum_boundary_internal_diag", "gpu_rans_momentum_boundary_relaxation_coefficients", "gpu_rans_add_scalar_field", "gpu_rans_add_vector_field", "gpu_rans_zero_tensor_field", "gpu_rans_momentum_gauss_grad_u_internal", "gpu_rans_momentum_gauss_grad_u_boundary", "gpu_rans_momentum_gauss_grad_u_finish", "gpu_rans_momentum_internal_dev_tau_source", "gpu_rans_momentum_wall_dev_tau_source", "gpu_rans_momentum_linear_upwind_source", "gpu_rans_zero_scalar_field", "gpu_rans_momentum_diag_from_offdiag_by_cell", "gpu_rans_momentum_offdiag_abs_accumulate_by_cell", "gpu_rans_momentum_equation_relaxation", "gpu_rans_momentum_pressure_gradient_source", "gpu_rans_momentum_pressure_boundary_source"], + "pressure_assembly": ["gpu_rans_pressure_inputs", "gpu_copy_scalar", "gpu_rans_momentum_h1_finish", "gpu_rans_consistent_rAtU", "gpu_rans_momentum_hbyA_by_cell", "gpu_rans_momentum_hbyA_finish", "gpu_rans_surface_flux_from_cells", "gpu_rans_consistent_phiHbyA_correction", "gpu_rans_pressure_assembly", "gpu_rans_pressure_laplacian_coefficients", "gpu_rans_pressure_source_from_flux", "gpu_rans_pressure_source_from_boundary_flux"], + "linear_solve_results": ["gpu_extract_vector_component", "gpu_scatter_vector_component", "gpu_ldu_matvec_vector_asymmetric_diag", "gpu_ldu_matvec_vector_asymmetric_face_accumulate", "gpu_bicgstab_initialize_vector", "gpu_bicgstab_dot_vector", "gpu_bicgstab_update_direction_vector", "gpu_bicgstab_precondition_vector", "gpu_dilu_apply_vector_asymmetric_faces", "gpu_dilu_forward_level_vector", "gpu_dilu_backward_level_vector", "gpu_bicgstab_update_intermediate_vector_preconditioned", "gpu_bicgstab_update_solution_residual_vector_preconditioned", "gpu_vector_residual_squared", "gpu_rans_negate_scalar_field", "gpu_rans_pressure_mixed_solver_diag", "gpu_ldu_matvec_scalar_symmetric_diag", "gpu_ldu_matvec_scalar_symmetric_face_accumulate", "gpu_pcg_initialize_scalar", "gpu_pcg_initialize_residual_scalar", "gpu_dic_forward_level_scalar", "gpu_dic_backward_level_scalar", "gpu_cg_dot_scalar", "gpu_pcg_update_solution_residual_scalar", "gpu_pcg_update_solution_residual_only_scalar", "gpu_pcg_update_direction_scalar"], "final_correction": ["gpu_rans_pressure_flux_correction", "gpu_rans_final_correction", "gpu_rans_pressure_velocity_correction"], "turbulence_updates": ["gpu_rans_turbulence_update", "gpu_rans_omega_wall_update"], } diff --git a/python/src/foam_stepper/gpu/kernels.py b/python/src/foam_stepper/gpu/kernels.py index 2000d80..296c22a 100644 --- a/python/src/foam_stepper/gpu/kernels.py +++ b/python/src/foam_stepper/gpu/kernels.py @@ -27,10 +27,11 @@ def gpu_rans_momentum_diffusion_coefficients( owner: qd.types.NDArray[qd.i32, 1], neighbour: qd.types.NDArray[qd.i32, 1], nut_internal: qd.types.NDArray[qd.f64, 1], - laminar_nu: float, + laminar_nu: qd.f64, cell_centres: qd.types.NDArray[qd.f64, 2], sf: qd.types.NDArray[qd.f64, 2], mag_sf: qd.types.NDArray[qd.f64, 1], + face_weights: qd.types.NDArray[qd.f64, 1], diag: qd.types.NDArray[qd.f64, 1], upper: qd.types.NDArray[qd.f64, 1], lower: qd.types.NDArray[qd.f64, 1], @@ -46,7 +47,8 @@ def gpu_rans_momentum_diffusion_coefficients( projected_delta = -projected_delta if projected_delta < 1.0e-300: projected_delta = 1.0e-300 - effective_nu = laminar_nu + 0.5 * (nut_internal[owner_cell] + nut_internal[neighbour_cell]) + owner_weight = face_weights[face] + effective_nu = laminar_nu + owner_weight * nut_internal[owner_cell] + (1.0 - owner_weight) * nut_internal[neighbour_cell] coeff = effective_nu * mag_sf[face] * mag_sf[face] / projected_delta upper[face] = -coeff lower[face] = -coeff @@ -55,13 +57,32 @@ def gpu_rans_momentum_diffusion_coefficients( +@qd.kernel +def gpu_rans_momentum_diffusion_coefficients_from_face_coeff( + n_internal_faces: int, + owner: qd.types.NDArray[qd.i32, 1], + neighbour: qd.types.NDArray[qd.i32, 1], + diffusion_coeff: qd.types.NDArray[qd.f64, 1], + diag: qd.types.NDArray[qd.f64, 1], + upper: qd.types.NDArray[qd.f64, 1], + lower: qd.types.NDArray[qd.f64, 1], +) -> None: + for face in range(n_internal_faces): + coeff = diffusion_coeff[face] + upper[face] = -coeff + lower[face] = -coeff + qd.atomic_add(diag[owner[face]], coeff) + qd.atomic_add(diag[neighbour[face]], coeff) + + + @qd.kernel def gpu_rans_momentum_wall_diffusion_coefficients( n_boundary_faces: int, face_cells: qd.types.NDArray[qd.i32, 1], boundary_values: qd.types.NDArray[qd.f64, 2], nut_boundary: qd.types.NDArray[qd.f64, 1], - laminar_nu: float, + laminar_nu: qd.f64, cell_centres: qd.types.NDArray[qd.f64, 2], face_centres: qd.types.NDArray[qd.f64, 2], face_area_vectors: qd.types.NDArray[qd.f64, 2], @@ -91,6 +112,18 @@ def gpu_rans_momentum_wall_diffusion_coefficients( qd.atomic_add(source[cell, 2], coeff * boundary_values[boundary_face, 2]) +@qd.kernel +def gpu_rans_add_vector_field( + n_cells: int, + source: qd.types.NDArray[qd.f64, 2], + out: qd.types.NDArray[qd.f64, 2], +) -> None: + for cell in range(n_cells): + out[cell, 0] += source[cell, 0] + out[cell, 1] += source[cell, 1] + out[cell, 2] += source[cell, 2] + + @qd.kernel def gpu_rans_momentum_convection_coefficients( @@ -263,19 +296,21 @@ def gpu_rans_momentum_gauss_grad_u_internal( neighbour: qd.types.NDArray[qd.i32, 1], u_internal: qd.types.NDArray[qd.f64, 2], sf: qd.types.NDArray[qd.f64, 2], + face_weights: qd.types.NDArray[qd.f64, 1], grad_u: qd.types.NDArray[qd.f64, 3], ) -> None: for face in range(n_internal_faces): owner_cell = owner[face] neighbour_cell = neighbour[face] + owner_weight = face_weights[face] + neighbour_weight = 1.0 - owner_weight for component in range(3): - face_value = 0.5 * (u_internal[owner_cell, component] + u_internal[neighbour_cell, component]) + face_value = owner_weight * u_internal[owner_cell, component] + neighbour_weight * u_internal[neighbour_cell, component] for direction in range(3): flux_value = face_value * sf[face, direction] qd.atomic_add(grad_u[owner_cell, component, direction], flux_value) qd.atomic_add(grad_u[neighbour_cell, component, direction], -flux_value) - @qd.kernel def gpu_rans_momentum_gauss_grad_u_boundary( n_boundary_faces: int, @@ -305,6 +340,141 @@ def gpu_rans_momentum_gauss_grad_u_finish( grad_u[cell, component, direction] *= inv_volume +@qd.kernel +def gpu_rans_momentum_internal_dev_tau_source( + n_internal_faces: int, + owner: qd.types.NDArray[qd.i32, 1], + neighbour: qd.types.NDArray[qd.i32, 1], + nut_internal: qd.types.NDArray[qd.f64, 1], + laminar_nu: qd.f64, + cell_centres: qd.types.NDArray[qd.f64, 2], + sf: qd.types.NDArray[qd.f64, 2], + mag_sf: qd.types.NDArray[qd.f64, 1], + face_weights: qd.types.NDArray[qd.f64, 1], + grad_u: qd.types.NDArray[qd.f64, 3], + source: qd.types.NDArray[qd.f64, 2], +) -> None: + for face in range(n_internal_faces): + owner_cell = owner[face] + neighbour_cell = neighbour[face] + owner_weight = face_weights[face] + neighbour_weight = 1.0 - owner_weight + gamma = laminar_nu + owner_weight * nut_internal[owner_cell] + neighbour_weight * nut_internal[neighbour_cell] + trace_owner = grad_u[owner_cell, 0, 0] + grad_u[owner_cell, 1, 1] + grad_u[owner_cell, 2, 2] + trace_neighbour = grad_u[neighbour_cell, 0, 0] + grad_u[neighbour_cell, 1, 1] + grad_u[neighbour_cell, 2, 2] + dot0 = sf[face, 0] * 0.0 + dot1 = sf[face, 1] * 0.0 + dot2 = sf[face, 2] * 0.0 + for row in range(3): + dev_owner0 = grad_u[owner_cell, row, 0] + dev_owner1 = grad_u[owner_cell, row, 1] + dev_owner2 = grad_u[owner_cell, row, 2] + dev_neighbour0 = grad_u[neighbour_cell, row, 0] + dev_neighbour1 = grad_u[neighbour_cell, row, 1] + dev_neighbour2 = grad_u[neighbour_cell, row, 2] + if row == 0: + dev_owner0 -= (2.0 / 3.0) * trace_owner + dev_neighbour0 -= (2.0 / 3.0) * trace_neighbour + elif row == 1: + dev_owner1 -= (2.0 / 3.0) * trace_owner + dev_neighbour1 -= (2.0 / 3.0) * trace_neighbour + else: + dev_owner2 -= (2.0 / 3.0) * trace_owner + dev_neighbour2 -= (2.0 / 3.0) * trace_neighbour + dev0 = owner_weight * dev_owner0 + neighbour_weight * dev_neighbour0 + dev1 = owner_weight * dev_owner1 + neighbour_weight * dev_neighbour1 + dev2 = owner_weight * dev_owner2 + neighbour_weight * dev_neighbour2 + area = sf[face, row] + dot0 += area * dev0 + dot1 += area * dev1 + dot2 += area * dev2 + qd.atomic_add(source[owner_cell, 0], gamma * dot0) + qd.atomic_add(source[owner_cell, 1], gamma * dot1) + qd.atomic_add(source[owner_cell, 2], gamma * dot2) + qd.atomic_add(source[neighbour_cell, 0], -gamma * dot0) + qd.atomic_add(source[neighbour_cell, 1], -gamma * dot1) + qd.atomic_add(source[neighbour_cell, 2], -gamma * dot2) + + + +@qd.kernel +def gpu_rans_momentum_wall_dev_tau_source( + n_boundary_faces: int, + face_cells: qd.types.NDArray[qd.i32, 1], + boundary_values: qd.types.NDArray[qd.f64, 2], + nut_boundary: qd.types.NDArray[qd.f64, 1], + laminar_nu: qd.f64, + u_internal: qd.types.NDArray[qd.f64, 2], + cell_centres: qd.types.NDArray[qd.f64, 2], + face_centres: qd.types.NDArray[qd.f64, 2], + face_area_vectors: qd.types.NDArray[qd.f64, 2], + face_area_magnitudes: qd.types.NDArray[qd.f64, 1], + grad_u: qd.types.NDArray[qd.f64, 3], + source: qd.types.NDArray[qd.f64, 2], +) -> None: + for face in range(n_boundary_faces): + cell = face_cells[face] + mag_sf = face_area_magnitudes[face] + nx = face_area_vectors[face, 0] / mag_sf + ny = face_area_vectors[face, 1] / mag_sf + nz = face_area_vectors[face, 2] / mag_sf + dx = face_centres[face, 0] - cell_centres[cell, 0] + dy = face_centres[face, 1] - cell_centres[cell, 1] + dz = face_centres[face, 2] - cell_centres[cell, 2] + projected_delta = dx * nx + dy * ny + dz * nz + if projected_delta < 0.0: + projected_delta = -projected_delta + if projected_delta < 1.0e-300: + projected_delta = 1.0e-300 + delta_coeff = 1.0 / projected_delta + gamma = laminar_nu + nut_boundary[face] + corrected_grad00 = grad_u[cell, 0, 0] + corrected_grad01 = grad_u[cell, 0, 1] + corrected_grad02 = grad_u[cell, 0, 2] + corrected_grad10 = grad_u[cell, 1, 0] + corrected_grad11 = grad_u[cell, 1, 1] + corrected_grad12 = grad_u[cell, 1, 2] + corrected_grad20 = grad_u[cell, 2, 0] + corrected_grad21 = grad_u[cell, 2, 1] + corrected_grad22 = grad_u[cell, 2, 2] + + sn_grad = (boundary_values[face, 0] - u_internal[cell, 0]) * delta_coeff + normal_grad = nx * corrected_grad00 + ny * corrected_grad01 + nz * corrected_grad02 + correction = sn_grad - normal_grad + corrected_grad00 += correction * nx + corrected_grad01 += correction * ny + corrected_grad02 += correction * nz + + sn_grad = (boundary_values[face, 1] - u_internal[cell, 1]) * delta_coeff + normal_grad = nx * corrected_grad10 + ny * corrected_grad11 + nz * corrected_grad12 + correction = sn_grad - normal_grad + corrected_grad10 += correction * nx + corrected_grad11 += correction * ny + corrected_grad12 += correction * nz + + sn_grad = (boundary_values[face, 2] - u_internal[cell, 2]) * delta_coeff + normal_grad = nx * corrected_grad20 + ny * corrected_grad21 + nz * corrected_grad22 + correction = sn_grad - normal_grad + corrected_grad20 += correction * nx + corrected_grad21 += correction * ny + corrected_grad22 += correction * nz + + trace = corrected_grad00 + corrected_grad11 + corrected_grad22 + dev00 = corrected_grad00 - (2.0 / 3.0) * trace + dev01 = corrected_grad01 + dev02 = corrected_grad02 + dev10 = corrected_grad10 + dev11 = corrected_grad11 - (2.0 / 3.0) * trace + dev12 = corrected_grad12 + dev20 = corrected_grad20 + dev21 = corrected_grad21 + dev22 = corrected_grad22 - (2.0 / 3.0) * trace + qd.atomic_add(source[cell, 0], gamma * (face_area_vectors[face, 0] * dev00 + face_area_vectors[face, 1] * dev10 + face_area_vectors[face, 2] * dev20)) + qd.atomic_add(source[cell, 1], gamma * (face_area_vectors[face, 0] * dev01 + face_area_vectors[face, 1] * dev11 + face_area_vectors[face, 2] * dev21)) + qd.atomic_add(source[cell, 2], gamma * (face_area_vectors[face, 0] * dev02 + face_area_vectors[face, 1] * dev12 + face_area_vectors[face, 2] * dev22)) + + + @qd.kernel def gpu_rans_momentum_linear_upwind_source( n_internal_faces: int, @@ -366,11 +536,78 @@ def gpu_rans_momentum_offdiag_abs_accumulate( qd.atomic_add(offdiag_sum[neighbour[face]], lower_abs) + +@qd.kernel +def gpu_rans_momentum_offdiag_abs_accumulate_by_cell( + n_cells: int, + cell_offsets: qd.types.NDArray[qd.i32, 1], + cell_faces: qd.types.NDArray[qd.i32, 1], + cell_sides: qd.types.NDArray[qd.i32, 1], + upper: qd.types.NDArray[qd.f64, 1], + lower: qd.types.NDArray[qd.f64, 1], + offdiag_sum: qd.types.NDArray[qd.f64, 1], +) -> None: + for cell in range(n_cells): + total = upper[0] * 0.0 + for slot in range(cell_offsets[cell], cell_offsets[cell + 1]): + face = cell_faces[slot] + coeff = upper[face] + if cell_sides[slot] != 0: + coeff = lower[face] + if coeff < 0.0: + coeff = -coeff + total += coeff + offdiag_sum[cell] = total + + +@qd.kernel +def gpu_rans_momentum_h1_by_cell( + n_cells: int, + cell_offsets: qd.types.NDArray[qd.i32, 1], + cell_faces: qd.types.NDArray[qd.i32, 1], + cell_sides: qd.types.NDArray[qd.i32, 1], + upper: qd.types.NDArray[qd.f64, 1], + lower: qd.types.NDArray[qd.f64, 1], + H1: qd.types.NDArray[qd.f64, 1], +) -> None: + for cell in range(n_cells): + total = upper[0] * 0.0 + for slot in range(cell_offsets[cell], cell_offsets[cell + 1]): + face = cell_faces[slot] + if cell_sides[slot] == 0: + total -= upper[face] + else: + total -= lower[face] + H1[cell] = total + + +@qd.kernel +def gpu_rans_momentum_diag_from_offdiag_by_cell( + n_cells: int, + cell_offsets: qd.types.NDArray[qd.i32, 1], + cell_faces: qd.types.NDArray[qd.i32, 1], + cell_sides: qd.types.NDArray[qd.i32, 1], + upper: qd.types.NDArray[qd.f64, 1], + lower: qd.types.NDArray[qd.f64, 1], + diag: qd.types.NDArray[qd.f64, 1], +) -> None: + for cell in range(n_cells): + total = upper[0] * 0.0 + for slot in range(cell_offsets[cell], cell_offsets[cell + 1]): + face = cell_faces[slot] + if cell_sides[slot] == 0: + total -= upper[face] + else: + total -= lower[face] + diag[cell] = total + + + @qd.kernel def gpu_rans_momentum_equation_relaxation( n_cells: int, u_internal: qd.types.NDArray[qd.f64, 2], - alpha: float, + alpha: qd.f64, offdiag_sum: qd.types.NDArray[qd.f64, 1], boundary_relax_add: qd.types.NDArray[qd.f64, 1], boundary_relax_subtract: qd.types.NDArray[qd.f64, 1], @@ -441,6 +678,38 @@ def gpu_rans_momentum_hbyA_source( HbyA[cell, 2] = source[cell, 2] + +@qd.kernel +def gpu_rans_momentum_hbyA_by_cell( + n_cells: int, + cell_offsets: qd.types.NDArray[qd.i32, 1], + cell_faces: qd.types.NDArray[qd.i32, 1], + cell_sides: qd.types.NDArray[qd.i32, 1], + owner: qd.types.NDArray[qd.i32, 1], + neighbour: qd.types.NDArray[qd.i32, 1], + upper: qd.types.NDArray[qd.f64, 1], + lower: qd.types.NDArray[qd.f64, 1], + u_internal: qd.types.NDArray[qd.f64, 2], + source: qd.types.NDArray[qd.f64, 2], + HbyA: qd.types.NDArray[qd.f64, 2], +) -> None: + for cell in range(n_cells): + total0 = source[cell, 0] + total1 = source[cell, 1] + total2 = source[cell, 2] + for slot in range(cell_offsets[cell], cell_offsets[cell + 1]): + face = cell_faces[slot] + neighbour_cell = neighbour[face] + coeff = upper[face] + if cell_sides[slot] == 1: + neighbour_cell = owner[face] + coeff = lower[face] + total0 += -coeff * u_internal[neighbour_cell, 0] + total1 += -coeff * u_internal[neighbour_cell, 1] + total2 += -coeff * u_internal[neighbour_cell, 2] + HbyA[cell, 0] = total0 + HbyA[cell, 1] = total1 + HbyA[cell, 2] = total2 @qd.kernel def gpu_rans_momentum_hbyA_face_accumulate( n_internal_faces: int, @@ -524,11 +793,12 @@ def gpu_rans_consistent_rAtU( n_cells: int, rAU: qd.types.NDArray[qd.f64, 1], H1: qd.types.NDArray[qd.f64, 1], + floor_factor: qd.f64, rAtU: qd.types.NDArray[qd.f64, 1], ) -> None: for cell in range(n_cells): inv_rAU = 1.0 / rAU[cell] - floor = 0.1 * inv_rAU + floor = inv_rAU / floor_factor denominator = inv_rAU - H1[cell] if denominator < floor: denominator = floor @@ -542,15 +812,18 @@ def gpu_rans_surface_flux_from_cells( neighbour: qd.types.NDArray[qd.i32, 1], cell_vector: qd.types.NDArray[qd.f64, 2], sf: qd.types.NDArray[qd.f64, 2], + face_weights: qd.types.NDArray[qd.f64, 1], out: qd.types.NDArray[qd.f64, 1], ) -> None: for face in range(n_internal_faces): owner_cell = owner[face] neighbour_cell = neighbour[face] - out[face] = 0.5 * ( - (cell_vector[owner_cell, 0] + cell_vector[neighbour_cell, 0]) * sf[face, 0] - + (cell_vector[owner_cell, 1] + cell_vector[neighbour_cell, 1]) * sf[face, 1] - + (cell_vector[owner_cell, 2] + cell_vector[neighbour_cell, 2]) * sf[face, 2] + owner_weight = face_weights[face] + neighbour_weight = 1.0 - owner_weight + out[face] = ( + (owner_weight * cell_vector[owner_cell, 0] + neighbour_weight * cell_vector[neighbour_cell, 0]) * sf[face, 0] + + (owner_weight * cell_vector[owner_cell, 1] + neighbour_weight * cell_vector[neighbour_cell, 1]) * sf[face, 1] + + (owner_weight * cell_vector[owner_cell, 2] + neighbour_weight * cell_vector[neighbour_cell, 2]) * sf[face, 2] ) @@ -566,6 +839,7 @@ def gpu_rans_consistent_phiHbyA_correction( cell_centres: qd.types.NDArray[qd.f64, 2], sf: qd.types.NDArray[qd.f64, 2], mag_sf: qd.types.NDArray[qd.f64, 1], + face_weights: qd.types.NDArray[qd.f64, 1], phiHbyA: qd.types.NDArray[qd.f64, 1], ) -> None: for face in range(n_internal_faces): @@ -579,9 +853,11 @@ def gpu_rans_consistent_phiHbyA_correction( projected_delta = -projected_delta if projected_delta < 1.0e-300: projected_delta = 1.0e-300 - interpolated_delta = 0.5 * ( - (rAtU[owner_cell] - rAU[owner_cell]) - + (rAtU[neighbour_cell] - rAU[neighbour_cell]) + owner_weight = face_weights[face] + neighbour_weight = 1.0 - owner_weight + interpolated_delta = ( + owner_weight * (rAtU[owner_cell] - rAU[owner_cell]) + + neighbour_weight * (rAtU[neighbour_cell] - rAU[neighbour_cell]) ) pressure_jump = p_internal[neighbour_cell] - p_internal[owner_cell] phiHbyA[face] += interpolated_delta * pressure_jump * mag_sf[face] * mag_sf[face] / projected_delta @@ -608,6 +884,7 @@ def gpu_rans_pressure_laplacian_coefficients( cell_centres: qd.types.NDArray[qd.f64, 2], sf: qd.types.NDArray[qd.f64, 2], mag_sf: qd.types.NDArray[qd.f64, 1], + face_weights: qd.types.NDArray[qd.f64, 1], diag: qd.types.NDArray[qd.f64, 1], upper: qd.types.NDArray[qd.f64, 1], ) -> None: @@ -622,7 +899,7 @@ def gpu_rans_pressure_laplacian_coefficients( projected_delta = -projected_delta if projected_delta < 1.0e-300: projected_delta = 1.0e-300 - coeff = 0.5 * (rAtU[owner_cell] + rAtU[neighbour_cell]) * mag_sf[face] * mag_sf[face] / projected_delta + coeff = (face_weights[face] * rAtU[owner_cell] + (1.0 - face_weights[face]) * rAtU[neighbour_cell]) * mag_sf[face] * mag_sf[face] / projected_delta upper[face] = coeff qd.atomic_add(diag[owner_cell], -coeff) qd.atomic_add(diag[neighbour_cell], -coeff) @@ -652,22 +929,44 @@ def gpu_rans_pressure_source_from_boundary_flux( for face in range(n_boundary_faces): qd.atomic_add(source[face_cells[face]], phi_boundary[face]) - @qd.kernel -def gpu_rans_pressure_mixed_boundary_laplacian( +def gpu_rans_pressure_mixed_solver_diag( n_boundary_faces: int, face_cells: qd.types.NDArray[qd.i32, 1], - boundary_scales: qd.types.NDArray[qd.f64, 1], - boundary_values: qd.types.NDArray[qd.f64, 1], + u_boundary: qd.types.NDArray[qd.f64, 2], + face_centres: qd.types.NDArray[qd.f64, 2], + face_area_vectors: qd.types.NDArray[qd.f64, 2], + face_area_magnitudes: qd.types.NDArray[qd.f64, 1], + cell_centres: qd.types.NDArray[qd.f64, 2], rAtU: qd.types.NDArray[qd.f64, 1], - diag: qd.types.NDArray[qd.f64, 1], - source: qd.types.NDArray[qd.f64, 1], + solve_diag: qd.types.NDArray[qd.f64, 1], ) -> None: - for face in range(n_boundary_faces): - cell = face_cells[face] - coeff = rAtU[cell] * boundary_scales[face] - qd.atomic_add(diag[cell], -coeff) - qd.atomic_add(source[cell], -coeff * boundary_values[face]) + for worker in range(1): + for face in range(n_boundary_faces): + cell = face_cells[face] + ux: qd.f64 = u_boundary[face, 0] + uy: qd.f64 = u_boundary[face, 1] + uz: qd.f64 = u_boundary[face, 2] + area_x: qd.f64 = face_area_vectors[face, 0] + area_y: qd.f64 = face_area_vectors[face, 1] + area_z: qd.f64 = face_area_vectors[face, 2] + mag_area: qd.f64 = face_area_magnitudes[face] + mag_u: qd.f64 = (ux * ux + uy * uy + uz * uz) ** 0.5 + value_fraction: qd.f64 = 0.5 + if mag_u > 1.0e-300 and mag_area > 1.0e-300: + normal_dot_u: qd.f64 = (ux * area_x + uy * area_y + uz * area_z) / mag_area + value_fraction = 0.5 + 0.5 * normal_dot_u / mag_u + dx = face_centres[face, 0] - cell_centres[cell, 0] + dy = face_centres[face, 1] - cell_centres[cell, 1] + dz = face_centres[face, 2] - cell_centres[cell, 2] + projected_delta = dx * area_x + dy * area_y + dz * area_z + if projected_delta < 0.0: + projected_delta = -projected_delta + if projected_delta < 1.0e-300: + projected_delta = 1.0e-300 + solve_diag[cell] += value_fraction * rAtU[cell] * mag_area * mag_area / projected_delta + + @qd.kernel @@ -775,8 +1074,8 @@ def gpu_rans_omega_wall_update( n_wall_cells: int, wall_cells: qd.types.NDArray[qd.i32, 1], wall_distances: qd.types.NDArray[qd.f64, 1], - laminar_nu: float, - beta1: float, + laminar_nu: qd.f64, + beta1: qd.f64, omega_out: qd.types.NDArray[qd.f64, 1], ) -> None: for index in range(n_wall_cells): @@ -804,10 +1103,13 @@ __all__ = [ "gpu_rans_momentum_gauss_grad_u_internal", "gpu_rans_momentum_gauss_grad_u_boundary", "gpu_rans_momentum_gauss_grad_u_finish", + "gpu_rans_momentum_internal_dev_tau_source", + "gpu_rans_momentum_wall_dev_tau_source", "gpu_rans_momentum_linear_upwind_source", "gpu_rans_momentum_equation_relaxation", "gpu_rans_momentum_pressure_gradient_source", "gpu_rans_momentum_hbyA_source", + "gpu_rans_momentum_hbyA_by_cell", "gpu_rans_momentum_hbyA_face_accumulate", "gpu_rans_momentum_hbyA_finish", "gpu_rans_pressure_inputs", @@ -819,7 +1121,7 @@ __all__ = [ "gpu_rans_pressure_laplacian_coefficients", "gpu_rans_pressure_source_from_flux", "gpu_rans_pressure_source_from_boundary_flux", - "gpu_rans_pressure_mixed_boundary_laplacian", + "gpu_rans_pressure_mixed_solver_diag", "gpu_rans_negate_scalar_field", "gpu_rans_face_flux_copy", "gpu_rans_surface_flux_from_cells", diff --git a/python/src/foam_stepper/gpu/linear_solve.py b/python/src/foam_stepper/gpu/linear_solve.py index 1a8f51f..c65c3e4 100644 --- a/python/src/foam_stepper/gpu/linear_solve.py +++ b/python/src/foam_stepper/gpu/linear_solve.py @@ -98,6 +98,44 @@ def gpu_zero_scalar_accumulator(accumulator: qd.types.NDArray[qd.f64, 1]) -> Non accumulator[0] = 0.0 +@qd.kernel +def gpu_scalar_abs_sum( + n_cells: int, + field: qd.types.NDArray[qd.f64, 1], + out: qd.types.NDArray[qd.f64, 1], +) -> None: + for worker in range(1): + total: qd.f64 = 0.0 + for cell in range(n_cells): + value = field[cell] + if value < 0.0: + value = -value + total += value + out[0] = total + + +@qd.kernel +def gpu_vector_abs_sum( + n_cells: int, + field: qd.types.NDArray[qd.f64, 2], + out: qd.types.NDArray[qd.f64, 1], +) -> None: + for worker in range(1): + total: qd.f64 = 0.0 + for cell in range(n_cells): + value0 = field[cell, 0] + value1 = field[cell, 1] + value2 = field[cell, 2] + if value0 < 0.0: + value0 = -value0 + if value1 < 0.0: + value1 = -value1 + if value2 < 0.0: + value2 = -value2 + total += value0 + value1 + value2 + out[0] = total + + @qd.kernel def gpu_ldu_matvec_scalar_symmetric_diag( n_cells: int, @@ -166,14 +204,17 @@ def gpu_cg_dot_scalar( right: qd.types.NDArray[qd.f64, 1], out: qd.types.NDArray[qd.f64, 1], ) -> None: - for cell in range(n_cells): - qd.atomic_add(out[0], left[cell] * right[cell]) + for worker in range(1): + total: qd.f64 = 0.0 + for cell in range(n_cells): + total += left[cell] * right[cell] + out[0] = total @qd.kernel def gpu_cg_update_solution_residual_scalar( n_cells: int, - alpha: float, + alpha: qd.f64, solution: qd.types.NDArray[qd.f64, 1], direction: qd.types.NDArray[qd.f64, 1], residual: qd.types.NDArray[qd.f64, 1], @@ -191,7 +232,7 @@ def gpu_cg_update_solution_residual_scalar( @qd.kernel def gpu_cg_update_direction_scalar( n_cells: int, - beta: float, + beta: qd.f64, residual: qd.types.NDArray[qd.f64, 1], direction: qd.types.NDArray[qd.f64, 1], ) -> None: @@ -229,11 +270,27 @@ def gpu_pcg_initialize_scalar( qd.atomic_add(residual_squared[0], cell_residual * cell_residual) qd.atomic_add(preconditioned_dot[0], cell_residual * z) +@qd.kernel +def gpu_pcg_initialize_residual_scalar( + n_cells: int, + source: qd.types.NDArray[qd.f64, 1], + operator_current: qd.types.NDArray[qd.f64, 1], + initial: qd.types.NDArray[qd.f64, 1], + solution: qd.types.NDArray[qd.f64, 1], + residual: qd.types.NDArray[qd.f64, 1], + residual_squared: qd.types.NDArray[qd.f64, 1], +) -> None: + for cell in range(n_cells): + solution[cell] = initial[cell] + cell_residual = source[cell] - operator_current[cell] + residual[cell] = cell_residual + qd.atomic_add(residual_squared[0], cell_residual * cell_residual) + @qd.kernel def gpu_pcg_update_solution_residual_scalar( n_cells: int, - alpha: float, + alpha: qd.f64, diag: qd.types.NDArray[qd.f64, 1], solution: qd.types.NDArray[qd.f64, 1], direction: qd.types.NDArray[qd.f64, 1], @@ -259,11 +316,28 @@ def gpu_pcg_update_solution_residual_scalar( qd.atomic_add(residual_squared[0], next_residual * next_residual) qd.atomic_add(preconditioned_dot[0], next_residual * z) +@qd.kernel +def gpu_pcg_update_solution_residual_only_scalar( + n_cells: int, + alpha: qd.f64, + solution: qd.types.NDArray[qd.f64, 1], + direction: qd.types.NDArray[qd.f64, 1], + residual: qd.types.NDArray[qd.f64, 1], + operator_direction: qd.types.NDArray[qd.f64, 1], + residual_squared: qd.types.NDArray[qd.f64, 1], +) -> None: + for cell in range(n_cells): + next_solution = solution[cell] + alpha * direction[cell] + next_residual = residual[cell] - alpha * operator_direction[cell] + solution[cell] = next_solution + residual[cell] = next_residual + qd.atomic_add(residual_squared[0], next_residual * next_residual) + @qd.kernel def gpu_pcg_update_direction_scalar( n_cells: int, - beta: float, + beta: qd.f64, preconditioned_residual: qd.types.NDArray[qd.f64, 1], direction: qd.types.NDArray[qd.f64, 1], ) -> None: @@ -347,24 +421,49 @@ def gpu_ldu_pcg_scalar_symmetric_faces( residual_squared: Any, denominator: Any, preconditioner_diag: Any | None = None, + preconditioner_reciprocal_diag: Any | None = None, + dic_schedule: GpuLduLevelSchedule | None = None, *, iterations: int, residual_tolerance_squared: float = 0.0, + residual_tolerance: float | None = None, + min_iterations: int = 0, + normalization_factor: float | None = None, + label: str | None = None, ) -> dict[str, Any]: - """Run GPU CG for a symmetric per-face LDU matrix with a diagonal preconditioner.""" + """Run GPU CG for a symmetric per-face LDU matrix with an optional DIC preconditioner.""" precond_diag = diag if preconditioner_diag is None else preconditioner_diag + use_dic_levels = dic_schedule is not None and preconditioner_reciprocal_diag is not None gpu_ldu_matvec_scalar_symmetric_faces(n_cells, n_internal_faces, owner, neighbour, upper, diag, initial, operator_work) gpu_zero_scalar_accumulator(residual_squared) gpu_zero_scalar_accumulator(denominator) - gpu_pcg_initialize_scalar(n_cells, source, operator_work, precond_diag, initial, out, residual, preconditioned_residual, direction, residual_squared, denominator) + if use_dic_levels: + gpu_pcg_initialize_residual_scalar(n_cells, source, operator_work, initial, out, residual, residual_squared) + gpu_dic_apply_scalar_symmetric_levels(dic_schedule, owner, neighbour, upper, preconditioner_reciprocal_diag, residual, preconditioned_residual) + gpu_copy_scalar(n_cells, preconditioned_residual, direction) + gpu_cg_dot_scalar(n_cells, residual, preconditioned_residual, denominator) + else: + gpu_pcg_initialize_scalar(n_cells, source, operator_work, precond_diag, initial, out, residual, preconditioned_residual, direction, residual_squared, denominator) qd.sync() rr_value = float(np.asarray(residual_squared.to_numpy())[0]) rho_value = float(np.asarray(denominator.to_numpy())[0]) performed_iterations = 0 + residual_l1_accumulator = qd.ndarray(qd.f64, shape=(1,)) if normalization_factor is not None and normalization_factor > 0.0 else None + normalized_residual = None + initial_normalized_residual = None + if residual_l1_accumulator is not None: + gpu_zero_scalar_accumulator(residual_l1_accumulator) + gpu_scalar_abs_sum(n_cells, residual, residual_l1_accumulator) + qd.sync() + normalized_residual = float(np.asarray(residual_l1_accumulator.to_numpy())[0]) / float(normalization_factor) + initial_normalized_residual = normalized_residual for _ in range(iterations): - if rr_value <= residual_tolerance_squared: + if residual_l1_accumulator is not None: + if performed_iterations >= min_iterations and normalized_residual is not None and normalized_residual <= float(residual_tolerance or 0.0): + break + elif performed_iterations >= min_iterations and rr_value <= residual_tolerance_squared: break if not np.isfinite(rho_value) or abs(rho_value) <= 1.0e-300: break @@ -378,12 +477,29 @@ def gpu_ldu_pcg_scalar_symmetric_faces( alpha = rho_value / denominator_value gpu_zero_scalar_accumulator(residual_squared) gpu_zero_scalar_accumulator(denominator) - gpu_pcg_update_solution_residual_scalar(n_cells, alpha, precond_diag, out, direction, residual, operator_work, preconditioned_residual, residual_squared, denominator) + if use_dic_levels: + gpu_pcg_update_solution_residual_only_scalar(n_cells, alpha, out, direction, residual, operator_work, residual_squared) + gpu_dic_apply_scalar_symmetric_levels(dic_schedule, owner, neighbour, upper, preconditioner_reciprocal_diag, residual, preconditioned_residual) + gpu_cg_dot_scalar(n_cells, residual, preconditioned_residual, denominator) + else: + gpu_pcg_update_solution_residual_scalar(n_cells, alpha, precond_diag, out, direction, residual, operator_work, preconditioned_residual, residual_squared, denominator) qd.sync() next_rr_value = float(np.asarray(residual_squared.to_numpy())[0]) next_rho_value = float(np.asarray(denominator.to_numpy())[0]) + next_normalized_residual = None + if residual_l1_accumulator is not None: + gpu_zero_scalar_accumulator(residual_l1_accumulator) + gpu_scalar_abs_sum(n_cells, residual, residual_l1_accumulator) + qd.sync() + next_normalized_residual = float(np.asarray(residual_l1_accumulator.to_numpy())[0]) / float(normalization_factor) performed_iterations += 1 - if next_rr_value <= residual_tolerance_squared: + if residual_l1_accumulator is not None: + if performed_iterations >= min_iterations and next_normalized_residual is not None and next_normalized_residual <= float(residual_tolerance or 0.0): + rr_value = next_rr_value + rho_value = next_rho_value + normalized_residual = next_normalized_residual + break + elif performed_iterations >= min_iterations and next_rr_value <= residual_tolerance_squared: rr_value = next_rr_value rho_value = next_rho_value break @@ -391,15 +507,31 @@ def gpu_ldu_pcg_scalar_symmetric_faces( gpu_pcg_update_direction_scalar(n_cells, beta, preconditioned_residual, direction) rr_value = next_rr_value rho_value = next_rho_value + normalized_residual = next_normalized_residual - return { + preconditioner_name = "diagonal_jacobi" if preconditioner_diag is None else "dic_reciprocal_diagonal" + if use_dic_levels: + preconditioner_name = "dic_level_scheduled" + report: dict[str, Any] = { "iterations": performed_iterations, "final_residual_squared": rr_value, "final_preconditioned_dot": rho_value, - "converged": rr_value <= residual_tolerance_squared, + "converged": ( + normalized_residual is not None and normalized_residual <= float(residual_tolerance or 0.0) + ) + if residual_l1_accumulator is not None + else rr_value <= residual_tolerance_squared, "residual_tolerance_squared": residual_tolerance_squared, - "preconditioner": "diagonal_jacobi" if preconditioner_diag is None else "dic_reciprocal_diagonal", + "residual_tolerance": residual_tolerance, + "initial_normalized_residual": initial_normalized_residual, + "final_normalized_residual": normalized_residual, + "normalization_factor": normalization_factor, + "min_iterations": min_iterations, + "preconditioner": preconditioner_name, } + if use_dic_levels and dic_schedule is not None: + report["preconditioner_schedule"] = dict(dic_schedule.metadata) + return report def gpu_ldu_jacobi_scalar_symmetric_faces( n_cells: int, n_internal_faces: int, @@ -466,6 +598,31 @@ def gpu_copy_vector( out[cell, 2] = source[cell, 2] +@qd.kernel +def gpu_extract_vector_component( + n_cells: int, + component: int, + source: qd.types.NDArray[qd.f64, 2], + out: qd.types.NDArray[qd.f64, 2], +) -> None: + for cell in range(n_cells): + out[cell, 0] = 0.0 + out[cell, 1] = 0.0 + out[cell, 2] = 0.0 + out[cell, component] = source[cell, component] + + +@qd.kernel +def gpu_scatter_vector_component( + n_cells: int, + component: int, + source: qd.types.NDArray[qd.f64, 2], + out: qd.types.NDArray[qd.f64, 2], +) -> None: + for cell in range(n_cells): + out[cell, component] = source[cell, component] + + @qd.kernel def gpu_ldu_matvec_vector_asymmetric_diag( @@ -561,11 +718,11 @@ def gpu_bicgstab_dot_vector( right: qd.types.NDArray[qd.f64, 2], out: qd.types.NDArray[qd.f64, 1], ) -> None: - for cell in range(n_cells): - qd.atomic_add( - out[0], - left[cell, 0] * right[cell, 0] + left[cell, 1] * right[cell, 1] + left[cell, 2] * right[cell, 2], - ) + for worker in range(1): + total: qd.f64 = 0.0 + for cell in range(n_cells): + total += left[cell, 0] * right[cell, 0] + left[cell, 1] * right[cell, 1] + left[cell, 2] * right[cell, 2] + out[0] = total @@ -573,8 +730,8 @@ def gpu_bicgstab_dot_vector( @qd.kernel def gpu_bicgstab_update_direction_vector( n_cells: int, - beta: float, - omega: float, + beta: qd.f64, + omega: qd.f64, residual: qd.types.NDArray[qd.f64, 2], direction: qd.types.NDArray[qd.f64, 2], operator_direction: qd.types.NDArray[qd.f64, 2], @@ -588,7 +745,7 @@ def gpu_bicgstab_update_direction_vector( @qd.kernel def gpu_bicgstab_update_intermediate_vector( n_cells: int, - alpha: float, + alpha: qd.f64, solution: qd.types.NDArray[qd.f64, 2], direction: qd.types.NDArray[qd.f64, 2], residual: qd.types.NDArray[qd.f64, 2], @@ -615,7 +772,7 @@ def gpu_bicgstab_update_intermediate_vector( @qd.kernel def gpu_bicgstab_update_solution_residual_vector( n_cells: int, - omega: float, + omega: qd.f64, solution: qd.types.NDArray[qd.f64, 2], intermediate: qd.types.NDArray[qd.f64, 2], operator_intermediate: qd.types.NDArray[qd.f64, 2], @@ -802,10 +959,98 @@ def gpu_dilu_apply_vector_asymmetric_levels( reciprocal_diag, out, ) + + +@qd.kernel +def gpu_dic_forward_level_scalar( + level: int, + level_offsets: qd.types.NDArray[qd.i32, 1], + level_cells: qd.types.NDArray[qd.i32, 1], + owner: qd.types.NDArray[qd.i32, 1], + incoming_offsets: qd.types.NDArray[qd.i32, 1], + incoming_faces: qd.types.NDArray[qd.i32, 1], + upper: qd.types.NDArray[qd.f64, 1], + reciprocal_diag: qd.types.NDArray[qd.f64, 1], + source: qd.types.NDArray[qd.f64, 1], + out: qd.types.NDArray[qd.f64, 1], +) -> None: + for index in range(level_offsets[level], level_offsets[level + 1]): + cell = level_cells[index] + scale = reciprocal_diag[cell] + value = scale * source[cell] + for face_slot in range(incoming_offsets[cell], incoming_offsets[cell + 1]): + face = incoming_faces[face_slot] + owner_cell = owner[face] + value -= scale * upper[face] * out[owner_cell] + out[cell] = value + + +@qd.kernel +def gpu_dic_backward_level_scalar( + level: int, + level_offsets: qd.types.NDArray[qd.i32, 1], + level_cells: qd.types.NDArray[qd.i32, 1], + neighbour: qd.types.NDArray[qd.i32, 1], + outgoing_offsets: qd.types.NDArray[qd.i32, 1], + outgoing_faces: qd.types.NDArray[qd.i32, 1], + upper: qd.types.NDArray[qd.f64, 1], + reciprocal_diag: qd.types.NDArray[qd.f64, 1], + out: qd.types.NDArray[qd.f64, 1], +) -> None: + for index in range(level_offsets[level], level_offsets[level + 1]): + cell = level_cells[index] + scale = reciprocal_diag[cell] + value = out[cell] + begin = outgoing_offsets[cell] + end = outgoing_offsets[cell + 1] + for reverse_slot in range(end - begin): + face = outgoing_faces[end - 1 - reverse_slot] + neighbour_cell = neighbour[face] + value -= scale * upper[face] * out[neighbour_cell] + out[cell] = value + + +def gpu_dic_apply_scalar_symmetric_levels( + schedule: GpuLduLevelSchedule, + owner: Any, + neighbour: Any, + upper: Any, + reciprocal_diag: Any, + source: Any, + out: Any, +) -> None: + """Apply OpenFOAM DIC using parallel cell work within each dependency level.""" + + for level in range(schedule.n_levels): + gpu_dic_forward_level_scalar( + level, + schedule.level_offsets, + schedule.level_cells, + owner, + schedule.incoming_offsets, + schedule.incoming_faces, + upper, + reciprocal_diag, + source, + out, + ) + for reverse_level in range(schedule.n_levels): + level = schedule.n_levels - 1 - reverse_level + gpu_dic_backward_level_scalar( + level, + schedule.level_offsets, + schedule.level_cells, + neighbour, + schedule.outgoing_offsets, + schedule.outgoing_faces, + upper, + reciprocal_diag, + out, + ) @qd.kernel def gpu_bicgstab_update_intermediate_vector_preconditioned( n_cells: int, - alpha: float, + alpha: qd.f64, solution: qd.types.NDArray[qd.f64, 2], preconditioned_direction: qd.types.NDArray[qd.f64, 2], residual: qd.types.NDArray[qd.f64, 2], @@ -832,7 +1077,7 @@ def gpu_bicgstab_update_intermediate_vector_preconditioned( @qd.kernel def gpu_bicgstab_update_solution_residual_vector_preconditioned( n_cells: int, - omega: float, + omega: qd.f64, solution: qd.types.NDArray[qd.f64, 2], preconditioned_intermediate: qd.types.NDArray[qd.f64, 2], intermediate: qd.types.NDArray[qd.f64, 2], @@ -1193,11 +1438,16 @@ def gpu_ldu_pbicgstab_vector_asymmetric_faces( omega_numerator: Any, omega_denominator: Any, preconditioner_diag: Any | None = None, + preconditioner_reciprocal_diag: Any | None = None, + dilu_schedule: GpuLduLevelSchedule | None = None, *, iterations: int, residual_tolerance_squared: float = 0.0, + residual_tolerance: float | None = None, + min_iterations: int = 0, + normalization_factor: float | None = None, ) -> dict[str, Any]: - """Run GPU BiCGStab for an asymmetric vector LDU matrix with a diagonal preconditioner.""" + """Run GPU BiCGStab for an asymmetric vector LDU matrix with an optional DILU preconditioner.""" gpu_ldu_matvec_vector_asymmetric_faces( n_cells, @@ -1230,9 +1480,22 @@ def gpu_ldu_pbicgstab_vector_asymmetric_faces( omega = 1.0 performed_iterations = 0 precond_diag = diag if preconditioner_diag is None else preconditioner_diag + use_dilu_levels = dilu_schedule is not None and preconditioner_reciprocal_diag is not None + residual_l1_accumulator = qd.ndarray(qd.f64, shape=(1,)) if normalization_factor is not None and normalization_factor > 0.0 else None + normalized_residual = None + initial_normalized_residual = None + if residual_l1_accumulator is not None: + gpu_zero_scalar_accumulator(residual_l1_accumulator) + gpu_vector_abs_sum(n_cells, residual, residual_l1_accumulator) + qd.sync() + normalized_residual = float(np.asarray(residual_l1_accumulator.to_numpy())[0]) / float(normalization_factor) + initial_normalized_residual = normalized_residual for _ in range(iterations): - if residual_value <= residual_tolerance_squared: + if residual_l1_accumulator is not None: + if performed_iterations >= min_iterations and normalized_residual is not None and normalized_residual <= float(residual_tolerance or 0.0): + break + elif performed_iterations >= min_iterations and residual_value <= residual_tolerance_squared: break gpu_zero_scalar_accumulator(rho) @@ -1250,7 +1513,19 @@ def gpu_ldu_pbicgstab_vector_asymmetric_faces( beta = (rho_new / rho_old) * (alpha / omega) gpu_bicgstab_update_direction_vector(n_cells, beta, omega, residual, direction, operator_direction) - gpu_bicgstab_precondition_vector(n_cells, precond_diag, direction, intermediate) + if use_dilu_levels: + gpu_dilu_apply_vector_asymmetric_levels( + dilu_schedule, + owner, + neighbour, + upper, + lower, + preconditioner_reciprocal_diag, + direction, + intermediate, + ) + else: + gpu_bicgstab_precondition_vector(n_cells, precond_diag, direction, intermediate) gpu_ldu_matvec_vector_asymmetric_faces( n_cells, n_internal_faces, @@ -1284,12 +1559,35 @@ def gpu_ldu_pbicgstab_vector_asymmetric_faces( ) qd.sync() intermediate_residual = float(np.asarray(residual_squared.to_numpy())[0]) + intermediate_normalized_residual = None + if residual_l1_accumulator is not None: + gpu_zero_scalar_accumulator(residual_l1_accumulator) + gpu_vector_abs_sum(n_cells, intermediate, residual_l1_accumulator) + qd.sync() + intermediate_normalized_residual = float(np.asarray(residual_l1_accumulator.to_numpy())[0]) / float(normalization_factor) performed_iterations += 1 - if intermediate_residual <= residual_tolerance_squared: + if residual_l1_accumulator is not None: + if performed_iterations >= min_iterations and intermediate_normalized_residual is not None and intermediate_normalized_residual <= float(residual_tolerance or 0.0): + residual_value = intermediate_residual + normalized_residual = intermediate_normalized_residual + break + elif performed_iterations >= min_iterations and intermediate_residual <= residual_tolerance_squared: residual_value = intermediate_residual break - gpu_bicgstab_precondition_vector(n_cells, precond_diag, intermediate, residual) + if use_dilu_levels: + gpu_dilu_apply_vector_asymmetric_levels( + dilu_schedule, + owner, + neighbour, + upper, + lower, + preconditioner_reciprocal_diag, + intermediate, + residual, + ) + else: + gpu_bicgstab_precondition_vector(n_cells, precond_diag, intermediate, residual) gpu_ldu_matvec_vector_asymmetric_faces(n_cells, n_internal_faces, owner, neighbour, upper, lower, diag, residual, operator_intermediate) gpu_zero_scalar_accumulator(omega_numerator) gpu_zero_scalar_accumulator(omega_denominator) @@ -1315,14 +1613,137 @@ def gpu_ldu_pbicgstab_vector_asymmetric_faces( ) qd.sync() residual_value = float(np.asarray(residual_squared.to_numpy())[0]) + if residual_l1_accumulator is not None: + gpu_zero_scalar_accumulator(residual_l1_accumulator) + gpu_vector_abs_sum(n_cells, residual, residual_l1_accumulator) + qd.sync() + normalized_residual = float(np.asarray(residual_l1_accumulator.to_numpy())[0]) / float(normalization_factor) rho_old = rho_new - return { + preconditioner_name = "diagonal_jacobi" if preconditioner_diag is None else "dilu_reciprocal_diagonal" + if use_dilu_levels: + preconditioner_name = "dilu_level_scheduled" + report: dict[str, Any] = { "iterations": performed_iterations, "final_residual_squared": residual_value, - "converged": residual_value <= residual_tolerance_squared, + "converged": ( + normalized_residual is not None and normalized_residual <= float(residual_tolerance or 0.0) + ) + if residual_l1_accumulator is not None + else residual_value <= residual_tolerance_squared, "residual_tolerance_squared": residual_tolerance_squared, - "preconditioner": "diagonal_jacobi" if preconditioner_diag is None else "dilu_reciprocal_diagonal", + "residual_tolerance": residual_tolerance, + "initial_normalized_residual": initial_normalized_residual, + "final_normalized_residual": normalized_residual, + "normalization_factor": normalization_factor, + "min_iterations": min_iterations, + "preconditioner": preconditioner_name, + } + if use_dilu_levels and dilu_schedule is not None: + report["preconditioner_schedule"] = dict(dilu_schedule.metadata) + return report + + +def gpu_ldu_pbicgstab_vector_asymmetric_components( + n_cells: int, + n_internal_faces: int, + owner: Any, + neighbour: Any, + upper: Any, + lower: Any, + diag: Any, + source: Any, + initial: Any, + out: Any, + residual: Any, + shadow: Any, + direction: Any, + operator_direction: Any, + intermediate: Any, + operator_intermediate: Any, + residual_squared: Any, + rho: Any, + denominator: Any, + omega_numerator: Any, + omega_denominator: Any, + preconditioner_diag: Any | None = None, + preconditioner_reciprocal_diag: Any | None = None, + dilu_schedule: GpuLduLevelSchedule | None = None, + *, + iterations: int, + residual_tolerance_squared: float = 0.0, + residual_tolerance: float | None = None, + min_iterations: int = 0, + normalization_factors: Any | None = None, +) -> dict[str, Any]: + """Run OpenFOAM-style component-wise PBiCGStab for a vector LDU matrix.""" + + component_source = qd.ndarray(qd.f64, shape=(n_cells, 3)) + component_initial = qd.ndarray(qd.f64, shape=(n_cells, 3)) + component_out = qd.ndarray(qd.f64, shape=(n_cells, 3)) + component_residual = qd.ndarray(qd.f64, shape=(n_cells, 3)) + component_shadow = qd.ndarray(qd.f64, shape=(n_cells, 3)) + component_direction = qd.ndarray(qd.f64, shape=(n_cells, 3)) + component_operator_direction = qd.ndarray(qd.f64, shape=(n_cells, 3)) + component_intermediate = qd.ndarray(qd.f64, shape=(n_cells, 3)) + component_operator_intermediate = qd.ndarray(qd.f64, shape=(n_cells, 3)) + component_reports: list[dict[str, Any]] = [] + for component in range(3): + gpu_extract_vector_component(n_cells, component, source, component_source) + gpu_extract_vector_component(n_cells, component, initial, component_initial) + normalization_factor = None + if normalization_factors is not None: + normalization_factor = float(normalization_factors[component]) + report = gpu_ldu_pbicgstab_vector_asymmetric_faces( + n_cells, + n_internal_faces, + owner, + neighbour, + upper, + lower, + diag, + component_source, + component_initial, + component_out, + component_residual, + component_shadow, + component_direction, + component_operator_direction, + component_intermediate, + component_operator_intermediate, + residual_squared, + rho, + denominator, + omega_numerator, + omega_denominator, + preconditioner_diag=preconditioner_diag, + preconditioner_reciprocal_diag=preconditioner_reciprocal_diag, + dilu_schedule=dilu_schedule, + iterations=iterations, + residual_tolerance_squared=residual_tolerance_squared, + residual_tolerance=residual_tolerance, + min_iterations=min_iterations, + normalization_factor=normalization_factor, + ) + component_reports.append(report) + gpu_scatter_vector_component(n_cells, component, component_out, out) + gpu_scatter_vector_component(n_cells, component, component_residual, residual) + qd.sync() + final_residual_squared = sum(float(report.get("final_residual_squared", 0.0) or 0.0) for report in component_reports) + return { + "iterations": max(int(report.get("iterations", 0) or 0) for report in component_reports), + "component_iterations": [int(report.get("iterations", 0) or 0) for report in component_reports], + "final_residual_squared": final_residual_squared, + "component_final_residual_squared": [float(report.get("final_residual_squared", 0.0) or 0.0) for report in component_reports], + "converged": all(bool(report.get("converged")) for report in component_reports), + "residual_tolerance_squared": residual_tolerance_squared, + "residual_tolerance": residual_tolerance, + "component_initial_normalized_residual": [report.get("initial_normalized_residual") for report in component_reports], + "component_final_normalized_residual": [report.get("final_normalized_residual") for report in component_reports], + "component_normalization_factor": [report.get("normalization_factor") for report in component_reports], + "min_iterations": min_iterations, + "preconditioner": component_reports[0].get("preconditioner") if component_reports else None, + "preconditioner_schedule": component_reports[0].get("preconditioner_schedule") if component_reports else None, } def gpu_ldu_jacobi_vector_symmetric_faces( @@ -1422,6 +1843,25 @@ def build_losort_addr(n_cells: int, neighbour: np.ndarray) -> np.ndarray: cursor[int(neighbour_cell)] += 1 return np.ascontiguousarray(losort, dtype=np.int32) +def _build_cell_face_offsets(n_cells: int, face_cells: np.ndarray) -> tuple[np.ndarray, np.ndarray]: + cells = np.asarray(face_cells, dtype=np.int32).reshape(-1) + if cells.size == 0: + return np.zeros(n_cells + 1, dtype=np.int32), np.zeros(0, dtype=np.int32) + if int(cells.min()) < 0 or int(cells.max()) >= n_cells: + raise ValueError("face-cell addresses exceed cell range") + + counts = np.bincount(cells, minlength=n_cells).astype(np.int32, copy=False) + offsets = np.empty(n_cells + 1, dtype=np.int32) + offsets[0] = 0 + np.cumsum(counts, out=offsets[1:]) + faces = np.empty(int(offsets[-1]), dtype=np.int32) + cursor = offsets[:-1].copy() + for face, cell in enumerate(cells): + slot = int(cursor[int(cell)]) + faces[slot] = int(face) + cursor[int(cell)] += 1 + return offsets, np.ascontiguousarray(faces, dtype=np.int32) + def _gpu_i32(values: np.ndarray) -> Any: host = np.ascontiguousarray(values.astype(np.int32, copy=False)) @@ -1440,6 +1880,74 @@ def _gpu_f64(values: np.ndarray) -> Any: gpu.from_numpy(host) return gpu +def build_ldu_level_schedule(n_cells: int, owner: np.ndarray, neighbour: np.ndarray, *, name: str) -> GpuLduLevelSchedule: + """Build a GPU level schedule for DILU/DIC triangular preconditioner sweeps.""" + + owner_i32 = np.asarray(owner, dtype=np.int32).reshape(-1) + neighbour_i32 = np.asarray(neighbour, dtype=np.int32).reshape(-1) + if owner_i32.shape != neighbour_i32.shape: + raise ValueError(f"owner/neighbour shape mismatch: {owner_i32.shape} != {neighbour_i32.shape}") + if owner_i32.size and (int(owner_i32.min()) < 0 or int(neighbour_i32.min()) < 0 or int(max(owner_i32.max(), neighbour_i32.max())) >= n_cells): + raise ValueError("owner/neighbour addresses exceed cell range") + + incoming_offsets, incoming_faces = _build_cell_face_offsets(n_cells, neighbour_i32) + outgoing_offsets, outgoing_faces = _build_cell_face_offsets(n_cells, owner_i32) + in_degree = np.diff(incoming_offsets).astype(np.int32, copy=True) + levels = np.zeros(n_cells, dtype=np.int32) + ready = [int(cell) for cell in np.flatnonzero(in_degree == 0)] + cursor = 0 + visited = 0 + while cursor < len(ready): + cell = ready[cursor] + cursor += 1 + visited += 1 + next_level = int(levels[cell]) + 1 + for slot in range(int(outgoing_offsets[cell]), int(outgoing_offsets[cell + 1])): + face = int(outgoing_faces[slot]) + neighbour_cell = int(neighbour_i32[face]) + if levels[neighbour_cell] < next_level: + levels[neighbour_cell] = next_level + in_degree[neighbour_cell] -= 1 + if in_degree[neighbour_cell] == 0: + ready.append(neighbour_cell) + if visited != n_cells: + raise ValueError("owner/neighbour graph contains a dependency cycle") + + n_levels = int(levels.max()) + 1 if n_cells else 0 + level_counts = np.bincount(levels, minlength=n_levels).astype(np.int32, copy=False) if n_levels else np.zeros(0, dtype=np.int32) + level_offsets = np.empty(n_levels + 1, dtype=np.int32) + level_offsets[0] = 0 + np.cumsum(level_counts, out=level_offsets[1:]) + level_cells = np.empty(n_cells, dtype=np.int32) + level_cursor = level_offsets[:-1].copy() + for cell, level in enumerate(levels): + slot = int(level_cursor[int(level)]) + level_cells[slot] = int(cell) + level_cursor[int(level)] += 1 + + non_empty_level_sizes = level_counts[level_counts > 0] + metadata = { + "name": name, + "format": "ldu_dependency_levels", + "n_cells": int(n_cells), + "n_internal_faces": int(owner_i32.size), + "n_levels": int(n_levels), + "min_level_size": int(non_empty_level_sizes.min()) if non_empty_level_sizes.size else 0, + "max_level_size": int(non_empty_level_sizes.max()) if non_empty_level_sizes.size else 0, + "mean_level_size": float(non_empty_level_sizes.mean()) if non_empty_level_sizes.size else 0.0, + "gpu_backed": True, + } + return GpuLduLevelSchedule( + level_offsets=_gpu_i32(level_offsets), + level_cells=_gpu_i32(np.ascontiguousarray(level_cells, dtype=np.int32)), + incoming_offsets=_gpu_i32(incoming_offsets), + incoming_faces=_gpu_i32(incoming_faces), + outgoing_offsets=_gpu_i32(outgoing_offsets), + outgoing_faces=_gpu_i32(outgoing_faces), + n_levels=n_levels, + metadata=metadata, + ) + def ldu_csr_to_gpu(offsets: np.ndarray, columns: np.ndarray, coefficients: np.ndarray, *, name: str) -> GpuLduCsr: return GpuLduCsr( @@ -1468,8 +1976,10 @@ def empty_ldu_csr_gpu(n_cells: int, name: str) -> GpuLduCsr: __all__ = [ "GpuLduCsr", + "GpuLduLevelSchedule", "build_ldu_csr", "build_losort_addr", + "build_ldu_level_schedule", "empty_ldu_csr_gpu", "gpu_bicgstab_dot_vector", "gpu_bicgstab_initialize_vector", @@ -1481,8 +1991,17 @@ __all__ = [ "gpu_bicgstab_update_solution_residual_vector_preconditioned", "gpu_copy_scalar", "gpu_copy_vector", + "gpu_extract_vector_component", + "gpu_scatter_vector_component", "gpu_dilu_apply_vector_asymmetric_faces", + "gpu_dilu_apply_vector_asymmetric_levels", + "gpu_dilu_forward_level_vector", + "gpu_dilu_backward_level_vector", + "gpu_dic_apply_scalar_symmetric_levels", + "gpu_dic_forward_level_scalar", + "gpu_dic_backward_level_scalar", "gpu_ldu_bicgstab_vector_asymmetric_faces", + "gpu_ldu_pbicgstab_vector_asymmetric_components", "gpu_ldu_pbicgstab_vector_asymmetric_faces", "gpu_ldu_cg_scalar_symmetric_faces", "gpu_ldu_pcg_scalar_symmetric_faces", @@ -1490,8 +2009,10 @@ __all__ = [ "gpu_ldu_jacobi_scalar_symmetric_accumulate", "gpu_ldu_jacobi_scalar_symmetric_faces", "gpu_pcg_initialize_scalar", + "gpu_pcg_initialize_residual_scalar", "gpu_pcg_update_direction_scalar", "gpu_pcg_update_solution_residual_scalar", + "gpu_pcg_update_solution_residual_only_scalar", "gpu_ldu_jacobi_vector", "gpu_ldu_jacobi_vector_asymmetric_accumulate", "gpu_ldu_jacobi_vector_asymmetric_faces",