diff --git a/examples/2D_probe_rerun/README.md b/examples/2D_probe_rerun/README.md new file mode 100644 index 000000000..4d2c1ea41 --- /dev/null +++ b/examples/2D_probe_rerun/README.md @@ -0,0 +1,57 @@ +# Probe files across a re-run + +`s_open_probe_files` appends whenever `D/probe*_prim.dat` already exists: + +```fortran +if (file_exist) then + open (..., STATUS='old', POSITION='append') +``` + +That is correct when a run is being **continued**. It is wrong when one is being **started over**, which is +what happens every time a case is re-run in place after a parameter change. The second run's rows land on top +of the first's, nothing in the file marks the join, and the time column simply resets partway down. A reader +sees one monotonic series and is silently wrong. + +There is no warning, no header and no separator. The two runs need not even share a grid — a case whose +resolution changed between runs produces a file whose first half was recorded at different probe locations. + +## Reproducing + +``` +./mfc.sh run examples/2D_probe_rerun/case.py -n 1 # 20 steps +./mfc.sh run examples/2D_probe_rerun/case.py -n 1 # same again, from scratch +wc -l D/probe1_prim.dat +``` + +| | rows in `D/probe1_prim.dat` | +| --- | --- | +| before | **40** — two runs of 20, spliced | +| after | **20** | + +And the time column, read straight through, before the fix: + +``` +0.028977 +0.030682 +0.032386 <- end of run 1 +0.000000 <- run 2 starts, time goes backwards +0.001705 +0.003409 +``` + +## The fix + +Append only when continuing: `t_step_start > 0`, or `n_start > 0` under `cfl_dt`. A fresh start replaces the +file, which is what every other output MFC writes already does. `s_open_com_files` had the same pattern and +gets the same treatment. + +## Why it matters beyond tidiness + +This produced three separate wrong numbers in one project before it was noticed. The worst was a jet whose +probe files held a t = 15 run of 4,962 rows followed by a t = 40 run of 13,233 — **on different grids**. +Read together they manufactured a velocity drop of 0.99 U_j in a single sample, which was time running +backwards at the seam and was diagnosed as a physical instability first. + +A related symptom is louder and easier to spot: if the probe output format changes between runs, the column +count changes partway down the file and `numpy.loadtxt` refuses it outright ("the number of columns changed +from 11 to 18"). That one at least announces itself. The time reset does not. diff --git a/examples/2D_probe_rerun/case.py b/examples/2D_probe_rerun/case.py new file mode 100644 index 000000000..a3dc1a564 --- /dev/null +++ b/examples/2D_probe_rerun/case.py @@ -0,0 +1,83 @@ +#!/usr/bin/env python3 +"""A case with probes, meant to be run twice in the same directory. + +MFC appends to D/probe*_prim.dat when the file already exists. That is right when a run is being continued +and wrong when one is being started over: the second run's rows land on top of the first's, nothing in the +file marks the join, and the time column simply resets partway down. A reader sees one monotonic series. + +Run it, run it again, and count the rows. See README.md. +""" + +import json +import os + +Re, Ma, gamma, U, rho = 40.0, 0.1, 1.4, 1.0, 1.0 +P = rho * U**2 / (gamma * Ma**2) +L = 2.0 +N = int(os.environ.get("N", 64)) +NSTEP = int(os.environ.get("NSTEP", 20)) +dx = 2 * L / N + +case = { + "run_time_info": "T", + "parallel_io": "T", + "prim_vars_wrt": "T", + "format": "silo", + "precision": "double", + "x_domain%beg": -L, + "x_domain%end": L, + "y_domain%beg": -L, + "y_domain%end": L, + "m": N - 1, + "n": N - 1, + "p": 0, + "cyl_coord": "F", + "dt": 0.3 * dx / (U + U / Ma), + "t_step_start": 0, + "t_step_stop": NSTEP, + "t_step_save": NSTEP, + "num_patches": 1, + "num_fluids": 1, + "model_eqns": "5eq", + "alt_soundspeed": "F", + "mpp_lim": "F", + "mixture_err": "T", + "time_stepper": "rk3", + "weno_order": 5, + "weno_eps": 1.0e-10, + "weno_Re_flux": "T", + "weno_avg": "T", + "avg_state": "arithmetic", + "mapped_weno": "T", + "null_weights": "F", + "mp_weno": "F", + "riemann_solver": "hllc", + "low_Mach": 2, + "wave_speeds": "direct", + "viscous": "T", + "fd_order": 4, + "patch_icpp(1)%geometry": 3, + "patch_icpp(1)%x_centroid": 0.0, + "patch_icpp(1)%y_centroid": 0.0, + "patch_icpp(1)%length_x": 2 * L, + "patch_icpp(1)%length_y": 2 * L, + "patch_icpp(1)%vel(1)": U, + "patch_icpp(1)%vel(2)": 0.0, + "patch_icpp(1)%pres": P, + "patch_icpp(1)%alpha_rho(1)": rho, + "patch_icpp(1)%alpha(1)": 1.0, + "fluid_pp(1)%gamma": 1.0 / (gamma - 1.0), + "fluid_pp(1)%eos": "ideal_gas", + "fluid_pp(1)%Re(1)": Re, + "bc_x%beg": -3, + "bc_x%end": -3, + "bc_y%beg": -3, + "bc_y%end": -3, + "probe_wrt": "T", + "num_probes": 2, + "probe(1)%x": 0.0, + "probe(1)%y": 0.0, + "probe(2)%x": 0.5, + "probe(2)%y": 0.5, +} +print(json.dumps(case, indent=4)) diff --git a/src/simulation/m_data_output.fpp b/src/simulation/m_data_output.fpp index 57f7eaf9f..43741e27c 100644 --- a/src/simulation/m_data_output.fpp +++ b/src/simulation/m_data_output.fpp @@ -113,11 +113,18 @@ contains character(len=path_len + 3*name_len) :: file_path !< Relative path to the CoM file in the case directory integer :: i !< Generic loop iterator + logical :: fresh_start + + fresh_start = (t_step_start == 0) .or. (cfl_dt .and. n_start == 0) do i = 1, num_fluids write (file_path, '(A,I0,A)') '/fluid', i, '_com.dat' file_path = trim(case_dir) // trim(file_path) - open (i + 120, file=trim(file_path), form='formatted', position='append', status='unknown') + if (fresh_start) then + open (i + 120, file=trim(file_path), form='formatted', status='replace') + else + open (i + 120, file=trim(file_path), form='formatted', position='append', status='unknown') + end if if (n == 0) then write (i + 120, '(A)') ' Non-Dimensional Time ' // ' Total Mass ' // ' x-loc ' // ' Total Volume ' else if (p == 0) then @@ -139,6 +146,9 @@ contains character(LEN=path_len + 3*name_len) :: file_path !< Relative path to the probe data file in the case directory integer :: i !< Generic loop iterator logical :: file_exist + logical :: fresh_start + + fresh_start = (t_step_start == 0) .or. (cfl_dt .and. n_start == 0) do i = 1, num_probes write (file_path, '(A,I0,A)') '/D/probe', i, '_prim.dat' @@ -146,10 +156,14 @@ contains inquire (file=trim(file_path), exist=file_exist) - if (file_exist) then + ! Append only when continuing a run. A fresh start that appends splices the previous run's rows + ! onto this one's, and nothing in the file marks the join: the time column simply resets partway + ! down, and the two runs need not even share a grid. Readers see one monotonic series and are + ! silently wrong. + if (file_exist .and. .not. fresh_start) then open (i + 30, FILE=trim(file_path), form='formatted', STATUS='old', POSITION='append') else - open (i + 30, FILE=trim(file_path), form='formatted', STATUS='unknown') + open (i + 30, FILE=trim(file_path), form='formatted', STATUS='replace') end if end do diff --git a/toolchain/mfc/test/cases.py b/toolchain/mfc/test/cases.py index fff3b4552..d48cae472 100644 --- a/toolchain/mfc/test/cases.py +++ b/toolchain/mfc/test/cases.py @@ -3167,6 +3167,8 @@ def foreach_example(): # the transverse momentum drifts past the 1e-3 Example tolerance across compilers # (nvhpc passes; Intel and CCE disagree by ~2e-3 absolute). No single golden is portable. "2D_hybrid_slab", + # The bug it shows needs two consecutive runs in one directory; a single Example run cannot see it. + "2D_probe_rerun", ] if path in casesToSkip: continue