hypre_obj_SolveS Subroutine

    impure subroutine hypre_obj_SolveS(this,sol)
      !> Solves the system Ax=b and return the solution, Struct interface.
      implicit none
      class(hypre_obj),     intent(inout) :: this                              !! Hypre machinery
      type(eulerian_obj_r), intent(inout) :: sol                               !! Solution vector
      ! Work variables
      integer :: ierr
      integer :: nval
      integer :: i,j,k,m
      real(wp), allocatable :: values(:)

      ! Number of grid points
      nval=(this%block%hi(3)-this%block%lo(3)+1) &
          *(this%block%hi(2)-this%block%lo(2)+1) &
          *(this%block%hi(1)-this%block%lo(1)+1)

      select case (this%solver_name)
      case (HYPRE_SOL_S_PCG_NONE)
        call HYPRE_StructPCGSolve(this%solver,this%matrix,this%rhs,this%sol, ierr)
        call HYPRE_StructPCGGetNumIterations(this%solver,this%it,ierr)
        call HYPRE_StructPCGGetFinalRelative(this%solver,this%rel,ierr)
      case (HYPRE_SOL_S_SMG_NONE)
        call HYPRE_StructSMGSolve(this%solver,this%matrix,this%rhs,this%sol, ierr)
        call HYPRE_StructSMGGetNumIterations(this%solver,this%it,ierr)
        call HYPRE_StructSMGGetFinalRelative(this%solver,this%rel,ierr)
      case (HYPRE_SOL_S_PFMG_NONE)
        call HYPRE_StructPFMGSolve(this%solver,this%matrix,this%rhs,this%sol, ierr)
        call HYPRE_StructPFMGGetNumIteration(this%solver,this%it,ierr)
        call HYPRE_StructPFMGGetFinalRelativ(this%solver,this%rel,ierr)
      case default
        call this%parallel%Stop("HYPRE: Unknown solver")
      end select

      allocate(values(nval))

      ! Get the values
      call HYPRE_StructVectorGetBoxValues(this%sol,this%block%lo(1:this%dim),&
           this%block%hi(1:this%dim),values(1:nval),ierr)
      ! Loop over cell following Hypre convention (from bottom left, to top right)
      m=1
      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)
            sol%cell(i,j,k)=values(m)
            m=m+1
          end do
        end do
      end do

      deallocate(values)

      return
    end subroutine hypre_obj_SolveS