CI probe: in-situ check of the alt_soundspeed bulk moduli

Author: sbryngelson

src/simulation/m_rhs.fpp

@@ -121,6 +121,7 @@ module m_rhs
     real(wp), allocatable, dimension(:,:,:,:) :: qL_rsx_vf, qR_rsx_vf
     real(wp), allocatable, dimension(:,:,:,:) :: dqL_rsx_vf, dqR_rsx_vf
     $:GPU_DECLARE(create='[blkmod1, blkmod2, alpha1, alpha2, Kterm]')
+    logical :: probe_first = .true.
     $:GPU_DECLARE(create='[qL_rsx_vf, qR_rsx_vf]')
     $:GPU_DECLARE(create='[dqL_rsx_vf, dqR_rsx_vf]')
 
@@ -1057,6 +1058,8 @@ contains
         real(wp) :: inv_ds, flux_face1, flux_face2
         real(wp) :: advected_qty_val, pressure_val, velocity_val
         real(wp) :: G1_eff, G2_eff
+        real(wp) :: probe_ref
+        integer :: probe_nbad
 
         G1_eff = 0._wp
         G2_eff = 0._wp
@@ -1096,6 +1099,35 @@ contains
                 end do
             end do
             $:END_GPU_PARALLEL_LOOP()
+            ! CI probe: the same bulk moduli recomputed on the host from the same fields, on the first call only.
+            if (probe_first) then
+                probe_first = .false.
+                $:GPU_UPDATE(host='[blkmod1, blkmod2]')
+                $:GPU_UPDATE(host='[q_prim_vf%vf(eqn_idx%E)%sf, q_prim_vf%vf(eqn_idx%adv%beg)%sf, q_prim_vf%vf(eqn_idx%adv%end)%sf]')
+                $:GPU_UPDATE(host='[q_prim_vf%vf(eqn_idx%cont%beg)%sf, q_prim_vf%vf(eqn_idx%cont%end)%sf]')
+                probe_nbad = 0
+                do q_loop = 0, p
+                    do l_loop = 0, n
+                        do k_loop = 0, m
+                            call s_phase_bulk_modulus(q_prim_vf%vf(eqn_idx%E)%sf(k_loop, l_loop, q_loop), &
+                                                      & q_prim_vf%vf(eqn_idx%adv%beg)%sf(k_loop, l_loop, q_loop), &
+                                                      & q_prim_vf%vf(eqn_idx%cont%beg)%sf(k_loop, l_loop, q_loop), 1, probe_ref)
+                            probe_ref = probe_ref + (4._wp/3._wp)*G1_eff
+                            if (.not. abs(blkmod1(k_loop, l_loop, q_loop) - probe_ref) <= 1.e-12_wp*abs(probe_ref)) then
+                                probe_nbad = probe_nbad + 1
+                                if (probe_nbad == 1 .and. proc_rank == 0) print '(A, 3I6, 2ES24.15)', &
+                                    & ' eos probe rhs blkmod1 first bad cell ', k_loop, l_loop, q_loop, blkmod1(k_loop, l_loop, &
+                                    & q_loop), probe_ref
+                            end if
+                        end do
+                    end do
+                end do
+                if (proc_rank == 0) print '(A, I8, A, I8)', ' eos probe rhs blkmod1 mismatches ', probe_nbad, ' of ', &
+                    & (m + 1)*(n + 1)*(p + 1)
+                if (proc_rank == 0) print '(A, 4ES24.15)', ' eos probe rhs fields at cell 0: E adv cont blkmod1 ', &
+                    & q_prim_vf%vf(eqn_idx%E)%sf(0, 0, 0), q_prim_vf%vf(eqn_idx%adv%beg)%sf(0, 0, 0), &
+                    & q_prim_vf%vf(eqn_idx%cont%beg)%sf(0, 0, 0), blkmod1(0, 0, 0)
+            end if
         end if
 
         select case (idir)