cdifs_obj_ComputeSolidVF Subroutine

    subroutine cdifs_obj_ComputeSolidVF(this)
      !> Computes the solid volume fraction (ibVF) by solving a Poisson
      ! equation.
      ! This subroutine computes the IB normals fields, applies approriate
      ! boundary conditions, and takes the divergence, which forms the
      ! right-hand side of the poisson equation.
      ! It then uses a HYPRE solver to solve this equation.
      ! In case where the Poisson operator is singular and using resolved
      ! particles, this subroutine will apply a correction that sets the solid
      ! volume fraction at particle centroids.
      implicit none
      class(cdifs_obj), intent(inout) :: this                                  !! CDIFS solver
      ! Work variables
      integer :: i,j,k
      real(wp):: buffr
      real(wp):: vol
      real(wp):: intRHS
      real(wp):: dV

      ! Update normals
      this%ibN(1) = 0.0_wp; this%ibN(2) = 0.0_wp; this%ibN(3) = 0.0_wp
      if (this%use_RP) call this%RP%UpdateNormals(this%ibN)
      if (this%use_IB) call this%IB%UpdateNormals(this%ibN)
      call this%bcs%UpdateBoundary(this%ibN)

      ! Build right hand side of the equation and compute its integral.
      ! If no Dirichlet conditions are imposed, this needs to be close
      ! to machine precision
      intRHS = 0.0_wp
      vol    = 0.0_wp
      do k=this%block%lo(3),this%block%hi(3)
        do j=this%block%lo(2),this%block%hi(2)
          do i=this%block%lo(1),this%block%hi(1)
            this%rhs%cell(i,j,k) = - dot_product(this%op%c_d1dx1(:,i),this%ibN(1)%cell(i-this%op%st+1:i+this%op%st,j,k)) &
                                   - dot_product(this%op%c_d1dx2(:,j),this%ibN(2)%cell(i,j-this%op%st+1:j+this%op%st,k)) &
                                   - dot_product(this%op%c_d1dx3(:,k),this%ibN(3)%cell(i,j,k-this%op%st+1:k+this%op%st))

            dV     = this%block%axis(1)%dxm(i)*this%block%axis(2)%dxm(j)*this%block%axis(3)%dxm(k)
            vol    = vol + dV
            intRHS = intRHS + this%rhs%cell(i,j,k)*dV
          end do
        end do
      end do
      call this%parallel%Sum(vol,   buffr); vol    = buffr
      call this%parallel%Sum(intRHS,buffr); intRHS = buffr/vol

      ! Add Dirichlet contributions to RHS
      call this%op%ApplyLaplacianDC(this%rhs,this%bcs,'ibVF')

      ! Send RHS to HYPRE
      call this%VFSolver%SetRHS(this%rhs)

      ! Send RHS to HYPRE
      ! Solve Poisson equation with HYPRE
      call this%VFSolver%Solve(this%ibVF)

      !if (this%parallel%RankIsRoot()) then
      !  print*, 'iter:',     this%timer%iter
      !  print*, 'time:',     this%timer%time
      !  print*, 'Rel:',      this%VFSolver%rel
      !  print*, 'intRHS:',   intRHS
      !  print*, 'HYPRE it:', this%VFSolver%it
      !end if

      if (this%singularVF.and.this%use_RP) &
        this%ibVF%cell = this%ibVF%cell + (1.0_wp - this%RP%GetCentroidVF(this%ibVF))

      call this%bcs%UpdateBoundary(this%ibVF)

      return
    end subroutine cdifs_obj_ComputeSolidVF