grans_obj_PrepareSolverOperators Subroutine

    subroutine grans_obj_PrepareSolverOperators(this)
      !> Prepares operators used by GRANS.
      ! This subroutine also builds the volume fraction Poisson operator
      ! with the specified boundary conditions.
      use leapDiffOp
      implicit none
      class(grans_obj), intent(inout) :: this                                  !! GRANS solver
      ! Work variables
      integer               :: scheme_order
      character(str64)      :: hypre_solver
      integer               :: hypre_max_it
      real(wp)              :: hypre_res_tol
      real(wp), allocatable :: GRAD(:,:,:,:,:)
      real(wp), allocatable :: DIV (:,:,:,:,:)
      type(extent_obj)      :: extents
      integer               :: bcdir
      integer               :: side
      integer               :: slo(3)
      integer               :: shi(3)
      integer               :: i,j,k,m,n

      ! Get information from input
      call this%parser%Get("Scheme order",       scheme_order,  default = 2            )
      call this%parser%Get('VF solver',          hypre_solver , default = "IJ-AMG-NONE")
      call this%parser%Get('VF Max Iterm',       hypre_max_it , default = 30           )
      call this%parser%Get('VF Convergence Tol', hypre_res_tol, default = 1.0e-10_wp   )

      ! Initialize basic differential operators
      call this%op%Initialize(this%block,this%parallel,order=scheme_order)

      ! Intialize volume fraction solver
      call this%VFSolver%Initialize(this%block,this%parallel)
      call this%VFSolver%SelectSolver(hypre_solver,MaxTol=hypre_res_tol,MaxIt=hypre_max_it)

      ! Use differential operators to build the laplacian inside HYPRE
      associate (lo => this%block%lo, hi => this%block%hi, st => this%op%st, bcs => this%bcs)
        allocate(GRAD(-st:st-1,lo(1)-st+1:hi(1)+st,lo(2)-st+1:hi(2)+st,lo(3)-st+1:hi(3)+st,3))
        allocate(DIV (-st+1:st,lo(1)     :hi(1)   ,lo(2)     :hi(2)   ,lo(3)     :hi(3)   ,3))
        GRAD = 0.0_wp
        DIV  = 0.0_wp

        ! Define GRAD
        do k=lo(3)-st+1,hi(3)+st
          do j=lo(2)-st+1,hi(2)+st
            do i=lo(1)-st+1,hi(1)+st
              do m=-st,st-1
                GRAD(m,i,j,k,1) = GRAD(m,i,j,k,1) + this%op%c_d1dx1m(m,i)
                GRAD(m,i,j,k,2) = GRAD(m,i,j,k,2) + this%op%c_d1dx2m(m,j)
                GRAD(m,i,j,k,3) = GRAD(m,i,j,k,3) + this%op%c_d1dx3m(m,k)
              end do
            end do
          end do
        end do

        ! Apply BC
        do n=1,bcs%count
          ! Check whether we have a boundary condition on this region
          if (bcs%CheckBCExists(bcs%region(n)%name, 'ibVF')) then

            select case ( bcs%GetBCType(bcs%region(n)%name, 'ibVF') )
            case (BC_NEUMANN, BC_SYMMETRY)
                ! Get direction and side of BC
                call bcs%GetSideDirByRegion(bcs%region(n)%name,side,bcdir)

                ! Get extents
                extents = bcs%GetExtents(bcs%region(n)%name)

                slo = extents%lo
                shi = extents%hi

                select case (side)
                case (BC_LEFT)
                  slo(bcdir) = slo(bcdir) - min(st-1,this%block%ngc)
                case (BC_RIGHT)
                  slo(bcdir) = slo(bcdir) + 1
                  shi(bcdir) = shi(bcdir) + min(st,this%block%ngc)
                end select

                do k=slo(3),shi(3)
                  do j=slo(2),shi(2)
                    do i=slo(1),shi(1)
                      GRAD(:,i,j,k,bcdir) = 0.0_wp
                    end do
                  end do
                end do

            case default
              ! Do nothing
            end select
          end if
        end do

        ! Define DIV
        do k=lo(3),hi(3)
          do j=lo(2),hi(2)
            do i=lo(1),hi(1)
              do n=-st+1,st
                DIV(n,i,j,k,1) = DIV(n,i,j,k,1) + this%op%c_d1dx1(n,i)
                DIV(n,i,j,k,2) = DIV(n,i,j,k,2) + this%op%c_d1dx2(n,j)
                DIV(n,i,j,k,3) = DIV(n,i,j,k,3) + this%op%c_d1dx3(n,k)
              end do
            end do
          end do
        end do

      end associate

      this%VFSolver%stm = 2*this%op%st-1
      associate (lo => this%block%lo, hi => this%block%hi, stm =>this%VFSolver%stm)
        allocate(this%VFSolver%mat(lo(1):hi(1),lo(2):hi(2),lo(3):hi(3),3,-stm:stm),source=0.0_wp)
      end associate

      associate (lo => this%block%lo, hi => this%block%hi, st => this%op%st)
        do k=lo(3),hi(3)
          do j=lo(2),hi(2)
            do i=lo(1),hi(1)
              do n= -st+1,st
                do m=-st,st-1
                 this%VFSolver%mat(i,j,k,1,n+m) = this%VFSolver%mat(i,j,k,1,n+m) + DIV(n,i,j,k,1)*GRAD(m,i+n,j,k,1)
                 this%VFSolver%mat(i,j,k,2,n+m) = this%VFSolver%mat(i,j,k,2,n+m) + DIV(n,i,j,k,2)*GRAD(m,i,j+n,k,2)
                 this%VFSolver%mat(i,j,k,3,n+m) = this%VFSolver%mat(i,j,k,3,n+m) + DIV(n,i,j,k,3)*GRAD(m,i,j,k+n,3)
                end do
              end do

            end do
          end do
        end do
      end associate

      deallocate(GRAD,DIV)

      call this%VFSolver%Setup()

      return
    end subroutine grans_obj_PrepareSolverOperators