marker_set_ComputeSolidVolFrac Subroutine

    impure subroutine marker_set_ComputeSolidVolFrac (this,VF,solver,MaxTol,MaxIt,RelaxType,Rel,It,intRHS)
      !> Computes the solid volume fraction on the mesh.
      use leapDiffOp
      implicit none
      class(marker_set),    intent(inout) :: this                              !! A collection of tessellation elements
      type(eulerian_obj_r), intent(inout) :: VF                                !! Volume fraction
      character(len=*),     intent(in)    :: solver                            !! Name of solver to be used
      real(wp),             intent(in)    :: MaxTol                            !! Maximum relative tolerance
      integer,              intent(in)    :: MaxIt                             !! Maximum number of subiterations
      integer,              intent(in),    &
                                 optional :: RelaxType                         !! Relaxation type
      real(wp),             intent(out),   &
                                 optional :: Rel                               !! Relative error at end of solve
      integer,              intent(out),   &
                                 optional :: It                                !! Number of iterations performed
      real(wp),             intent(out),   &
                                 optional :: intRHS                            !! Magnitude of RHS
      ! Work variables
      type(hypre_obj)      :: hypre
      type(op_obj)         :: op
      type(eulerian_obj_r) :: N1,  N2,  N3
      integer              :: i,j,k
      real(wp)             :: dx,dy,dz,Vol,buffr

      ! Initialize HYPRE
      ! -------------------------- !
      call hypre%Initialize(this%block,this%parallel)

      ! Chose solver to use
      if (present(RelaxType)) then
        call hypre%SelectSolver(solver,MaxTol=MaxTol,MaxIt=MaxIt,RelaxType=RelaxType)
      else
        call hypre%SelectSolver(solver, MaxTol=MaxTol,MaxIt=MaxIt)
      end if

      ! Setup the solver
      call hypre%Setup()

      ! Setting up right hand side
      ! -------------------------- !
      call N1%Initialize('N1',this%block,this%parallel,0)
      call N2%Initialize('N2',this%block,this%parallel,0)
      call N3%Initialize('N3',this%block,this%parallel,0)

      ! Filter normals on mesh
      call this%Filter('N1',N1)
      call this%Filter('N2',N2)
      call this%Filter('N3',N3)

      ! Compute the divergence
      call op%Initialize(this%block,this%parallel)
      VF=op%div(VF%name,N1,N2,N3)
      VF%cell = - VF%cell

      call N1%Finalize()
      call N2%Finalize()
      call N3%Finalize()

      ! Setting up right hand side
      ! -------------------------- !
      if (present(intRHS)) then
        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)
              dx = this%block%x(i+1)-this%block%x(i)
              dy = this%block%y(j+1)-this%block%y(j)
              dz = this%block%z(k+1)-this%block%z(k)
              intRHS = intRHS + VF%cell(i,j,k)* dx*dy*dz
              Vol= Vol + dx*dy*dz
            end do
          end do
        end do
        call this%parallel%Sum(Vol,    buffr); Vol   = buffr
        call this%parallel%Sum(intRHS, buffr); intRHS= buffr/Vol
      end if

      call hypre%SetRHS(VF)

      ! Solve Poisson equation
      ! -------------------------- !
      call hypre%Solve(VF)

      ! Update ghostcells
      call VF%UpdateGhostCells

      if (present(Rel)) Rel=hypre%rel
      if (present(It))  It=hypre%it

      ! Throw warning if appropriate
      if (hypre%rel > MaxTol .and. this%parallel%RankIsRoot()) &
            write(*,fmt='(a,ES12.5)') "Warning: Volume Fraction solver did not converge -- Rel. Res.=",hypre%rel

      ! Finalize
      call hypre%Finalize()
      call op%Finalize()

      return
    end subroutine marker_set_ComputeSolidVolFrac