diff --git a/fem/src/SolverUtils.F90 b/fem/src/SolverUtils.F90 index 62be482eaa..532a269355 100644 --- a/fem/src/SolverUtils.F90 +++ b/fem/src/SolverUtils.F90 @@ -15913,14 +15913,25 @@ END SUBROUTINE BlockSolveExt Nmode = 0 20 CALL ConstraintModesDriver( A, x, b, Solver, .TRUE., Nmode, LinModes, FirstLoop = FirstLoop ) - ! Here activate constraint solve only if constraints are not treated as blocks - IF( BlockMode .AND. RestrictionMode ) THEN - CALL Warn(Caller,'Matrix is restricted and block matrix, giving precedence to block nature!') - END IF IF( BlockMode ) THEN CALL Info(Caller,'Solving linear system with block strategy',Level=10) - CALL BlockSolveExt( A, x, bb, Solver ) + ! Here activate constraint solve only if constraints are not treated as blocks + IF( RestrictionMode ) THEN + BLOCK + TYPE(Matrix_t), POINTER :: Acoll + Acoll => AllocateMatrix() + Acoll % FORMAT = MATRIX_LIST + CALL Warn(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 FreeMatrix(Acoll) + Acoll => NULL() + END BLOCK + ELSE + CALL BlockSolveExt( A, x, bb, Solver ) + END IF ELSE IF ( RestrictionMode ) THEN CALL Info(Caller,'Solving linear system with linear restrictions!',Level=10) IF( ListGetLogical( Params,'Save Constraint Matrix',Found ) ) THEN @@ -18515,6 +18526,342 @@ END SUBROUTINE ChangeToHarmonicSystem !------------------------------------------------------------------------------ + +SUBROUTINE EliminateLinearRestriction( StiffMatrix, ForceVector, RestMatrix, & + CollectionMatrix, Solver, CopyStiffMatrix ) + TYPE(Matrix_t) :: StiffMatrix + REAL(KIND=dp) :: ForceVector(:) + TYPE(Matrix_t), POINTER :: RestMatrix + TYPE(Matrix_t) :: CollectionMatrix + TYPE(Solver_t) :: Solver + LOGICAL, OPTIONAL :: CopyStiffMatrix + + INTEGER :: m,n,i,j,k,l,ix,p,q,Loop + INTEGER, ALLOCATABLE, TARGET :: SlavePerm(:),MasterPerm(:),SlaveIPerm(:),MasterIPerm(:) + REAL(KIND=dp), ALLOCATABLE, TARGET :: SlaveDiag(:), MasterDiag(:), DiagDiag(:) + INTEGER, POINTER :: UsePerm(:), UseIPerm(:) + REAL(KIND=dp), POINTER :: UseDiag(:) + REAL(KIND=dp) :: scl + REAL(KIND=dp), POINTER :: TVals(:), Vals(:) + REAL(KIND=dp), POINTER :: CollectionVector(:), RestVector(:) + TYPE(ListMatrix_t), POINTER :: Lmat(:) + TYPE(Matrix_t), POINTER :: Xmat, Tmat + TYPE(ListMatrixEntry_t), POINTER :: cTmp + LOGICAL :: Found, EliminateSlave, EliminateFromMaster, UseTranspose + TYPE(ValueList_t), POINTER :: Params + CHARACTER(*), PARAMETER :: Caller = 'EliminateLinearRestriction' + + + CALL Info(Caller,'Eliminating Constraints from CollectionMatrix',Level=12) + + Params => Solver % Values + + + EliminateSlave = ListGetLogical( Params, 'Eliminate Slave',Found ) + EliminateFromMaster = ListGetLogical( Params, 'Eliminate From Master',Found ) + + UseTranspose = ListGetLogical(Params, 'Use Transpose values', Found) + IF( UseTranspose ) THEN + CALL Info(Caller,'Using transpose values in elimination',Level=15) + END IF + + + n = StiffMatrix % NumberOfRows + m = RestMatrix % NumberOfRows + + RestVector => NULL() + IF(ASSOCIATED(RestMatrix)) RestVector => RestMatrix % RHS + + CollectionVector => CollectionMatrix % RHS + + ! We may optionally ask that the stiffness matrix is copied to the base. + IF( PRESENT(CopyStiffMatrix)) THEN + IF(CopyStiffMatrix) THEN + DO i=StiffMatrix % NumberOfRows,1,-1 + DO j=StiffMatrix % Rows(i+1)-1,StiffMatrix % Rows(i),-1 + CALL AddToMatrixElement( CollectionMatrix, & + i, StiffMatrix % Cols(j), StiffMatrix % Values(j) ) + END DO + CollectionVector(i) = CollectionVector(i) + ForceVector(i) + END DO + END IF + END IF + + + ALLOCATE(SlaveDiag(m),MasterDiag(m),SlavePerm(n),MasterPerm(n),& + SlaveIPerm(m),MasterIPerm(m),DiagDiag(m)) + SlavePerm = 0; SlaveIPerm = 0; + MasterPerm = 0; MasterIPerm = 0 + SlaveDiag = 0.0_dp; MasterDiag = 0.0_dp + DiagDiag = 0.0_dp + + Tvals => RestMatrix % TValues + IF (.NOT.ASSOCIATED(Tvals)) Tvals => RestMatrix % Values + + ! Extract diagonal entries for constraints: + !------------------------------------------ + CALL Info(Caller,'Extracting diagonal entries for constraints',Level=15) + + DO i=1, RestMatrix % NumberOfRows + m = RestMatrix % InvPerm(i) + + IF( m == 0 ) THEN + PRINT *,'InvPerm is zero:',ParEnv % MyPe, i + CYCLE + END IF + + m = MOD(m-1,n) + 1 + SlavePerm(m) = i + SlaveIperm(i) = m + + DO j=RestMatrix % Rows(i), RestMatrix % Rows(i+1)-1 + k = RestMatrix % Cols(j) + IF(k>n) THEN + DiagDiag(i) = Tvals(j) + CYCLE + END IF + + IF( ABS( TVals(j) ) < TINY( 1.0_dp ) ) THEN + PRINT *,'Tvals too small',ParEnv % MyPe,j,i,k,RestMatrix % InvPerm(i),Tvals(j) + END IF + + IF(k == RestMatrix % InvPerm(i)) THEN + SlaveDiag(i) = Tvals(j) + ELSE + MasterDiag(i) = Tvals(j) + MasterPerm(k) = i + MasterIperm(i) = k + END IF + END DO + END DO + + IF(InfoActive(25)) THEN + PRINT *,'SlaveSum:',SUM(SlaveDiag) + PRINT *,'MasterSum:',SUM(MasterDiag) + PRINT *,'SlaveSum abs:',SUM(ABS(SlaveDiag)) + PRINT *,'MasterSum abs:',SUM(ABS(MasterDiag)) + END IF + + IF(EliminateFromMaster) THEN + CALL Info(Caller,'Eliminating from master',Level=15) + UsePerm => MasterPerm + UseDiag => MasterDiag + UseIPerm => MasterIPerm + ELSE + CALL Info(Caller,'Eliminating from slave',Level=15) + UsePerm => SlavePerm + UseDiag => SlaveDiag + UseIPerm => SlaveIPerm + END IF + + IF(UseTranspose) THEN + Vals => Tvals + ELSE + Vals => RestMatrix % Values + END IF + + + ! Replace elimination equations by the constraints (could done be as a postprocessing + ! step, if eq's totally eliminated from linsys.) + ! ---------------------------------------------------------------------------------- + CALL Info(Caller,'Deleting rows from equation to be eliminated',Level=15) + + Lmat => CollectionMatrix % ListMatrix + DO m=1,RestMatrix % NumberOfRows + i = UseIPerm(m) + CALL List_DeleteRow(Lmat, i, Keep=.TRUE.) + END DO + + CALL Info(Caller,'Copying rows from constraint matrix to eliminate dofs',Level=15) + DO m=1,RestMatrix % NumberOfRows + i = UseIPerm(m) + DO l=RestMatrix % Rows(m+1)-1, RestMatrix % Rows(m), -1 + j = RestMatrix % Cols(l) + + ! skip l-coefficient entries, handled separately afterwards: + ! -------------------------------------------------------- + IF(j > n) CYCLE + CALL List_AddToMatrixElement( Lmat, i, j, Vals(l) ) + END DO + CollectionVector(i) = RestVector(m) + END DO + + ! Eliminate slave dof cycles: + ! --------------------------- + Xmat => RestMatrix + Found = .TRUE. + Loop = 0 + DO WHILE(Found) + DO i=Xmat % NumberofRows,1,-1 + q = 0 + DO j = Xmat % Rows(i+1)-1, Xmat % Rows(i),-1 + k = Xmat % Cols(j) + IF(k>n) CYCLE + IF(UsePerm(k)>0 .AND. ABS(TVals(j))>AEPS) q=q+1 + END DO + IF(q>1) EXIT + END DO + Found = (q>1) + + Tmat => Xmat + IF(Found) THEN + Loop = Loop + 1 + CALL Info(Caller,'Recursive elimination round: '//I2S(Loop),Level=15) + + Tmat => AllocateMatrix() + Tmat % Format = MATRIX_LIST + + DO i=Xmat % NumberofRows,1,-1 + DO j = Xmat % Rows(i+1)-1, Xmat % Rows(i),-1 + k = Xmat % Cols(j) + IF ( ABS(Tvals(j))>AEPS ) & + CALL List_AddToMatrixElement(Tmat % ListMatrix, i, k, TVals(j)) + END DO + END DO + + DO m=1,Xmat % NumberOfRows + i = UseIPerm(m) + DO j=Xmat % Rows(m), Xmat % Rows(m+1)-1 + k = Xmat % Cols(j) + + ! The size of SlavePerm is often exceeded but I don't really undersrtand the operation... + ! so this is just a dirty fix. + IF( k > SIZE( SlavePerm ) ) CYCLE + + l = SlavePerm(k) + + IF(l>0 .AND. k/=i) THEN + IF(ABS(Tvals(j)) Tmat % Values + IF(.NOT.ASSOCIATED(Xmat,RestMatrix)) CALL FreeMatrix(Xmat) + END IF + Xmat => TMat + END DO + + ! Eliminate Lagrange Coefficients: + ! -------------------------------- + CALL Info(Caller,'Eliminating Largrange Coefficients',Level=15) + + DO m=1,Tmat % NumberOfRows + i = UseIPerm(m) + IF( ABS( UseDiag(m) ) < TINY( 1.0_dp ) ) THEN + PRINT *,'UseDiag too small:',m,ParEnv % MyPe,UseDiag(m) + CYCLE + END IF + + DO j=TMat % Rows(m), TMat % Rows(m+1)-1 + k = TMat % Cols(j) + IF(k<=n) THEN + IF(UsePerm(k)/=0) CYCLE + scl = -Tvals(j) / UseDiag(m) + ELSE + k = UseIPerm(k-n) + scl = -Tvals(j) / UseDiag(m) + END IF + + DO l=StiffMatrix % Rows(i+1)-1, StiffMatrix % Rows(i),-1 + CALL List_AddToMatrixElement( Lmat, k, & + StiffMatrix % Cols(l), scl * StiffMatrix % Values(l) ) + END DO + CollectionVector(k) = CollectionVector(k) + scl * ForceVector(i) + END DO + END DO + + IF ( .NOT.ASSOCIATED(Tmat, RestMatrix ) ) CALL FreeMatrix(Tmat) + + ! Eliminate slave dofs, using the constraint equations: + ! ----------------------------------------------------- + IF ( EliminateSlave ) THEN + CALL Info(Caller,'Eliminate slave dofs using constraint equations',Level=15) + + CALL List_ToCRSMatrix(CollectionMatrix) + Tmat => AllocateMatrix() + Tmat % Format = MATRIX_LIST + + DO i=1,StiffMatrix % NumberOfRows + IF(UsePerm(i)/=0) CYCLE + + DO m = CollectionMatrix % Rows(i), CollectionMatrix % Rows(i+1)-1 + j = SlavePerm(CollectionMatrix % Cols(m)) + + IF(j==0) THEN + CYCLE + END IF + IF( ABS( SlaveDiag(j) ) < TINY( 1.0_dp ) ) THEN + PRINT *,'SlaveDiag too small:',j,ParEnv % MyPe,SlaveDiag(j) + CYCLE + END IF + + scl = -CollectionMatrix % Values(m) / SlaveDiag(j) + CollectionMatrix % Values(m) = 0._dp + + ! ... and add replacement values: + ! ------------------------------- + k = UseIPerm(j) + DO p=CollectionMatrix % Rows(k+1)-1, CollectionMatrix % Rows(k), -1 + l = CollectionMatrix % Cols(p) + IF ( l /= SlaveIPerm(j) ) & + CALL List_AddToMatrixElement( Tmat % listmatrix, i, l, scl*CollectionMatrix % Values(p) ) + END DO + CollectionVector(i) = CollectionVector(i) + scl * CollectionVector(k) + END DO + END DO + + CALL List_ToListMatrix(CollectionMatrix) + Lmat => CollectionMatrix % ListMatrix + + CALL List_ToCRSMatrix(Tmat) + DO i=TMat % NumberOfRows,1,-1 + DO j=TMat % Rows(i+1)-1,TMat % Rows(i),-1 + CALL List_AddToMatrixElement( Lmat, i, TMat % cols(j), TMat % Values(j) ) + END DO + END DO + CALL FreeMatrix(Tmat) + END IF + + ! Optimize bandwidth, if needed: + ! ------------------------------ + IF(EliminateFromMaster) THEN + CALL Info(Caller,'Optimizing bandwidth after elimination',Level=15) + DO i=1,RestMatrix % NumberOfRows + j = SlaveIPerm(i) + k = MasterIPerm(i) + + Ctmp => Lmat(j) % Head + Lmat(j) % Head => Lmat(k) % Head + Lmat(k) % Head => Ctmp + + l = Lmat(j) % Degree + Lmat(j) % Degree = Lmat(k) % Degree + Lmat(k) % Degree = l + + scl = CollectionVector(j) + CollectionVector(j) = CollectionVector(k) + CollectionVector(k) = scl + END DO + END IF + + CALL Info(Caller,'Finished Eliminating Restrictions',Level=12) + +END SUBROUTINE EliminateLinearRestriction + + + !------------------------------------------------------------------------------ !> This subroutine will solve the system with some linear restriction. !> The restriction matrix is assumed to be in the ConstraintMatrix-field of @@ -18926,6 +19273,11 @@ RECURSIVE SUBROUTINE SolveWithLinearRestriction( StiffMatrix, ForceVector, & ! necessarily biorthogonal constraint equation test functions. !------------------------------------------------------------------------------ IF (ASSOCIATED(RestMatrix) .AND. EliminateConstraints) THEN + + IF( .NOT. ListGetLogical( Params,'Old Elimination', Found ) ) THEN + CALL EliminateLinearRestriction( StiffMatrix, ForceVector, RestMatrix, CollectionMatrix, Solver) + + ELSE CALL Info(Caller,'Eliminating Constraints from CollectionMatrix',Level=12) n = StiffMatrix % NumberOfRows @@ -18942,6 +19294,7 @@ RECURSIVE SUBROUTINE SolveWithLinearRestriction( StiffMatrix, ForceVector, & ! Extract diagonal entries for constraints: !------------------------------------------ CALL Info(Caller,'Extracting diagonal entries for constraints',Level=15) + DO i=1, RestMatrix % NumberOfRows m = RestMatrix % InvPerm(i) @@ -19202,6 +19555,8 @@ RECURSIVE SUBROUTINE SolveWithLinearRestriction( StiffMatrix, ForceVector, & END IF CALL Info(Caller,'Finished Adding ConstraintMatrix',Level=12) + END IF + END IF CALL Info(Caller,'Reverting CollectionMatrix back to CRS matrix',Level=10) diff --git a/fem/tests/DisContBoundaryMortarContElim/discont.sif b/fem/tests/DisContBoundaryMortarContElim/discont.sif index 5bbdb6a71b..eaee247184 100644 --- a/fem/tests/DisContBoundaryMortarContElim/discont.sif +++ b/fem/tests/DisContBoundaryMortarContElim/discont.sif @@ -17,7 +17,7 @@ Header End Simulation - Max Output Level = 5 + Max Output Level = 25 Coordinate System = "Cartesian" Coordinate Mapping(3) = 1 2 3 @@ -93,6 +93,9 @@ Solver 1 Nonlinear System Consistent Norm = True Linear System Abort Not Converged = False + + +! Old Elimination = Logical True End