Skip to content

Commit

Permalink
Fix to previous
Browse files Browse the repository at this point in the history
  • Loading branch information
raback committed Nov 26, 2024
1 parent 06c3689 commit aa67912
Showing 1 changed file with 26 additions and 29 deletions.
55 changes: 26 additions & 29 deletions fem/src/SolverUtils.F90
Original file line number Diff line number Diff line change
Expand Up @@ -15917,15 +15917,17 @@ END SUBROUTINE BlockSolveExt
IF( BlockMode ) THEN
CALL Info(Caller,'Solving linear system with block strategy',Level=10)
! Here activate constraint solve only if constraints are not treated as blocks
IF( RestrictionMode ) THEN
IF( RestrictionMode .AND. &
ListGetLogical( Params, 'Eliminate Linear Constraints', Found) ) THEN
BLOCK
TYPE(Matrix_t), POINTER :: Acoll
Acoll => AllocateMatrix()
Acoll % FORMAT = MATRIX_LIST
CALL Warn(Caller,'Eliminating constraints before going into block matrix!')
CALL Info(Caller,'Eliminating constraints before going into block matrix!')
CALL EliminateLinearRestriction( A, bb, A % ConstraintMatrix, Acoll, Solver, .TRUE. )
CALL List_ToCRSMatrix(Acoll)
CALL BlockSolveExt( Acoll, x, Acoll % rhs, Solver )
CALL Info(Caller,'Freeing collection matrix after solution',Level=10)
CALL FreeMatrix(Acoll)
Acoll => NULL()
END BLOCK
Expand Down Expand Up @@ -18527,8 +18529,14 @@ END SUBROUTINE ChangeToHarmonicSystem



!------------------------------------------------------------------------------
!> Eliminate linear restriction only using the ListMatrix structure.
!> This is placed in a separate routine such that it can be called
!> when solving with and without block matrix being active.
!------------------------------------------------------------------------------
SUBROUTINE EliminateLinearRestriction( StiffMatrix, ForceVector, RestMatrix, &
CollectionMatrix, Solver, CopyStiffMatrix )
IMPLICIT NONE
TYPE(Matrix_t) :: StiffMatrix
REAL(KIND=dp) :: ForceVector(:)
TYPE(Matrix_t), POINTER :: RestMatrix
Expand Down Expand Up @@ -18869,49 +18877,39 @@ END SUBROUTINE EliminateLinearRestriction
!> ConstraintMatrix.
!------------------------------------------------------------------------------
RECURSIVE SUBROUTINE SolveWithLinearRestriction( StiffMatrix, ForceVector, &
Solution, Norm, DOFs, Solver )
Solution, Norm, DOFs, Solver )
!------------------------------------------------------------------------------
IMPLICIT NONE
TYPE(Matrix_t), POINTER :: StiffMatrix !< Linear equation matrix information.
!< The restriction matrix is assumed to be in the EMatrix-field
REAL(KIND=dp),TARGET :: ForceVector(:) !< The right hand side of the linear equation
REAL(KIND=dp),TARGET :: Solution(:) !< Previous solution as input, new solution as output.
REAL(KIND=dp) :: Norm !< The L2 norm of the solution.
INTEGER :: DOFs !< Number of degrees of freedom of the equation.
TYPE(Solver_t), TARGET :: Solver !< Linear equation solver options.
!------------------------------------------------------------------------------
TYPE(Solver_t), POINTER :: SolverPointer
TYPE(Matrix_t), POINTER :: CollectionMatrix, RestMatrix, AddMatrix, &
RestMatrixTranspose, TMat, XMat
REAL(KIND=dp), POINTER CONTIG :: CollectionVector(:), RestVector(:),&
AddVector(:), Tvals(:), Vals(:)
TYPE(Matrix_t), POINTER :: CollectionMatrix, RestMatrix, AddMatrix, RestMatrixTranspose
REAL(KIND=dp), POINTER CONTIG :: CollectionVector(:), RestVector(:), AddVector(:)
REAL(KIND=dp), POINTER :: MultiplierValues(:), pSol(:),DiagScaling(:)
REAL(KIND=dp), ALLOCATABLE, TARGET :: CollectionSolution(:), TotValues(:)
REAL(KIND=dp), ALLOCATABLE, TARGET :: CollectionSolution(:)
INTEGER :: NumberOfRows, NumberOfValues, MultiplierDOFs, istat, NoEmptyRows
INTEGER :: i, j, k, l, m, n, p,q, ix, Loop, colj, nIter
TYPE(Variable_t), POINTER :: MultVar, iterV
REAL(KIND=dp) :: scl, rowsum, Relax, val
LOGICAL :: Found, ExportMultiplier, NotExplicit, Refactorize, EnforceDirichlet, &
NonEmptyRow, ComplexSystem, ConstraintScaling, UseTranspose, EliminateConstraints, &
SkipConstraints, ResidualMode
SAVE MultiplierValues, SolverPointer

TYPE(ListMatrix_t), POINTER :: cList
TYPE(ListMatrixEntry_t), POINTER :: cPtr, cPrev, cTmp

INTEGER, ALLOCATABLE, TARGET :: SlavePerm(:), SlaveIPerm(:), MasterPerm(:), MasterIPerm(:)
INTEGER, POINTER :: UsePerm(:), UseIPerm(:)
REAL(KIND=dp), POINTER :: UseDiag(:), svals(:)
TYPE(ListMatrix_t), POINTER :: Lmat(:)
LOGICAL :: EliminateFromMaster, EliminateSlave, Parallel, UseTreeGauge, &
NeedMassDampValues, DoOwnScaling
REAL(KIND=dp), ALLOCATABLE, TARGET :: SlaveDiag(:), MasterDiag(:), DiagDiag(:)
NonEmptyRow, ComplexSystem, ConstraintScaling, UseTranspose, EliminateConstraints, &
SkipConstraints, ResidualMode
INTEGER, POINTER :: UseIPerm(:)
REAL(KIND=dp), POINTER :: UseDiag(:)
LOGICAL :: Parallel, UseTreeGauge, NeedMassDampValues, DoOwnScaling
LOGICAL, ALLOCATABLE :: TrueDof(:)
INTEGER, ALLOCATABLE :: Iperm(:)
REAL(KIND=dp) :: t0,rt0,st,rst
CHARACTER(:), ALLOCATABLE :: str,MultiplierName
TYPE(ValueList_t), POINTER :: Params
CHARACTER(*), PARAMETER :: Caller = 'SolveWithLinearRestriction'

SAVE MultiplierValues, SolverPointer

!------------------------------------------------------------------------------
CALL Info( Caller, ' ', Level=12 )
Expand Down Expand Up @@ -19274,10 +19272,9 @@ RECURSIVE SUBROUTINE SolveWithLinearRestriction( StiffMatrix, ForceVector, &
!------------------------------------------------------------------------------
IF (ASSOCIATED(RestMatrix) .AND. EliminateConstraints) THEN

IF( .NOT. ListGetLogical( Params,'Old Elimination', Found ) ) THEN
CALL EliminateLinearRestriction( StiffMatrix, ForceVector, RestMatrix, CollectionMatrix, Solver)

ELSE
#if 1
CALL EliminateLinearRestriction( StiffMatrix, ForceVector, RestMatrix, CollectionMatrix, Solver)
#else
CALL Info(Caller,'Eliminating Constraints from CollectionMatrix',Level=12)

n = StiffMatrix % NumberOfRows
Expand Down Expand Up @@ -19555,9 +19552,9 @@ RECURSIVE SUBROUTINE SolveWithLinearRestriction( StiffMatrix, ForceVector, &
END IF

CALL Info(Caller,'Finished Adding ConstraintMatrix',Level=12)
END IF

#endif
END IF


CALL Info(Caller,'Reverting CollectionMatrix back to CRS matrix',Level=10)
IF(CollectionMatrix % FORMAT==MATRIX_LIST) THEN
Expand Down

0 comments on commit aa67912

Please sign in to comment.