scalar_face_value Function

public function scalar_face_value(mesh, phi, grad_phi, scheme_name, limiter_name, c, neighbor, fid, flux_out, boundary_value, is_dirichlet) result(phi_face)

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.

Arguments

Type IntentOptional Attributes Name
type(mesh_t), intent(in) :: mesh
real(kind=rk), intent(in) :: phi(:)
real(kind=rk), intent(in) :: grad_phi(:,:)
character(len=*), intent(in) :: scheme_name
character(len=*), intent(in) :: limiter_name
integer, intent(in) :: c
integer, intent(in) :: neighbor
integer, intent(in) :: fid
real(kind=rk), intent(in) :: flux_out
real(kind=rk), intent(in) :: boundary_value
logical, intent(in) :: is_dirichlet

Return Value real(kind=rk)


Source Code

   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