From e7064226de8fbd28cd934b4c65802a192d5f17b7 Mon Sep 17 00:00:00 2001 From: Adalberto Date: Mon, 5 Jan 2026 15:24:16 +0100 Subject: [PATCH] Example script for hpod --- scripts/2-rom/5-hpod/hpod_run.py | 321 ++++++++++++++++++ scripts/2-rom/5-hpod/inputs.json | 45 +++ .../2-rom/5-hpod/meshgrid_data/cylinder.py | 271 +++++++++++++++ scripts/2-rom/5-hpod/sem_run.py | 184 ++++++++++ scripts/2-rom/5-hpod/visualize.py | 149 ++++++++ 5 files changed, 970 insertions(+) create mode 100644 scripts/2-rom/5-hpod/hpod_run.py create mode 100644 scripts/2-rom/5-hpod/inputs.json create mode 100644 scripts/2-rom/5-hpod/meshgrid_data/cylinder.py create mode 100644 scripts/2-rom/5-hpod/sem_run.py create mode 100644 scripts/2-rom/5-hpod/visualize.py diff --git a/scripts/2-rom/5-hpod/hpod_run.py b/scripts/2-rom/5-hpod/hpod_run.py new file mode 100644 index 0000000..705af56 --- /dev/null +++ b/scripts/2-rom/5-hpod/hpod_run.py @@ -0,0 +1,321 @@ +# Import general modules +import sys +import os +os.environ["OMP_NUM_THREADS"] = "1" +os.environ["OPENBLAS_NUM_THREADS"] = "1" +os.environ["MKL_NUM_THREADS"] = "1" +os.environ["VECLIB_MAXIMUM_THREADS"] = "1" +os.environ["NUMEXPR_NUM_THREADS"] = "1" + +import json +import numpy as np +import matplotlib.pyplot as plt + +def _mirror_right(a, axis=0): + """ + Mirror 'a' about its last sample along `axis` (one-sided). + """ + a = np.asarray(a) + N = a.shape[axis] + if N < 2: + raise ValueError("Need at least 2 samples along axis to mirror.") + sl = [slice(None)] * a.ndim + sl[axis] = slice(N-2, 0, -1) # N-2, ..., 1 + a_m = a[tuple(sl)] + return np.concatenate([a, a_m], axis=axis) + +def _hilbert_fft(a, axis=0): + """ + Analytic signal via FFT along `axis`. + """ + a = np.asarray(a) + M = a.shape[axis] + A = np.fft.fft(a, axis=axis) + + # Build the kernel in fourier space + h = np.zeros(M, dtype=float) + if M % 2 == 0: + h[0] = 1.0 + h[M//2] = 1.0 + h[1:M//2] = 2.0 + else: + h[0] = 1.0 + h[1:(M+1)//2] = 2.0 + + shape = [1] * a.ndim + shape[axis] = M + H = h.reshape(shape) + + # Perform the convolution (multiplication in fourier space) + z = np.fft.ifft(A * H, axis=axis) + return z + +def hilbert_axis0_mirror_right(a, axis=0, return_analytic=False): + """ + 1) Extend by mirroring around the last point (right side) + 2) Hilbert (via FFT) along `axis` + 3) Return only the first N samples (original length) + + If return_analytic=False: returns Hilbert transform (real, same dtype as input float) + If return_analytic=True: returns analytic signal (complex) + """ + a = np.asarray(a) + N = a.shape[axis] + + a_ext = _mirror_right(a, axis=axis) # length 2N-2 + z_ext = _hilbert_fft(a_ext, axis=axis) # analytic signal on extended + + # slice back to original length + sl = [slice(None)] * a.ndim + sl[axis] = slice(0, N) + z = z_ext[tuple(sl)] + + if return_analytic: + return z + return np.imag(z) + +if __name__ == "__main__": + + # Import MPI + from mpi4py import MPI #equivalent to the use of MPI_init() in C + # Import IO helper functions + from pysemtools.io.utils import IoPathData, get_fld_from_ndarray + # Import modules for reading and writing + from pysemtools.io.ppymech import pynekread, pynekwrite + # Data types + from pysemtools.datatypes import FieldRegistry, Mesh, Coef + # POD + from pysemtools.rom import POD, IoHelp + # Interpolation + from pysemtools.interpolation import Probes + from pysemtools.io.wrappers import read_data, write_data + + # Split communicator for MPI - MPMD + worldcomm = MPI.COMM_WORLD + worldrank = worldcomm.Get_rank() + worldsize = worldcomm.Get_size() + col = 1 + comm = worldcomm.Split(col,worldrank) + rank = comm.Get_rank() + size = comm.Get_size() + + #========================================= + # Define inputs + #========================================= + + # Open input file to see path + f = open ("inputs.json", "r") + params_file = json.loads(f.read()) + f.close() + + # Read the POD inputs + pod_number_of_snapshots = params_file["number_of_snapshots"] + pod_fields = params_file["fields"] + number_of_pod_fields = len(pod_fields) + pod_batch_size = params_file["batch_size"] + pod_keep_modes = params_file["keep_modes"] + pod_write_modes = params_file["write_modes"] + hpod_flag = params_file["hpod"] + dtype_string = params_file.get("dtype", "double") + print(dtype_string) + backend = params_file.get("backend", "numpy") + if dtype_string == "single": + dtype = np.float32 + else: + dtype = np.float64 + if not hpod_flag: + pod_dtype = dtype + else: + pod_dtype = np.complex128 if dtype == np.float64 else np.complex64 + + # Read the interpolation inputs + query_points_fname = params_file["interpolation"]["query_points_fname"] + interpolated_fields_output_fname = params_file["interpolation"]["output_sequence_fname"] + distributed_axis = params_file["interpolation"]["distributed_axis"] + parallel_io = params_file["interpolation"]["parallel_io"] + interpolation_settings = params_file["interpolation"].get('interpolation_settings', {}) + interpolation_settings['point_interpolator_type'] = interpolation_settings.get('point_interpolator_type', 'multiple_point_legendre_numpy') + + # Start time + start_time = MPI.Wtime() + + #========================================= + # Get the mesh + #========================================= + + # Read the data paths from the input file + mesh_data = IoPathData(params_file["IO"]["mesh_data"]) + field_data = IoPathData(params_file["IO"]["field_data"]) + + # Initialize the mesh file + path = mesh_data.dataPath + casename = mesh_data.casename + index = mesh_data.index + fname = path+casename+'0.f'+str(index).zfill(5) + msh = Mesh(comm) + pynekread(fname, comm, data_dtype=dtype, msh = msh) + + # Read the meshgrid data + print("Reading query points from: ", query_points_fname) + meshgrid = read_data(comm, fname=query_points_fname, keys=["x", "y", "z", "mass"], parallel_io=parallel_io, distributed_axis=distributed_axis, dtype=dtype) + + #========================================= + # Initialize the interpolation + #========================================= + + xyz = [meshgrid["x"].flatten(), meshgrid["y"].flatten() , meshgrid["z"].flatten() ] + xyz = np.ascontiguousarray(np.array(xyz).T) + + # Initialize the probe object + probes = Probes(comm, probes=xyz, msh = msh, output_fname=interpolated_fields_output_fname, **interpolation_settings) + + #========================================= + # Initialize the POD + #========================================= + + # Instance the POD object + pod = POD(comm, number_of_modes_to_update = pod_keep_modes, global_updates = True, auto_expand = False, bckend = backend) + + # Initialize coef to get the mass matrix + bm = meshgrid["mass"] + + # Instance io helper that will serve as buffer for the snapshots + ioh = IoHelp(comm, number_of_fields = number_of_pod_fields, batch_size = pod_batch_size, field_size = bm.size, field_data_type=pod_dtype, mass_matrix_data_type=pod_dtype) + + # Put the mass matrix in the appropiate format (long 1d array) + mass_list = [] + for i in range(0, number_of_pod_fields): + if not hpod_flag: + mass_list.append(np.copy(np.sqrt(bm))) + else: + mass_list.append(np.sqrt(bm)+np.sqrt(bm)*0j) + ioh.copy_fieldlist_to_xi(mass_list) + ioh.bm1sqrt[:,:] = np.copy(ioh.xi[:,:]) + + #========================================= + # Perform the streaming of data + #========================================= + + j = 0 + while j < pod_number_of_snapshots: + + # Recieve the data from fortran + path = field_data.dataPath + casename = field_data.casename + index = field_data.index + fname=path+casename+'0.f'+str(index + j).zfill(5) + fld = FieldRegistry(comm) + pynekread(fname, comm, data_dtype=dtype, fld = fld) + + # Interpolate + field_list = [fld.registry[key] for key in pod_fields] + field_names = pod_fields + probes.interpolate_from_field_list(fld.t, field_list, comm, write_data=False, field_names=field_names) + + # Reshape the data + output_data = {} + for i, key in enumerate(field_names): + output_data[key] = probes.interpolated_fields[:, i+1].reshape(bm.shape).astype(dtype) + print(output_data['u'].shape, output_data['v'].dtype) + + output_data_h = {} + if hpod_flag: + # Hilber transform + for i, key in enumerate(field_names): + output_data_h[key] = hilbert_axis0_mirror_right(output_data[key], axis=0, return_analytic=True).astype(pod_dtype) + output_data = output_data_h + print(output_data['u'].shape, output_data['v'].dtype) + + # Put the snapshot data into a column array + ioh.copy_fieldlist_to_xi([output_data[key] for key in field_names]) + + # Load the column array into the buffer + ioh.load_buffer(scale_snapshot = True) + + # Update POD modes + if ioh.update_from_buffer: + pod.update(comm, buff = ioh.buff[:,:(ioh.buffer_index)]) + + j += 1 + + #========================================= + # Perform post-stream operations + #========================================= + + # Check if there is information in the buffer that should be taken in case the loop exit without flushing + if ioh.buffer_index > ioh.buffer_max_index: + ioh.log.write("info","All snapshots where properly included in the updates") + else: + ioh.log.write("warning","Last loaded snapshot to buffer was: "+repr(ioh.buffer_index-1)) + ioh.log.write("warning","The buffer updates when it is full to position: "+repr(ioh.buffer_max_index)) + ioh.log.write("warning","Data must be updated now to not lose anything, Performing an update with data in buffer ") + pod.update(comm, buff = ioh.buff[:,:(ioh.buffer_index)]) + + # Scale back the modes + pod.scale_modes(comm, bm1sqrt = ioh.bm1sqrt, op = "div") + + # Rotate local modes back to global, This only enters in effect if global_update = false + pod.rotate_local_modes_to_global(comm) + + #========================================= + # Write out the modes + #========================================= + + # Write the data out + output_index = 1 + for j in range(0, pod_write_modes): + + if (j+1) < pod.u_1t.shape[1]: + + ## Split the snapshots into the proper fields + output_data = {} + field_list1d = ioh.split_narray_to_1dfields(pod.u_1t[:,j]) + output_data = {key: field.reshape(bm.shape) for key, field in zip(field_names, field_list1d)} + print(output_data['u'].shape, output_data['u'].dtype) + + if not hpod_flag: + #write the data + prefix = interpolated_fields_output_fname.split(".")[0] + suffix = interpolated_fields_output_fname.split(".")[1] + out_fname = f"{prefix}{str(output_index).zfill(5)}.{suffix}" + write_data(comm, fname=out_fname, data_dict = output_data, parallel_io=parallel_io, distributed_axis=distributed_axis, msh = [meshgrid[key] for key in ["x", "y", "z"]], write_mesh = True, uniform_shape = True) + output_index += 1 + # Clear the fields + else: + #write the data + + # Real part + prefix = interpolated_fields_output_fname.split(".")[0] + suffix = interpolated_fields_output_fname.split(".")[1] + out_fname = f"{prefix}{str(output_index).zfill(5)}.{suffix}" + _output_data = {key: np.real(output_data[key]) for key in field_names} + write_data(comm, fname=out_fname, data_dict = _output_data, parallel_io=parallel_io, distributed_axis=distributed_axis, msh = [meshgrid[key] for key in ["x", "y", "z"]], write_mesh = True, uniform_shape = True) + output_index += 1 + + # Imaginary part + prefix = interpolated_fields_output_fname.split(".")[0] + suffix = interpolated_fields_output_fname.split(".")[1] + out_fname = f"{prefix}{str(output_index).zfill(5)}.{suffix}" + _output_data = {key: np.imag(output_data[key]) for key in field_names} + write_data(comm, fname=out_fname, data_dict = _output_data, parallel_io=parallel_io, distributed_axis=distributed_axis, msh = [meshgrid[key] for key in ["x", "y", "z"]], write_mesh = True, uniform_shape = True) + output_index += 1 + + fld.clear() + + #========================================= + # Write out singular values and right + # singular vectors + #========================================= + + # Write the singular values and vectors + if comm.Get_rank() == 0: + np.save("singular_values", pod.d_1t) + print("Wrote signular values") + np.save("right_singular_vectors", pod.vt_1t) + print("Wrote right signular values") + + # End time + end_time = MPI.Wtime() + # Print the time + if comm.Get_rank() == 0: + print("Time to complete: ", end_time - start_time) diff --git a/scripts/2-rom/5-hpod/inputs.json b/scripts/2-rom/5-hpod/inputs.json new file mode 100644 index 0000000..be1fecc --- /dev/null +++ b/scripts/2-rom/5-hpod/inputs.json @@ -0,0 +1,45 @@ +{ +"Version": 1.0, +"number_of_snapshots": 20, +"batch_size" : 20, +"keep_modes" : 20, +"write_modes" : 20, +"fields" : ["u","v"], +"dtype" : "double", +"backend" : "numpy", +"hpod" : true, +"IO":{ + "mesh_data": { + "dataPath" : "./temp_data/", + "casename" : "field", + "first_index" : 0 + }, + "field_data": { + "dataPath" : "./temp_data/", + "casename" : "field", + "first_index" : 1500 + } +}, +"interpolation":{ + "query_points_fname" : "./meshgrid_data/points.hdf5", + "output_sequence_fname" : "field.hdf5", + "parallel_io" : true, + "distributed_axis" : 1, + "interpolation_settings": { + "find_points_max_iter": 10, + "max_pts":256, + "find_points_iterative": [ + false, + 1 + ], + "find_points_tol": 1e-7, + "write_coords": false, + "global_tree_type": "rank_bbox", + "global_tree_nbins": 100000, + "find_points_comm_pattern": "rma", + "local_data_structure": "kdtree", + "point_interpolator_type": "multiple_point_legendre_numpy", + "use_oriented_bbox": true + } +} +} diff --git a/scripts/2-rom/5-hpod/meshgrid_data/cylinder.py b/scripts/2-rom/5-hpod/meshgrid_data/cylinder.py new file mode 100644 index 0000000..dd514e1 --- /dev/null +++ b/scripts/2-rom/5-hpod/meshgrid_data/cylinder.py @@ -0,0 +1,271 @@ +import numpy as np +import matplotlib.pyplot as plt +import sys +import h5py + +def trap_weights_1d(x): + """Trapezoidal-rule weights for possibly nonuniform 1D grid x (size n).""" + x = np.asarray(x) + n = x.size + w = np.zeros(n, dtype=float) + if n == 1: + w[0] = 0.0 + return w + dx = np.diff(x) + w[0] = 0.5 * dx[0] + w[-1] = 0.5 * dx[-1] + if n > 2: + w[1:-1] = 0.5 * (dx[:-1] + dx[1:]) + return w + +def smoothstep(t): + """C^1 smoothstep from 0 to 1 on t∈[0,1].""" + t = np.clip(t, 0.0, 1.0) + return 3*t*t - 2.0*t*t*t + +def plot_structured_grid(ax, X2, Y2, stride_i=3, stride_j=3, title=""): + nx, ny = X2.shape + for i in range(0, nx, stride_i): + ax.plot(X2[i, :], Y2[i, :], linewidth=0.6) + for j in range(0, ny, stride_j): + ax.plot(X2[:, j], Y2[:, j], linewidth=0.6) + ax.set_aspect("equal", adjustable="box") + ax.set_title(title) + ax.set_xlabel("x") + ax.set_ylabel("y") + +def mass_matrix(Yd, x_1d, y_1d, z_1d, yc=0.0): + """ + Diagonal mass weights w[i,j,k] + + Assumes mapping: + X = x (unchanged), Z = z (unchanged), Y = Yd(x,y,z) + so J = dY/dy (w.r.t computational y coordinate y_1d). + + """ + + # 1D trap weights (works for uniform/nonuniform) + wx = trap_weights_1d(x_1d) + wy = trap_weights_1d(y_1d) + wz = trap_weights_1d(z_1d) + if np.sum(wz) <= 1e-8: # Then I am using only a 2D slice + wz = np.ones_like(wz) + + # Base tensor-product weights + w = wx[:, None, None] * wy[None, :, None] * wz[None, None, :] + + # Split indices + lower = np.where(y_1d < yc)[0] + upper = np.where(y_1d > yc)[0] + equal = np.where(np.isclose(y_1d, yc))[0] + + # Lower half Jacobian: gradient only on that subarray + if lower.size >= 2: + J_lo = np.gradient(Yd[:, lower, :], y_1d[lower], axis=1) + w[:, lower, :] *= J_lo + elif lower.size == 1: + # too few points to differentiate; safest: zero it out + w[:, lower, :] *= 0.0 + + # Upper half Jacobian + if upper.size >= 2: + J_up = np.gradient(Yd[:, upper, :], y_1d[upper], axis=1) + w[:, upper, :] *= J_up + elif upper.size == 1: + w[:, upper, :] *= 0.0 + + # If you have a gridline exactly at yc, set its weight to 0 (it lies on the cut) + if equal.size: + w[:, equal, :] = 0.0 + + return w + +def create_cylinder( + X, Y, + x_bbox, y_bbox, + center=(0,0), + x_smoothing_npoints = 1, + R=1, + y_smooth_window = [1,5], +): + + # ======================================================= + # Obtain the required displacement to create the cylinder + # ======================================================= + + # Perform some checks + xc, yc = center + y_min, y_max = y_bbox + + # computational dy (minimum) + if Y.ndim == 2: + # Y is built from y_1d; use min dy from the first x-line + dy_comp = float(np.min(np.diff(Y[0, :]))) + else: + dy_comp = float(np.min(np.diff(Y[0, :, 0]))) + + if dy_comp <= 0: + raise ValueError("Computational y grid must be strictly increasing.") + + dx = X - xc # Translated coordinates to origin at center + disps = np.sqrt(np.clip(R*R - dx*dx, 0.0, None)) # This is the absolute value of y that every point should be displaced to create the cylinder + + # Mask what is in the top and bottom of the domain + mask_up = (Y >= yc) + mask_lo = ~mask_up + + # =============================================================================== + # Get a smooth step function in the y direction to avoid deforming the boundaries + # =============================================================================== + + y_1d = Y[0,:] + s_1d = np.zeros_like(y_1d) + mask_up_1d = (y_1d >= yc) + mask_lo_1d = ~mask_up_1d + y_up_1d = y_1d[mask_up_1d] + y_lo_1d = y_1d[mask_lo_1d] + + # Set the new coordinate system from 1 - 0 + smooth_window = y_smooth_window + w0, w1 = smooth_window + r = w1 - w0 # Range + s_up = np.ones_like(y_up_1d) + mid = (y_up_1d > w0) & (y_up_1d < w1) # Mask that marks where smoothing happens + s_up[y_up_1d >= w1] = 0.0 + s_up[mid] = 1.0 - (y_up_1d[mid] - w0)/r # map y to [0,1] inside the smoothing window and substract to ramp down + s_up = np.clip(s_up, 0.0, 1.0) # 1 close to cylinder, 0 far away, and linear in-between + + # Mirror the values on the lower half of the domain + s_1d[mask_up_1d] = s_up + s_1d[mask_lo_1d] = np.copy(np.flip(s_up)) + # Smooth step function + smooth_y = smoothstep(s_1d) + + # Apply the smoothing + disps = disps * smooth_y[None, :] + + # ======================================================== + # Smooth in x too, with as imple convolution sort of thing + # ======================================================== + + if x_smoothing_npoints > 1: + kernel_lenght = x_smoothing_npoints + kernel = np.ones(kernel_lenght) + disp_conv = np.zeros_like(disps) + for k in range(X.shape[2]): + for j in range(Y.shape[1]): + disp_conv[:, j, k] = np.convolve(disps[:, j, k], kernel, mode='same') / kernel_lenght + + for k in range(X.shape[2]): + for j in range(Y.shape[1]): + for i in range(X.shape[0]): + disps[i, j, k] = max(disp_conv[i, j, k], disps[i, j, k]) + #disps[i, j, k] = disp_conv[i, j, k] + + # ======================= + # Apply the displacements + # ======================= + # Perform the displacements + Yd = Y.copy() + Yd[mask_up] = Yd[mask_up] + disps[mask_up] + Yd[mask_lo] = Yd[mask_lo] - disps[mask_lo] + + + return X, Yd + +# =========================== + +# Domain +x_bbox = [-10, 30] +y_bbox = [-15, 15] +z_bbox = [2.5 , 2.5] +# Resolution +ddx = 0.1 +ddy = 0.1 +ddz = 1.0 +# Cylinder info +center = (0.0, 0.0) +R = 1.0 +x_kernel_support = 1 # This indicates the support in physical units of the convolution kernel used to smooth in the x direction +x_smoothing_npoints = int(x_kernel_support/ddx) +y_smooth_window = [1,5] # This indicates where the smooth step function is valid in the y direction + + + +if __name__ == "__main__": + + # Number of points (from resolution and domain) + nx = int((x_bbox[1] - x_bbox[0]) / ddx) + ny = int((y_bbox[1] - y_bbox[0]) / ddy) + if np.mod(nx, 2) == 0: + nx += 1 # make it odd + if np.mod(ny, 2) != 0: + ny += 1 # make it even + nz = 1 + ## 1D grids + x_1d = np.linspace(*x_bbox, nx) + y_1d = np.linspace(*y_bbox, ny) + z_1d = np.linspace(*z_bbox, nz) + + # 3D cube + X, Y, Z = np.meshgrid(x_1d, y_1d, z_1d, indexing="ij") # (nx, ny, nz) + + # Cylinder center (all through the z direction) + Xd, Yd = create_cylinder(X, Y, x_bbox, y_bbox, center=center, R=R, x_smoothing_npoints=x_smoothing_npoints, y_smooth_window=y_smooth_window) + Zd = Z + + if 1 == 1: + # --------------------------- + # Plot one z-slice (e.g. mid-plane) + # --------------------------- + k = nz // 2 # slice index + theta = np.linspace(0, 2*np.pi, 400) + cx = center[0] + R*np.cos(theta) + cy = center[1] + R*np.sin(theta) + + fig2, ax2 = plt.subplots() + plot_structured_grid(ax2, Xd[:, :, k], Yd[:, :, k], stride_i=1, stride_j=1, title="Deformed grid (z-slice)") + ax2.plot(cx, cy, linewidth=1.2) + plt.show() + + # --------------------------- + # Sanity checks + # --------------------------- + rmin = np.min(np.sqrt((Xd-center[0])**2 + (Yd-center[1])**2)) + print("min radius after deformation =", rmin) + print("X unchanged:", np.allclose(Xd, X)) + print("Shapes:", X.shape, Y.shape, Z.shape) + + # Generate mass matrix + B = mass_matrix(Yd, x_1d, y_1d, z_1d, yc=center[1]) + B = np.ascontiguousarray(B) + print("Sum of mass matrix:", np.sum(B)) + + cube = (x_bbox[1] - x_bbox[0]) * (y_bbox[1] - y_bbox[0]) * (z_bbox[1] - z_bbox[0]) + + if (z_bbox[1] - z_bbox[0]) > 1e-8: + vol_exact = cube - np.pi * R * R * (z_bbox[1] - z_bbox[0]) + print("Exact volume =", vol_exact) + else: + area_exact = (x_bbox[1] - x_bbox[0]) * (y_bbox[1] - y_bbox[0]) - np.pi * R * R + print("Exact area =", area_exact) + + if nz > 1: + raise ValueError("z_1d is just one value. If you want a 2d slice, then make sure to make nz = 1") + + print(Xd.shape, Yd.shape, Zd.shape, B.shape) + + fname = 'points.hdf5' + with h5py.File(fname, 'w') as f: + + # Create a header + f.attrs['nx'] = nx + f.attrs['ny'] = ny + f.attrs['nz'] = nz + #f.attrs['probe_list_key'] = 'xyz' + + # Include data sets + f.create_dataset('x', data=Xd) + f.create_dataset('y', data=Yd) + f.create_dataset('z', data=Zd) + f.create_dataset('mass', data=B) diff --git a/scripts/2-rom/5-hpod/sem_run.py b/scripts/2-rom/5-hpod/sem_run.py new file mode 100644 index 0000000..b9967f4 --- /dev/null +++ b/scripts/2-rom/5-hpod/sem_run.py @@ -0,0 +1,184 @@ +# Import general modules +import sys +import os +os.environ["OMP_NUM_THREADS"] = "1" +os.environ["OPENBLAS_NUM_THREADS"] = "1" +os.environ["MKL_NUM_THREADS"] = "1" +os.environ["VECLIB_MAXIMUM_THREADS"] = "1" +os.environ["NUMEXPR_NUM_THREADS"] = "1" + +import json +import numpy as np +import matplotlib.pyplot as plt + +# Import MPI +from mpi4py import MPI #equivalent to the use of MPI_init() in C + +# Split communicator for MPI - MPMD +worldcomm = MPI.COMM_WORLD +worldrank = worldcomm.Get_rank() +worldsize = worldcomm.Get_size() +col = 1 +comm = worldcomm.Split(col,worldrank) +rank = comm.Get_rank() +size = comm.Get_size() + +#========================================= +# Define inputs +#========================================= + +# Open input file to see path +f = open ("inputs.json", "r") +params_file = json.loads(f.read()) +f.close() + +# Read the POD inputs +pod_number_of_snapshots = params_file["number_of_snapshots"] +pod_fields = params_file["fields"] +number_of_pod_fields = len(pod_fields) +pod_batch_size = params_file["batch_size"] +pod_keep_modes = params_file["keep_modes"] +pod_write_modes = params_file["write_modes"] +dtype_string = params_file.get("dtype", "double") +backend = params_file.get("backend", "numpy") +if dtype_string == "single": + dtype = np.float32 +else: + dtype = np.float64 + +# Import IO helper functions +from pysemtools.io.utils import get_fld_from_ndarray, IoPathData +# Import modules for reading and writing +from pysemtools.io.ppymech.neksuite import pynekread, pynekwrite +# Import the data types +from pysemtools.datatypes.msh import Mesh +from pysemtools.datatypes.coef import Coef +from pysemtools.datatypes.field import Field, FieldRegistry +from pysemtools.datatypes.utils import create_hexadata_from_msh_fld +# Import types asociated with POD +from pysemtools.rom.pod import POD +from pysemtools.rom.io_help import IoHelp + +# Start time +start_time = MPI.Wtime() + +#========================================= +# Get the mesh +#========================================= + +# Read the data paths from the input file +mesh_data = IoPathData(params_file["IO"]["mesh_data"]) +field_data = IoPathData(params_file["IO"]["field_data"]) + +#========================================= +# Initialize the POD +#========================================= + +# Instance the POD object +pod = POD(comm, number_of_modes_to_update = pod_keep_modes, global_updates = True, auto_expand = False, bckend = backend) + +# Initialize the mesh file +path = mesh_data.dataPath +casename = mesh_data.casename +index = mesh_data.index +fname = path+casename+'0.f'+str(index).zfill(5) +msh = Mesh(comm) +pynekread(fname, comm, data_dtype=dtype, msh = msh) + +# Initialize coef to get the mass matrix +coef = Coef(msh, comm) +bm = coef.B + +# Instance io helper that will serve as buffer for the snapshots +ioh = IoHelp(comm, number_of_fields = number_of_pod_fields, batch_size = pod_batch_size, field_size = bm.size, field_data_type=dtype, mass_matrix_data_type=dtype) + +# Put the mass matrix in the appropiate format (long 1d array) +mass_list = [] +for i in range(0, number_of_pod_fields): + mass_list.append(np.copy(np.sqrt(bm))) +ioh.copy_fieldlist_to_xi(mass_list) +ioh.bm1sqrt[:,:] = np.copy(ioh.xi[:,:]) + +#========================================= +# Perform the streaming of data +#========================================= + +j = 0 +while j < pod_number_of_snapshots: + + # Recieve the data from fortran + path = field_data.dataPath + casename = field_data.casename + index = field_data.index + fname=path+casename+'0.f'+str(index + j).zfill(5) + fld = FieldRegistry(comm) + pynekread(fname, comm, data_dtype=dtype, fld = fld) + + # Put the snapshot data into a column array + ioh.copy_fieldlist_to_xi([fld.registry[key] for key in pod_fields]) + + # Load the column array into the buffer + ioh.load_buffer(scale_snapshot = True) + + # Update POD modes + if ioh.update_from_buffer: + pod.update(comm, buff = ioh.buff[:,:(ioh.buffer_index)]) + + j += 1 + +#========================================= +# Perform post-stream operations +#========================================= + +# Check if there is information in the buffer that should be taken in case the loop exit without flushing +if ioh.buffer_index > ioh.buffer_max_index: + ioh.log.write("info","All snapshots where properly included in the updates") +else: + ioh.log.write("warning","Last loaded snapshot to buffer was: "+repr(ioh.buffer_index-1)) + ioh.log.write("warning","The buffer updates when it is full to position: "+repr(ioh.buffer_max_index)) + ioh.log.write("warning","Data must be updated now to not lose anything, Performing an update with data in buffer ") + pod.update(comm, buff = ioh.buff[:,:(ioh.buffer_index)]) + +# Scale back the modes +pod.scale_modes(comm, bm1sqrt = ioh.bm1sqrt, op = "div") + +# Rotate local modes back to global, This only enters in effect if global_update = false +pod.rotate_local_modes_to_global(comm) + +#========================================= +# Write out the modes +#========================================= + +# Write the data out +for j in range(0, pod_write_modes): + + if (j+1) < pod.u_1t.shape[1]: + + ## Split the snapshots into the proper fields + field_list1d = ioh.split_narray_to_1dfields(pod.u_1t[:,j]) + u_mode = get_fld_from_ndarray(field_list1d[0], msh.lx, msh.ly, msh.lz, msh.nelv) + v_mode = get_fld_from_ndarray(field_list1d[1], msh.lx, msh.ly, msh.lz, msh.nelv) + + # write the data + fld = FieldRegistry(comm) + for i, key in enumerate(pod_fields): + fld.add_field(comm, field_name = key, field = get_fld_from_ndarray(field_list1d[i], msh.lx, msh.ly, msh.lz, msh.nelv), dtype = dtype) + pynekwrite(f"./modes0.f{str(j).zfill(5)}", comm=comm, msh=msh, fld=fld, wdsz=4, istep = j) + +#========================================= +# Write out singular values and right +# singular vectors +#========================================= + +# Write the singular values and vectors +if comm.Get_rank() == 0: + np.save("singular_values", pod.d_1t) + print("Wrote signular values") + np.save("right_singular_vectors", pod.vt_1t) + print("Wrote right signular values") + +# End time +end_time = MPI.Wtime() +# Print the time +if comm.Get_rank() == 0: + print("Time to complete: ", end_time - start_time) diff --git a/scripts/2-rom/5-hpod/visualize.py b/scripts/2-rom/5-hpod/visualize.py new file mode 100644 index 0000000..c9ceb49 --- /dev/null +++ b/scripts/2-rom/5-hpod/visualize.py @@ -0,0 +1,149 @@ +import numpy as np +import h5py +from pyevtk.hl import unstructuredGridToVTK +from pyevtk.vtk import VtkHexahedron, VtkQuad + +def ijk_to_id(i, j, k, ny, nz): + return i*(ny*nz) + j*nz + k + +def ij_to_id(i, j, ny): + return i*ny + j + +def build_cells_with_hole(x, y, z, center=(0.0, 0.0), R=1.0, yc=0.0, tol=0.0): + """ + Build unstructured cell connectivity from a structured point grid: + - if nz > 1: hexahedra (3D) + - if nz == 1: quads (2D slice) + + Removes: + 1) cells that "bridge" across the cut y=yc (some vertices below and some above) + 2) cells whose center is inside cylinder r < R+tol + """ + nx, ny, nz = x.shape + xc, yc0 = center + + if nz == 1: + # ---- 2D quads ---- + conn_list = [] + for i in range(nx - 1): + for j in range(ny - 1): + v0 = ij_to_id(i, j, ny) + v1 = ij_to_id(i+1, j, ny) + v2 = ij_to_id(i+1, j+1, ny) + v3 = ij_to_id(i, j+1, ny) + + xverts = np.array([x[i, j, 0], x[i+1, j, 0], x[i+1, j+1, 0], x[i, j+1, 0]]) + yverts = np.array([y[i, j, 0], y[i+1, j, 0], y[i+1, j+1, 0], y[i, j+1, 0]]) + + # remove bridge cells across cut + if (yverts.min() < yc) and (yverts.max() > yc) and (xverts.min() < xc - R) and (xverts.max() > xc + R): + continue + + # remove inside-cylinder cells via cell center + xcen = 0.25 * (x[i, j, 0] + x[i+1, j, 0] + x[i+1, j+1, 0] + x[i, j+1, 0]) + ycen = 0.25 * (y[i, j, 0] + y[i+1, j, 0] + y[i+1, j+1, 0] + y[i, j+1, 0]) + r = np.sqrt((xcen - xc)**2 + (ycen - yc0)**2) + if r < (R + tol): + continue + + conn_list.append([v0, v1, v2, v3]) + + conn = np.asarray(conn_list, dtype=np.int64).ravel() + ncells = len(conn_list) + offsets = (4 * np.arange(1, ncells + 1, dtype=np.int64)) + celltypes = np.full(ncells, VtkQuad.tid, dtype=np.uint8) + return conn, offsets, celltypes + + else: + # ---- 3D hexes ---- + conn_list = [] + for i in range(nx - 1): + for j in range(ny - 1): + for k in range(nz - 1): + v0 = ijk_to_id(i, j, k, ny, nz) + v1 = ijk_to_id(i+1, j, k, ny, nz) + v2 = ijk_to_id(i+1, j+1, k, ny, nz) + v3 = ijk_to_id(i, j+1, k, ny, nz) + v4 = ijk_to_id(i, j, k+1, ny, nz) + v5 = ijk_to_id(i+1, j, k+1, ny, nz) + v6 = ijk_to_id(i+1, j+1, k+1, ny, nz) + v7 = ijk_to_id(i, j+1, k+1, ny, nz) + + yverts = np.array([ + y[i, j, k], y[i+1, j, k], y[i+1, j+1, k], y[i, j+1, k], + y[i, j, k+1], y[i+1, j, k+1], y[i+1, j+1, k+1], y[i, j+1, k+1], + ]) + + # remove bridge cells across cut + if (yverts.min() < yc) and (yverts.max() > yc): + continue + + # remove inside-cylinder cells via cell center (x,y only) + xcen = 0.125 * ( + x[i, j, k] + x[i+1, j, k] + x[i+1, j+1, k] + x[i, j+1, k] + + x[i, j, k+1] + x[i+1, j, k+1] + x[i+1, j+1, k+1] + x[i, j+1, k+1] + ) + ycen = 0.125 * ( + y[i, j, k] + y[i+1, j, k] + y[i+1, j+1, k] + y[i, j+1, k] + + y[i, j, k+1] + y[i+1, j, k+1] + y[i+1, j+1, k+1] + y[i, j+1, k+1] + ) + r = np.sqrt((xcen - xc)**2 + (ycen - yc0)**2) + if r < (R + tol): + continue + + conn_list.append([v0, v1, v2, v3, v4, v5, v6, v7]) + + conn = np.asarray(conn_list, dtype=np.int64).ravel() + ncells = len(conn_list) + offsets = (8 * np.arange(1, ncells + 1, dtype=np.int64)) + celltypes = np.full(ncells, VtkHexahedron.tid, dtype=np.uint8) + return conn, offsets, celltypes + + +# ----------------------- +# ----------------------- +mesh_fname = "./meshgrid_data/points.hdf5" +with h5py.File(mesh_fname, 'r') as f: + x = f["x"][:] # (nx, ny, nz) + y = f["y"][:] + z = f["z"][:] + +center = (0.0, 0.0) +R = 1.0 +yc = 0.0 + +nx, ny, nz = x.shape + +# Flatten points: for nz==1, flatten (nx,ny) slice; for nz>1 flatten full 3D +if nz == 1: + x1 = np.ascontiguousarray(x[:, :, 0].ravel(order="C")) + y1 = np.ascontiguousarray(y[:, :, 0].ravel(order="C")) + z1 = np.ascontiguousarray(z[:, :, 0].ravel(order="C")) +else: + x1 = np.ascontiguousarray(x.ravel(order="C")) + y1 = np.ascontiguousarray(y.ravel(order="C")) + z1 = np.ascontiguousarray(z.ravel(order="C")) + +# Build connectivity once (reuse per timestep) +conn, offsets, celltypes = build_cells_with_hole( + x, y, z, center=center, R=R, yc=yc, tol=0.0 +) + +for it in range(0, 20): + field_fname = "field" + str(it+1).zfill(5) + ".hdf5" + with h5py.File(field_fname, 'r') as f: + field_dict = {key: np.ascontiguousarray(f[key][:]) for key in f.keys()} + + if nz == 1: + pointData = {k: np.ascontiguousarray(v[:, :, 0].ravel(order="C")) for k, v in field_dict.items()} + else: + pointData = {k: np.ascontiguousarray(v.ravel(order="C")) for k, v in field_dict.items()} + + unstructuredGridToVTK( + "field" + str(it+1).zfill(5), + x1, y1, z1, + connectivity=conn, + offsets=offsets, + cell_types=celltypes, + pointData=pointData + )