From 762a915af10db056d2fd7639297bd7fb1482647b Mon Sep 17 00:00:00 2001 From: Philippe Blain Date: Tue, 12 Jul 2022 16:38:06 -0400 Subject: [PATCH] ice_dyn_vp: do not skip halo updates in 'pgmres' under 'bfbflag' The 'pgmres' subroutine implements a separate GMRES solver and is used as a preconditioner for the FGMRES linear solver. Since it is only a preconditioner, it was decided to skip the halo updates after computing the matrix-vector product (in 'matvec'), for efficiency. This leads to non-reproducibility since the content of the non-updated halos depend on the block / MPI distribution. Add the required halo updates, but only perform them when we are explicitely asking for bit-for-bit global sums, i.e. when 'bfbflag' is set to something else than 'not'. Adjust the interfaces of 'pgmres' and 'precondition' (from which 'pgmres' is called) to accept 'halo_info_mask', since it is needed for masked updates. Closes https://github.com/CICE-Consortium/CICE/issues/518 --- cicecore/cicedynB/dynamics/ice_dyn_vp.F90 | 35 +++++++++++++++++++++-- 1 file changed, 33 insertions(+), 2 deletions(-) diff --git a/cicecore/cicedynB/dynamics/ice_dyn_vp.F90 b/cicecore/cicedynB/dynamics/ice_dyn_vp.F90 index 4e0740e15..b62f8aca5 100644 --- a/cicecore/cicedynB/dynamics/ice_dyn_vp.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_vp.F90 @@ -3016,6 +3016,7 @@ subroutine fgmres (zetax2 , etax2 , & call precondition(zetax2 , etax2 , & Cb , vrel , & umassdti , & + halo_info_mask, & arnoldi_basis_x(:,:,:,initer), & arnoldi_basis_y(:,:,:,initer), & diagx , diagy , & @@ -3212,6 +3213,7 @@ end subroutine fgmres subroutine pgmres (zetax2 , etax2 , & Cb , vrel , & umassdti , & + halo_info_mask, & bx , by , & diagx , diagy , & tolerance, maxinner, & @@ -3219,6 +3221,11 @@ subroutine pgmres (zetax2 , etax2 , & solx , soly , & nbiter) + use ice_boundary, only: ice_HaloUpdate + use ice_domain, only: maskhalo_dyn, halo_info + use ice_fileunits, only: bfbflag + use ice_timers, only: ice_timer_start, ice_timer_stop, timer_bound + real (kind=dbl_kind), dimension(nx_block,ny_block,max_blocks,4), intent(in) :: & zetax2 , & ! zetax2 = 2*zeta (bulk viscosity) etax2 ! etax2 = 2*eta (shear viscosity) @@ -3228,6 +3235,9 @@ subroutine pgmres (zetax2 , etax2 , & Cb , & ! seabed stress coefficient umassdti ! mass of U-cell/dte (kg/m^2 s) + type (ice_halo), intent(in) :: & + halo_info_mask ! ghost cell update info for masked halo + real (kind=dbl_kind), dimension(nx_block, ny_block, max_blocks), intent(in) :: & bx , & ! Right hand side of the linear system (x components) by ! Right hand side of the linear system (y components) @@ -3397,14 +3407,29 @@ subroutine pgmres (zetax2 , etax2 , & call precondition(zetax2 , etax2 , & Cb , vrel , & umassdti , & + halo_info_mask, & arnoldi_basis_x(:,:,:,initer), & arnoldi_basis_y(:,:,:,initer), & diagx , diagy , & precond_type, & workspace_x , workspace_y) - ! NOTE: halo updates for (workspace_x, workspace_y) - ! are skipped here for efficiency since this is just a preconditioner + ! Update workspace with boundary values + ! NOTE: skipped for efficiency since this is just a preconditioner + ! unless bfbflag is active + if (bfbflag /= 'off') then + call stack_fields(workspace_x, workspace_y, fld2) + call ice_timer_start(timer_bound) + if (maskhalo_dyn) then + call ice_HaloUpdate (fld2, halo_info_mask, & + field_loc_NEcorner, field_type_vector) + else + call ice_HaloUpdate (fld2, halo_info, & + field_loc_NEcorner, field_type_vector) + endif + call ice_timer_stop(timer_bound) + call unstack_fields(fld2, workspace_x, workspace_y) + endif !$OMP PARALLEL DO PRIVATE(iblk) do iblk = 1, nblocks @@ -3528,6 +3553,7 @@ subroutine pgmres (zetax2 , etax2 , & call precondition(zetax2 , etax2 , & Cb , vrel , & umassdti , & + halo_info_mask, & workspace_x , workspace_y, & diagx , diagy , & precond_type, & @@ -3594,6 +3620,7 @@ end subroutine pgmres subroutine precondition(zetax2 , etax2, & Cb , vrel , & umassdti , & + halo_info_mask, & vx , vy , & diagx , diagy, & precond_type, & @@ -3608,6 +3635,9 @@ subroutine precondition(zetax2 , etax2, & Cb , & ! seabed stress coefficient umassdti ! mass of U-cell/dte (kg/m^2 s) + type (ice_halo), intent(in) :: & + halo_info_mask ! ghost cell update info for masked halo + real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks), intent(in) :: & vx , & ! input vector (x components) vy ! input vector (y components) @@ -3669,6 +3699,7 @@ subroutine precondition(zetax2 , etax2, & call pgmres (zetax2, etax2 , & Cb , vrel , & umassdti , & + halo_info_mask , & vx , vy , & diagx , diagy , & tolerance, maxinner, &