From 1fa35eea9492acf67a7017b765d36076a17a39b9 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Thu, 24 Apr 2025 16:08:54 -0400 Subject: [PATCH] (*)Fix northern Flather OBCs in southern halos Added additional logic to avoid segmentation faults from applying Flather open boundary conditions near the opposite side of the computational domain (e.g., applying a Northern Flather boundary condition on the southern edge of the domain). This situation does not arise when open boundaries are only found on the edges of a domain, but in parallel cases with sloped OBCs (such as the ESMG-configs rotated seamount test case) this atypical situation can arise. The key point to the specific fix is to note that when such a case arises, the OBC velocity in the outer halo point does not subsequently impact the solution on that PE either because its effects are masked out by land or OBC points at the neighboring orthogonal velocities, and then they are overwritten by a halo update (perhaps several steps later). All answers are bitwise identical in all cases that have been tested and were working before, and the ESMG/rotated_seamount test case is now running correctly (and giving identical answers) for a variety of PE counts. --- src/core/MOM_barotropic.F90 | 65 +++++++++++++++++++++++++++++-------- 1 file changed, 51 insertions(+), 14 deletions(-) diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90 index dd34f722c3..f8bc982d18 100644 --- a/src/core/MOM_barotropic.F90 +++ b/src/core/MOM_barotropic.F90 @@ -3686,7 +3686,11 @@ subroutine apply_u_velocity_OBCs(ubt, uhbt, ubt_trans, eta, SpV_avg, ubt_old, BT elseif (BT_OBC%u_OBC_type(I,j) == FLATHER_OBC) then ! Eastern Flather OBC cfl = dtbt * BT_OBC%Cg_u(I,j) * G%IdxCu(I,j) ! CFL u_inlet = cfl*ubt_old(I-1,j) + (1.0-cfl)*ubt_old(I,j) ! Valid for cfl<1 - if (GV%Boussinesq) then + if (I <= MS%isdw) then + ! Do not apply an Eastern Flather OBC at the western halo points on a PE, as doing so would + ! create a segmentation fault and this velocity will not propagate through to the next iteration. + ssh_in = BT_OBC%SSH_outer_u(I,j) + elseif (GV%Boussinesq) then ssh_in = GV%H_to_Z*(eta(i,j) + (0.5-cfl)*(eta(i,j)-eta(i-1,j))) ! internal else ssh_1 = GV%H_to_RZ * eta(i,j) * SpV_avg(i,j) - (CS%bathyT(i,j) + G%Z_ref) @@ -3739,7 +3743,11 @@ subroutine apply_u_velocity_OBCs(ubt, uhbt, ubt_trans, eta, SpV_avg, ubt_old, BT elseif (BT_OBC%u_OBC_type(I,j) == -FLATHER_OBC) then ! Western Flather OBC cfl = dtbt * BT_OBC%Cg_u(I,j) * G%IdxCu(I,j) ! CFL u_inlet = cfl*ubt_old(I+1,j) + (1.0-cfl)*ubt_old(I,j) ! Valid for cfl<1 - if (GV%Boussinesq) then + if (I >= MS%iedw-1) then + ! Do not apply a Western Flather OBC at the eastern halo points on a PE, as doing so would + ! create a segmentation fault and this velocity will not propagate through to the next iteration. + ssh_in = BT_OBC%SSH_outer_u(I,j) + elseif (GV%Boussinesq) then ssh_in = GV%H_to_Z*(eta(i+1,j) + (0.5-cfl)*(eta(i+1,j)-eta(i+2,j))) ! internal else ssh_1 = GV%H_to_RZ * eta(i+1,j) * SpV_avg(i+1,j) - (CS%bathyT(i+1,j) + G%Z_ref) @@ -3872,7 +3880,11 @@ subroutine apply_v_velocity_OBCs(vbt, vhbt, vbt_trans, eta, SpV_avg, vbt_old, BT elseif (BT_OBC%v_OBC_type(i,J) == FLATHER_OBC) then ! Northern Flather OBC cfl = dtbt * BT_OBC%Cg_v(i,J) * G%IdyCv(i,J) ! CFL v_inlet = cfl*vbt_old(i,J-1) + (1.0-cfl)*vbt_old(i,J) ! Valid for cfl<1 - if (GV%Boussinesq) then + if (J <= MS%jsdw) then + ! Do not apply a Northern Flather OBC at the southern halo points on a PE, as doing so would + ! create a segmentation fault and this velocity will not propagate through to the next iteration. + ssh_in = BT_OBC%SSH_outer_v(i,J) + elseif (GV%Boussinesq) then ssh_in = GV%H_to_Z*(eta(i,j) + (0.5-cfl)*(eta(i,j)-eta(i,j-1))) ! internal else ssh_1 = GV%H_to_RZ * eta(i,j) * SpV_avg(i,j) - (CS%bathyT(i,j) + G%Z_ref) @@ -3926,7 +3938,11 @@ subroutine apply_v_velocity_OBCs(vbt, vhbt, vbt_trans, eta, SpV_avg, vbt_old, BT elseif (BT_OBC%v_OBC_type(i,J) == -FLATHER_OBC) then ! Southern Flather OBC cfl = dtbt * BT_OBC%Cg_v(i,J) * G%IdyCv(i,J) ! CFL v_inlet = cfl*vbt_old(i,J+1) + (1.0-cfl)*vbt_old(i,J) ! Valid for cfl <1 - if (GV%Boussinesq) then + if (J >= MS%jedw-1) then + ! Do not apply a Southern Flather OBC at the northern halo points on a PE, as doing so would + ! create a segmentation fault and this velocity will not propagate through to the next iteration. + ssh_in = BT_OBC%SSH_outer_v(i,J) + elseif (GV%Boussinesq) then ssh_in = GV%H_to_Z*(eta(i,j+1) + (0.5-cfl)*(eta(i,j+1)-eta(i,j+2))) ! internal else ssh_1 = GV%H_to_RZ * eta(i,j+1) * SpV_avg(i,j+1) - (CS%bathyT(i,j+1) + G%Z_ref) @@ -3987,7 +4003,7 @@ subroutine initialize_BT_OBC(OBC, BT_OBC, G, CS) real :: OBC_sign ! A sign encoding the direction of the OBC being used at a point [nondim] real :: OBC_type ! A real copy of the integer encoding the type of OBC being used at a point [nondim] integer :: i, j, isdw, iedw, jsdw, jedw - integer :: l_seg + integer :: l_seg, Flather_OBC_in_halo isdw = CS%isdw ; iedw = CS%iedw ; jsdw = CS%jsdw ; jedw = CS%jedw @@ -4037,27 +4053,48 @@ subroutine initialize_BT_OBC(OBC, BT_OBC, G, CS) BT_OBC%is_v_N_obc = iedw + 1 ; BT_OBC%ie_v_N_obc = isdw - 1 BT_OBC%Js_v_N_obc = jedw + 1 ; BT_OBC%Je_v_N_obc = jsdw - 2 + Flather_OBC_in_halo = 0 do j=jsdw,jedw ; do I=isdw-1,iedw BT_OBC%u_OBC_type(I,j) = nint(u_OBC(I,j)) if (BT_OBC%u_OBC_type(I,j) < 0) then ! This point has OBC_DIRECTION_W. - BT_OBC%Is_u_W_obc = min(I, BT_OBC%Is_u_W_obc) ; BT_OBC%Ie_u_W_obc = max(I, BT_OBC%Ie_u_W_obc) - BT_OBC%js_u_W_obc = min(j, BT_OBC%js_u_W_obc) ; BT_OBC%je_u_W_obc = max(j, BT_OBC%je_u_W_obc) + if ((BT_OBC%u_OBC_type(I,j) == -FLATHER_OBC) .and. (I >= iedw-1)) then + ! There is no need to specify the OBC at this point, but the stencil might need to be increased. + Flather_OBC_in_halo = 1 + else + BT_OBC%Is_u_W_obc = min(I, BT_OBC%Is_u_W_obc) ; BT_OBC%Ie_u_W_obc = max(I, BT_OBC%Ie_u_W_obc) + BT_OBC%js_u_W_obc = min(j, BT_OBC%js_u_W_obc) ; BT_OBC%je_u_W_obc = max(j, BT_OBC%je_u_W_obc) + endif endif if (BT_OBC%u_OBC_type(I,j) > 0) then ! This point has OBC_DIRECTION_E. - BT_OBC%Is_u_E_obc = min(I, BT_OBC%Is_u_E_obc) ; BT_OBC%Ie_u_E_obc = max(I, BT_OBC%Ie_u_E_obc) - BT_OBC%js_u_E_obc = min(j, BT_OBC%js_u_E_obc) ; BT_OBC%je_u_E_obc = max(j, BT_OBC%je_u_E_obc) + if ((BT_OBC%u_OBC_type(I,j) == FLATHER_OBC) .and. (I <= isdw)) then + ! There is no need to specify the OBC at this point, but the stencil might need to be increased. + Flather_OBC_in_halo = 1 + else + BT_OBC%Is_u_E_obc = min(I, BT_OBC%Is_u_E_obc) ; BT_OBC%Ie_u_E_obc = max(I, BT_OBC%Ie_u_E_obc) + BT_OBC%js_u_E_obc = min(j, BT_OBC%js_u_E_obc) ; BT_OBC%je_u_E_obc = max(j, BT_OBC%je_u_E_obc) + endif endif enddo ; enddo do J=jsdw-1,jedw ; do i=isdw,iedw BT_OBC%v_OBC_type(i,J) = nint(v_OBC(i,J)) - if (BT_OBC%v_OBC_type(i,J) < 0) then ! This point has OBC_DIRECTION_S. - BT_OBC%is_v_S_obc = min(i, BT_OBC%is_v_S_obc) ; BT_OBC%ie_v_S_obc = max(i, BT_OBC%ie_v_S_obc) - BT_OBC%Js_v_S_obc = min(J, BT_OBC%Js_v_S_obc) ; BT_OBC%Je_v_S_obc = max(J, BT_OBC%Je_v_S_obc) + if (BT_OBC%v_OBC_type(i,J) < 0) then ! This point has OBC_DIRECTION_S. + if ((BT_OBC%v_OBC_type(i,J) == -FLATHER_OBC) .and. (J >= jedw-1)) then + ! There is no need to specify the OBC at this point, but the stencil might need to be increased. + Flather_OBC_in_halo = 1 + else + BT_OBC%is_v_S_obc = min(i, BT_OBC%is_v_S_obc) ; BT_OBC%ie_v_S_obc = max(i, BT_OBC%ie_v_S_obc) + BT_OBC%Js_v_S_obc = min(J, BT_OBC%Js_v_S_obc) ; BT_OBC%Je_v_S_obc = max(J, BT_OBC%Je_v_S_obc) + endif endif if (BT_OBC%v_OBC_type(i,J) > 0) then ! This point has OBC_DIRECTION_N. - BT_OBC%is_v_N_obc = min(i, BT_OBC%is_v_N_obc) ; BT_OBC%ie_v_N_obc = max(i, BT_OBC%ie_v_N_obc) - BT_OBC%Js_v_N_obc = min(J, BT_OBC%Js_v_N_obc) ; BT_OBC%Je_v_N_obc = max(J, BT_OBC%Je_v_N_obc) + if ((BT_OBC%v_OBC_type(i,J) == FLATHER_OBC) .and. (J <= jsdw)) then + ! There is no need to specify the OBC at this point, but the stencil might need to be increased. + Flather_OBC_in_halo = 1 + else + BT_OBC%is_v_N_obc = min(i, BT_OBC%is_v_N_obc) ; BT_OBC%ie_v_N_obc = max(i, BT_OBC%ie_v_N_obc) + BT_OBC%Js_v_N_obc = min(J, BT_OBC%Js_v_N_obc) ; BT_OBC%Je_v_N_obc = max(J, BT_OBC%Je_v_N_obc) + endif endif enddo ; enddo