ensure_pressure_operator_cache Subroutine

private subroutine ensure_pressure_operator_cache(mesh, flow, bc)

Pre-computes Laplacian coefficients for the Poisson operator.

Arguments

Type IntentOptional Attributes Name
type(mesh_t), intent(in) :: mesh
type(flow_mpi_t), intent(in) :: flow
type(bc_set_t), intent(in) :: bc

Calls

proc~~ensure_pressure_operator_cache~~CallsGraph proc~ensure_pressure_operator_cache mod_flow_projection::ensure_pressure_operator_cache mpi_allreduce mpi_allreduce proc~ensure_pressure_operator_cache->mpi_allreduce proc~boundary_pressure_type mod_bc::boundary_pressure_type proc~ensure_pressure_operator_cache->proc~boundary_pressure_type proc~face_effective_neighbor mod_bc::face_effective_neighbor proc~ensure_pressure_operator_cache->proc~face_effective_neighbor proc~face_normal_distance mod_flow_projection::face_normal_distance proc~ensure_pressure_operator_cache->proc~face_normal_distance proc~fatal_error mod_kinds::fatal_error proc~ensure_pressure_operator_cache->proc~fatal_error proc~finalize_flow_projection_workspace mod_flow_projection::finalize_flow_projection_workspace proc~ensure_pressure_operator_cache->proc~finalize_flow_projection_workspace proc~patch_type_for_face mod_bc::patch_type_for_face proc~ensure_pressure_operator_cache->proc~patch_type_for_face proc~is_periodic_face mod_bc::is_periodic_face proc~face_effective_neighbor->proc~is_periodic_face proc~face_normal_distance->proc~fatal_error proc~face_normal_distance->proc~patch_type_for_face proc~outward_normal mod_flow_projection::outward_normal proc~face_normal_distance->proc~outward_normal proc~is_periodic_face->proc~patch_type_for_face

Called by

proc~~ensure_pressure_operator_cache~~CalledByGraph proc~ensure_pressure_operator_cache mod_flow_projection::ensure_pressure_operator_cache proc~advance_projection_step mod_flow_projection::advance_projection_step proc~advance_projection_step->proc~ensure_pressure_operator_cache proc~solve_pressure_correction mod_flow_projection::solve_pressure_correction proc~advance_projection_step->proc~solve_pressure_correction proc~pressure_matvec mod_flow_projection::pressure_matvec proc~pressure_matvec->proc~ensure_pressure_operator_cache proc~solve_pressure_correction->proc~ensure_pressure_operator_cache proc~solve_pressure_correction->proc~pressure_matvec program~lowmach_react_hex lowmach_react_hex program~lowmach_react_hex->proc~advance_projection_step

Source Code

   subroutine ensure_pressure_operator_cache(mesh, flow, bc)
      type(mesh_t), intent(in) :: mesh
      type(flow_mpi_t), intent(in) :: flow
      type(bc_set_t), intent(in) :: bc
      integer :: ierr, c, lf, f, nb
      real(rk) :: dist, coeff
      integer :: local_dirichlet_flag, global_dirichlet_flag

      if (pressure_cache%initialized) then
         if (pressure_cache%ncells == mesh%ncells) return
         call finalize_flow_projection_workspace()
      end if

      pressure_cache%ncells = mesh%ncells
      pressure_cache%max_faces = size(mesh%cell_faces, 1)

      allocate(pressure_cache%nb(pressure_cache%max_faces, mesh%ncells))
      allocate(pressure_cache%coeff(pressure_cache%max_faces, mesh%ncells))
      allocate(pressure_cache%diag(mesh%ncells))

      pressure_cache%nb = 0
      pressure_cache%coeff = zero
      pressure_cache%diag = zero
      pressure_cache%has_neumann_outlet = .false.

      local_dirichlet_flag = 0

      do c = 1, mesh%ncells
         do lf = 1, mesh%ncell_faces(c)
            f = mesh%cell_faces(lf, c)
            nb = face_effective_neighbor(mesh, bc, f, c)

            if (nb <= 0) then
               if (patch_type_for_face(mesh, bc, f) == bc_neumann) then
                  pressure_cache%has_neumann_outlet = .true.
               end if

               if (boundary_pressure_type(mesh, bc, f) == bc_dirichlet) then
                  dist = face_normal_distance(mesh, bc, f, c, 0)
                  coeff = mesh%faces(f)%area / dist
                  pressure_cache%diag(c) = pressure_cache%diag(c) + coeff
                  local_dirichlet_flag = 1
               end if
               cycle
            end if

            dist = face_normal_distance(mesh, bc, f, c, nb)
            coeff = mesh%faces(f)%area / dist

            pressure_cache%nb(lf, c) = nb
            pressure_cache%coeff(lf, c) = coeff
            pressure_cache%diag(c) = pressure_cache%diag(c) + coeff
         end do

         if (pressure_cache%diag(c) <= tiny_safe) then
            call fatal_error('flow', 'non-positive cached pressure diagonal')
         end if
      end do

      call MPI_Allreduce(local_dirichlet_flag, global_dirichlet_flag, 1, MPI_INTEGER, MPI_MAX, flow%comm, ierr)
      pressure_cache%has_dirichlet_pressure = (global_dirichlet_flag > 0)

      pressure_cache%initialized = .true.
   end subroutine ensure_pressure_operator_cache