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

@@ -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

src/Solution/ParallelSolution.f90

@@ -183,24 +183,30 @@ 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
     integer :: ierr
 
     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

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)