test: implement analytical synthetic grid verification test suite

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")

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) * <cos(y)> * <cos(z)>
+    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) * <cos(y)> * <cos(z)>
+    expected["avg_d2y_dx"] = -np.cos(x) * cos_y_avg * cos_z_avg
+    
+    # Fluctuation averages should be exactly 0 (u' = 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 * y> / <y>
+    # 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()

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}