From 7f542864eff0d15c77010d22decc73ddb27e47fc Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Fri, 20 Feb 2026 22:17:01 -0500 Subject: [PATCH 01/11] Fix z-boundary backward FD stencil using fields(2) instead of fields(3) The z-direction upper-boundary backward difference at iz_s%end uses fields(2) (y-component) instead of fields(3) (z-component) in the third term, corrupting the divergence in all 3D simulations using this finite difference routine. Co-Authored-By: Claude Opus 4.6 --- src/common/m_finite_differences.fpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/common/m_finite_differences.fpp b/src/common/m_finite_differences.fpp index 09dba0a903..e44b6905c0 100644 --- a/src/common/m_finite_differences.fpp +++ b/src/common/m_finite_differences.fpp @@ -50,7 +50,7 @@ contains if (z == iz_s%beg) then divergence = divergence + (-3._wp*fields(3)%sf(x, y, z) + 4._wp*fields(3)%sf(x, y, z + 1) - fields(3)%sf(x, y, z + 2))/(z_cc(z + 2) - z_cc(z)) else if (z == iz_s%end) then - divergence = divergence + (+3._wp*fields(3)%sf(x, y, z) - 4._wp*fields(3)%sf(x, y, z - 1) + fields(2)%sf(x, y, z - 2))/(z_cc(z) - z_cc(z - 2)) + divergence = divergence + (+3._wp*fields(3)%sf(x, y, z) - 4._wp*fields(3)%sf(x, y, z - 1) + fields(3)%sf(x, y, z - 2))/(z_cc(z) - z_cc(z - 2)) else divergence = divergence + (fields(3)%sf(x, y, z + 1) - fields(3)%sf(x, y, z - 1))/(z_cc(z + 1) - z_cc(z - 1)) end if From 0af98aeb5ec2f0ee68270bcb1d54dabee9f21a8d Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Fri, 20 Feb 2026 22:17:33 -0500 Subject: [PATCH 02/11] Fix bc_x%ve3 never MPI-broadcast (duplicate ve2 instead) The bc_x velocity-end broadcast list has bc_x%ve2 duplicated where bc_x%ve3 should be. bc_y and bc_z rows are correct. Non-root ranks get uninitialized bc_x%ve3 in multi-rank 3D runs with x-boundary velocity BCs. Co-Authored-By: Claude Opus 4.6 --- src/simulation/m_mpi_proxy.fpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/simulation/m_mpi_proxy.fpp b/src/simulation/m_mpi_proxy.fpp index be737782be..95f61fc7d7 100644 --- a/src/simulation/m_mpi_proxy.fpp +++ b/src/simulation/m_mpi_proxy.fpp @@ -145,7 +145,7 @@ contains #:for VAR in [ 'dt','weno_eps','teno_CT','pref','rhoref','R0ref','Web','Ca', 'sigma', & & 'Re_inv', 'poly_sigma', 'palpha_eps', 'ptgalpha_eps', 'pi_fac', & - & 'bc_x%vb1','bc_x%vb2','bc_x%vb3','bc_x%ve1','bc_x%ve2','bc_x%ve2', & + & 'bc_x%vb1','bc_x%vb2','bc_x%vb3','bc_x%ve1','bc_x%ve2','bc_x%ve3', & & 'bc_y%vb1','bc_y%vb2','bc_y%vb3','bc_y%ve1','bc_y%ve2','bc_y%ve3', & & 'bc_z%vb1','bc_z%vb2','bc_z%vb3','bc_z%ve1','bc_z%ve2','bc_z%ve3', & & 'bc_x%pres_in','bc_x%pres_out','bc_y%pres_in','bc_y%pres_out', 'bc_z%pres_in','bc_z%pres_out', & From 85102be01143fe27e0b7efe006421a43a362af9a Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Fri, 20 Feb 2026 22:18:19 -0500 Subject: [PATCH 03/11] Fix fluid_rho broadcast using MPI_LOGICAL instead of mpi_p fluid_rho is a real(wp) array but is broadcast with MPI_LOGICAL type, silently corrupting reference densities via bit reinterpretation on non-root ranks. Co-Authored-By: Claude Opus 4.6 --- src/pre_process/m_mpi_proxy.fpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/pre_process/m_mpi_proxy.fpp b/src/pre_process/m_mpi_proxy.fpp index 30ef061689..cbfac0571b 100644 --- a/src/pre_process/m_mpi_proxy.fpp +++ b/src/pre_process/m_mpi_proxy.fpp @@ -58,7 +58,7 @@ contains & 'igr', 'down_sample', 'simplex_perturb','fft_wrt', 'hyper_cleaning' ] call MPI_BCAST(${VAR}$, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr) #:endfor - call MPI_BCAST(fluid_rho(1), num_fluids_max, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr) + call MPI_BCAST(fluid_rho(1), num_fluids_max, mpi_p, 0, MPI_COMM_WORLD, ierr) #:for VAR in [ 'x_domain%beg', 'x_domain%end', 'y_domain%beg', & & 'y_domain%end', 'z_domain%beg', 'z_domain%end', 'a_x', 'a_y', & From abdabcfde0fa2a4b9386dff3504350a5b47b346c Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Fri, 20 Feb 2026 22:20:18 -0500 Subject: [PATCH 04/11] Fix loc_violations used uninitialized in MPI_Allreduce loc_violations is never set to 0 before the conditional that may or may not assign it. Non-violating ranks sum garbage in the reduction. Co-Authored-By: Claude Opus 4.6 --- src/pre_process/m_data_output.fpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/pre_process/m_data_output.fpp b/src/pre_process/m_data_output.fpp index bac4fdc038..e057234fbd 100644 --- a/src/pre_process/m_data_output.fpp +++ b/src/pre_process/m_data_output.fpp @@ -476,7 +476,7 @@ contains ! Generic loop iterators integer :: i, j, k, l - real(wp) :: loc_violations, glb_violations + real(wp) :: loc_violations = 0._wp, glb_violations ! Downsample variables integer :: m_ds, n_ds, p_ds From 4051479a7afeea4259f74bb0d78591448283e87c Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Fri, 20 Feb 2026 23:48:21 -0500 Subject: [PATCH 05/11] Fix loc_violations gaining implicit SAVE attribute Initializing a local variable in its declaration gives it the SAVE attribute in Fortran, meaning it would not reset to zero on subsequent calls. Move the initialization to an executable assignment so the variable is properly zeroed each time the subroutine is entered. Co-Authored-By: Claude Sonnet 4.6 --- src/pre_process/m_data_output.fpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/pre_process/m_data_output.fpp b/src/pre_process/m_data_output.fpp index e057234fbd..1311a105b7 100644 --- a/src/pre_process/m_data_output.fpp +++ b/src/pre_process/m_data_output.fpp @@ -476,13 +476,15 @@ contains ! Generic loop iterators integer :: i, j, k, l - real(wp) :: loc_violations = 0._wp, glb_violations + real(wp) :: loc_violations, glb_violations ! Downsample variables integer :: m_ds, n_ds, p_ds integer :: m_glb_ds, n_glb_ds, p_glb_ds integer :: m_glb_save, n_glb_save, p_glb_save ! Size of array being saved + loc_violations = 0._wp + if (down_sample) then if ((mod(m + 1, 3) > 0) .or. (mod(n + 1, 3) > 0) .or. (mod(p + 1, 3) > 0)) then loc_violations = 1._wp From c275404796371d720a6e489cc46607da41a2b943 Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Sat, 21 Feb 2026 17:56:10 -0500 Subject: [PATCH 06/11] Fix GRCBC subsonic inflow using wrong L index In the GRCBC subsonic inflow loop (do i = 2, momxb), L(2) was hardcoded instead of L(i), causing only the second wave amplitude to be updated rather than each wave amplitude in the loop. Co-Authored-By: Claude Opus 4.6 --- src/simulation/m_cbc.fpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/simulation/m_cbc.fpp b/src/simulation/m_cbc.fpp index f9d3f3a929..f5f9cb4c30 100644 --- a/src/simulation/m_cbc.fpp +++ b/src/simulation/m_cbc.fpp @@ -930,7 +930,7 @@ contains if (bc_${XYZ}$%grcbc_in) then $:GPU_LOOP(parallelism='[seq]') do i = 2, momxb - L(2) = c**3._wp*Ma*(alpha_rho(i - 1) - alpha_rho_in(i - 1, ${CBC_DIR}$))/Del_in(${CBC_DIR}$) - c*Ma*(pres - pres_in(${CBC_DIR}$))/Del_in(${CBC_DIR}$) + L(i) = c**3._wp*Ma*(alpha_rho(i - 1) - alpha_rho_in(i - 1, ${CBC_DIR}$))/Del_in(${CBC_DIR}$) - c*Ma*(pres - pres_in(${CBC_DIR}$))/Del_in(${CBC_DIR}$) end do if (n > 0) then L(momxb + 1) = c*Ma*(vel(dir_idx(2)) - vel_in(${CBC_DIR}$, dir_idx(2)))/Del_in(${CBC_DIR}$) From 6eea7eae9b98daf98d396e45c4b9d01d547db611 Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Sat, 21 Feb 2026 17:57:11 -0500 Subject: [PATCH 07/11] Fix G_K exponential degradation in damage model The damage factor was applied inside the stress component loop, causing G_K (and G) to be multiplied by the damage factor on every iteration. With N stress components, the effective shear modulus was reduced by damage^N instead of damage^1. Move the damage application before the loop so it is applied exactly once per cell. Co-Authored-By: Claude Opus 4.6 --- src/common/m_variables_conversion.fpp | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/common/m_variables_conversion.fpp b/src/common/m_variables_conversion.fpp index 4f4aca7f05..9c4e72258f 100644 --- a/src/common/m_variables_conversion.fpp +++ b/src/common/m_variables_conversion.fpp @@ -852,11 +852,11 @@ contains end if if (hypoelasticity) then + if (cont_damage) G_K = G_K*max((1._wp - qK_cons_vf(damage_idx)%sf(j, k, l)), 0._wp) $:GPU_LOOP(parallelism='[seq]') do i = strxb, strxe ! subtracting elastic contribution for pressure calculation if (G_K > verysmall) then - if (cont_damage) G_K = G_K*max((1._wp - qK_cons_vf(damage_idx)%sf(j, k, l)), 0._wp) qK_prim_vf(E_idx)%sf(j, k, l) = qK_prim_vf(E_idx)%sf(j, k, l) - & ((qK_prim_vf(i)%sf(j, k, l)**2._wp)/(4._wp*G_K))/gamma_K ! Double for shear stresses @@ -1123,11 +1123,10 @@ contains end if if (hypoelasticity) then + if (cont_damage) G = G*max((1._wp - q_prim_vf(damage_idx)%sf(j, k, l)), 0._wp) do i = strxb, strxe ! adding elastic contribution if (G > verysmall) then - if (cont_damage) G = G*max((1._wp - q_prim_vf(damage_idx)%sf(j, k, l)), 0._wp) - q_cons_vf(E_idx)%sf(j, k, l) = q_cons_vf(E_idx)%sf(j, k, l) + & (q_prim_vf(i)%sf(j, k, l)**2._wp)/(4._wp*G) ! Double for shear stresses From e1790e6e99a84b565b196399cca28c99889c0007 Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Sat, 21 Feb 2026 17:57:29 -0500 Subject: [PATCH 08/11] Fix domain decomposition overwriting muscl_order The IGR conditional block unconditionally reset recon_order to weno_order in its else branch, overwriting the muscl_order that was correctly set by the recon_type check above. Remove the else branch so the original recon_order is preserved when IGR is inactive. Co-Authored-By: Claude Opus 4.6 --- src/common/m_mpi_common.fpp | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/common/m_mpi_common.fpp b/src/common/m_mpi_common.fpp index de3ed8150a..197b5b912c 100644 --- a/src/common/m_mpi_common.fpp +++ b/src/common/m_mpi_common.fpp @@ -1153,8 +1153,6 @@ contains if (igr) then recon_order = igr_order - else - recon_order = weno_order end if ! 3D Cartesian Processor Topology From 358519a9bb6144a0e8b8b52494b38df92fcd9b59 Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Sat, 21 Feb 2026 17:57:51 -0500 Subject: [PATCH 09/11] Fix NaN check reading wrong index in MPI unpack The NaN diagnostic check used q_comm(i)%sf(j, k, l) but the value was unpacked into q_comm(i)%sf(j + unpack_offset, k, l). This meant the check was reading a stale or unrelated cell instead of the just- received value. Co-Authored-By: Claude Opus 4.6 --- src/common/m_mpi_common.fpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/common/m_mpi_common.fpp b/src/common/m_mpi_common.fpp index 197b5b912c..000203386a 100644 --- a/src/common/m_mpi_common.fpp +++ b/src/common/m_mpi_common.fpp @@ -936,7 +936,7 @@ contains (j + buff_size*((k + 1) + (n + 1)*l)) q_comm(i)%sf(j + unpack_offset, k, l) = real(buff_recv(r), kind=stp) #if defined(__INTEL_COMPILER) - if (ieee_is_nan(q_comm(i)%sf(j, k, l))) then + if (ieee_is_nan(q_comm(i)%sf(j + unpack_offset, k, l))) then print *, "Error", j, k, l, i error stop "NaN(s) in recv" end if @@ -991,7 +991,7 @@ contains ((k + buff_size) + buff_size*l)) q_comm(i)%sf(j, k + unpack_offset, l) = real(buff_recv(r), kind=stp) #if defined(__INTEL_COMPILER) - if (ieee_is_nan(q_comm(i)%sf(j, k, l))) then + if (ieee_is_nan(q_comm(i)%sf(j, k + unpack_offset, l))) then print *, "Error", j, k, l, i error stop "NaN(s) in recv" end if @@ -1050,7 +1050,7 @@ contains (l + buff_size))) q_comm(i)%sf(j, k, l + unpack_offset) = real(buff_recv(r), kind=stp) #if defined(__INTEL_COMPILER) - if (ieee_is_nan(q_comm(i)%sf(j, k, l))) then + if (ieee_is_nan(q_comm(i)%sf(j, k, l + unpack_offset))) then print *, "Error", j, k, l, i error stop "NaN(s) in recv" end if From 8efac1dd136098f83af4816f03ab685ab3f7215f Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Mon, 23 Feb 2026 21:33:19 -0500 Subject: [PATCH 10/11] Replace error stop with s_mpi_abort in NaN diagnostic CLAUDE.md Critical Rules forbid error stop. Use s_mpi_abort for proper MPI-aware termination in the Intel NaN detection blocks. Co-Authored-By: Claude Opus 4.6 --- src/common/m_mpi_common.fpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/common/m_mpi_common.fpp b/src/common/m_mpi_common.fpp index 000203386a..8768298306 100644 --- a/src/common/m_mpi_common.fpp +++ b/src/common/m_mpi_common.fpp @@ -938,7 +938,7 @@ contains #if defined(__INTEL_COMPILER) if (ieee_is_nan(q_comm(i)%sf(j + unpack_offset, k, l))) then print *, "Error", j, k, l, i - error stop "NaN(s) in recv" + call s_mpi_abort("NaN(s) in recv") end if #endif end do @@ -993,7 +993,7 @@ contains #if defined(__INTEL_COMPILER) if (ieee_is_nan(q_comm(i)%sf(j, k + unpack_offset, l))) then print *, "Error", j, k, l, i - error stop "NaN(s) in recv" + call s_mpi_abort("NaN(s) in recv") end if #endif end do @@ -1052,7 +1052,7 @@ contains #if defined(__INTEL_COMPILER) if (ieee_is_nan(q_comm(i)%sf(j, k, l + unpack_offset))) then print *, "Error", j, k, l, i - error stop "NaN(s) in recv" + call s_mpi_abort("NaN(s) in recv") end if #endif end do From 944c6adbed9834e10c0a89721a9e07cee30118c9 Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Tue, 24 Feb 2026 16:04:14 -0500 Subject: [PATCH 11/11] Replace error stop with s_mpi_abort in NaN timestep output check MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit CLAUDE.md forbids error stop; use s_mpi_abort for proper MPI-aware termination. The print * before the abort is retained since it logs the exact coordinates (j,k,l), variable index, rank, and timestep — information that doesn't fit in the abort message string. Co-Authored-By: Claude Sonnet 4.6 --- src/simulation/m_start_up.fpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/simulation/m_start_up.fpp b/src/simulation/m_start_up.fpp index d8d0c87417..925ad55db0 100644 --- a/src/simulation/m_start_up.fpp +++ b/src/simulation/m_start_up.fpp @@ -965,7 +965,7 @@ contains do j = 0, m if (ieee_is_nan(real(q_cons_ts(stor)%vf(i)%sf(j, k, l), kind=wp))) then print *, "NaN(s) in timestep output.", j, k, l, i, proc_rank, t_step, m, n, p - error stop "NaN(s) in timestep output." + call s_mpi_abort("NaN(s) in timestep output.") end if end do end do