From 15bcdadad701d17316bdfc60a4ef7cf542bd00e6 Mon Sep 17 00:00:00 2001 From: alex-huth Date: Tue, 30 Apr 2024 12:11:06 -0400 Subject: [PATCH] For iKID and calving from shelves: fixed some parallelization issues related to bergs that are on the shared edges of multiple grid cells --- src/ice_shelf_calving.F90 | 34 +++++++++++-------- src/icebergs.F90 | 1 + src/icebergs_framework.F90 | 68 +++++++++++++++++++++++--------------- 3 files changed, 63 insertions(+), 40 deletions(-) diff --git a/src/ice_shelf_calving.F90 b/src/ice_shelf_calving.F90 index 5276246..acdd526 100644 --- a/src/ice_shelf_calving.F90 +++ b/src/ice_shelf_calving.F90 @@ -11,8 +11,9 @@ module ice_shelf_tabular_calving use ice_bergs_framework, only : add_new_berg_to_list use ice_bergs_framework, only : spread_variable_across_cells, sum_up_spread_fields use ice_bergs_framework, only : hexagon_into_quadrants_using_triangles, Rearth +use ice_bergs_framework, only : square_into_quadrants_using_triangles use ice_bergs_framework, only : initialize_iceberg_bonds, count_bonds -use ice_bergs_framework, only : find_cell_wide, pos_within_cell, generate_id +use ice_bergs_framework, only : find_cell, pos_within_cell, generate_id use ice_bergs_framework, only : debug, footloose, connect_all_bonds, delete_all_bonds use ice_bergs_framework, only : update_halo_calved_tabular_icebergs, assign_n_bonds,transfer_mts_bergs use fms_mod, only : error_mesg, FATAL, WARNING, stderr @@ -552,12 +553,11 @@ subroutine ice_shelf_to_bonded_bergs(bergs, TC) real :: minlon, maxlon, minlat, maxlat real :: minlon0, maxlon0, minlat0, maxlat0 real :: minx, miny, cos_lat_ref, dlonscale - integer :: stderrunit, i, k, nbonds, berg_count + integer :: stderrunit, i, k, nbonds logical :: check_bond_quality !The number of tabular icebergs to initialize on the current PE bcount = TC%berg_pe_count - berg_count = 0 !Note: account for the calving mask to be between 0 and 1 !Initialize bergs over all cells with mask>0. Eliminate a berg if its groundfrac is greater than some threshold @@ -612,10 +612,10 @@ subroutine ice_shelf_to_bonded_bergs(bergs, TC) endif !grid bounds - maxlon0=maxval(grd%lon(grd%isc-1:grd%iec, grd%jsc-1:grd%jec)) - maxlat0=maxval(grd%lat(grd%isc-1:grd%iec, grd%jsc-1:grd%jec)) - minlon0=minval(grd%lon(grd%isc-1:grd%iec, grd%jsc-1:grd%jec)) - minlat0=minval(grd%lat(grd%isc-1:grd%iec, grd%jsc-1:grd%jec)) + maxlon0=maxval(grd%lon( grd%iec, grd%jsc-1:grd%jec)) + maxlat0=maxval(grd%lat(grd%isc-1:grd%iec, grd%jec)) + minlon0=minval(grd%lon(grd%isc-1 , grd%jsc-1:grd%jec)) + minlat0=minval(grd%lat(grd%isc-1:grd%iec, grd%jsc-1 )) !initialize the particles for each new tabular iceberg do i = 1,bcount @@ -681,11 +681,13 @@ subroutine ice_shelf_to_bonded_bergs(bergs, TC) lat = minlat k=0 - berg_count=0 - do while (lat<=maxlat) + !Start calving bergs on the current PE, keeping in mind that bergs that fall on the eastern and northern + !boundary of the PE should actually be included in the PEs to the east or north, respectively, instead of + !the current PE + do while (lat<=maxlat .and. latminlat0) then + if (lat>=minlat0) then if (grd%grid_is_latlon) dlon = dlonscale/cos(lat*(pi/180.)) @@ -695,9 +697,8 @@ subroutine ice_shelf_to_bonded_bergs(bergs, TC) lon=minlon endif - do while (lon<=maxlon) - if (lon>minlon0) then - berg_count=berg_count+1 + do while (lon<=maxlon .and. lon=minlon0) then !Calve particles that overlap the mask. Save the overlapping area of !each particle with neighboring cells. Their thickness and scaling !will be determined below in new_tabular_bergs_thickness_and_pressure. @@ -1195,8 +1196,9 @@ subroutine begin_calving_tabular_iceberg_from_shelf(bergs, grd, lon, lat, calve_ ! allocations_done=.false. - lres=find_cell_wide(grd, lon, lat, i, j) + lres=find_cell(grd, lon, lat, i, j) + if (.not. lres) return !Assume that the calve mask spans over the cells that shelf associated with the berg may advect into !Get rid of bergs that clearly do not overlap the calve mask if (all(frac_shelf(i-1:i+1,j-1:j+1)==0) .and. all(calve_mask(i-1:i+1,j-1:j+1)==0)) return @@ -1210,6 +1212,10 @@ subroutine begin_calving_tabular_iceberg_from_shelf(bergs, grd, lon, lat, calve_ write(stderrunit,*) 'KID, calve_icebergs: something went very wrong!',i,j,xi,yj call error_mesg('KID, calve_icebergs', 'berg xi,yj is not correct!', FATAL) endif + + !Ignore bergs on the N and E boundry of the PE, as they will be included in the PEs to the N or E, respectively + if ((i==grd%iec .and. xi==1) .or. (j==grd%jec .and. yj==1)) return + ! if (grd%msk(i,j)<0.5) then ! write(stderrunit,*) 'KID, calve_icebergs: WARNING!!! Iceberg born in land cell',i,j,newberg%lon,newberg%lat ! if (debug) call error_mesg('KID, calve_icebergs', 'Iceberg born in Land Cell!', FATAL) diff --git a/src/icebergs.F90 b/src/icebergs.F90 index 99ab056..2f0f8fb 100644 --- a/src/icebergs.F90 +++ b/src/icebergs.F90 @@ -55,6 +55,7 @@ module ice_bergs use ice_bergs_framework, only: break_bonds_on_sub_steps, initialize_iceberg_bonds use ice_bergs_framework, only: short_step_mts_grounding, radius_based_drag use ice_bergs_framework, only: hexagon_into_quadrants_using_triangles +use ice_bergs_framework, only: square_into_quadrants_using_triangles use ice_bergs_framework, only: sum_up_spread_fields, sum_up_spread_fields, Area_of_triangle use ice_bergs_framework, only: point_in_triangle, point_in_interval, point_is_on_the_line use ice_bergs_framework, only: convert_from_grid_to_meters, convert_from_meters_to_grid diff --git a/src/icebergs_framework.F90 b/src/icebergs_framework.F90 index ba01803..17e414d 100644 --- a/src/icebergs_framework.F90 +++ b/src/icebergs_framework.F90 @@ -108,6 +108,7 @@ module ice_bergs_framework public set_constant_interaction_length_and_width public break_bonds_on_sub_steps, skip_first_outer_mts_step, no_frac_first_ts public hexagon_into_quadrants_using_triangles, Area_of_triangle, point_in_triangle +public square_into_quadrants_using_triangles public point_in_interval, point_is_on_the_line public ij_component_of_id, spread_variable_across_cells, sum_up_spread_fields public initialize_iceberg_bonds, find_orientation_using_iceberg_bonds @@ -2570,7 +2571,8 @@ subroutine mts_pack_in_dir(bergs, nbergs_to_send, dir) integer :: ics,ice,jcs,jce integer :: inhs,inhe,jnhs,jnhe real :: pfix !0.0 .and. grd%lon(grd%iec,grd%jec).eq.grd%maxlon_c) pfix=grd%Lx !lat/lon range of halo cells - rhc(1)=grd%lon(is-1,js); rhc(2)=grd%lat(is-1,js) !min corner - rhc(3)=grd%lon(grd%iec,je); rhc(4)=grd%lat(grd%iec,je) !max corner + ! rhc(1)=grd%lon(is-1,js); rhc(2)=grd%lat(is-1,js) !min corner + ! rhc(3)=grd%lon(grd%iec,je); rhc(4)=grd%lat(grd%iec,je) !max corner + rhc(1)=is; rhc(2)=js+1 + rhc(3)=grd%iec; rhc(4)=je !range for non-halo cells (for contact with conglom) inhs=grd%isd; inhe=is-1; jnhs=js; jnhe=je !range for cells w/in contact dist of edge of comp domain: @@ -2602,8 +2606,10 @@ subroutine mts_pack_in_dir(bergs, nbergs_to_send, dir) case ("w") is=grd%isd; ie=grd%isc+halo_width-1; js=grd%jsd; je=grd%jed if (grd%Lx>0.0 .and. grd%lon(grd%isc-1,grd%jsc-1).eq.grd%minlon_c) pfix=-grd%Lx - rhc(1)=grd%lon(grd%isc-1,js); rhc(2)=grd%lat(grd%isc-1,js) - rhc(3)=grd%lon(ie,je); rhc(4)=grd%lat(ie,je) + ! rhc(1)=grd%lon(grd%isc-1,js); rhc(2)=grd%lat(grd%isc-1,js) + ! rhc(3)=grd%lon(ie,je); rhc(4)=grd%lat(ie,je) + rhc(1)=grd%isc; rhc(2)=js+1 + rhc(3)=ie; rhc(4)=je inhs=ie+1; inhe=grd%ied; jnhs=js; jnhe=je ice=min(grd%isc+nc_x,grd%ied) if (ice>ie) then @@ -2612,8 +2618,10 @@ subroutine mts_pack_in_dir(bergs, nbergs_to_send, dir) endif case ("n") is=grd%isd; ie=grd%ied; js=grd%jec-halo_width+2; je=grd%jed - rhc(1)=grd%lon(is,js-1); rhc(2)=grd%lat(is,js-1) - rhc(3)=grd%lon(ie,grd%jec); rhc(4)=grd%lat(ie,grd%jec) + ! rhc(1)=grd%lon(is,js-1); rhc(2)=grd%lat(is,js-1) + ! rhc(3)=grd%lon(ie,grd%jec); rhc(4)=grd%lat(ie,grd%jec) + rhc(1)=is+1; rhc(2)=js + rhc(3)=ie; rhc(4)=grd%jec inhs=is; inhe=ie; jnhs=grd%jsd; jnhe=js-1 jcs=max(grd%jec-nc_y,grd%jsd+1) if (jcsje) then @@ -2729,7 +2739,8 @@ recursive subroutine mts_mark_and_pack_halo_and_congloms(bergs, berg, dir, nberg character(len=1) :: dir !< north,south,east,or west integer :: nbergs_to_send !< Number of bergs to send real :: pfix !rhc(3).or.berg%lat>rhc(4)) then + ! if (berg%lonrhc(3).or.berg%lat>rhc(4)) then + if (berg%inerhc(3).or.berg%jne>rhc(4)) then berg%halo_berg=max(2.,current_halo_id) !outside neighboring PE halo else berg%halo_berg=1 !inside neighboring PE halo @@ -3934,12 +3946,12 @@ subroutine unpack_berg_from_buffer2(bergs, buff, n, grd, force_append, max_bonds if (bergs%mts .and. abs(localberg%halo_berg)>1) then temp_lon = localberg%lon; temp_lat = localberg%lat - if (temp_lon > grd%lon(grd%ied,grd%jed)) then + if (temp_lon >= grd%lon(grd%ied,grd%jed)) then temp_lon = grd%lonc(grd%ied,grd%jed) elseif (temp_lon < grd%lon(grd%isd,grd%jsd)) then temp_lon = grd%lonc(grd%isd+1,grd%jsd+1) endif - if (temp_lat > grd%lat(grd%ied,grd%jed)) then + if (temp_lat >= grd%lat(grd%ied,grd%jed)) then temp_lat = grd%latc(grd%ied,grd%jed) elseif (temp_lat < grd%lat(grd%isd,grd%jsd)) then temp_lat = grd%latc(grd%isd+1,grd%jsd+1) @@ -3999,7 +4011,10 @@ subroutine unpack_berg_from_buffer2(bergs, buff, n, grd, force_append, max_bonds else write(stderrunit,'("KID, unpack_berg_from_buffer pe=(",i3,a,2i4,a,2f8.2)')& & mpp_pe(),') Failed to find i,j=',localberg%ine,localberg%jne,' for lon,lat=',localberg%lon,localberg%lat - write(stderrunit,*) localberg%halo_berg +! if (localberg%id==4294972731) then + lres=find_cell_wide(grd, localberg%lon, localberg%lat, localberg%ine, localberg%jne, explain=.true.) +! endif + write(stderrunit,*) localberg%id, localberg%halo_berg write(stderrunit,*) localberg%lon,localberg%lat write(stderrunit,*) localberg%uvel,localberg%vvel write(stderrunit,*) localberg%axn,localberg%ayn !Alon @@ -6489,13 +6504,14 @@ logical function find_cell(grd, x, y, oi, oj) end function find_cell !> Scans each all grid cells until is_point_in_cell() is true (includes halos) -logical function find_cell_wide(grd, x, y, oi, oj) +logical function find_cell_wide(grd, x, y, oi, oj, explain) ! Arguments type(icebergs_gridded), intent(in) :: grd !< Container for gridded fields real, intent(in) :: x !< Longitude of position real, intent(in) :: y !< Latitude of position integer, intent(out) :: oi !< i-index of cell containing position or -999 integer, intent(out) :: oj !< j-index of cell containing position or -999 +logical, intent(in), optional :: explain !< If true, print debugging ! Local variables integer :: i,j @@ -6504,7 +6520,7 @@ logical function find_cell_wide(grd, x, y, oi, oj) oi=floor((x-grd%lon(grd%isd,grd%jsd))/(grd%lon(grd%isd+1,grd%jsd+1)-grd%lon(grd%isd,grd%jsd)))+grd%isd+1 oj=floor((y-grd%lat(grd%isd,grd%jsd))/(grd%lat(grd%isd+1,grd%jsd+1)-grd%lat(grd%isd,grd%jsd)))+grd%jsd+1 if (.not. (oi-1.lt.grd%isd.or.oi.gt.grd%ied.or.oj-1.lt.grd%jsd.or.oj.gt.grd%jed)) then - if (is_point_in_cell(grd, x, y, oi, oj)) then + if (is_point_in_cell(grd, x, y, oi, oj, explain=explain)) then find_cell_wide=.true.; return endif endif @@ -6512,7 +6528,7 @@ logical function find_cell_wide(grd, x, y, oi, oj) oi=-999; oj=-999 do j=grd%jsd+1,grd%jed; do i=grd%isd+1,grd%ied - if (is_point_in_cell(grd, x, y, i, j)) then + if (is_point_in_cell(grd, x, y, i, j, explain=explain)) then oi=i; oj=j; find_cell_wide=.true. return endif @@ -6639,15 +6655,15 @@ logical function sum_sign_dot_prod4(x0, y0, x1, y1, x2, y2, x3, y3, x, y, Lx, ex xx3= apply_modulo_around_point(x3,x0,Lx) - l0=(xx-xx0)*(y1-y0)-(y-y0)*(xx1-xx0) - l1=(xx-xx1)*(y2-y1)-(y-y1)*(xx2-xx1) - l2=(xx-xx2)*(y3-y2)-(y-y2)*(xx3-xx2) - l3=(xx-xx3)*(y0-y3)-(y-y3)*(xx0-xx3) + l0=(xx-xx0)*(y1-y0)-(y-y0)*(xx1-xx0) !S segment + l1=(xx-xx1)*(y2-y1)-(y-y1)*(xx2-xx1) !E segment + l2=(xx-xx2)*(y3-y2)-(y-y2)*(xx3-xx2) !N segment + l3=(xx-xx3)*(y0-y3)-(y-y3)*(xx0-xx3) !W segment - !We use an asymmetry between South and East line boundaries and North and East + !We use an asymmetry between South and West line boundaries and North and East !to avoid icebergs appearing to two cells (half values used for debugging) - !This is intended to make the South and East boundaries be part of the - !cell, while the North and West are not part of the cell. + !This is intended to make the South and West boundaries be part of the + !cell, while the North and East are not part of the cell. p0=sign(1., l0); if (l0.eq.0.) p0=-0.5 p1=sign(1., l1); if (l1.eq.0.) p1=0.5 p2=sign(1., l2); if (l2.eq.0.) p2=0.5 @@ -6983,7 +6999,7 @@ end subroutine calc_xiyj !> Returns true if non-dimensional position xi,yj is in unit interval !! -!! Includes South and East boundaries, and excludes North and West. +!! Includes South and West boundaries, and excludes North and East. !! \todo Double check definition of is_point_within_xi_yj_bounds() logical function is_point_within_xi_yj_bounds(xi,yj) ! Arguments @@ -8590,7 +8606,7 @@ subroutine Square_into_quadrants_using_triangles(x0, y0, L, theta, Area_square , C1x= L2 ; C1y= L2 !Corner 1 (top right) C2x=-L2 ; C2y= L2 !Corner 2 (top left) C3x=-L2 ; C3y=-L2 !Corner 3 (bottom left) - C1x= L2 ; C1y=-L2 !Corner 4 (bottom right) + C4x= L2 ; C4y=-L2 !Corner 4 (bottom right) ! Finding positions of corners call rotate_and_translate(C1x,C1y,theta,x0,y0)