merge: viscous AMR (SP11) + fine-ghost-coordinate fix into the PR branch

Author: sbryngelson

docs/documentation/case.md

@@ -676,7 +676,7 @@ To restart the simulation from $k$-th time step, see @ref running "Restarting Ca
 | `sfc_partition_wrt`     | Logical | Report SFC-weighted load-balance partition |
 | `rank_time_wrt`         | Logical | Report per-rank RHS compute-time imbalance (max/mean) |
 | `load_balance`          | Logical | (Experimental/diagnostic) Weighted static Cartesian decomposition at init (requires `parallel_io = T`, >1 rank). Measured gain is small on CPU (~5%) and can be slower on GPU due to the occupancy floor; equal decomposition is near-optimal for uniform-cost workloads. |
-| `amr`                   | Logical | (Experimental) Enable block-structured AMR: a 2:1 refined level-1 patch with gradient-based dynamic regrid, optional dt/2 subcycling, and conservative coupling with refluxing. Requires WENO reconstruction, SSP-RK3, model_eqns=2; num_fluids > 1 requires mpp_lim. |
+| `amr`                   | Logical | (Experimental) Enable block-structured AMR: a 2:1 refined level-1 patch with gradient-based dynamic regrid, optional dt/2 subcycling, and conservative coupling with refluxing. Requires WENO reconstruction, SSP-RK3, model_eqns=2; num_fluids > 1 requires mpp_lim; supports physical viscosity. |
 | `amr_patch_beg(i)`      | Integer | Refined-patch start cell index in direction $i$ (level-0 index space) |
 | `amr_patch_end(i)`      | Integer | Refined-patch end cell index in direction $i$ (level-0 index space) |
 | `amr_regrid_int`        | Integer | Steps between AMR regrid events (0 = static patch) |
@@ -787,7 +787,17 @@ Multiple fluids (`num_fluids > 1`) are supported and additionally require `mpp_l
 whose volume-fraction clamp+renormalize maintains coarse/fine alpha consistency; the
 per-fluid masses are refluxed exactly, and volume fractions are prolonged with a
 sum-preserving closure (fine-level volume fractions sum to one by construction).
-It is incompatible with viscosity, surface tension, bubble models, phase-change (relax),
+Physical viscosity (`viscous = T`) is supported: the viscous stress/work travels through
+the momentum- and energy-equation source fluxes, which are captured into the same
+coarse–fine flux registers as the advective fluxes, so the interface is refluxed against
+the matched *total* (advective + viscous) flux and energy — including viscous work — is
+conserved. Fine-ghost velocity gradients at the coarse–fine boundary are taken from the
+conservative-linear prolongation of the coarse state (no special gradient reconstruction);
+that interface inconsistency is bounded and conservation is enforced by the flux-register
+matching. The density-gradient regrid tagger does not sense shear or boundary layers well,
+so viscous features may need a static or generously buffered patch (error-estimator taggers
+are future work).
+It is incompatible with surface tension, bubble models, phase-change (relax),
 immersed boundaries, IGR, cylindrical coordinates, MHD, chemistry, `hybrid_weno`,
 `hybrid_riemann`, and `acoustic_source`.
 Multi-rank runs are supported: the fine level mirrors the base decomposition (each rank

src/simulation/m_amr.fpp

@@ -637,6 +637,42 @@ contains
             z_cc(0:amr_fine%p) = amr_fine%z_cc(0:amr_fine%p)
             dz(0:amr_fine%p) = amr_fine%dz(0:amr_fine%p)
         end if
+        ! Extend the fine grid into the ghost shell (s_build_level_coords only fills the interior 0:m).
+        ! The viscous velocity gradient divides by ghost cell-center spacings (x_cc(j+-1) - x_cc(j)), so at
+        ! a rank seam the fine subdomain-edge cell would otherwise use a stale coarse spacing and diverge
+        ! from the np=1 (fully-owned) result. Uniform 2:1 continuation of the (locally-uniform) base region.
+        block
+            integer  :: jg
+            real(wp) :: sp
+            sp = amr_fine%dx(amr_fine%m)
+            do jg = amr_fine%m + 1, amr_fine%m + buff_size
+                dx(jg) = sp; x_cb(jg) = x_cb(jg - 1) + sp; x_cc(jg) = x_cc(jg - 1) + sp
+            end do
+            sp = amr_fine%dx(0)
+            do jg = -1, -buff_size, -1
+                dx(jg) = sp; x_cb(jg - 1) = x_cb(jg) - sp; x_cc(jg) = x_cc(jg + 1) - sp
+            end do
+            if (n_glb > 0) then
+                sp = amr_fine%dy(amr_fine%n)
+                do jg = amr_fine%n + 1, amr_fine%n + buff_size
+                    dy(jg) = sp; y_cb(jg) = y_cb(jg - 1) + sp; y_cc(jg) = y_cc(jg - 1) + sp
+                end do
+                sp = amr_fine%dy(0)
+                do jg = -1, -buff_size, -1
+                    dy(jg) = sp; y_cb(jg - 1) = y_cb(jg) - sp; y_cc(jg) = y_cc(jg + 1) - sp
+                end do
+            end if
+            if (p_glb > 0) then
+                sp = amr_fine%dz(amr_fine%p)
+                do jg = amr_fine%p + 1, amr_fine%p + buff_size
+                    dz(jg) = sp; z_cb(jg) = z_cb(jg - 1) + sp; z_cc(jg) = z_cc(jg - 1) + sp
+                end do
+                sp = amr_fine%dz(0)
+                do jg = -1, -buff_size, -1
+                    dz(jg) = sp; z_cb(jg - 1) = z_cb(jg) - sp; z_cc(jg) = z_cc(jg + 1) - sp
+                end do
+            end if
+        end block
         ! sync the swapped extents/bounds/coordinates to the device: RHS kernels read the
         ! device copies of these GPU_DECLARE'd globals (stale coarse bounds = OOB kernels)
         call s_amr_sync_grid_state_to_device()

src/simulation/m_amr_registers.fpp

@@ -18,6 +18,13 @@
 !! remainder (the +alpha*d(u_star)/dx compression term m_rhs assembles from flux_src_n = u_star) is deliberately NOT captured:
 !! alpha is genuinely non-conservative, so forcing flux-matching on u_star would be wrong; coarse/fine volume-fraction
 !! consistency is instead maintained by mpp_lim's clamp+renormalize (required by the checker for amr with num_fluids > 1).
+!!
+!! Viscous (SP11): the viscous stress/work face fluxes travel through flux_src_n for the momentum and energy equations
+!! (m_rhs s_compute_additional_physics_rhs: rhs += (flux_src_n(j-1) - flux_src_n(j))/dx, identical face indexing and sign to the
+!! advective flux_n). They are captured into the SAME registers (added on top of the advective flux for mom..E) so the c/f reflux
+!! matches the TOTAL advective+viscous flux; energy conservation therefore includes the viscous work. Fine-ghost velocity gradients
+!! at the c/f boundary come from the conservative-linear cons prolongation (no special gradient reconstruction) - like the alpha
+!! K-term, that inconsistency is bounded, and conservation is enforced by the flux-register matching.
 module m_amr_registers
 
     use m_derived_types
@@ -112,10 +119,11 @@ contains
     !! call (amr_in_fine_advance false, coarse globals) fills creg at the patch boundary faces; fine call (flag true, globals
     !! swapped to the fine patch) fills freg at fine faces -1 and m/n/p. creg uses relative 0-based transverse; freg uses 0-based
     !! fine.
-    impure subroutine s_amr_capture_boundary_flux(id, flux_dir, stage)
+    impure subroutine s_amr_capture_boundary_flux(id, flux_dir, flux_src, stage)
 
         integer, intent(in)            :: id
         type(vector_field), intent(in) :: flux_dir
+        type(vector_field), intent(in) :: flux_src
         integer, intent(in)            :: stage
         integer                        :: eq, t1, t2, jlo, jhi, t1_hi, t2_hi, o1, o2
         integer                        :: sidx(3), ext(3)
@@ -176,6 +184,30 @@ contains
                 end do
             end do
             $:END_GPU_PARALLEL_LOOP()
+            ! total-flux matching: add the viscous momentum/energy face fluxes (flux_src) into the same fine
+            ! registers so the c/f reflux sees advective+viscous. Base coef/accum are applied above; always
+            ! accumulate here. Inviscid path skips this entirely (registers stay byte-identical).
+            if (viscous) then
+                $:GPU_PARALLEL_LOOP(collapse=3)
+                do t2 = 0, t2_hi
+                    do t1 = 0, t1_hi
+                        do eq = eqn_idx%mom%beg, eqn_idx%E
+                            select case (id)
+                            case (1)
+                                freg(1)%lo(eq, t1, t2) = freg(1)%lo(eq, t1, t2) + coef*real(flux_src%vf(eq)%sf(jlo, t1, t2), wp)
+                                freg(1)%hi(eq, t1, t2) = freg(1)%hi(eq, t1, t2) + coef*real(flux_src%vf(eq)%sf(jhi, t1, t2), wp)
+                            case (2)
+                                freg(2)%lo(eq, t1, t2) = freg(2)%lo(eq, t1, t2) + coef*real(flux_src%vf(eq)%sf(t1, jlo, t2), wp)
+                                freg(2)%hi(eq, t1, t2) = freg(2)%hi(eq, t1, t2) + coef*real(flux_src%vf(eq)%sf(t1, jhi, t2), wp)
+                            case (3)
+                                freg(3)%lo(eq, t1, t2) = freg(3)%lo(eq, t1, t2) + coef*real(flux_src%vf(eq)%sf(t1, t2, jlo), wp)
+                                freg(3)%hi(eq, t1, t2) = freg(3)%hi(eq, t1, t2) + coef*real(flux_src%vf(eq)%sf(t1, t2, jhi), wp)
+                            end select
+                        end do
+                    end do
+                end do
+                $:END_GPU_PARALLEL_LOOP()
+            end if
         else
             ! coarse branch: a face's capture runs on the rank owning the coarse cells just OUTSIDE it (its
             ! flux_n covers that face; at a rank-interior face the same rank also holds the inside cells).
@@ -258,6 +290,35 @@ contains
                 end do
             end do
             $:END_GPU_PARALLEL_LOOP()
+            ! total-flux matching (coarse side): add viscous momentum/energy face fluxes into creg, same
+            ! face gating and transverse offsets as the base capture; always accumulate.
+            if (viscous) then
+                $:GPU_PARALLEL_LOOP(collapse=3)
+                do t2 = 0, t2_hi
+                    do t1 = 0, t1_hi
+                        do eq = eqn_idx%mom%beg, eqn_idx%E
+                            select case (id)
+                            case (1)
+                                if (cap_lo) creg(1)%lo(eq, t1, t2) = creg(1)%lo(eq, t1, t2) + coef*real(flux_src%vf(eq)%sf(jlo, &
+                                    & o1 + t1, o2 + t2), wp)
+                                if (cap_hi) creg(1)%hi(eq, t1, t2) = creg(1)%hi(eq, t1, t2) + coef*real(flux_src%vf(eq)%sf(jhi, &
+                                    & o1 + t1, o2 + t2), wp)
+                            case (2)
+                                if (cap_lo) creg(2)%lo(eq, t1, t2) = creg(2)%lo(eq, t1, &
+                                    & t2) + coef*real(flux_src%vf(eq)%sf(o1 + t1, jlo, o2 + t2), wp)
+                                if (cap_hi) creg(2)%hi(eq, t1, t2) = creg(2)%hi(eq, t1, &
+                                    & t2) + coef*real(flux_src%vf(eq)%sf(o1 + t1, jhi, o2 + t2), wp)
+                            case (3)
+                                if (cap_lo) creg(3)%lo(eq, t1, t2) = creg(3)%lo(eq, t1, &
+                                    & t2) + coef*real(flux_src%vf(eq)%sf(o1 + t1, o2 + t2, jlo), wp)
+                                if (cap_hi) creg(3)%hi(eq, t1, t2) = creg(3)%hi(eq, t1, &
+                                    & t2) + coef*real(flux_src%vf(eq)%sf(o1 + t1, o2 + t2, jhi), wp)
+                            end select
+                        end do
+                    end do
+                end do
+                $:END_GPU_PARALLEL_LOOP()
+            end if
         end if
 
     end subroutine s_amr_capture_boundary_flux

src/simulation/m_checker.fpp

@@ -95,8 +95,8 @@ contains
             @:PROHIBIT(model_eqns /= 2, "amr requires model_eqns = 2 (5-equation)")
             @:PROHIBIT(num_fluids > 1 .and. .not. mpp_lim, &
                        & "amr with num_fluids > 1 requires mpp_lim (its volume-fraction clamp+renormalize maintains coarse/fine alpha consistency)")
-            @:PROHIBIT(viscous .or. surface_tension .or. hypoelasticity .or. hyperelasticity .or. mhd .or. chemistry, &
-                       & "amr does not support viscous/elastic/surface-tension/MHD/chemistry")
+            @:PROHIBIT(surface_tension .or. hypoelasticity .or. hyperelasticity .or. mhd .or. chemistry, &
+                       & "amr does not support elastic/surface-tension/MHD/chemistry")
             @:PROHIBIT(bubbles_euler .or. bubbles_lagrange .or. qbmm .or. relax .or. ib .or. igr .or. cyl_coord, &
                        & "amr does not support bubbles/phase-change/IB/IGR/cylindrical")
             @:PROHIBIT(active_box, "amr is incompatible with active_box (unvalidated combination)")

src/simulation/m_rhs.fpp

@@ -781,7 +781,7 @@ contains
                 ! AMR refluxing: record the c/f boundary-face fluxes exactly as the assembly above
                 ! used them (after s_cbc may have modified flux_n inside the advection source call);
                 ! must run before the next direction's sweep reuses the aliased flux_n storage
-                if (amr) call s_amr_capture_boundary_flux(id, flux_n(id), stage)
+                if (amr) call s_amr_capture_boundary_flux(id, flux_n(id), flux_src_n(id), stage)
 
                 ! RHS additions for hypoelasticity
                 call nvtxStartRange("RHS-HYPOELASTICITY")