Merge pull request #610 from ElmerCSC/HypreUpdate

Hypre update

fem/src/CRSMatrix.F90

@@ -1805,7 +1805,11 @@ FUNCTION CRS_Transpose( A ) RESULT(B)
        CALL Info('CRS_Transpose','Creating a transpose of matrix',Level=20)
 
        B => AllocateMatrix()
-
+
+       IF(.NOT. ASSOCIATED(A) ) THEN
+         CALL Fatal('CRS_Transpose','Matrix not associated!')
+       END IF
+
        na = A % NumberOfRows
        IF( na == 0 ) THEN
          B % NumberOfRows = 0
@@ -1824,23 +1828,23 @@ FUNCTION CRS_Transpose( A ) RESULT(B)
          ALLOCATE( B % Diag(nb) )       
          B % Diag = 0
        END IF
-
+
+       ! Count how many hits there are in A for each column
        Row = 0       
        DO i = 1, NVals
          Row( A % Cols(i) ) = Row( A % Cols(i) ) + 1
        END DO
-
+
+       ! For transpose the column hits are row hits
        B % Rows = 0
        B % Rows(1) = 1
        DO i = 1, nB
          B % Rows(i+1) = B % Rows(i) + Row(i)
        END DO
        B % Cols = 0
 
-       DO i = 1, nB
-         Row(i) = B % Rows(i)
-       END DO
-
+       ! Location of 1st entry in each row
+       Row(1:nB) = B % Rows(1:nB)
 
        DO i = 1, nA
 
@@ -1852,7 +1856,7 @@ FUNCTION CRS_Transpose( A ) RESULT(B)
              B % Values( Row(k) ) = A % Values(j)
              Row(k) = Row(k) + 1
            ELSE
-             WRITE( Message, * ) 'Trying to access non-existent column', i,k,j
+             WRITE( Message, * ) 'Trying to access column beyond allocation: ', i,k,j
              CALL Error( 'CRS_Transpose', Message )
              RETURN
            END IF

fem/src/ElementUtils.F90

@@ -299,6 +299,7 @@ END SUBROUTINE SolveTrilinos4
 #ifdef HAVE_HYPRE
      IF (Matrix % Hypre /= 0) THEN
        CALL SolveHypre4(Matrix % Hypre)
+       Matrix % Hypre = 0
      END IF
 #endif
 

fem/src/Lists.F90

@@ -4307,36 +4307,41 @@ END FUNCTION ListGetConstReal
 !> Returns a scalar real value, that may depend on other scalar values such as 
 !> time or timestep size etc.
 !------------------------------------------------------------------------------
-  RECURSIVE FUNCTION ListGetCReal( List, Name, Found, minv, maxv, UnfoundFatal) RESULT(s)
+  RECURSIVE FUNCTION ListGetCReal( List, Name, Found, minv, maxv, UnfoundFatal, DefValue ) RESULT(s)
 !------------------------------------------------------------------------------
      TYPE(ValueList_t), POINTER :: List
      CHARACTER(LEN=*) :: Name
      REAL(KIND=dp), OPTIONAL :: minv,maxv
      LOGICAL, OPTIONAL :: Found,UnfoundFatal
      INTEGER, TARGET :: Dnodes(1)
      INTEGER, POINTER :: NodeIndexes(:)
-
+     REAL(KIND=dp), OPTIONAL :: DefValue
+
      REAL(KIND=dp) :: s
      REAL(KIND=dp) :: x(1)
      TYPE(Element_t), POINTER :: Element
-
+     LOGICAL :: LFound 
+
      INTEGER :: n, istat
 
-     IF ( PRESENT( Found ) ) Found = .FALSE.
-
+     LFound = .FALSE.
      NodeIndexes => Dnodes
      n = 1
      NodeIndexes(n) = 1
 
      x = 0.0_dp
      IF ( ASSOCIATED(List % head) ) THEN
-        IF ( PRESENT( Found ) ) THEN
-           x(1:n) = ListGetReal( List, Name, n, NodeIndexes, Found, minv=minv, maxv=maxv, UnfoundFatal=UnfoundFatal )
-        ELSE
-           x(1:n) = ListGetReal( List, Name, n, NodeIndexes, minv=minv, maxv=maxv, UnfoundFatal=UnfoundFatal)
-        END IF
+       x(1:n) = ListGetReal( List, Name, n, NodeIndexes, LFound, minv=minv, maxv=maxv, &
+           UnfoundFatal=UnfoundFatal )
      END IF
      s = x(1)
+
+     IF( PRESENT( DefValue ) ) THEN
+       IF(.NOT. LFound ) s = DefValue
+     END IF
+
+     IF ( PRESENT( Found ) ) Found = LFound
+
 !------------------------------------------------------------------------------
   END FUNCTION ListGetCReal
 !------------------------------------------------------------------------------

fem/src/Multigrid.F90

@@ -918,6 +918,8 @@ END SUBROUTINE BlockSolveExt
 !------------------------------------------------------------------------------
        tt = CPUTime()
 
+       CALL Info('PMGSolve','Solving multigrid Level '//I2S(Level),Level=20)
+
 !
 !      Initialize:
 !      -----------
@@ -936,7 +938,6 @@ END SUBROUTINE BlockSolveExt
 !      If at lowest level, solve directly:
 !      -----------------------------------
        IF ( Level <= 1 ) THEN
-
           CALL ListPushNamespace('mglowest:')
 
           CALL ListAddLogical( Params,'mglowest: Linear System Free Factorization', .FALSE. )
@@ -955,10 +956,21 @@ END SUBROUTINE BlockSolveExt
             IF ( LIter ) LowestSolver='iterative'
           END IF
 
+          ! This fixes an issue that comes from having different calling convention for
+          ! serial Hypre.
+          IF(.NOT. Parallel) THEN
+            IF(ListGetLogical(Params, 'Linear System Use Hypre', Found) &
+                .AND. LowestSolver == 'iterative') LowestSolver = 'hypre'
+          END IF
+
           CALL Info('PMGSolve','Starting lowest linear solver using: '//TRIM(LowestSolver),Level=10 )
 
           SELECT CASE(LowestSolver)
 
+            ! This one only for serial!
+          CASE('hypre')
+            CALL SolveHypre( Matrix1, Solution, ForceVector, Solver )
+
           CASE('block')
             CALL BlockSolveExt( Matrix1, Solution, ForceVector, Solver )
 

fem/src/ParallelUtils.F90

@@ -121,8 +121,15 @@ SUBROUTINE ParallelInitMatrix( Solver, Matrix, inPerm )
        Mesh => Solver % Mesh
        DOFs = Solver % Variable % DOFs
 
-       Perm => Solver % Variable % Perm
-       IF(PRESENT(inPerm)) Perm=>InPerm
+       IF(PRESENT(inPerm)) THEN
+         Perm => InPerm
+       ELSE
+         Perm => Solver % Variable % Perm
+       END IF
+       IF(.NOT. ASSOCIATED(Perm)) THEN
+         CALL Fatal('ParallelInitMatrix','Cannot initialize matrix without Perm vector!')
+       END IF
+
 
        n = SIZE(Perm)
        k = n*DOFs + Matrix % ExtraDOFs
@@ -155,7 +162,7 @@ SUBROUTINE ParallelInitMatrix( Solver, Matrix, inPerm )
        Matrix % INVPerm = 0
        DO i=1,SIZE(Matrix % Perm)
           IF ( Matrix % Perm(i) /= 0 ) THEN
-             Matrix % INVPerm(Matrix % Perm(i)) = i
+             Matrix % InvPerm(Matrix % Perm(i)) = i
           END IF
        END DO