diff --git a/fem/src/MeshUtils.F90 b/fem/src/MeshUtils.F90 index 07d11b3309..87fda1121e 100644 --- a/fem/src/MeshUtils.F90 +++ b/fem/src/MeshUtils.F90 @@ -2970,7 +2970,10 @@ SUBROUTINE PrepareMesh( Model, Mesh, Parallel, Def_Dofs, mySolver ) CALL NonNodalElements() IF( Parallel ) THEN + CALL Info(Caller,'Generating parallel communications for the non-nodal mesh',Level=20) + CALL ResetTimer('ParallelNonNodal') CALL ParallelNonNodalElements() + CALL CheckTimer('ParallelNonNodal',Level=7,Delete=.TRUE.) END IF CALL EnlargeCoordinates( Mesh ) diff --git a/fem/src/ModelDescription.F90 b/fem/src/ModelDescription.F90 index c798893f78..13a183f164 100644 --- a/fem/src/ModelDescription.F90 +++ b/fem/src/ModelDescription.F90 @@ -553,6 +553,18 @@ RECURSIVE SUBROUTINE LoadInputFile( Model, InFileUnit, FileName, & ELSE IF ( Section == 'run' ) THEN IF ( PRESENT(runc) ) runc=.TRUE. EXIT + ELSE IF ( Section == 'stop' ) THEN + CALL Warn(Caller,'Encountered "STOP" in sif, rest will be ignored!') + EXIT + ELSE IF( Section == '/*' ) THEN + CALL Info(Caller,'Starting comment section!') + DO WHILE( ReadAndTrim( InFileUnit, Section, Echo ) ) + IF ( Section == '*/' ) THEN + CALL Info(Caller,'Finished comment section!') + EXIT + END IF + END DO + CYCLE END IF FreeNames = ( CheckAbort <= 0 ) diff --git a/fem/src/SParIterSolver.F90 b/fem/src/SParIterSolver.F90 index bff3a848db..ba7f82c03a 100644 --- a/fem/src/SParIterSolver.F90 +++ b/fem/src/SParIterSolver.F90 @@ -2055,6 +2055,10 @@ END SUBROUTINE SolveTrilinos4 ! !------------------------------------------------------------------ + CALL Info(Caller,'Copying Matrix values into SplittedMatrix',Level=20) + CALL ResetTimer('SplittedMatrix') + + GT => SplittedMatrix % GlueTable DO i = 1, SourceMatrix % NumberOfRows GRow = ParallelInfo % GlobalDOFs(i) @@ -2154,12 +2158,14 @@ END SUBROUTINE SolveTrilinos4 END DO CALL GlueFinalize( SourceMatrix, SplittedMatrix, ParallelInfo ) + CALL CheckTimer('SplittedMatrix',Level=7,Delete=.TRUE.) + !------------------------------------------------------------------ - ! ! Call the actual solver routine (based on older design) - ! !------------------------------------------------------------------ + CALL Info(Caller,'Going into actual parallel solution',Level=20) + CALL Solve( SourceMatrix, SParMatrixDesc % SplittedMatrix, & ParallelInfo, RHSVec, XVec, Solver, Errinfo ) @@ -2436,9 +2442,7 @@ SUBROUTINE Solve( SourceMatrix, SplittedMatrix, ParallelInfo, & PIGpntr => GlobalData !---------------------------------------------------------------------- - ! ! Initialize Right-Hand-Side - ! !---------------------------------------------------------------------- ALLOCATE(TmpRHSVec(SplittedMatrix % InsideMatrix % NumberOfRows)) TmpRHSVec = 0 @@ -2495,12 +2499,10 @@ SUBROUTINE Solve( SourceMatrix, SplittedMatrix, ParallelInfo, & GlobalMatrix % Ematrix => SourceMatrix GlobalMatrix % COMPLEX = SourceMatrix % COMPLEX + !---------------------------------------------------------------------- - ! ! Set up the preconditioner - ! !---------------------------------------------------------------------- - IF (SplittedMatrix % InsideMatrix % NumberOFRows>0) THEN ! IF (SplittedMatrix % InsideMatrix % Diag(1)==0) THEN DO i = 1, SplittedMatrix % InsideMatrix % NumberOfRows @@ -2533,11 +2535,8 @@ SUBROUTINE Solve( SourceMatrix, SplittedMatrix, ParallelInfo, & !---------------------------------------------------------------------- - ! ! Call the main iterator routine - ! !---------------------------------------------------------------------- - CM => SourceMatrix % ConstraintMatrix IF (ASSOCIATED(CM)) THEN ALLOCATE(SPerm(SourceMatrix % NumberOfRows)); SPerm=0 @@ -2576,9 +2575,7 @@ SUBROUTINE Solve( SourceMatrix, SplittedMatrix, ParallelInfo, & GlobalMatrix => SaveMatrix !---------------------------------------------------------------------- - ! ! Collect the result - ! !---------------------------------------------------------------------- ALLOCATE( VecEPerNB( ParEnv % PEs ) ) VecEPerNB = 0 @@ -2609,12 +2606,10 @@ SUBROUTINE Solve( SourceMatrix, SplittedMatrix, ParallelInfo, & END DO CALL ExchangeResult( SourceMatrix,SplittedMatrix,ParallelInfo,XVec ) + !---------------------------------------------------------------------- - ! ! Clean the work space - ! !---------------------------------------------------------------------- - DEALLOCATE( TmpXVec, TmpRHSVec, VecEPerNB ) !---------------------------------------------------------------------- END SUBROUTINE Solve diff --git a/fem/src/SolverUtils.F90 b/fem/src/SolverUtils.F90 index c55b0fcd73..78865f901e 100644 --- a/fem/src/SolverUtils.F90 +++ b/fem/src/SolverUtils.F90 @@ -10703,6 +10703,12 @@ SUBROUTINE ComputeChange(Solver,SteadyState,nsize,values,values0,Matrix,RHS) Parallel = Solver % Parallel + IF(PRESENT(Matrix)) THEN + A => Matrix + ELSE + A => Solver % Matrix + END IF + IF(SteadyState) THEN Skip = ListGetLogical( SolverParams,'Skip Compute Steady State Change',Stat) IF( Skip ) THEN @@ -10768,7 +10774,7 @@ SUBROUTINE ComputeChange(Solver,SteadyState,nsize,values,values0,Matrix,RHS) END IF ResidualMode = ListGetLogical( SolverParams,'Linear System Residual Mode',Stat) - + ConvergenceType = ListGetString(SolverParams,& 'Nonlinear System Convergence Measure',Stat) IF(.NOT. stat) ConvergenceType = 'norm' @@ -10826,7 +10832,7 @@ SUBROUTINE ComputeChange(Solver,SteadyState,nsize,values,values0,Matrix,RHS) END IF IF( SkipConstraints ) THEN - n = MIN( n, Solver % Matrix % NumberOfRows ) + n = MIN( n, A % NumberOfRows ) END IF ! If requested (for p-elements) only use the dofs associated to nodes. @@ -10926,16 +10932,10 @@ SUBROUTINE ComputeChange(Solver,SteadyState,nsize,values,values0,Matrix,RHS) ! x is solution of A(x0)x=b(x0) thus residual should really be r=b(x)-A(x)x ! Instead we use r=b(x0)-A(x0)x0 which unfortunately is one step behind. !-------------------------------------------------------------------------- - IF(PRESENT(Matrix)) THEN - A => Matrix - ELSE - A => Solver % Matrix - END IF - IF(PRESENT(RHS)) THEN b => RHS ELSE - b => Solver % Matrix % rhs + b => A % rhs END IF ALLOCATE(r(n)) @@ -10982,8 +10982,7 @@ SUBROUTINE ComputeChange(Solver,SteadyState,nsize,values,values0,Matrix,RHS) ! Here the true linear system residual r=b(x0)-A(x0)x is computed. ! This option is useful for certain special solvers. !-------------------------------------------------------------------------- - A => Solver % Matrix - b => Solver % Matrix % rhs + b => A % rhs IF (Parallel) THEN @@ -12982,24 +12981,21 @@ SUBROUTINE ScaleLinearSystemDiagonal() REAL(KIND=dp), POINTER :: Diag(:) TYPE(Matrix_t), POINTER :: CM - IF( Parallel ) THEN - CALL Info('ScaleLinearSystem','Scaling diagonal entries to unity in parallel',Level=10) - ELSE - CALL Info('ScaleLinearSystem','Scaling diagonal entries to unity in serial',Level=10) - END IF - + A % ScalingMethod = 1 IF( PRESENT( DiagScaling ) ) THEN CALL Info('ScaleLinearSystem','Reusing existing > DiagScaling < vector',Level=12) Diag => DiagScaling ELSE - CALL Info('ScaleLinearSystem','Computing > DiagScaling < vector',Level=12) IF(.NOT. ASSOCIATED(A % DiagScaling)) THEN + CALL Info('ScaleLinearSystem','Creating > DiagScaling < vector of size '//I2S(n),Level=10) ALLOCATE( A % DiagScaling(n) ) + ELSE + CALL Info('ScaleLinearSystem','Recomputing > DiagScaling < vector of size '//I2S(n),Level=12) END IF Diag => A % DiagScaling - Diag = 0._dp + Diag(1:n) = 0._dp IF ( ComplexMatrix ) THEN CALL Info('ScaleLinearSystem','Assuming complex matrix while scaling',Level=20) @@ -13033,7 +13029,10 @@ SUBROUTINE ScaleLinearSystemDiagonal() !$OMP END PARALLEL DO END IF - IF ( Parallel ) CALL ParallelSumVector(A, Diag) + IF ( Parallel ) THEN + CALL Info('ScaleLinearSystem','Performing parallel summation of > DiagScaling < vector',Level=20) + CALL ParallelSumVector(A, Diag) + END IF IF ( ComplexMatrix ) THEN !$OMP PARALLEL DO & @@ -13090,7 +13089,10 @@ SUBROUTINE ScaleLinearSystemDiagonal() ! Optionally we may just create the diag and leave the scaling undone !-------------------------------------------------------------------- IF( PRESENT( ApplyScaling ) ) THEN - IF(.NOT. ApplyScaling ) RETURN + IF(.NOT. ApplyScaling ) THEN + CALL Info('ScaleLinearSystem','Application of scaling skipped!',Level=20) + RETURN + END IF END IF CALL Info('ScaleLinearSystem','Scaling matrix values',Level=20) @@ -13297,7 +13299,7 @@ SUBROUTINE ScaleLinearSystemConstant() DoRHS = .TRUE. IF (PRESENT(RhsScaling)) DoRHS = RhsScaling IF (DoRHS) THEN - bsum = SUM( ABS( b ) ) + bsum = SUM( ABS( b(1:n) ) ) nSum = n IF ( Parallel ) THEN @@ -13635,8 +13637,12 @@ SUBROUTINE BackScaleLinearSystemDiagonal() END IF A % RhsScaling=1._dp - DEALLOCATE(A % DiagScaling); A % DiagScaling=>NULL() + IF(.NOT. PRESENT(DiagScaling) ) THEN + DEALLOCATE(A % DiagScaling) + A % DiagScaling=>NULL() + END IF + END SUBROUTINE BackScaleLinearSystemDiagonal @@ -14434,7 +14440,7 @@ RECURSIVE SUBROUTINE SolveLinearSystem( A, b, & CHARACTER(:), ALLOCATABLE :: Method, Prec, SaveSlot INTEGER(KIND=AddrInt) :: Proc REAL(KIND=dp), ALLOCATABLE, TARGET :: Px(:), & - TempVector(:), TempRHS(:), NonlinVals(:) + TempRHS(:), NonlinVals(:) REAL(KIND=dp), POINTER :: Diag(:) REAL(KIND=dp) :: s,Relaxation,Beta,Gamma,bnorm,Energy,xn,bn TYPE(ValueList_t), POINTER :: Params @@ -14766,45 +14772,8 @@ END SUBROUTINE BlockSolveExt IF( NoSolve ) GOTO 110 END IF - ! Sometimes the r.h.s. may abruptly diminish in value resulting to significant - ! convergence issues or it may be that the system scales linearly with the source. - ! This flag tries to improve on the initial guess of the linear solvers, and may - ! sometimes even result to the exact solution. IF( ListGetLogical( Params,'Linear System Normalize Guess',GotIt ) ) THEN - CALL Info(Caller,'Normalizing initial guess!',Level=30) - - ALLOCATE( TempVector(A % NumberOfRows) ) - - IF ( Parallel ) THEN - IF( .NOT. ALLOCATED( TempRHS ) ) THEN - ALLOCATE( TempRHS(A % NumberOfRows) ); TempRHS=0._dp - END IF - - Tempvector = 0._dp - TempRHS(1:n) = b(1:n) - CALL ParallelInitSolve( A, x, TempRHS, Tempvector ) - - MP => ParallelMatrix(A,mx,mb,mr) - mn = MP % NumberOfRows - - TempVector = 0._dp - CALL ParallelMatrixVector( A, mx, TempVector ) - - bn = ParallelDot( mn, TempVector, mb ) - xn = ParallelDot( mn, TempVector, TempVector ) - DEALLOCATE( TempRHS ) - ELSE - CALL MatrixVectorMultiply( A, x, TempVector ) - xn = SUM( TempVector(1:n)**2 ) - bn = SUM( TempVector(1:n) * b(1:n) ) - END IF - - IF( xn > TINY( xn ) ) THEN - x(1:n) = x(1:n) * ( bn / xn ) - WRITE( Message,'(A,ES12.3)') 'Linear System Normalizing Factor: ',bn/xn - CALL Info(Caller,Message,Level=6) - END IF - DEALLOCATE( TempVector ) + CALL NormalizeInitialGuess() END IF IF( ListGetLogical( Params,'Linear System Nullify Guess',GotIt ) ) THEN @@ -14828,8 +14797,7 @@ END SUBROUTINE BlockSolveExt CALL ListAddNewInteger( Params,'Vanka Mode',i) END IF CALL VankaCreate(A,Solver) - END IF - IF ( Prec=='circuit' ) THEN + ELSE IF ( Prec=='circuit' ) THEN CALL CircuitPrecCreate(A,Solver) END IF CALL CheckTimer("Prec0-"//TRIM(Prec),Level=8,Delete=.TRUE.) @@ -15022,6 +14990,59 @@ END SUBROUTINE BlockSolveExt END IF END IF + + CONTAINS + + + ! Sometimes the r.h.s. may abruptly diminish in value resulting to significant + ! convergence issues or it may be that the system scales linearly with the source. + ! This flag tries to improve on the initial guess of the linear solvers, and may + ! sometimes even result to the exact solution. + !-------------------------------------------------------------------------------- + SUBROUTINE NormalizeInitialGuess() + REAL(KIND=dp) :: xn, bn + REAL(KIND=dp), ALLOCATABLE, TARGET :: TempVector(:) + + + CALL Info(Caller,'Normalizing initial guess!',Level=30) + + ALLOCATE( TempVector(A % NumberOfRows) ) + + IF ( Parallel ) THEN + IF( .NOT. ALLOCATED( TempRHS ) ) THEN + ALLOCATE( TempRHS(A % NumberOfRows) ); TempRHS=0._dp + END IF + + Tempvector = 0._dp + TempRHS(1:n) = b(1:n) + CALL ParallelInitSolve( A, x, TempRHS, Tempvector ) + + MP => ParallelMatrix(A,mx,mb,mr) + mn = MP % NumberOfRows + + TempVector = 0._dp + CALL ParallelMatrixVector( A, mx, TempVector ) + + bn = ParallelDot( mn, TempVector, mb ) + xn = ParallelDot( mn, TempVector, TempVector ) + DEALLOCATE( TempRHS ) + ELSE + CALL MatrixVectorMultiply( A, x, TempVector ) + xn = SUM( TempVector(1:n)**2 ) + bn = SUM( TempVector(1:n) * b(1:n) ) + END IF + + IF( xn > TINY( xn ) ) THEN + x(1:n) = x(1:n) * ( bn / xn ) + WRITE( Message,'(A,ES12.3)') 'Linear System Normalizing Factor: ',bn/xn + CALL Info(Caller,Message,Level=6) + END IF + DEALLOCATE( TempVector ) + + END SUBROUTINE NormalizeInitialGuess + + + !------------------------------------------------------------------------------ END SUBROUTINE SolveLinearSystem !------------------------------------------------------------------------------ @@ -15804,8 +15825,10 @@ END SUBROUTINE BlockSolveExt n = A % NumberOfRows - ResidualMode = ListGetLogical( Params,'Linear System Residual Mode',Found ) - + RestrictionMode = HaveRestrictionMatrix( A ) + + ResidualMode = ListGetLogical( Params,'Linear System Residual Mode',Found ) + BlockMode = ListGetLogical( Params,'Linear System Block Mode',Found ) !------------------------------------------------------------------------------ @@ -15813,7 +15836,7 @@ END SUBROUTINE BlockSolveExt ! work properly with the Dirichlet elimination. !------------------------------------------------------------------------------ NeedPrevSol = ResidualMode - + IF(.NOT. NeedPrevSol ) THEN Relaxation = ListGetCReal( Params, & 'Nonlinear System Relaxation Factor', Found ) @@ -15868,8 +15891,6 @@ END SUBROUTINE BlockSolveExt bb => b END IF - RestrictionMode = HaveRestrictionMatrix( A ) - FirstLoop = .TRUE. Nmode = 0 20 CALL ConstraintModesDriver( A, x, b, Solver, .TRUE., Nmode, LinModes, FirstLoop = FirstLoop ) @@ -18481,7 +18502,6 @@ END SUBROUTINE ChangeToHarmonicSystem !> The restriction matrix is assumed to be in the ConstraintMatrix-field of !> the StiffMatrix. The restriction vector is the RHS-field of the !> ConstraintMatrix. -!> NOTE: Only serial solver implemented so far ... !------------------------------------------------------------------------------ RECURSIVE SUBROUTINE SolveWithLinearRestriction( StiffMatrix, ForceVector, & Solution, Norm, DOFs, Solver ) @@ -18500,15 +18520,15 @@ RECURSIVE SUBROUTINE SolveWithLinearRestriction( StiffMatrix, ForceVector, & RestMatrixTranspose, TMat, XMat REAL(KIND=dp), POINTER CONTIG :: CollectionVector(:), RestVector(:),& AddVector(:), Tvals(:), Vals(:) - REAL(KIND=dp), POINTER :: MultiplierValues(:), pSol(:) + REAL(KIND=dp), POINTER :: MultiplierValues(:), pSol(:),DiagScaling(:) REAL(KIND=dp), ALLOCATABLE, TARGET :: CollectionSolution(:), TotValues(:) INTEGER :: NumberOfRows, NumberOfValues, MultiplierDOFs, istat, NoEmptyRows - INTEGER :: i, j, k, l, m, n, p,q, ix, Loop, colj - TYPE(Variable_t), POINTER :: MultVar + 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, EliminateDiscont, & NonEmptyRow, ComplexSystem, ConstraintScaling, UseTranspose, EliminateConstraints, & - SkipConstraints + SkipConstraints, ResidualMode SAVE MultiplierValues, SolverPointer TYPE(ListMatrix_t), POINTER :: cList @@ -18518,29 +18538,35 @@ RECURSIVE SUBROUTINE SolveWithLinearRestriction( StiffMatrix, ForceVector, & INTEGER, POINTER :: UsePerm(:), UseIPerm(:) REAL(KIND=dp), POINTER :: UseDiag(:), svals(:) TYPE(ListMatrix_t), POINTER :: Lmat(:) - LOGICAL :: EliminateFromMaster, EliminateSlave, Parallel, UseTreeGauge, NeedMassDampValues + LOGICAL :: EliminateFromMaster, EliminateSlave, Parallel, UseTreeGauge, & + NeedMassDampValues, DoOwnScaling REAL(KIND=dp), ALLOCATABLE, TARGET :: SlaveDiag(:), MasterDiag(:), DiagDiag(:) 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' - !------------------------------------------------------------------------------ CALL Info( Caller, ' ', Level=12 ) SolverPointer => Solver + Params => Solver % Values t0 = CPUTime() rt0 = RealTime() Parallel = Solver % Parallel - - NotExplicit = ListGetLogical(Solver % Values,'No Explicit Constrained Matrix',Found) + + ResidualMode = ListGetLogical( Params,'Restriction System Residual Mode',Found ) + iterV => VariableGet(Solver % Mesh % Variables,'nonlin iter') + nIter = NINT(iterV % Values(1)) + + NotExplicit = ListGetLogical(Params,'No Explicit Constrained Matrix',Found) IF(.NOT. Found) NotExplicit=.FALSE. - NeedMassDampValues = ListGetLogical( Solver % Values, 'Eigen Analysis', Found ) + NeedMassDampValues = ListGetLogical(Params, 'Eigen Analysis', Found ) RestMatrix => NULL() IF(.NOT.NotExplicit) RestMatrix => StiffMatrix % ConstraintMatrix @@ -18554,9 +18580,11 @@ RECURSIVE SUBROUTINE SolveWithLinearRestriction( StiffMatrix, ForceVector, & NumberOfRows = StiffMatrix % NumberOfRows CollectionMatrix => StiffMatrix % CollectionMatrix - Refactorize = ListGetLogical(Solver % Values,'Linear System Refactorize',Found) - IF(.NOT.Found) Refactorize = .TRUE. - + Refactorize = ListGetLogical(Params,'Linear System Refactorize',Found) + IF(.NOT.Found) THEN + Refactorize = .NOT. ( ResidualMode .AND. nIter > 1) + END IF + IF(ASSOCIATED(CollectionMatrix)) THEN IF(Refactorize.AND..NOT.NotExplicit) THEN CALL Info( Caller,'Freeing previous collection matrix structures',Level=10) @@ -18573,7 +18601,7 @@ RECURSIVE SUBROUTINE SolveWithLinearRestriction( StiffMatrix, ForceVector, & ELSE DEALLOCATE(CollectionMatrix % RHS) CollectionMatrix % Values = 0.0_dp - + IF(NeedMassDampValues) THEN IF(ASSOCIATED(CollectionMatrix % MassValues)) CollectionMatrix % MassValues = 0.0_dp IF(ASSOCIATED(CollectionMatrix % DampValues)) CollectionMatrix % DampValues = 0.0_dp @@ -18583,7 +18611,7 @@ RECURSIVE SUBROUTINE SolveWithLinearRestriction( StiffMatrix, ForceVector, & NumberOfRows = StiffMatrix % NumberOfRows IF(ASSOCIATED(AddMatrix)) NumberOfRows = MAX(NumberOfRows,AddMatrix % NumberOfRows) - EliminateConstraints = ListGetLogical( Solver % Values, 'Eliminate Linear Constraints', Found) + EliminateConstraints = ListGetLogical( Params, 'Eliminate Linear Constraints', Found) IF(ASSOCIATED(RestMatrix)) THEN IF(.NOT.EliminateConstraints) NumberOfRows = NumberOFRows + RestMatrix % NumberOfRows END IF @@ -18597,78 +18625,78 @@ RECURSIVE SUBROUTINE SolveWithLinearRestriction( StiffMatrix, ForceVector, & CollectionSolution = 0.0_dp ComplexSystem = StiffMatrix % COMPLEX .OR. & - ListGetLogical( Solver % Values,'Linear System Complex', Found ) + ListGetLogical(Params,'Linear System Complex', Found ) !------------------------------------------------------------------------------ ! If multiplier should be exported, allocate memory and export the variable. !------------------------------------------------------------------------------ - ExportMultiplier = ListGetLogical( Solver % Values, 'Export Lagrange Multiplier', Found ) + ExportMultiplier = ListGetLogical(Params, 'Export Lagrange Multiplier', Found ) IF ( ExportMultiplier ) THEN - MultiplierName = LagrangeMultiplierName( Solver ) - MultVar => VariableGet(Solver % Mesh % Variables, MultiplierName) - j = 0 - IF(ASSOCIATED(RestMatrix)) j = RestMatrix % NumberofRows - IF(ASSOCIATED(AddMatrix)) j = j+MAX(0,AddMatrix % NumberofRows-StiffMatrix % NumberOfRows) - - IF ( .NOT. ASSOCIATED(MultVar) ) THEN - CALL Info(Caller,'Creating variable for Lagrange multiplier',Level=8) - ALLOCATE( MultiplierValues(j), STAT=istat ) - IF ( istat /= 0 ) CALL Fatal(Caller,'Memory allocation error.') - - MultiplierValues = 0.0_dp - IF( ComplexSystem ) THEN - CALL VariableAddVector(Solver % Mesh % Variables, Solver % Mesh, SolverPointer, & - MultiplierName, 2, MultiplierValues) - ELSE - CALL VariableAdd(Solver % Mesh % Variables, Solver % Mesh, SolverPointer, & - MultiplierName, 1, MultiplierValues) - END IF - MultVar => VariableGet(Solver % Mesh % Variables, MultiplierName) - END IF + MultiplierName = LagrangeMultiplierName( Solver ) + MultVar => VariableGet(Solver % Mesh % Variables, MultiplierName) + j = 0 + IF(ASSOCIATED(RestMatrix)) j = RestMatrix % NumberofRows + IF(ASSOCIATED(AddMatrix)) j = j+MAX(0,AddMatrix % NumberofRows-StiffMatrix % NumberOfRows) + + IF ( .NOT. ASSOCIATED(MultVar) ) THEN + CALL Info(Caller,'Creating variable for Lagrange multiplier',Level=8) + ALLOCATE( MultiplierValues(j), STAT=istat ) + IF ( istat /= 0 ) CALL Fatal(Caller,'Memory allocation error.') + + MultiplierValues = 0.0_dp + IF( ComplexSystem ) THEN + CALL VariableAddVector(Solver % Mesh % Variables, Solver % Mesh, SolverPointer, & + MultiplierName, 2, MultiplierValues) + ELSE + CALL VariableAdd(Solver % Mesh % Variables, Solver % Mesh, SolverPointer, & + MultiplierName, 1, MultiplierValues) + END IF + MultVar => VariableGet(Solver % Mesh % Variables, MultiplierName) + END IF - IF( InfoActive( 20 ) ) THEN - CALL VectorValuesRange(MultVar % Values,SIZE(MultVar % Values),TRIM(MultVar % Name)) - END IF - - MultiplierValues => MultVar % Values - - IF (j > SIZE(MultiplierValues)) THEN - CALL Info(Caller,'Increasing Lagrange multiplier size to: '//I2S(j),Level=8) - ALLOCATE(MultiplierValues(j)); MultiplierValues=0._dp - MultiplierValues(1:SIZE(MultVar % Values)) = MultVar % Values - - ! If the Lagrange variable includes history also change its size. - IF( ASSOCIATED( MultVar % PrevValues ) ) THEN - MultVar % Values = MultVar % PrevValues(:,1) - DEALLOCATE( MultVar % PrevValues ) - ALLOCATE( MultVar % PrevValues(j,1) ) - MultVar % PrevValues = 0.0_dp - MultVar % PrevValues(:,1) = MultVar % Values - END IF + IF( InfoActive( 20 ) ) THEN + CALL VectorValuesRange(MultVar % Values,SIZE(MultVar % Values),TRIM(MultVar % Name)) + END IF - DEALLOCATE(MultVar % Values) - MultVar % Values => MultiplierValues - END IF + MultiplierValues => MultVar % Values - IF( InfoActive(25) ) THEN - CALL VectorValuesRange(MultVar % values,SIZE(MultVar % values),'MultVar') - END IF + IF (j > SIZE(MultiplierValues)) THEN + CALL Info(Caller,'Increasing Lagrange multiplier size to: '//I2S(j),Level=8) + ALLOCATE(MultiplierValues(j)); MultiplierValues=0._dp + MultiplierValues(1:SIZE(MultVar % Values)) = MultVar % Values + + ! If the Lagrange variable includes history also change its size. + IF( ASSOCIATED( MultVar % PrevValues ) ) THEN + MultVar % Values = MultVar % PrevValues(:,1) + DEALLOCATE( MultVar % PrevValues ) + ALLOCATE( MultVar % PrevValues(j,1) ) + MultVar % PrevValues = 0.0_dp + MultVar % PrevValues(:,1) = MultVar % Values + END IF + + DEALLOCATE(MultVar % Values) + MultVar % Values => MultiplierValues + END IF + + IF( InfoActive(25) ) THEN + CALL VectorValuesRange(MultVar % values,SIZE(MultVar % values),'MultVar') + END IF ELSE - MultiplierValues => NULL() + MultiplierValues => NULL() END IF - UseTreeGauge = ListGetlogical( Solver % Values, 'Use Tree Gauge', Found ) + UseTreeGauge = ListGetlogical(Params, 'Use Tree Gauge', Found ) !------------------------------------------------------------------------------ ! Put the RestMatrix to lower part of CollectionMatrix !------------------------------------------------------------------------------ - EnforceDirichlet = ListGetLogical( Solver % Values, 'Enforce Exact Dirichlet BCs',Found) + EnforceDirichlet = ListGetLogical(Params, 'Enforce Exact Dirichlet BCs',Found) IF(.NOT.Found) EnforceDirichlet = .TRUE. EnforceDirichlet = EnforceDirichlet .AND. ALLOCATED(StiffMatrix % ConstrainedDOF) - UseTranspose = ListGetLogical( Solver % Values, 'Use Transpose values', Found) + UseTranspose = ListGetLogical(Params, 'Use Transpose values', Found) IF( UseTranspose ) THEN CALL Info(Caller,'Using transpose values in elimination',Level=15) END IF @@ -18684,9 +18712,9 @@ RECURSIVE SUBROUTINE SolveWithLinearRestriction( StiffMatrix, ForceVector, & //I2S(RestMatrix % NumberOfRows)//' / '//I2S(SIZE(RestMatrix % Values)),Level=12) NoEmptyRows = 0 - ConstraintScaling = ListGetLogical(Solver % Values, 'Constraint Scaling',Found) + ConstraintScaling = ListGetLogical(Params, 'Constraint Scaling',Found) IF(ConstraintScaling) THEN - rowsum = ListGetConstReal( Solver % Values, 'Constraint Scale', Found) + rowsum = ListGetConstReal(Params, 'Constraint Scale', Found) IF(Found) RestMatrix % Values = RestMatrix % Values * rowsum END IF @@ -19213,9 +19241,10 @@ RECURSIVE SUBROUTINE SolveWithLinearRestriction( StiffMatrix, ForceVector, & END IF CALL Info(Caller,'Reverting CollectionMatrix back to CRS matrix',Level=10) - IF(CollectionMatrix % FORMAT==MATRIX_LIST) & - CALL List_toCRSMatrix(CollectionMatrix) - + IF(CollectionMatrix % FORMAT==MATRIX_LIST) THEN + CALL List_toCRSMatrix(CollectionMatrix) + END IF + ! CRS-format matrix needed here IF ( NeedMassDampValues ) THEN ! Doesn't work with constraints, "AddMatrix" only !! CALL CopyMassDampValues(CollectionMatrix, StiffMatrix, AddMatrix) @@ -19232,7 +19261,7 @@ RECURSIVE SUBROUTINE SolveWithLinearRestriction( StiffMatrix, ForceVector, & i = StiffMatrix % NumberOfRows+1 j = SIZE(CollectionSolution) - IF( j>= i) THEN + IF( j >= i) THEN CollectionSolution(i:j) = 0._dp IF(ExportMultiplier) CollectionSolution(i:j) = MultiplierValues(1:j-i+1) END IF @@ -19245,7 +19274,6 @@ RECURSIVE SUBROUTINE SolveWithLinearRestriction( StiffMatrix, ForceVector, & CollectionMatrix % ExtraDOFs = CollectionMatrix % NumberOfRows - & StiffMatrix % NumberOfRows - CollectionMatrix % ParallelDOFs = 0 IF(ASSOCIATED(AddMatrix)) & CollectionMatrix % ParallelDOFs = MAX(AddMatrix % NumberOfRows - & @@ -19259,45 +19287,6 @@ RECURSIVE SUBROUTINE SolveWithLinearRestriction( StiffMatrix, ForceVector, & ! Collectionmatrix % Complex = StiffMatrix % Complex - ! We may want to skip the constraints for norm if we use certain other options - SkipConstraints = ListGetLogical( Solver % values, & - 'Nonlinear System Convergence Without Constraints',Found ) - IF(.NOT. Found ) THEN - SkipConstraints = ListGetLogical( Solver % values, 'Linear System Residual Mode',Found ) - IF( SkipConstraints ) THEN - CALL Info(Caller,'Linear system residual mode must skip constraints',Level=10) - ELSE - SkipConstraints = ListGetLogical( Solver % values, 'NonLinear System Consistent Norm',Found ) - IF( SkipConstraints ) THEN - CALL Info(Caller,'Nonlinear system consistent norm must skip constraints',Level=10) - END IF - END IF - str = ListGetString( Solver % values, 'NonLinear System Convergence Measure',Found ) - IF( str == 'solution' ) THEN - SkipConstraints = .TRUE. - CALL Info(Caller,& - 'Nonlinear system convergence measure == "solution" must skip constraints',Level=10) - END IF - IF( SkipConstraints ) THEN - CALL Info(Caller,'Enforcing convergence without constraints to True',Level=10) - CALL ListAddLogical( Solver % Values, & - 'Nonlinear System Convergence Without Constraints',.TRUE.) - END IF - END IF - - !------------------------------------------------------------------------------ - ! Look at the nonlinear system previous values again, not taking the constrained - ! system into account... - !------------------------------------------------------------------------------ - Found = ASSOCIATED(Solver % Variable % NonlinValues) - IF( Found .AND. .NOT. SkipConstraints ) THEN - k = CollectionMatrix % NumberOfRows - IF ( SIZE(Solver % Variable % NonlinValues) /= k) THEN - DEALLOCATE(Solver % Variable % NonlinValues) - ALLOCATE(Solver % Variable % NonlinValues(k)) - END IF - Solver % Variable % NonlinValues(1:k) = CollectionSolution(1:k) - END IF CollectionMatrix % Comm = StiffMatrix % Comm @@ -19316,19 +19305,86 @@ RECURSIVE SUBROUTINE SolveWithLinearRestriction( StiffMatrix, ForceVector, & CALL VectorValuesRange(CollectionMatrix % Values,SIZE(CollectionMatrix % Values),'A') CALL VectorValuesRange(CollectionMatrix % rhs,SIZE(CollectionMatrix % rhs),'b') END IF - - CALL Info(Caller,'Now going for the coupled linear system',Level=10) + + IF( ResidualMode ) THEN + BLOCK + REAL(KIND=dp), POINTER :: Res(:) + ! If residual mode is requested make change of variables: + ! Ax=b -> Adx = b-Ax0 = r + IF( niter > 1 ) THEN + CALL Info(Caller,'Changing the equation to residual based mode',Level=10) + ALLOCATE( Res(SIZE(CollectionSolution)) ) + CALL LinearSystemResidual( CollectionMatrix, CollectionVector, CollectionSolution, res ) + CollectionVector = Res + CollectionSolution = 0.0_dp + DEALLOCATE(Res) + END IF + END BLOCK + END IF + + ! We may want to skip ComputeChange including the constraints if we use certain other options + SkipConstraints = ResidualMode .OR. & + ListGetLogical( Params, 'Nonlinear System Convergence Without Constraints',Found ) .OR. & + ListGetLogical( Params, 'NonLinear System Consistent Norm',Found ) + str = ListGetString( Params, 'NonLinear System Convergence Measure',Found ) + IF( str == 'solution' ) THEN + SkipConstraints = .TRUE. + CALL Info(Caller,& + 'Nonlinear system convergence measure == "solution" must skip constraints',Level=10) + END IF + IF( SkipConstraints ) THEN + CALL ListAddLogical( Params,'Skip Compute Nonlinear Change',.TRUE.) + CALL ListAddLogical( Params,'Skip Advance Nonlinear iter',.TRUE.) + END IF + DoOwnScaling = ListGetLogical( Params,'Linear System Scaling',Found) + IF(.NOT. Found) DoOwnScaling = .TRUE. + IF(.NOT. ResidualMode) DoOwnScaling = .FALSE. + IF(DoOwnScaling) THEN + CALL Info(Caller,'Performing special scaling with constraints',Level=10) + DiagScaling => CollectionMatrix % DiagScaling + IF(Niter == 1 ) THEN + IF(.NOT. ASSOCIATED(DiagScaling) ) THEN + ALLOCATE( DiagScaling(SIZE(CollectionVector))) + CollectionMatrix % DiagScaling => DiagScaling + END IF + + ! Should we scale only part or the full matrix? + IF(.FALSE.) THEN + DiagScaling = 1.0_dp + StiffMatrix % DiagScaling => DiagScaling + ! Just build the scaling matrix using only the original stiffness matrix. + CALL ScaleLinearSystem(Solver,StiffMatrix,ApplyScaling=.FALSE.) + CollectionMatrix % ScalingMethod = StiffMatrix % ScalingMethod + StiffMatrix % DiagScaling => NULL() + ELSE + CALL ScaleLinearSystem(Solver,CollectionMatrix,ApplyScaling=.FALSE.) + END IF + END IF + + CALL ScaleLinearSystem(Solver,CollectionMatrix,CollectionVector,& + CollectionSolution,DiagScaling=CollectionMatrix % DiagScaling) + CALL ListAddLogical( Params,'Linear System Skip Scaling',.TRUE. ) + END IF + + CALL Info(Caller,'Now solving the linear system with constraints!',Level=10) Collectionmatrix % DGMatrix = StiffMatrix % DGMatrix CALL SolveLinearSystem( CollectionMatrix, CollectionVector, & CollectionSolution, Norm, DOFs, Solver, StiffMatrix ) + + IF(DoOwnScaling) THEN + CALL BackScaleLinearSystem( Solver,CollectionMatrix,CollectionVector,& + CollectionSolution,CollectionMatrix % DiagScaling) + CALL ListAddLogical( Params,'Linear System Skip Scaling',.FALSE. ) + END IF + !------------------------------------------------------------------------------- ! For restricted systems study the norm without some block components. ! For example, excluding gauge constraints may give valuable information ! of the real accuracy of the unconstrained system. Currently just for info. !------------------------------------------------------------------------------- - IF( ListGetLogical( Solver % Values,'Restricted System Norm',Found ) ) THEN + IF( ListGetLogical( Params,'Restricted System Norm',Found ) ) THEN ALLOCATE( TrueDof( CollectionMatrix % NumberOfRows ) ) TrueDof = .TRUE. @@ -19338,14 +19394,14 @@ RECURSIVE SUBROUTINE SolveWithLinearRestriction( StiffMatrix, ForceVector, & WRITE( Message,'(A,ES13.6)') 'Residual norm of the original system:',Norm CALL Info(Caller,Message, Level = 5 ) - IF( ListGetLogical( Solver % Values,'Restricted System Norm Skip Nodes',Found ) ) THEN + IF( ListGetLogical( Params,'Restricted System Norm Skip Nodes',Found ) ) THEN i = 1 j = MAXVAL( Solver % Variable % Perm(1:Solver % Mesh % NumberOfNodes) ) CALL Info(Caller,'Skipping nodal dof range: '//I2S(i)//'-'//I2S(j),Level=8) TrueDof(i:j) = .FALSE. END IF - IF( ListGetLogical( Solver % Values,'Restricted System Norm Skip Constraints',Found ) ) THEN + IF( ListGetLogical( Params,'Restricted System Norm Skip Constraints',Found ) ) THEN i = StiffMatrix % NumberOfRows + 1 j = CollectionMatrix % NumberOfRows CALL Info(Caller,'Skipping constraints dof range: '& @@ -19368,11 +19424,14 @@ RECURSIVE SUBROUTINE SolveWithLinearRestriction( StiffMatrix, ForceVector, & !------------------------------------------------------------------------------ CALL Info(Caller,'Picking solution from collection solution',Level=10) - Solution = 0.0_dp - i = 1 j = StiffMatrix % NumberOfRows - Solution(i:j) = CollectionSolution(i:j) + IF( ResidualMode .AND. nIter > 1) THEN + Solution(1:j) = Solution(1:j) + CollectionSolution(1:j) + ELSE + Solution(1:j) = CollectionSolution(1:j) + END IF + IF ( ExportMultiplier ) THEN CALL Info(Caller,'Separating Lagrange multiplier from collection solution',Level=10) @@ -19385,6 +19444,10 @@ RECURSIVE SUBROUTINE SolveWithLinearRestriction( StiffMatrix, ForceVector, & IF(ASSOCIATED(RestMatrix) .AND. EliminateConstraints) THEN ! Compute eliminated l-coefficient values: ! --------------------------------------- + IF( ResidualMode ) THEN + CALL Fatal(Caller,'Elimination not possible with ResidualMode!') + END IF + MultiplierValues = 0.0_dp DO i=1,RestMatrix % NumberOfRows scl = 1._dp / UseDiag(i) @@ -19392,19 +19455,29 @@ RECURSIVE SUBROUTINE SolveWithLinearRestriction( StiffMatrix, ForceVector, & MultiplierValues(i) = scl * ForceVector(m) DO j=StiffMatrix % Rows(m), StiffMatrix % Rows(m+1)-1 MultiplierValues(i) = MultiplierValues(i) - & - scl * StiffMatrix % Values(j) * Solution(StiffMatrix % Cols(j)) + scl * StiffMatrix % Values(j) * Solution(StiffMatrix % Cols(j)) END DO END DO ELSE - Relax = ListGetCReal( Solver % Values,'Lagrange Multiplier Relaxation Factor', Found ) - IF( Found ) THEN - MultiplierValues(1:j) = (1-Relax) * MultiplierValues(1:j) + & - Relax * CollectionSolution(i+1:i+j) - ELSE - MultiplierValues(1:j) = CollectionSolution(i+1:i+j) + Relax = ListGetCReal( Params,'Lagrange Multiplier Relaxation Factor', Found ) + IF( ResidualMode .AND. nIter > 1 ) THEN + IF( Found ) THEN + MultiplierValues(1:j) = MultiplierValues(1:j) + & + Relax * CollectionSolution(i+1:i+j) + ELSE + MultiplierValues(1:j) = MultiplierValues(1:j) + CollectionSolution(i+1:i+j) + END IF + ELSE + IF( Found ) THEN + MultiplierValues(1:j) = (1-Relax) * MultiplierValues(1:j) + & + Relax * CollectionSolution(i+1:i+j) + ELSE + MultiplierValues(1:j) = CollectionSolution(i+1:i+j) + END IF END IF END IF + IF(EliminateConstraints .AND. EliminateDiscont) THEN IF (EliminateFromMaster) THEN CALL totv(StiffMatrix,MultiplierValues,MasterIPerm) @@ -19416,10 +19489,17 @@ RECURSIVE SUBROUTINE SolveWithLinearRestriction( StiffMatrix, ForceVector, & !------------------------------------------------------------------------------ + IF( SkipConstraints ) THEN + CALL ListAddLogical( Params,'Skip Advance Nonlinear iter',.FALSE.) + CALL ListAddLogical( Params,'Skip Compute Nonlinear Change',.FALSE.) + CALL ComputeChange(Solver,.FALSE.,StiffMatrix % NumberOfRows,Matrix=StiffMatrix,Rhs=ForceVector) + END IF + StiffMatrix % CollectionMatrix => CollectionMatrix DEALLOCATE(CollectionSolution) CollectionMatrix % ConstraintMatrix => NULL() + CALL Info( Caller, 'All done', Level=10 ) CONTAINS diff --git a/fem/tests/FluxIntegralBCInduction/case.sif b/fem/tests/FluxIntegralBCInduction/case.sif index 0c38f2c2f9..cadae070ea 100755 --- a/fem/tests/FluxIntegralBCInduction/case.sif +++ b/fem/tests/FluxIntegralBCInduction/case.sif @@ -61,8 +61,8 @@ Solver 1 Exec Solver = Always Nonlinear System Max Iterations = 1 - Linear System Solver = Direct !Iterative - Linear System Direct Method = MUMPS + Linear System Solver = Direct + Linear System Direct Method = umfpack Linear System Scaling = False Optimize Bandwidth = False