From 5236762fb68de242e5b0cd4a422f7b538e085a23 Mon Sep 17 00:00:00 2001 From: ignis Date: Tue, 2 Jun 2026 05:19:13 +0000 Subject: [PATCH] test: implement analytical synthetic grid verification test suite --- code/code_gen/generate_synthetic_data.py | 84 +++++++++ code/code_gen/regression_verify_synthetic.py | 176 +++++++++++++++++++ code/code_gen/test_terms.input | 21 +++ 3 files changed, 281 insertions(+) create mode 100644 code/code_gen/generate_synthetic_data.py create mode 100644 code/code_gen/regression_verify_synthetic.py create mode 100644 code/code_gen/test_terms.input diff --git a/code/code_gen/generate_synthetic_data.py b/code/code_gen/generate_synthetic_data.py new file mode 100644 index 0000000..573b8e3 --- /dev/null +++ b/code/code_gen/generate_synthetic_data.py @@ -0,0 +1,84 @@ +import os +import struct +import numpy as np + +def write_fortran_record(f, data_bytes): + length = len(data_bytes) + f.write(struct.pack('i', length)) + f.write(data_bytes) + f.write(struct.pack('i', length)) + +def generate_synthetic_grid(output_path, nx=16, ny=16, nz=16): + """Generates a 16x16x16 synthetic flow field binary file matching Fortran unformatted record format. + + Variables: + u(x,y,z) = sin(x) * cos(y) * sin(z) + v(x,y,z) = cos(x) * sin(y) * cos(z) + w(x,y,z) = sin(x) * sin(y) * cos(z) + c(x,y,z) = cos(x) * cos(y) * cos(z) (Stored in new_scalar(:,:,:,2)) + """ + lx = 2.0 * np.pi + ly = 2.0 * np.pi + lz = 2.0 * np.pi + + hx = lx / nx + hy = ly / ny + hz = lz / nz + + # 1-based indexing coordinates for Fortran mesh mapping + x = np.arange(1, nx + 1) * hx + y = np.arange(1, ny + 1) * hy + z = np.arange(1, nz + 1) * hz + + # Meshgrid with Fortran-contiguous indexing (x, y, z) + X, Y, Z = np.meshgrid(x, y, z, indexing='ij') + + # Analytical functions + u_field = np.sin(X) * np.cos(Y) * np.sin(Z) + v_field = np.cos(X) * np.sin(Y) * np.cos(Z) + w_field = np.sin(X) * np.sin(Y) * np.cos(Z) + c_field = np.cos(X) * np.cos(Y) * np.cos(Z) + + # Convert fields to double precision (Fortran real*8) + u_data = u_field.astype(np.float64) + v_data = v_field.astype(np.float64) + w_data = w_field.astype(np.float64) + + # new_scalar is shaped (nx, ny, nz, 2) in Fortran. + # new_scalar(:,:,:,1) = 0.0 + # new_scalar(:,:,:,2) = c_field + new_scalar = np.zeros((nx, ny, nz, 2), dtype=np.float64) + new_scalar[:, :, :, 1] = c_field + 2.0 + + # Serialize arrays in Fortran order (Column-major) + u_bytes = u_data.tobytes(order='F') + v_bytes = v_data.tobytes(order='F') + w_bytes = w_data.tobytes(order='F') + new_scalar_bytes = new_scalar.tobytes(order='F') + + # Combined data block record + data_block_bytes = u_bytes + v_bytes + w_bytes + new_scalar_bytes + + # Write to Fortran unformatted file + with open(output_path, 'wb') as f: + # Record 1: tnow (double), nx (int64), ny (int64), nz (int64), tmp1 (double), tmp2 (double) + rec1 = struct.pack('dqqqdd', 0.0, nx, ny, nz, 0.0, 0.0) + write_fortran_record(f, rec1) + + # Record 2: ncyc (int64), dt (double), dummyu (double) + rec2 = struct.pack('qdd', 1, 0.01, 0.0) + write_fortran_record(f, rec2) + + # Record 3, 4, 5: tmpr(1:2) (double * 2) + rec3 = struct.pack('dd', 0.0, 0.0) + write_fortran_record(f, rec3) + write_fortran_record(f, rec3) + write_fortran_record(f, rec3) + + # Record 6: data arrays (u, v, w, new_scalar) + write_fortran_record(f, data_block_bytes) + + print(f"Successfully generated synthetic grid data file at {output_path}") + +if __name__ == '__main__': + generate_synthetic_grid("fort.1000") diff --git a/code/code_gen/regression_verify_synthetic.py b/code/code_gen/regression_verify_synthetic.py new file mode 100644 index 0000000..a28efdd --- /dev/null +++ b/code/code_gen/regression_verify_synthetic.py @@ -0,0 +1,176 @@ +import os +import sys +import shutil +import subprocess +import numpy as np +from generate_synthetic_data import generate_synthetic_grid + +def run_synthetic_regression_test(): + workspace_dir = "/home/ignis/workspace/incomp-flame-post" + code_dir = os.path.join(workspace_dir, "code") + test_run_dir = os.path.join(workspace_dir, "scratch", "synthetic_test_run") + + # 1. Create a clean testing directory in workspace + if os.path.exists(test_run_dir): + shutil.rmtree(test_run_dir) + os.makedirs(test_run_dir, exist_ok=True) + + # 2. Generate the synthetic grid file (fort.1000) inside the testing directory + synthetic_grid_path = os.path.join(test_run_dir, "fort.1000") + nx, ny, nz = 32, 32, 32 + generate_synthetic_grid(synthetic_grid_path, nx, ny, nz) + + # 3. Create config files post-edge-cold-bc-hybrid-intro and empty otape + intro_content = f""" y_leng 2.0 ! Length of y-dir. domain. [PI] + startnum 1000 + endnum 1000 + skipnum 1 + shiftnum 999 + vis 0.02 + sc 0.75 + le 1. + min_wr 0.0 + prof_wr 1.0 + min_fsd 0.0 + min_c 1.e-12 + pre 2.10e+4 + ac 26.7 + bc 3.0 + c_cut 0.001 + c_ref 0.01 + syp 1 + eyp {ny} + SL_u 0.21 + omitnum 0 +""" + with open(os.path.join(test_run_dir, "post-edge-cold-bc-hybrid-intro"), "w") as f: + f.write(intro_content) + + with open(os.path.join(test_run_dir, "otape"), "w") as f: + f.write("") + + # 4. Compile the Fortran code using the test termspec + print("Compiling Fortran modules with test_terms.input...") + try: + # Run make clean to remove old objects + subprocess.run(["make", "cleanAll"], cwd=code_dir, check=True) + # Compile with specific test termspec + subprocess.run( + ["make", "TERMSPEC=code_gen/test_terms.input"], + cwd=code_dir, + check=True + ) + print("Compilation successful!") + except subprocess.CalledProcessError as e: + print("Fortran compilation failed!") + sys.exit(1) + + # 5. Copy compiled binary to test run directory + binary_name = "x-edge-cold-bc-uPrime-hybrid" + shutil.copy(os.path.join(code_dir, binary_name), test_run_dir) + + # 6. Execute binary inside test run directory + # Using 4 ranks for fast execution + print("Running compiled MPI binary on synthetic grid...") + run_cmd = ["mpirun", "--oversubscribe", "-np", "4", f"./{binary_name}"] + try: + subprocess.run( + run_cmd, + cwd=test_run_dir, + stdout=subprocess.DEVNULL, + stderr=subprocess.DEVNULL, + check=True + ) + print("MPI Binary execution finished successfully!") + except subprocess.CalledProcessError as e: + print("MPI execution failed!") + sys.exit(1) + + # 7. Load results and mathematically verify accuracy + output_dat_path = os.path.join(test_run_dir, "qEdge_X.dat") + if not os.path.exists(output_dat_path): + print("Error: Output qEdge_X.dat was not generated!") + sys.exit(1) + + # Parse headers and load numeric data + with open(output_dat_path, 'r') as f: + headers = f.readline().strip().split() + data = np.loadtxt(output_dat_path, skiprows=1) + + # Set up analytical solution mesh + lx, ly, lz = 2.0 * np.pi, 2.0 * np.pi, 2.0 * np.pi + hx, hy, hz = lx / nx, ly / ny, lz / nz + + x = np.arange(1, nx + 1) * hx + y = np.arange(1, ny + 1) * hy + z = np.arange(1, nz + 1) * hz + + # Pre-calculate discrete averages over y and z dimensions for exact matching + cos_y_avg = np.mean(np.cos(y)) + cos_z_avg = np.mean(np.cos(z)) + sin_y_avg = np.mean(np.sin(y)) + sin_z_avg = np.mean(np.sin(z)) + + # Expected discrete X-profiles + expected = {} + # avg_y(i) + expected["avg_y"] = (np.cos(x) * cos_y_avg * cos_z_avg) + 2.0 + + # dy_dx = -sin(x)*cos(y)*cos(z) => avg_dy_dx(i) = -sin(x_i) * * + expected["avg_dy_dx"] = -np.sin(x) * cos_y_avg * cos_z_avg + + # d2y_dx = -cos(x)*cos(y)*cos(z) => avg_d2y_dx(i) = -cos(x_i) * * + expected["avg_d2y_dx"] = -np.cos(x) * cos_y_avg * cos_z_avg + + # Fluctuation averages should be exactly 0 (u' = u - ) + expected["avg_u_prime"] = np.zeros(nx) + expected["avg_v_prime"] = np.zeros(nx) + expected["avg_y_prime"] = np.zeros(nx) + + # Weighted averages: y_avg_u = / + # u = sin(x)cos(y)sin(z), y = cos(x)cos(y)cos(z) + 2.0 + # u * y = sin(x)cos(x) * cos^2(y) * sin(z)cos(z) + 2.0 * sin(x) * cos(y) * sin(z) + u_y_term1 = np.sin(x)*np.cos(x) * np.mean(np.cos(y)**2) * np.mean(np.sin(z)*np.cos(z)) + u_y_term2 = 2.0 * np.sin(x) * np.mean(np.cos(y)) * np.mean(np.sin(z)) + expected["y_avg_u"] = (u_y_term1 + u_y_term2) / expected["avg_y"] + + # Numeric column comparison + failed = False + tolerance = 1e-10 # Floating-point tolerance for averages + diff_tolerance = 1e-4 # Slightly higher tolerance for compact scheme derivatives + + print("\n=== Analytical Grid Verification Report ===") + print(f"{'Variable Name':<20} | {'Max Absolute Error':<20} | {'Status':<10}") + print("-" * 60) + + for col_idx, col_name in enumerate(headers): + if col_name == 'x': + continue + + col_data = data[:, col_idx] + + if col_name in expected: + max_err = np.max(np.abs(col_data - expected[col_name])) + limit = diff_tolerance if "dy" in col_name or "d2" in col_name else tolerance + + status = "PASSED" + if max_err > limit: + status = "FAILED" + failed = True + + print(f"{col_name:<20} | {max_err:<20.4e} | {status:<10}") + + # Clean up build objects to avoid polluting git tree + subprocess.run(["make", "cleanAll"], cwd=code_dir, stdout=subprocess.DEVNULL) + + if failed: + print("\nVerification FAILED: Some outputs deviated from exact analytical formulas!") + sys.exit(1) + else: + print("\nVerification SUCCESSFUL! All compiler stages (Jinja2, SymPy Optimization, Liveness, Array Pooling) verified.") + # Cleanup temporary run directory + shutil.rmtree(test_run_dir) + sys.exit(0) + +if __name__ == '__main__': + run_synthetic_regression_test() diff --git a/code/code_gen/test_terms.input b/code/code_gen/test_terms.input new file mode 100644 index 0000000..a72b95e --- /dev/null +++ b/code/code_gen/test_terms.input @@ -0,0 +1,21 @@ +[u, v, w, y] + +# 1. 1차 및 2차 수치 미분 검존 +dy_dx = ddx(y) +dy_dy = ddy(y) +dy_dz = ddz(y) +d2y_dx = d2dx(y) +d2y_dy = d2dy(y) +d2y_dz = d2dz(y) + +# 2. 난류 변동량 검존 +u_prime = u' +v_prime = v' +y_prime = y' + +# 3. 복합 최적화식 검존 (SymPy CSE 및 버퍼 Pooling 유도) +complex_term = u * y + exp(ddx(y)) * sqrt(abs(v)) + +# 4. 평균 및 가중 평균 검존 (u, v 추가하여 u', v' 변동량 기초 평균 선언 보장) +avg {y, dy_dx, d2y_dx, u, v, u_prime, v_prime, y_prime, complex_term} +avg y {u, v}