From e6202c0a79d70b978599e92986cc109bca40d38d Mon Sep 17 00:00:00 2001 From: mjr-deltares <45555666+mjr-deltares@users.noreply.github.com> Date: Wed, 8 May 2024 19:50:55 +0200 Subject: [PATCH] fix(par): fix backtracking and residual norm calculation for certain cases in parallel (#1780) - update petsc matrix before multiplication - backtracking now properly synchronized --- src/Solution/NumericalSolution.f90 | 90 ++++++++++++++++++---------- src/Solution/ParallelSolution.f90 | 16 +++-- src/Utilities/Matrix/PetscMatrix.F90 | 7 ++- 3 files changed, 74 insertions(+), 39 deletions(-) diff --git a/src/Solution/NumericalSolution.f90 b/src/Solution/NumericalSolution.f90 index 7405cba2af7..07629dc4d0c 100644 --- a/src/Solution/NumericalSolution.f90 +++ b/src/Solution/NumericalSolution.f90 @@ -168,6 +168,8 @@ module NumericalSolutionModule procedure :: sln_calc_ptc procedure :: sln_underrelax procedure :: sln_backtracking_xupdate + procedure :: get_backtracking_flag + procedure :: apply_backtracking ! private procedure, private :: sln_connect @@ -179,7 +181,7 @@ module NumericalSolutionModule procedure, private :: sln_calcdx procedure, private :: sln_calc_residual procedure, private :: sln_l2norm - procedure, private :: sln_outer_check + procedure, private :: sln_get_dxmax procedure, private :: sln_get_loc procedure, private :: sln_get_nodeu procedure, private :: allocate_scalars @@ -1617,7 +1619,7 @@ subroutine solve(this, kiter) !------------------------------------------------------- ! ! -- check convergence of solution - call this%sln_outer_check(this%hncg(kiter), this%lrch(1, kiter)) + call this%sln_get_dxmax(this%hncg(kiter), this%lrch(1, kiter)) this%icnvg = 0 if (this%sln_has_converged(this%hncg(kiter))) then this%icnvg = 1 @@ -1752,7 +1754,7 @@ subroutine solve(this, kiter) this%icnvg = 1 ! ! -- reset outer dependent-variable change and location for output - call this%sln_outer_check(this%hncg(kiter), this%lrch(1, kiter)) + call this%sln_get_dxmax(this%hncg(kiter), this%lrch(1, kiter)) ! ! -- write revised dependent-variable change data after ! newton under-relaxation @@ -2790,41 +2792,65 @@ end subroutine sln_backtracking !! update exceeds the dependent variable closure criteria. !! !< - subroutine sln_backtracking_xupdate(this, btflag) + subroutine sln_backtracking_xupdate(this, bt_flag) ! -- dummy variables class(NumericalSolutionType), intent(inout) :: this !< NumericalSolutionType instance - integer(I4B), intent(inout) :: btflag !< backtracking flag (1) backtracking performed (0) backtracking not performed - ! -- local variables + integer(I4B), intent(inout) :: bt_flag !< backtracking flag (1) backtracking performed (0) backtracking not performed + + bt_flag = this%get_backtracking_flag() + + ! perform backtracking if ... + if (bt_flag > 0) then + call this%apply_backtracking() + end if + + end subroutine sln_backtracking_xupdate + + !> @brief Check if backtracking should be applied for this solution, + !< returns 1: yes, 0: no + function get_backtracking_flag(this) result(bt_flag) + class(NumericalSolutionType) :: this !< NumericalSolutionType instance + integer(I4B) :: bt_flag !< backtracking flag (1) backtracking performed (0) backtracking not performed + ! local + integer(I4B) :: n + real(DP) :: dx + real(DP) :: dx_abs + real(DP) :: dx_abs_max + + ! default is off + bt_flag = 0 + + ! find max. change + dx_abs_max = 0.0 + do n = 1, this%neq + if (this%active(n) < 1) cycle + dx = this%x(n) - this%xtemp(n) + dx_abs = abs(dx) + if (dx_abs > dx_abs_max) dx_abs_max = dx_abs + end do + + ! if backtracking, set flag + if (this%breduc * dx_abs_max >= this%dvclose) then + bt_flag = 1 + end if + + end function get_backtracking_flag + + !> @brief Update x with backtracking + !< + subroutine apply_backtracking(this) + class(NumericalSolutionType) :: this !< NumericalSolutionType instance + ! local integer(I4B) :: n real(DP) :: delx - real(DP) :: absdelx - real(DP) :: chmax - ! - ! -- initialize dummy variables - btflag = 0 - ! - ! -- no backtracking if maximum change is less than closure so return - chmax = 0.0 + do n = 1, this%neq if (this%active(n) < 1) cycle delx = this%breduc * (this%x(n) - this%xtemp(n)) - absdelx = abs(delx) - if (absdelx > chmax) chmax = absdelx + this%x(n) = this%xtemp(n) + delx end do - ! - ! -- perform backtracking if free of constraints and set counter and flag - if (chmax >= this%dvclose) then - btflag = 1 - do n = 1, this%neq - if (this%active(n) < 1) cycle - delx = this%breduc * (this%x(n) - this%xtemp(n)) - this%x(n) = this%xtemp(n) + delx - end do - end if - ! - ! -- return - return - end subroutine sln_backtracking_xupdate + + end subroutine !> @ brief Calculate the solution L-2 norm for all !! active cells using @@ -3116,7 +3142,7 @@ end subroutine sln_underrelax !! Picard iteration. !! !< - subroutine sln_outer_check(this, hncg, lrch) + subroutine sln_get_dxmax(this, hncg, lrch) ! -- dummy variables class(NumericalSolutionType), intent(inout) :: this !< NumericalSolutionType instance real(DP), intent(inout) :: hncg !< maximum dependent-variable change @@ -3150,7 +3176,7 @@ subroutine sln_outer_check(this, hncg, lrch) ! ! -- return return - end subroutine sln_outer_check + end subroutine sln_get_dxmax function sln_has_converged(this, max_dvc) result(has_converged) class(NumericalSolutionType) :: this !< NumericalSolutionType instance diff --git a/src/Solution/ParallelSolution.f90 b/src/Solution/ParallelSolution.f90 index 09c00392ef9..22d5353d522 100644 --- a/src/Solution/ParallelSolution.f90 +++ b/src/Solution/ParallelSolution.f90 @@ -183,10 +183,10 @@ end subroutine par_underrelax !> @brief synchronize backtracking flag over processes !< - subroutine par_backtracking_xupdate(this, btflag) + subroutine par_backtracking_xupdate(this, bt_flag) ! -- dummy variables class(ParallelSolutionType), intent(inout) :: this !< ParallelSolutionType instance - integer(I4B), intent(inout) :: btflag !< global backtracking flag (1) backtracking performed (0) backtracking not performed + integer(I4B), intent(inout) :: bt_flag !< global backtracking flag (1) backtracking performed (0) backtracking not performed ! -- local variables integer(I4B) :: btflag_local type(MpiWorldType), pointer :: mpi_world @@ -194,13 +194,19 @@ subroutine par_backtracking_xupdate(this, btflag) mpi_world => get_mpi_world() - btflag_local = 0 - call this%NumericalSolutionType%sln_backtracking_xupdate(btflag_local) + ! get local bt flag + btflag_local = this%NumericalSolutionType%get_backtracking_flag() - call MPI_Allreduce(btflag_local, btflag, 1, MPI_INTEGER, & + ! reduce into global decision (if any, then all) + call MPI_Allreduce(btflag_local, bt_flag, 1, MPI_INTEGER, & MPI_MAX, mpi_world%comm, ierr) call CHECK_MPI(ierr) + ! perform backtracking if ... + if (bt_flag > 0) then + call this%NumericalSolutionType%apply_backtracking() + end if + end subroutine par_backtracking_xupdate end module ParallelSolutionModule diff --git a/src/Utilities/Matrix/PetscMatrix.F90 b/src/Utilities/Matrix/PetscMatrix.F90 index c6e6654637d..440e9b71f96 100644 --- a/src/Utilities/Matrix/PetscMatrix.F90 +++ b/src/Utilities/Matrix/PetscMatrix.F90 @@ -12,8 +12,8 @@ module PetscMatrixModule private type, public, extends(MatrixBaseType) :: PetscMatrixType - Mat :: mat - ! offset in the global matrix + Mat :: mat !< the PETSc matrix object, NOTE: update() should be called before using this, + !! in case the matrix CSR array has changed!!! integer(I4B) :: nrow !< number of rows in this portion of the global matrix integer(I4B) :: ncol !< number of columns in the matrix integer(I4B) :: nnz !< number of nonzeros in the matrix @@ -446,6 +446,9 @@ subroutine pm_multiply(this, vec_x, vec_y) y => vec_y end select + ! copy data into petsc object + call this%update() + ! and multiply call MatMult(this%mat, x%vec_impl, y%vec_impl, ierr) CHKERRQ(ierr)