Append to probe and CoM files only when continuing a run

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.

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

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,17 +146,24 @@ 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'
             file_path = trim(case_dir) // trim(file_path)
 
             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
 

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