feat(amr): IGR solver support with restriction-only coarse/fine coupling - the fine block runs its own fixed-iteration sigma solve seeded and Dirichlet-bounded by the converged coarse sigma (piecewise-constant parent injection over the full buffered range; per-iteration BC/halo populate skipped under amr_in_fine_advance; coarse jac/jac_old warm-start state bounced across the fine advance), fine RK update uses the IGR coefficient form (dt embedded in the rhs), flux-register capture/apply return early under igr (no face fluxes exposed by the fused IGR kernels - seam conservation is truncation-order, documented), amr_subcycle gated. Validated: free-stream through block EXACTLY preserved, AMR-vs-reference at resolution scale (rho 1.3e-4 rel-L2, static and dynamic regrid), transverse seam artifact 3e-4 absolute from the coarse/fine sigma jump (truncation-level, documented), conservation drift at the reference's own outflow scale. Goldens 6C20B752 (static) 660FFBFE (dynamic regrid)

Author: sbryngelson

docs/documentation/amr.md

@@ -217,7 +217,7 @@ a diagnostic message for unsupported combinations.
 | QBMM bubbles (non-polytropic) | Supported (static block, no subcycle) | Each block carries its own `pb`/`mv` quadrature side-state: prolonged piecewise-constant (realizability), advanced with the block's own rhs scratch, restricted back with the moments; dynamic regrid and `amr_subcycle` remain gated |
 | Lagrangian bubbles (`bubbles_lagrange = T`) | Supported (cloud excluded from blocks) | Two-way coupling lives on the coarse grid: regrid suppresses tags and clips candidate boxes around the cloud's padded bbox (positions + `mapCells` smearing + stencil + drift margin, recomputed collectively each regrid), the fine advance skips the EL hooks, EL volume fractions prolong WITHOUT the sum-to-one closure (their sum is the local liquid fraction), and a per-stage guard aborts if the cloud reaches an active block |
 | Immersed boundaries (`ib = T`; one or more non-STL bodies, static or prescribed-motion `moving_ibm=1`) | Supported | Per-block fine-grid IB markers/ghost points, rebuilt each fine substage at the body's sub-time position for a moving body; non-conservative ghost-cell forcing at the body; with dynamic regrid candidate boxes expand to fully contain each body at its live position plus margin, the fine IB state rebuilds after every regrid, and a per-substage guard aborts if a moving body reaches its block boundary between regrids; force-driven (`moving_ibm=2`)/STL gated; a body spanning a rank seam is rejected at startup |
-| IGR solver | **Not supported** | Unvalidated with the fine-level advance |
+| IGR solver (`igr = T`) | Supported (restriction-only coupling) | The fine block runs its own fixed-iteration sigma solve, seeded and Dirichlet-bounded by the converged coarse sigma (frozen ghost ring; the per-iteration BC populate is skipped); the coarse warm-start state is saved/restored across the fine advance. The Berger-Colella reflux is NOT captured from the fused IGR flux kernels, so seam conservation is truncation-order rather than exact; free-stream preservation is exact. `amr_subcycle` is gated |
 | 2D axisymmetric (`cyl_coord = T`, `p = 0`) | Supported | Geometric sources read the live (swapped) grid; the axis half-width cell's per-cell WENO coefficients are recomputed for each block on swap/restore; blocks stay `buff_size` off the axis (domain-edge clamp); the axis-singularity viscous treatment runs on the coarse pass only |
 | 3D cylindrical (`cyl_coord = T`, `p > 0`) | **Not supported** | The per-stage azimuthal Fourier filter is a global operation incompatible with the block-local fine advance |
 | Grid stretching (`stretch_x[y,z] = T`) | Supported | Fine ghost-shell coordinates extend by exact parent-cell bisection, and the spacing-dependent WENO coefficients are recomputed for the active grid on every block swap/restore (`amr_weno_coef_recompute`, armed automatically when the grid is nonuniform); prolongation stays conservative but its slope estimate is first-order on nonuniform parents |

docs/documentation/case.md

@@ -855,7 +855,11 @@ lives on the coarse grid, regrid suppresses tags and clips boxes around the clou
 bounding box, EL volume fractions prolong without the sum-to-one closure (their sum is the
 local liquid fraction, not 1), and a per-stage guard aborts if the cloud reaches an active
 block (reduce `amr_regrid_int` or increase `amr_buf`).
-AMR is incompatible with surface tension, IGR, 3D cylindrical
+The IGR solver is supported with restriction-only coarse/fine coupling: the fine block runs
+its own fixed-iteration sigma solve seeded and Dirichlet-bounded by the converged coarse
+sigma; seam conservation is truncation-order (no reflux capture from the fused IGR flux
+kernels), free-stream preservation is exact, and `amr_subcycle` is gated under IGR.
+AMR is incompatible with surface tension, 3D cylindrical
 coordinates (2D axisymmetric IS supported), MHD, hyperelasticity, Riemann-extrapolation
 boundaries (bc = -4), and `active_box`. `hybrid_weno`/`hybrid_riemann` are supported: each
 level recomputes the smoothness sensor over its own (swapped) bounds every RHS call.

src/simulation/m_amr.fpp

@@ -30,6 +30,7 @@ module m_amr
     use m_weno, only: s_compute_weno_coefficients
     use m_acoustic_src, only: acoustic_supp_lo, acoustic_supp_hi
     use m_bubbles_EL, only: s_lag_cloud_bbox_local
+    use m_igr, only: jac, jac_old
 
     implicit none
 
@@ -93,6 +94,12 @@ module m_amr
     type(int_bounds_info) :: sw_idwint(3), sw_idwbuff(3)
     logical               :: sw_acoustic_source
 
+    !> IGR sigma-state bounce (igr only): the fine solve reuses the module jac/jac_old arrays at fine indices (the extent guard
+    !! keeps fine bounds inside), so the coarse contents - jac_old is the Jacobi warm start persisting across steps - are saved here
+    !! across the fine advance.
+    real(wp), allocatable :: sw_jac(:,:,:), sw_jac_old(:,:,:)
+    $:GPU_DECLARE(create='[sw_jac, sw_jac_old]')
+
     !> Lagrangian bubble-cloud exclusion support: padded global coarse-index bbox of the cloud (positions + mapCells smearing +
     !! stencil headroom [+ drift margin at regrid]). Blocks and regrid boxes stay clear of it: a bubble inside a block loses two-way
     !! coupling (the fine advance skips the EL hooks and the coarse result under the block is discarded by restriction). Recomputed
@@ -222,6 +229,10 @@ contains
             allocate (sw_z_cc(-buff_size:p + buff_size))
             allocate (sw_dz(-buff_size:p + buff_size))
         end if
+        if (igr) then
+            @:ALLOCATE(sw_jac(idwbuff(1)%beg:idwbuff(1)%end, idwbuff(2)%beg:idwbuff(2)%end, idwbuff(3)%beg:idwbuff(3)%end))
+            @:ALLOCATE(sw_jac_old(idwbuff(1)%beg:idwbuff(1)%end, idwbuff(2)%beg:idwbuff(2)%end, idwbuff(3)%beg:idwbuff(3)%end))
+        end if
 
         ! Grid uniformity policy. The two spacing-uniformity consumers are both handled exactly:
         ! the fine-block ghost-shell coordinates extend by exact parent-cell bisection (reads
@@ -244,7 +255,7 @@ contains
                 amr_weno_coef_recompute = .true.
             end if
         end if
-        if (weno_order == 1) amr_weno_coef_recompute = .false.  ! order 1 has no grid-dependent coefficients
+        if (weno_order == 1 .or. igr) amr_weno_coef_recompute = .false.  ! order 1 / IGR: no grid-dependent WENO coefficients
 
         ! preallocate every slot at max buffered extents (each sized exactly as the single legacy block): coords (host) +
         ! fields (device-resident @:ALLOCATE so s_compute_rhs kernels can run on the fine block; q_cons/q_prim/rhs get the
@@ -1247,6 +1258,13 @@ contains
         ! must be rebuilt for the block's own uniform grid (no-op flag on uniform grids)
         if (amr_weno_coef_recompute) call s_amr_recompute_weno_coefs()
 
+        ! IGR: save the coarse sigma state and seed the fine solve. jac holds THIS stage's
+        ! converged coarse sigma (the coarse RHS ran first), so its parent values are both the
+        ! best initial guess and the frozen Dirichlet ghost data for the block-local Jacobi
+        ! solve (the per-iteration BC/halo populate is skipped under amr_in_fine_advance).
+        ! Piecewise-constant parent injection over the full buffered fine range.
+        if (igr) call s_amr_igr_swap_sigma()
+
     end subroutine s_amr_swap_to_fine
 
     !> Restore the global grid state saved by s_amr_swap_to_fine.
@@ -1265,9 +1283,65 @@ contains
         call s_amr_sync_grid_state_to_device()
         if (hypoelasticity) call s_hypoelastic_update_fd_coeffs()
         if (amr_weno_coef_recompute) call s_amr_recompute_weno_coefs()
+        if (igr) call s_amr_igr_restore_sigma()
 
     end subroutine s_amr_restore_coarse
 
+    !> Save the coarse jac/jac_old and seed the (already swapped-in) fine block's sigma state by piecewise-constant parent injection
+    !! from the saved coarse sigma: interior = initial guess, ghost shell = frozen Dirichlet coupling data for the block-local
+    !! iterative solve.
+    impure subroutine s_amr_igr_swap_sigma()
+
+        integer :: j, k, l, ci, cj, ck
+
+        $:GPU_PARALLEL_LOOP(collapse=3, private='[j, k, l]')
+        do l = sw_idwbuff(3)%beg, sw_idwbuff(3)%end
+            do k = sw_idwbuff(2)%beg, sw_idwbuff(2)%end
+                do j = sw_idwbuff(1)%beg, sw_idwbuff(1)%end
+                    sw_jac(j, k, l) = jac(j, k, l)
+                    sw_jac_old(j, k, l) = jac_old(j, k, l)
+                end do
+            end do
+        end do
+        $:END_GPU_PARALLEL_LOOP()
+        $:GPU_PARALLEL_LOOP(collapse=3, private='[j, k, l, ci, cj, ck]')
+        do l = idwbuff(3)%beg, idwbuff(3)%end
+            do k = idwbuff(2)%beg, idwbuff(2)%end
+                do j = idwbuff(1)%beg, idwbuff(1)%end
+                    ci = amr_isect_lo(1) + floor(real(j, wp)/2._wp) - start_idx(1)
+                    cj = 0; ck = 0
+                    if (n_glb > 0) cj = amr_isect_lo(2) + floor(real(k, wp)/2._wp) - start_idx(2)
+                    if (p_glb > 0) ck = amr_isect_lo(3) + floor(real(l, wp)/2._wp) - start_idx(3)
+                    ci = min(max(ci, sw_idwbuff(1)%beg), sw_idwbuff(1)%end)
+                    cj = min(max(cj, sw_idwbuff(2)%beg), sw_idwbuff(2)%end)
+                    ck = min(max(ck, sw_idwbuff(3)%beg), sw_idwbuff(3)%end)
+                    jac(j, k, l) = sw_jac(ci, cj, ck)
+                    jac_old(j, k, l) = sw_jac(ci, cj, ck)
+                end do
+            end do
+        end do
+        $:END_GPU_PARALLEL_LOOP()
+
+    end subroutine s_amr_igr_swap_sigma
+
+    !> Restore the coarse jac/jac_old saved by s_amr_igr_swap_sigma (bounds already restored).
+    impure subroutine s_amr_igr_restore_sigma()
+
+        integer :: j, k, l
+
+        $:GPU_PARALLEL_LOOP(collapse=3, private='[j, k, l]')
+        do l = idwbuff(3)%beg, idwbuff(3)%end
+            do k = idwbuff(2)%beg, idwbuff(2)%end
+                do j = idwbuff(1)%beg, idwbuff(1)%end
+                    jac(j, k, l) = sw_jac(j, k, l)
+                    jac_old(j, k, l) = sw_jac_old(j, k, l)
+                end do
+            end do
+        end do
+        $:END_GPU_PARALLEL_LOOP()
+
+    end subroutine s_amr_igr_restore_sigma
+
     !> Recompute the WENO reconstruction coefficient arrays from the CURRENT grid globals (the fine block's after a swap, the coarse
     !! grid's after a restore). s_compute_weno_coefficients reads the live cell-boundary arrays, refreshes uniform_grid, and pushes
     !! its own device updates; the coefficient arrays were sized to the coarse local ranges at init, which the fine ranges never
@@ -1642,9 +1716,10 @@ contains
         call s_amr_restore_coarse()
         amr_in_fine_advance = .false.
 
-        ! RK stage update (device kernel; mirror of the coarse non-IGR form)
+        ! RK stage update (device kernel; mirror of the coarse form - under IGR the rhs already
+        ! embeds dt, matching the coarse igr update, so the dt factor is 1)
         call s_amr_fine_rk_update(amr_slots(amr_cur)%q_cons, amr_slots(amr_cur)%q_cons_stor, amr_slots(amr_cur)%rhs, coefs(1), &
-                                  & coefs(2), coefs(3), coefs(4), dt)
+                                  & coefs(2), coefs(3), coefs(4), merge(1._wp, dt, igr))
         if (qbmm .and. .not. polytropic) call s_amr_fine_rk_update_pbmv(coefs(1), coefs(2), coefs(3), coefs(4), dt)
         ! 6-equation model: per-stage pressure relaxation on the block (before IB correct, coarse order)
         if (model_eqns == model_eqns_6eq .and. (.not. relax)) call s_amr_pressure_relax_fine()
@@ -2990,6 +3065,10 @@ contains
             if (allocated(sw_x_cb)) deallocate (sw_x_cb, sw_x_cc, sw_dx)
             if (allocated(sw_y_cb)) deallocate (sw_y_cb, sw_y_cc, sw_dy)
             if (allocated(sw_z_cb)) deallocate (sw_z_cb, sw_z_cc, sw_dz)
+            if (igr) then
+                @:DEALLOCATE(sw_jac)
+                @:DEALLOCATE(sw_jac_old)
+            end if
 
         end subroutine s_finalize_amr_module
 

src/simulation/m_amr_registers.fpp

@@ -145,6 +145,7 @@ contains
         logical                        :: accum
 
         if (.not. amr) return
+        if (igr) return  ! stage-1 IGR coupling is restriction-only: the fused IGR flux kernels do not expose face fluxes to capture
         if (amr_in_fine_advance .and. .not. amr_rank_owns_block) return
         islot = amr_cur  ! working block slot (local => captured by value in the device kernels below)
         ! flux data was just written by device kernels; the face reads below run as device kernels too
@@ -471,6 +472,7 @@ contains
         real(wp)                                               :: fblo, fbhi, mlo, mhi
 
         if (.not. amr) return
+        if (igr) return  ! stage-1 IGR: restriction-only coupling (no captured fluxes)
         islot = amr_cur  ! working block slot (local => captured by value in the device kernels below)
         ! per-face participation: each face's correction runs on the rank owning its OUTSIDE cell layer
         ! (all faces at np=1); freg slices from a rank-boundary face's fine side arrive via
@@ -593,6 +595,7 @@ contains
         integer :: d, eq, t1, t2, t1_hi, t2_hi, islot
 
         if (.not. amr) return
+        if (igr) return  ! stage-1 IGR: restriction-only coupling (no captured fluxes)
         if (.not. amr_rank_owns_block) return
         islot = amr_cur  ! working block slot (local => captured by value in the device kernels below)
         do d = 1, 3
@@ -626,6 +629,7 @@ contains
         real(wp)                                               :: fblo, fbhi, mlo, mhi, dtl
 
         if (.not. amr) return
+        if (igr) return  ! stage-1 IGR: restriction-only coupling (no captured fluxes)
         islot = amr_cur  ! working block slot (local => captured by value in the device kernels below)
         ! per-face participation and index conventions: see s_amr_apply_reflux
         call s_amr_reflux_face_flags(sidx, ext, own_lo, own_hi)

src/simulation/m_checker.fpp

@@ -93,7 +93,7 @@ contains
                    & "hybrid_riemann (cheap central/Rusanov flux) is incompatible with the low_Mach correction")
 
         if (amr) then
-            @:PROHIBIT(recon_type /= recon_type_weno, "amr requires WENO reconstruction")
+            @:PROHIBIT((.not. igr) .and. recon_type /= recon_type_weno, "amr requires WENO reconstruction (or the IGR solver)")
             @:PROHIBIT(time_stepper /= time_stepper_rk3, "amr requires time_stepper = 3 (SSP-RK3)")
             ! 6-equation support: the internal-energy equations prolong/restrict on the generic
             ! conservative path and the per-stage pressure relaxation (cell-local) also runs on
@@ -113,7 +113,11 @@ contains
             ! hypoelasticity is supported: stress components prolong via the generic conservative-linear
             ! path and the swap/restore recomputes the spacing-dependent FD coefficients per grid
             @:PROHIBIT(hyperelasticity .or. mhd, "amr does not support hyperelasticity/MHD")
-            @:PROHIBIT(igr, "amr does not support the IGR solver")
+            ! IGR is supported with restriction-only coarse/fine coupling (stage 1): the fine
+            ! block runs its own fixed-iteration sigma solve seeded and Dirichlet-bounded by the
+            ! converged coarse sigma; the Berger-Colella reflux is not yet captured from the
+            ! fused IGR flux kernels, so seam conservation is truncation-order, not exact
+            @:PROHIBIT(igr .and. amr_subcycle, "amr_subcycle with the IGR solver is not yet supported (lockstep only)")
             ! Lagrangian bubbles are supported with the cloud EXCLUDED from fine blocks (two-way
             ! coupling lives on the coarse grid): regrid suppresses tags and clips boxes around
             ! the cloud's padded bbox, and a per-stage guard aborts if the cloud reaches a block

src/simulation/m_igr.fpp

@@ -18,7 +18,7 @@ module m_igr
     implicit none
 
     private; public :: s_initialize_igr_module, s_igr_iterative_solve, s_igr_riemann_solver, s_igr_sigma_x, s_igr_flux_add, &
-        & s_finalize_igr_module
+        & s_finalize_igr_module, jac, jac_old
 
     !> @cond
 #ifdef __NVCOMPILER_GPU_UNIFIED_MEM
@@ -300,7 +300,10 @@ contains
             end do
             $:END_GPU_PARALLEL_LOOP()
 
-            call s_populate_F_igr_buffers(bc_type, jac_sf)
+            ! AMR fine advance: the block's jac ghost shell holds frozen Dirichlet data prolonged
+            ! from the converged coarse sigma (seeded by s_amr_igr_swap_sigma); the physical-BC/halo
+            ! populate would write wrong (coarse-indexed) data on the swapped block grid
+            if (.not. amr_in_fine_advance) call s_populate_F_igr_buffers(bc_type, jac_sf)
 
             if (igr_iter_solver == 1) then  ! Jacobi iteration
                 $:GPU_PARALLEL_LOOP(private='[j, k, l]', collapse=3)