Skip to content

Commit

Permalink
Hopefully final fix for AMS new interface. Added a pure AMS solver te…
Browse files Browse the repository at this point in the history
…st case.
  • Loading branch information
raback committed Nov 19, 2024
1 parent c4780c4 commit d903b4d
Show file tree
Hide file tree
Showing 11 changed files with 3,657 additions and 336 deletions.
77 changes: 42 additions & 35 deletions fem/src/SParIterSolver.F90
Original file line number Diff line number Diff line change
Expand Up @@ -942,10 +942,14 @@ SUBROUTINE SetHypreParameters(Params, hypremethod, Ilun, hypre_dppara, hypre_int
DO i=0,1

j = 0
IF(i==0) THEN
IterativeMethod = ListGetString( Params,'Linear System Iterative Method' )
IF(i==0) THEN
hypre_sol = ListGetInteger( Params,'Linear System Method Hypre Index',Found )
IF(Found ) CYCLE
Method => IterativeMethod
IterativeMethod = ListGetString( Params,'Linear System Iterative Method' )
ELSE
hypre_pre = ListGetInteger( Params,'Linear System Preconditioning Hypre Index',Found )
IF(Found ) CYCLE
PrecMethod = ListGetString(Params,'Linear System Preconditioning',Found)
IF(.NOT. Found) THEN
CALL Info(Caller,'"Linear System Preconditioning" not given, assuming "none"')
Expand Down Expand Up @@ -986,27 +990,27 @@ SUBROUTINE SetHypreParameters(Params, hypremethod, Ilun, hypre_dppara, hypre_int
END IF

IF(i==0) THEN
CALL Info(Caller,'Using HYPRE iterative method: '//TRIM(IterativeMethod),Level=7)
hypre_sol = j
ELSE
CALL Info(Caller,'Using HYPRE preconditioning method: '//TRIM(PrecMethod),Level=7)
hypre_pre = j
END IF
END DO

! Some methods (Krylov methods) can not act as preconditioners!
CALL Info(Caller,'Using HYPRE preconditioning method: '//TRIM(PrecMethod),Level=7)
IF( ANY( hypre_pre == [6,7,8,9,10,11] ) ) THEN
CALL Fatal(Caller,'Invalid preconditioner for Hypre: '//TRIM(PrecMethod))
END IF

! Some methods cannot act as solvers!
CALL Info(Caller,'Using HYPRE iterative method: '//TRIM(IterativeMethod),Level=7)
IF( ANY( hypre_sol == [0,4,5] ) ) THEN
CALL Fatal(Caller,'Invalid solver for Hypre: '//TRIM(IterativeMethod))
END IF

! We map the precondtioner + solver to one figure.
hypremethod = 100 * hypre_sol + hypre_pre
CALL Info(Caller,'Hypre method index: '//I2S(hypremethod),Level=8)
CALL Info(Caller,'Hypre method index: '//I2S(hypremethod),Level=6)

DO i=0,1
IF(i==0) THEN
Expand Down Expand Up @@ -1058,8 +1062,15 @@ SUBROUTINE SetHypreParameters(Params, hypremethod, Ilun, hypre_dppara, hypre_int
'BoomerAMG Strong Threshold', Found, DefValue = 0.25_dp)

CASE(2) ! AMS
hypre_intpara(1) = ListGetInteger( Params,&
'AMS Cycle Type', Found, DefValue = 1 )
! The numbering follows the old convention.
hypre_intpara(6) = ListGetInteger( Params,&
'AMS Smoothing Option 1', Found, DefValue = 2 )

hypre_intpara(3) = ListGetInteger( Params,&
'AMS Smoothing Option 2', Found, DefValue = 1 )

hypre_intpara(7) = ListGetInteger( Params,&
'AMS Cycle Type', Found, DefValue = 1)


CASE(3) ! ILU
Expand Down Expand Up @@ -1725,7 +1736,7 @@ SUBROUTINE SParIterSolver( SourceMatrix, ParallelInfo, XVec, &
TYPE (GlueTableT), POINTER :: GT

CHARACTER(:), ALLOCATABLE :: XmlFile
REAL(KIND=dp) :: hypre_dppara(10) = 0
REAL(KIND=dp) :: hypre_dppara(10) = 0, TOL
INTEGER :: ILUn, BILU, Rounds, buf(2), src, status(MPI_STATUS_SIZE), ssz,nob, &
hypremethod, gind, hypre_intpara(20) = 0

Expand Down Expand Up @@ -2171,7 +2182,7 @@ SUBROUTINE CreateHypreAMS(nrows,rows,cols,vals,n,grows,gcols,gvals, &
hypre_dppara(10), Gvals(*), xx_d(*), yy_d(*), zz_d(*)
INTEGER(KIND=C_INTPTR_T) :: hypreContainer
END SUBROUTINE CreateHypreAMS

SUBROUTINE UpdateHypre(TOL, hypremethod, hypreContainer) BIND(C,name="updatehypre")
USE, INTRINSIC :: iso_c_binding
REAL(KIND=c_double) :: TOL
Expand Down Expand Up @@ -2207,8 +2218,6 @@ END SUBROUTINE UpdateHypre
CALL ContinuousNumbering(ParallelInfo,Matrix % Perm,APerm,Owner)
! Newer hypre libraries require zero based indexing
Aperm = Aperm-1

CALL SParIterActiveBarrier()
ELSE
n = Matrix % NumberOfRows
ALLOCATE( Owner(n), Aperm(n) )
Expand All @@ -2217,7 +2226,10 @@ END SUBROUTINE UpdateHypre
END DO
Owner = 1
END IF


CALL SParIterActiveBarrier()


!------------------------------------------------------------
verbosity = ListGetInteger( CurrentModel % Simulation,'Max Output Level',Found )
IF( .NOT. Found ) verbosity = 10
Expand Down Expand Up @@ -2255,15 +2267,6 @@ END SUBROUTINE UpdateHypre
IF( DoAMS ) THEN
CALL PrepareHypreAMS()
nnd = Solver % Mesh % NumberOfNodes
IF( ListGetLogical(Solver % Values,'Hypre test',Found ) ) THEN
PRINT *,'Using old hypre for debugging!'
CALL CreateHYPREAMS( Matrix % NumberOfRows, Rows, Cols, Vals, &
nnd,GM % Rows,GM % Cols,GM % Values,Aperm,Aperm,Aperm,Owner, &
Bperm,NodeOwner,Xvec,RHSvec,ParEnv % myPE, ILUn, Rounds,TOL, &
xx_d,yy_d,zz_d,hypremethod,hypre_intpara, hypre_dppara,verbosity, &
Matrix % Hypre, Matrix % Comm)
RETURN
END IF
CALL CreateHYPREAMS( Matrix % NumberOfRows, Rows, Cols, Vals, &
nnd,GM % Rows,GM % Cols,GM % Values,Aperm,Aperm,Aperm,Owner, &
Bperm,NodeOwner,Xvec,RHSvec,ParEnv % myPE, ILUn, Rounds,TOL, &
Expand Down Expand Up @@ -2309,8 +2312,12 @@ END SUBROUTINE UpdateHypre
CONTAINS

SUBROUTINE PrepareHypreAMS()

TYPE(Mesh_t), POINTER :: Mesh

Mesh => Solver % Mesh

nnd = Solver % Mesh % NumberOfNodes
nnd = Mesh % NumberOfNodes
ALLOCATE( NodeOwner(nnd), NodePerm(nnd), Bperm(nnd))

DO i=1,nnd
Expand All @@ -2319,8 +2326,8 @@ SUBROUTINE PrepareHypreAMS()

IF( Parallel ) THEN
NodeOwner = 0
CALL ContinuousNumbering( Solver % Mesh % ParallelInfo, &
NodePerm, BPerm, NodeOwner, nnd, Solver % Mesh)
CALL ContinuousNumbering( Mesh % ParallelInfo, &
NodePerm, BPerm, NodeOwner, nnd, Mesh)
bPerm = bPerm-1
ELSE
NodeOwner = 1
Expand All @@ -2332,24 +2339,24 @@ SUBROUTINE PrepareHypreAMS()
GM => AllocateMatrix()
GM % FORMAT = MATRIX_LIST

DO i=Solver % Mesh % NumberofEdges,1,-1
ind = Solver % Mesh % Edges(i) % NodeIndexes
IF (Solver % Mesh % ParallelInfo % GlobalDOFs(ind(1))> &
Solver % Mesh % ParallelInfo % GlobalDOFs(ind(2))) THEN
DO i=Mesh % NumberofEdges,1,-1
ind = Mesh % Edges(i) % NodeIndexes
IF (Mesh % ParallelInfo % GlobalDOFs(ind(1))> &
Mesh % ParallelInfo % GlobalDOFs(ind(2))) THEN
k=ind(1); ind(1)=ind(2);ind(2)=k
END IF
CALL List_AddToMatrixElement(gm%listmatrix,i,ind(1),-1._dp)
CALL List_AddToMatrixElement(gm%listmatrix,i,ind(2), 1._dp)
CALL List_AddToMatrixElement(gm % listmatrix,i,ind(1),-1._dp)
CALL List_AddToMatrixElement(gm % listmatrix,i,ind(2), 1._dp)
END DO
CALL List_tocrsMatrix(gm)

nnd = Solver % Mesh % NumberOfEdges
nnd = Mesh % NumberOfEdges
ALLOCATE(xx_d(nnd),yy_d(nnd),zz_d(nnd) )
CALL CRS_MatrixVectorMultiply(gm,Solver % Mesh % Nodes % x,xx_d)
CALL CRS_MatrixVectorMultiply(gm,Solver % Mesh % Nodes % y,yy_d)
CALL CRS_MatrixVectorMultiply(gm,Solver % Mesh % Nodes % z,zz_d)
CALL CRS_MatrixVectorMultiply(gm,Mesh % Nodes % x,xx_d)
CALL CRS_MatrixVectorMultiply(gm,Mesh % Nodes % y,yy_d)
CALL CRS_MatrixVectorMultiply(gm,Mesh % Nodes % z,zz_d)

nnd = Solver % Mesh % NumberOfNodes
nnd = Mesh % NumberOfNodes
END SUBROUTINE PrepareHypreAMS


Expand Down
Loading

0 comments on commit d903b4d

Please sign in to comment.