add core changes to probe gen

Author: sbryngelson

src/post_process/m_data_output.fpp

@@ -725,8 +725,8 @@ contains
 
             if (p > 0) then
                 err = DBMKOPTLIST(2, optlist)
-                err = DBADDIOPT(optlist, DBOPT_LO_OFFSET, lo_offset)
-                err = DBADDIOPT(optlist, DBOPT_HI_OFFSET, hi_offset)
+                err = DBADDIAOPT(optlist, DBOPT_LO_OFFSET, 3, lo_offset)
+                err = DBADDIAOPT(optlist, DBOPT_HI_OFFSET, 3, hi_offset)
                 if (grid_geometry == 3) then
                     err = DBPUTQM(dbfile, 'rectilinear_grid', 16, &
                                   'x', 1, 'y', 1, 'z', 1, &
@@ -743,8 +743,8 @@ contains
                 err = DBFREEOPTLIST(optlist)
             elseif (n > 0) then
                 err = DBMKOPTLIST(2, optlist)
-                err = DBADDIOPT(optlist, DBOPT_LO_OFFSET, lo_offset)
-                err = DBADDIOPT(optlist, DBOPT_HI_OFFSET, hi_offset)
+                err = DBADDIAOPT(optlist, DBOPT_LO_OFFSET, 2, lo_offset)
+                err = DBADDIAOPT(optlist, DBOPT_HI_OFFSET, 2, hi_offset)
                 err = DBPUTQM(dbfile, 'rectilinear_grid', 16, &
                               'x', 1, 'y', 1, 'z', 1, &
                               x_cb, y_cb, DB_F77NULL, dims, 2, &

src/simulation/m_data_output.fpp

@@ -363,6 +363,7 @@ contains
             write (3, *) ! new line
 
             if (.not. f_approx_equal(icfl_max_glb, icfl_max_glb)) then
+                print *, 'icfl', icfl_max_glb
                 call s_mpi_abort('ICFL is NaN. Exiting.')
             elseif (icfl_max_glb > 1._wp) then
                 print *, 'icfl', icfl_max_glb
@@ -1129,6 +1130,8 @@ contains
         real(wp) :: G_local
         real(wp) :: dyn_p, T
         real(wp) :: damage_state
+        character(LEN=15) :: FMT_glb
+        character(len=30) :: FMT
 
         integer :: i, j, k, l, s, d !< Generic loop iterator
 
@@ -1145,6 +1148,13 @@ contains
 
         T = dflt_T_guess
 
+        ! Set format string based on precision
+        if (precision == 1) then
+            FMT_glb = 'F28.7'
+        else
+            FMT_glb = 'F28.16'
+        end if
+
         ! Non-dimensional time calculation
         if (time_stepper == 23) then
             nondim_time = mytime
@@ -1578,28 +1588,35 @@ contains
                             R(1), &
                             Rdot(1)
                     else if (elasticity) then
-                        write (i + 30, '(6X,F12.6,F24.8,F24.8,F24.8,F24.8,'// &
-                               'F24.8,F24.8,F24.8)') &
+                        FMT = '(6X,F12.6,8'//FMT_glb//')'
+                        write (i + 30, FMT) &
                             nondim_time, &
                             rho, &
                             vel(1), &
                             vel(2), &
+                            vel(3), &
                             pres, &
                             tau_e(1), &
                             tau_e(2), &
                             tau_e(3)
                     else
-                        write (i + 30, '(6X,F12.6,F24.8,F24.8,F24.8)') &
+                        FMT = '(6X,F12.6,11'//FMT_glb//')'
+                        write (i + 30, FMT) &
                             nondim_time, &
                             rho, &
                             vel(1), &
-                            pres
-                        print *, 'time =', nondim_time, 'rho =', rho, 'pres =', pres
+                            vel(2), &
+                            vel(3), &
+                            pres, & ! Out of tolerance
+                            gamma, &
+                            pi_inf, &
+                            qv, &
+                            c, &
+                            accel
                     end if
                 else
-                    write (i + 30, '(6X,F12.6,F24.8,F24.8,F24.8,F24.8,'// &
-                           'F24.8,F24.8,F24.8,F24.8,F24.8,'// &
-                           'F24.8)') &
+                    FMT = '(6X,F12.6,11'//FMT_glb//')'
+                    write (i + 30, FMT) &
                         nondim_time, &
                         rho, &
                         vel(1), &
@@ -1608,9 +1625,7 @@ contains
                         pres, &
                         gamma, &
                         pi_inf, &
-                        qv, &
-                        c, &
-                        accel
+                        qv
                 end if
             end if
         end do

src/simulation/m_derived_variables.fpp

@@ -168,9 +168,7 @@ contains
             $:GPU_UPDATE(host='[accel_mag]')
 
             call s_derive_center_of_mass(q_prim_ts2(2)%vf, c_mass)
-
             call s_write_probe_files(t_step, q_cons_ts(1)%vf, accel_mag)
-
             call s_write_com_files(t_step, c_mass)
         end if
 

toolchain/mfc/packer/pack.py

@@ -136,6 +136,10 @@ def _extract_doubles(s: str) -> list:
                 lines = content.splitlines()
                 content = "\n".join(lines[1:])  # Skip the first line
                 doubles = _extract_doubles(content)
+            elif 'probe' in short_filepath:
+                # Probes are constructed as <time> <value1> <value2> ...
+                # Just discard the first column
+                doubles = [float(e) for line in content.strip().splitlines() if line.strip() for e in _extract_doubles(line)[1:]]
             else:
                 # Every line is <x> <y> <z> <value> (<y> and <z> are optional). So the
                 # number of dimensions is the number of doubles in the first line minus 1.

toolchain/mfc/test/case.py

@@ -230,8 +230,36 @@ def create_directory(self):
         mods['omega_wrt(2)'] = 'T'
         mods['omega_wrt(3)'] = 'T'
 else:
-    mods['parallel_io']   = 'F'
+    mods['parallel_io'] = 'F'
     mods['prim_vars_wrt'] = 'F'
+    if ('fd_order' not in case) and ('fd_order' not in mods):
+        mods['fd_order'] = 1
+    if ('probe_wrt' not in case):
+        mods['probe_wrt'] = 'T'
+        if case['n'] == 0:
+            mods['num_probes'] = 2
+            dx = (case['x_domain%end'] - case['x_domain%beg']) / 3
+            for i in range(2):
+                mods[f'probe({{i+1}})%x'] = case['x_domain%beg'] + (i+1) * dx
+        elif case['p'] == 0:
+            mods['num_probes'] = 4
+            dx = (case['x_domain%end'] - case['x_domain%beg']) / 3
+            dy = (case['y_domain%end'] - case['y_domain%beg']) / 3
+            for i in range(2):
+                for j in range(2):
+                    mods[f'probe({{i*2+j+1}})%x'] = case['x_domain%beg'] + (i+1) * dx
+                    mods[f'probe({{i*2+j+1}})%y'] = case['y_domain%beg'] + (j+1) * dy
+        else:
+            mods['num_probes'] = 8
+            dx = (case['x_domain%end'] - case['x_domain%beg']) / 3
+            dy = (case['y_domain%end'] - case['y_domain%beg']) / 3
+            dz = (case['z_domain%end'] - case['z_domain%beg']) / 3
+            for i in range(2):
+                for j in range(2):
+                    for k in range(2):
+                        mods[f'probe({{i*4+j*2+k+1}})%x'] = case['x_domain%beg'] + (i+1) * dx
+                        mods[f'probe({{i*4+j*2+k+1}})%y'] = case['y_domain%beg'] + (j+1) * dy
+                        mods[f'probe({{i*4+j*2+k+1}})%z'] = case['z_domain%beg'] + (k+1) * dz
 
 print(json.dumps({{**case, **mods}}))
 """)