mod_reconstruction.f90 – LowMachReact-Hex

Source Code

!> Reusable scalar face reconstruction helpers.
!!
!! The routines here intentionally avoid equation-specific state.  Species,
!! enthalpy, and future transported scalars can share the same menu of face
!! interpolation choices:
!! - upwind
!! - central
!! - bounded_central / bounded_linear
!! - muscl / limited_linear / linear_upwind / higher_order_bounded
!!
!! The high-order branch uses a compact cell-gradient reconstruction from the
!! upwind cell and clamps the reconstructed face value to local cell-neighbor
!! bounds.  This is not a full FCT implementation, but it gives a robust,
!! conservative, bounded upgrade path over first-order upwind without widening
!! the MPI stencil.
module mod_reconstruction
   use mod_precision, only : rk, zero, tiny_safe, fatal_error, lowercase
   use mod_mesh, only : mesh_t
   implicit none

   private

   public :: compute_scalar_gradients
   public :: scalar_face_value
   public :: validate_advection_scheme
   public :: validate_limiter_name

contains

   !> Compute compact cell-centered gradients from immediately adjacent cells.
   !! Boundary faces without an effective/periodic neighbor are skipped; boundary
   !! values are still enforced by scalar_face_value at the boundary face.
   subroutine compute_scalar_gradients(mesh, phi, grad_phi)
      type(mesh_t), intent(in) :: mesh
      real(rk), intent(in) :: phi(:)
      real(rk), intent(out) :: grad_phi(:,:)

      integer :: c, lf, fid, nb, ncontrib
      real(rk) :: dx(3), denom

      if (size(phi) < mesh%ncells) then
         call fatal_error('reconstruction', 'field length is smaller than mesh%ncells')
      end if
      if (size(grad_phi, 1) /= 3 .or. size(grad_phi, 2) < mesh%ncells) then
         call fatal_error('reconstruction', 'gradient array must have shape (3,ncells)')
      end if

      grad_phi = zero

      do c = 1, mesh%ncells
         ncontrib = 0
         do lf = 1, mesh%ncell_faces(c)
            fid = mesh%cell_faces(lf, c)
            nb = adjacent_cell(mesh, fid, c)
            if (nb <= 0) cycle

            dx = mesh%cells(nb)%center - mesh%cells(c)%center
            denom = max(dot_product(dx, dx), tiny_safe)
            grad_phi(:, c) = grad_phi(:, c) + (phi(nb) - phi(c)) * dx / denom
            ncontrib = ncontrib + 1
         end do

         if (ncontrib > 0) grad_phi(:, c) = grad_phi(:, c) / real(ncontrib, rk)
      end do
   end subroutine compute_scalar_gradients


   !> Reconstruct a scalar value at a face using the requested scheme.
   !! flux_out is oriented outward from cell c.  neighbor is the effective
   !! neighbor seen from c (0 for non-periodic boundary).  For Dirichlet
   !! boundaries, boundary_value is used as the exterior state.
   function scalar_face_value(mesh, phi, grad_phi, scheme_name, limiter_name, c, neighbor, fid, &
                              flux_out, boundary_value, is_dirichlet) result(phi_face)
      type(mesh_t), intent(in) :: mesh
      real(rk), intent(in) :: phi(:)
      real(rk), intent(in) :: grad_phi(:,:)
      character(len=*), intent(in) :: scheme_name
      character(len=*), intent(in) :: limiter_name
      integer, intent(in) :: c, neighbor, fid
      real(rk), intent(in) :: flux_out
      real(rk), intent(in) :: boundary_value
      logical, intent(in) :: is_dirichlet
      real(rk) :: phi_face

      character(len=len(scheme_name)) :: scheme
      character(len=len(limiter_name)) :: limiter
      real(rk) :: phi_cell, phi_other, candidate
      integer :: upwind_cell
      logical :: clamp_to_bounds

      if (c < 1 .or. c > mesh%ncells) call fatal_error('reconstruction', 'invalid current cell index')
      if (size(phi) < mesh%ncells) call fatal_error('reconstruction', 'field length is smaller than mesh%ncells')
      if (size(grad_phi, 1) /= 3 .or. size(grad_phi, 2) < mesh%ncells) then
         call fatal_error('reconstruction', 'gradient array must have shape (3,ncells)')
      end if

      scheme = trim(lowercase(scheme_name))
      limiter = trim(lowercase(limiter_name))
      if (len_trim(scheme) == 0) scheme = 'upwind'
      if (len_trim(limiter) == 0) limiter = 'barth_jespersen'

      phi_face = zero
      phi_cell = phi(c)
      if (neighbor > 0) then
         phi_other = phi(neighbor)
      else if (is_dirichlet) then
         phi_other = boundary_value
      else
         phi_other = phi_cell
      end if

      select case (trim(scheme))
      case ('upwind', 'first_order_upwind', 'first-order-upwind')
         if (flux_out >= zero) then
            phi_face = phi_cell
         else
            phi_face = phi_other
         end if

      case ('central', 'central_difference', 'central-difference')
         phi_face = 0.5_rk * (phi_cell + phi_other)

      case ('bounded_central', 'bounded-central', 'bounded_linear', 'bounded-linear', 'deferred_central', 'deferred-correction')
         candidate = 0.5_rk * (phi_cell + phi_other)
         phi_face = clamp_value(candidate, min(phi_cell, phi_other), max(phi_cell, phi_other))

      case ('muscl', 'limited_linear', 'limited-linear', 'linear_upwind', 'linear-upwind', 'higher_order_bounded', 'higher-order-bounded')
         if (flux_out >= zero) then
            upwind_cell = c
            if (neighbor > 0) then
               select case (trim(limiter))
               case ('none')
                  phi_face = phi(upwind_cell) + dot_product(grad_phi(:, upwind_cell), &
                                                            mesh%faces(fid)%center - mesh%cells(upwind_cell)%center)
               case ('barth_jespersen', 'barth-jespersen')
                  block
                     real(rk) :: alpha
                     alpha = compute_bj_coefficient(mesh, phi, grad_phi, upwind_cell)
                     phi_face = phi(upwind_cell) + alpha * dot_product(grad_phi(:, upwind_cell), &
                                                                       mesh%faces(fid)%center - mesh%cells(upwind_cell)%center)
                  end block
               case ('minmod', 'vanleer', 'van_leer', 'superbee', 'mc', 'monotonized_central')
                  block
                     real(rk) :: grad_corr, corr
                     grad_corr = dot_product(grad_phi(:, upwind_cell), mesh%faces(fid)%center - mesh%cells(upwind_cell)%center)
                     corr = compute_tvd_limited_correction(limiter, phi(upwind_cell), phi(neighbor), grad_corr)
                     phi_face = phi(upwind_cell) + corr
                  end block
               case default
                  candidate = phi(upwind_cell) + dot_product(grad_phi(:, upwind_cell), &
                                                             mesh%faces(fid)%center - mesh%cells(upwind_cell)%center)
                  phi_face = bounded_to_local_neighbors(mesh, phi, candidate, upwind_cell, boundary_value, &
                                                        neighbor <= 0 .and. is_dirichlet)
               end select
            else
               ! Boundary face outflow
               candidate = phi(upwind_cell) + dot_product(grad_phi(:, upwind_cell), &
                                                          mesh%faces(fid)%center - mesh%cells(upwind_cell)%center)
               phi_face = bounded_to_local_neighbors(mesh, phi, candidate, upwind_cell, boundary_value, &
                                                     neighbor <= 0 .and. is_dirichlet)
            end if
         else
            if (neighbor > 0) then
               upwind_cell = neighbor
               select case (trim(limiter))
               case ('none')
                  phi_face = phi(upwind_cell) + dot_product(grad_phi(:, upwind_cell), &
                                                            mesh%faces(fid)%center - mesh%cells(upwind_cell)%center)
               case ('barth_jespersen', 'barth-jespersen')
                  block
                     real(rk) :: alpha
                     alpha = compute_bj_coefficient(mesh, phi, grad_phi, upwind_cell)
                     phi_face = phi(upwind_cell) + alpha * dot_product(grad_phi(:, upwind_cell), &
                                                                       mesh%faces(fid)%center - mesh%cells(upwind_cell)%center)
                  end block
               case ('minmod', 'vanleer', 'van_leer', 'superbee', 'mc', 'monotonized_central')
                  block
                     real(rk) :: grad_corr, corr
                     grad_corr = dot_product(grad_phi(:, upwind_cell), mesh%faces(fid)%center - mesh%cells(upwind_cell)%center)
                     corr = compute_tvd_limited_correction(limiter, phi(upwind_cell), phi(c), grad_corr)
                     phi_face = phi(upwind_cell) + corr
                  end block
               case default
                  candidate = phi(upwind_cell) + dot_product(grad_phi(:, upwind_cell), &
                                                             mesh%faces(fid)%center - mesh%cells(upwind_cell)%center)
                  phi_face = bounded_to_local_neighbors(mesh, phi, candidate, upwind_cell, zero, .false.)
               end select
            else
               ! Inflow from physical boundary
               phi_face = phi_other
            end if
         end if

      case ('quick', 'hlpa', 'quick_limited', 'quick-limited')
         if (flux_out >= zero) then
            upwind_cell = c
            if (neighbor > 0) then
               candidate = phi_cell + 0.375_rk * (phi_other - phi_cell) - &
                           0.125_rk * dot_product(grad_phi(:, c), mesh%cells(neighbor)%center - mesh%cells(c)%center)
            else
               ! Outflow boundary face: use boundary extrapolation
               candidate = phi_cell + 0.375_rk * (phi_other - phi_cell) - &
                           0.25_rk * dot_product(grad_phi(:, c), mesh%faces(fid)%center - mesh%cells(c)%center)
            end if

            clamp_to_bounds = (trim(scheme) == 'hlpa' .or. trim(scheme) == 'quick_limited' .or. &
                               trim(scheme) == 'quick-limited' .or. trim(limiter) /= 'none')
            if (clamp_to_bounds) then
               phi_face = bounded_to_local_neighbors(mesh, phi, candidate, upwind_cell, boundary_value, &
                                                     neighbor <= 0 .and. is_dirichlet)
            else
               phi_face = candidate
            end if
         else
            if (neighbor > 0) then
               upwind_cell = neighbor
               candidate = phi_other + 0.375_rk * (phi_cell - phi_other) - &
                           0.125_rk * dot_product(grad_phi(:, neighbor), mesh%cells(c)%center - mesh%cells(neighbor)%center)

               clamp_to_bounds = (trim(scheme) == 'hlpa' .or. trim(scheme) == 'quick_limited' .or. &
                                  trim(scheme) == 'quick-limited' .or. trim(limiter) /= 'none')
               if (clamp_to_bounds) then
                  phi_face = bounded_to_local_neighbors(mesh, phi, candidate, upwind_cell, zero, .false.)
               else
                  phi_face = candidate
               end if
            else
               ! Inflow from physical boundary: the boundary value is the highest quality available
               phi_face = phi_other
            end if
         end if

      case default
         call fatal_error('reconstruction', 'unknown scalar advection scheme: '//trim(scheme))
      end select
   end function scalar_face_value


   !> Validate scalar/momentum advection names at input time.
   subroutine validate_advection_scheme(scheme_name, scope, allow_scalar_high_order)
      character(len=*), intent(in) :: scheme_name
      character(len=*), intent(in) :: scope
      logical, intent(in) :: allow_scalar_high_order
      character(len=len(scheme_name)) :: scheme

      scheme = trim(lowercase(scheme_name))
      if (len_trim(scheme) == 0) scheme = 'upwind'

      select case (trim(scheme))
      case ('upwind', 'first_order_upwind', 'first-order-upwind', &
            'central', 'central_difference', 'central-difference')
         return
      case ('bounded_central', 'bounded-central', 'bounded_linear', 'bounded-linear', &
            'deferred_central', 'deferred-correction', &
            'muscl', 'limited_linear', 'limited-linear', 'linear_upwind', 'linear-upwind', &
            'higher_order_bounded', 'higher-order-bounded', &
            'quick', 'hlpa', 'quick_limited', 'quick-limited')
         if (allow_scalar_high_order) return
      case default
         continue
      end select

      if (allow_scalar_high_order) then
         call fatal_error('input', trim(scope)//' must be one of: upwind, central, bounded_central, bounded_linear, muscl, limited_linear, quick, hlpa')
      else
         call fatal_error('input', trim(scope)//' must be one of: upwind, central')
      end if
   end subroutine validate_advection_scheme


   !> Validate limiter names.  The current bounded implementation maps these
   !! limiters to the same local-bounds clamp; the menu is preserved in input so
   !! future limiter-specific formulas can be added without changing cases.
   subroutine validate_limiter_name(limiter_name)
      character(len=*), intent(in) :: limiter_name
      character(len=len(limiter_name)) :: limiter

      limiter = trim(lowercase(limiter_name))
      if (len_trim(limiter) == 0) limiter = 'barth_jespersen'

      select case (trim(limiter))
      case ('none', 'minmod', 'vanleer', 'van_leer', 'mc', 'monotonized_central', &
            'superbee', 'barth_jespersen', 'barth-jespersen', 'venkatakrishnan', 'hlpa')
         return
      case default
         call fatal_error('input', 'scalar_limiter must be one of: none, minmod, vanleer, mc, superbee, barth_jespersen, venkatakrishnan, hlpa')
      end select
   end subroutine validate_limiter_name


   integer function adjacent_cell(mesh, fid, c) result(nb)
      type(mesh_t), intent(in) :: mesh
      integer, intent(in) :: fid, c

      nb = 0
      if (mesh%faces(fid)%owner == c) then
         nb = mesh%faces(fid)%neighbor
         if (nb <= 0) nb = mesh%faces(fid)%periodic_neighbor
      else if (mesh%faces(fid)%neighbor == c .or. mesh%faces(fid)%periodic_neighbor == c) then
         nb = mesh%faces(fid)%owner
      end if
   end function adjacent_cell


   real(rk) function bounded_to_local_neighbors(mesh, phi, candidate, cell, boundary_value, include_boundary) result(out)
      type(mesh_t), intent(in) :: mesh
      real(rk), intent(in) :: phi(:)
      real(rk), intent(in) :: candidate
      integer, intent(in) :: cell
      real(rk), intent(in) :: boundary_value
      logical, intent(in) :: include_boundary

      integer :: lf, fid, nb
      real(rk) :: lo, hi

      lo = phi(cell)
      hi = phi(cell)
      do lf = 1, mesh%ncell_faces(cell)
         fid = mesh%cell_faces(lf, cell)
         nb = adjacent_cell(mesh, fid, cell)
         if (nb > 0) then
            lo = min(lo, phi(nb))
            hi = max(hi, phi(nb))
         end if
      end do
      if (include_boundary) then
         lo = min(lo, boundary_value)
         hi = max(hi, boundary_value)
      end if

      out = clamp_value(candidate, lo, hi)
   end function bounded_to_local_neighbors


   pure real(rk) function clamp_value(x, lo, hi) result(y)
      real(rk), intent(in) :: x, lo, hi

      y = min(max(x, lo), hi)
   end function clamp_value


   real(rk) function compute_bj_coefficient(mesh, phi, grad_phi, c) result(alpha)
      type(mesh_t), intent(in) :: mesh
      real(rk), intent(in) :: phi(:)
      real(rk), intent(in) :: grad_phi(:,:)
      integer, intent(in) :: c

      integer :: lf, fid, nb
      real(rk) :: phi_cell, lo, hi, phi_f, diff, val
      
      phi_cell = phi(c)
      lo = phi_cell
      hi = phi_cell
      
      ! 1. Find the min/max among neighbors
      do lf = 1, mesh%ncell_faces(c)
         fid = mesh%cell_faces(lf, c)
         nb = adjacent_cell(mesh, fid, c)
         if (nb > 0) then
            lo = min(lo, phi(nb))
            hi = max(hi, phi(nb))
         end if
      end do
      
      ! 2. Compute the BJ limiter over all faces of cell c
      alpha = 1.0_rk
      do lf = 1, mesh%ncell_faces(c)
         fid = mesh%cell_faces(lf, c)
         phi_f = phi_cell + dot_product(grad_phi(:, c), mesh%faces(fid)%center - mesh%cells(c)%center)
         diff = phi_f - phi_cell
         if (diff > 1.0e-14_rk) then
            val = min(1.0_rk, (hi - phi_cell) / diff)
            alpha = min(alpha, val)
         else if (diff < -1.0e-14_rk) then
            val = min(1.0_rk, (lo - phi_cell) / diff)
            alpha = min(alpha, val)
         end if
      end do
      
      alpha = max(0.0_rk, min(1.0_rk, alpha))
   end function compute_bj_coefficient


   real(rk) function compute_tvd_limited_correction(limiter, phi_cell, phi_nb, grad_corr) result(corr)
      character(len=*), intent(in) :: limiter
      real(rk), intent(in) :: phi_cell, phi_nb, grad_corr

      real(rk) :: theta, psi, eps, denom

      eps = 1.0e-14_rk
      denom = 2.0_rk * grad_corr
      
      if (abs(denom) < eps) then
         corr = 0.0_rk
         return
      end if
      
      theta = (phi_nb - phi_cell) / denom
      
      select case (trim(limiter))
      case ('minmod')
         psi = max(0.0_rk, min(1.0_rk, theta))
      case ('vanleer', 'van_leer')
         psi = (theta + abs(theta)) / (1.0_rk + abs(theta))
      case ('superbee')
         psi = max(0.0_rk, max(min(2.0_rk*theta, 1.0_rk), min(theta, 2.0_rk)))
      case ('mc', 'monotonized_central')
         psi = max(0.0_rk, min(2.0_rk, min(2.0_rk*theta, 0.5_rk*(1.0_rk + theta))))
      case default
         psi = 1.0_rk
      end select
      
      corr = psi * grad_corr
   end function compute_tvd_limited_correction

end module mod_reconstruction
mod_reconstruction.f90 Source File

Contents

Modules

Source Code

Source Code
