compute_scalar_gradients Subroutine

public subroutine compute_scalar_gradients(mesh, phi, grad_phi)

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.

Arguments

Type IntentOptional Attributes Name
type(mesh_t), intent(in) :: mesh
real(kind=rk), intent(in) :: phi(:)
real(kind=rk), intent(out) :: grad_phi(:,:)

Source Code

   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