From 6d4d4d4b05a18cfc867a148dc57cfa30872e67fc Mon Sep 17 00:00:00 2001 From: alex-huth Date: Tue, 27 Dec 2022 13:52:58 -0500 Subject: [PATCH] Added bond variable 'broken' to track if a bond is broken (1) or intact (0). Added namelist option 'use_broken_bonds_for_substep_contact', which if true, identifies a pair of particles as contacting during the sub-steps if and only if the bond between them is broken (according to the new 'broken' variable'). This is more computationally-efficient than searching for contact pairs in the current and neighboring grid cells. Both the new 'broken' variable and 'use_broken_bonds_for_substep_contact' option can only be used for iKID multiple time step runs in dem-mode (dem=.true.). --- src/icebergs.F90 | 99 +++++++++------- src/icebergs_framework.F90 | 167 +++++++++++++++++++-------- src/icebergs_io.F90 | 27 +++-- tests/collision_tests/.gitignore | 3 +- tests/collision_tests/input.nml | 161 ++++++++++++++++++++++++++ tests/dem_ground_frac_test/README | 9 +- tests/dem_ground_frac_test/input.nml | 17 ++- tests/footloose_tests/.gitignore | 3 +- 8 files changed, 379 insertions(+), 107 deletions(-) create mode 100644 tests/collision_tests/input.nml diff --git a/src/icebergs.F90 b/src/icebergs.F90 index 5910afb..d881b89 100644 --- a/src/icebergs.F90 +++ b/src/icebergs.F90 @@ -1174,6 +1174,17 @@ subroutine calculate_force_dem(bergs, berg, other_berg, current_bond, & if (current_bond%nstress>bergs%frac_thres_n .or. current_bond%sstress>bergs%frac_thres_t) then bergs%bond_break_detected=.true. + if (current_bond%broken.ne.1) then + current_bond%broken=1 + if (current_bond%nstress>bergs%frac_thres_n .and. current_bond%sstress>bergs%frac_thres_t) then + print *,'T and S break',berg%lat,current_bond%nstress,current_bond%sstress + elseif (current_bond%nstress>bergs%frac_thres_n) then + print *,'TENSILE break',berg%lat,current_bond%nstress,current_bond%sstress + else + print *,'SHEAR break',berg%lat,current_bond%nstress,current_bond%sstress + endif + endif + if (current_bond%nstress<0) then !if fail under shear with normal compression, then should still feel contact force damping_coef = bergs%dem_damping_coef*sqrt(bergs%dem_K_damp*M1*M2/(M1+M2)) @@ -1181,6 +1192,8 @@ subroutine calculate_force_dem(bergs, berg, other_berg, current_bond, & F_x = F_x + Fn_x; F_y = F_y + Fn_y !linear forces endif + if (bergs%use_broken_bonds_for_substep_contact) berg%n_bonds=berg%n_bonds-1 + if (savestress .and. save_bond_forces) then current_bond%T=0.; current_bond%T_d=0. if (current_bond%nstress<0) then @@ -1201,29 +1214,17 @@ subroutine calculate_force_dem(bergs, berg, other_berg, current_bond, & current_bond%other_bond%tangd1 = -current_bond%tangd1 current_bond%other_bond%tangd2 = -current_bond%tangd2 current_bond%other_bond%length = current_bond%length + current_bond%other_bond%broken=1 current_bond%other_bond%other_id=-current_bond%other_bond%other_id + if (bergs%use_broken_bonds_for_substep_contact) then + current_bond%other_berg%n_bonds=current_bond%other_berg%n_bonds-1 + endif endif return !Broken bond. Return without adding forces and torques endif - elseif (bergs%fracture_criterion=='computed') then - Rsig=current_bond%nstress/bergs%frac_thres_n - Rtau=current_bond%sstress/bergs%frac_thres_t - if ((0.5*Rsig+sqrt((0.5*Rsig)**2. + Rtau**2.)) > 1) then - bergs%bond_break_detected=.true. - if (save_bond_forces) current_bond%other_bond%other_id=-current_bond%other_bond%other_id - return - endif - elseif (bergs%fracture_criterion=='computed_s') then - Rsig=current_bond%nstress/bergs%frac_thres_n - Rtau=current_bond%sstress/bergs%frac_thres_t - if (sqrt(Rsig**2. + Rtau**2.) > 1) then - if (save_bond_forces) current_bond%other_bond%other_id=-current_bond%other_bond%other_id - bergs%bond_break_detected=.true. - return - endif else call error_mesg('KID, calculate_force_dem', & - 'break_bond_on_sub_steps currently only possible with fracture critera=computed or stress', FATAL) + 'break_bond_on_sub_steps currently only possible with fracture critera=stress', FATAL) endif endif @@ -1611,7 +1612,11 @@ subroutine accel_mts(bergs, berg, i, j, xi, yj, lat, uvel, vvel, uvel0, vvel0, d ! Determining number of bonds current_bond=>berg%first_bond do while (associated(current_bond)) ! loop over all bonds - N_bonds=N_bonds+1.0 + if (bergs%dem) then + if (current_bond%broken.ne.1) N_bonds=N_bonds+1.0 + else + N_bonds=N_bonds+1.0 + endif current_bond=>current_bond%next_bond enddo dragfrac = ((N_max-N_bonds)/N_max) @@ -2039,7 +2044,9 @@ subroutine accel_explicit_inner_mts(bergs, berg, i, j, xi, yj, lat, uvel, vvel, if (current_bond%other_id<0) then current_bond%other_id=-current_bond%other_id other_berg=>current_bond%other_berg - other_berg%id=-other_berg%id !mark, as to not repeat in contact search below + if (.not. bergs%use_broken_bonds_for_substep_contact) then + other_berg%id=-other_berg%id !mark, as to not repeat in contact search below + endif bond_matched=.true. F_x=F_x+current_bond%F_x; F_y=F_y+current_bond%F_y @@ -2055,15 +2062,22 @@ subroutine accel_explicit_inner_mts(bergs, berg, i, j, xi, yj, lat, uvel, vvel, call error_mesg('KID,bond interactions', 'Trying to do Bond interactions with unassosiated berg!' ,FATAL) else if (bergs%dem) then - call calculate_force_dem(bergs,berg,other_berg,current_bond,& - dt,F_x,F_y,T,Fd_x,Fd_y,T_d,savestress=.true.) - other_berg%id=-other_berg%id !mark, as to not repeat in contact search below + if (current_bond%broken.eq.1) then + call calculate_unbonded_same_conglom_dem_force(bergs, berg, other_berg, & + IA_x, IA_y, IAd_x, IAd_y, uvel0, vvel0, uvel0, vvel0) + else + call calculate_force_dem(bergs,berg,other_berg,current_bond,& + dt,F_x,F_y,T,Fd_x,Fd_y,T_d,savestress=.true.) + endif + if (.not. bergs%use_broken_bonds_for_substep_contact) then + other_berg%id=-other_berg%id !mark, as to not repeat in contact search below + endif else P_ia_11=0. ; P_ia_12=0. ; P_ia_21=0.; P_ia_22=0. P_ia_times_u_x=0. ; P_ia_times_u_y=0. call calculate_force(bergs, berg, other_berg, IA_x, IA_y, uvel0, vvel0, uvel0, vvel0, & P_ia_11, P_ia_12, P_ia_21, P_ia_22, P_ia_times_u_x, P_ia_times_u_y,bonded) - other_berg%id=-other_berg%id !mark, as to not repeat in contact search below + other_berg%id=-other_berg%id !mark, as to not repeat in contact search below (dem required for use_broken_bonds_for_substep_contact) IAd_x=IAd_x+P_ia_11*(other_berg%uvel_old-berg%uvel_old)+P_ia_12*(other_berg%vvel_old-berg%vvel_old) IAd_y=IAd_y+P_ia_12*(other_berg%uvel_old-berg%uvel_old)+P_ia_22*(other_berg%vvel_old-berg%vvel_old) endif @@ -2077,7 +2091,7 @@ subroutine accel_explicit_inner_mts(bergs, berg, i, j, xi, yj, lat, uvel, vvel, !for dem-mode, do not bother running the contact algorithm if the particle is within interior of conglomerate if (bergs%iceberg_bonds_on .and. bergs%dem) then - if (berg%n_bonds .eq. bergs%max_bonds) then + if ((berg%n_bonds .eq. bergs%max_bonds).or. bergs%use_broken_bonds_for_substep_contact) then run_contact=.false. else run_contact=.true. @@ -2107,14 +2121,15 @@ subroutine accel_explicit_inner_mts(bergs, berg, i, j, xi, yj, lat, uvel, vvel, enddo enddo;enddo endif - - !unmark the bonded bergs - current_bond=>berg%first_bond - do while (associated(current_bond)) - other_berg=>current_bond%other_berg - other_berg%id=abs(other_berg%id) - current_bond=>current_bond%next_bond - enddo + if (.not. bergs%use_broken_bonds_for_substep_contact) then + !unmark the bonded bergs + current_bond=>berg%first_bond + do while (associated(current_bond)) + other_berg=>current_bond%other_berg + other_berg%id=abs(other_berg%id) + current_bond=>current_bond%next_bond + enddo + endif endif if (bergs%dem) then @@ -2511,7 +2526,11 @@ subroutine accel(bergs, berg, i, j, xi, yj, lat, uvel, vvel, uvel0, vvel0, dt, r ! Determining number of bonds current_bond=>berg%first_bond do while (associated(current_bond)) ! loop over all bonds - N_bonds=N_bonds+1.0 + if (bergs%dem) then + if (current_bond%broken.ne.1) N_bonds=N_bonds+1.0 + else + N_bonds=N_bonds+1.0 + endif current_bond=>current_bond%next_bond enddo dragfrac = ((N_max-N_bonds)/N_max) @@ -3397,7 +3416,11 @@ subroutine thermodynamics(bergs) ! Determining number of bonds ! current_bond=>this%first_bond ! do while (associated(current_bond)) ! loop over all bonds + ! if (bergs%dem) then + ! if (current_bond%broken.ne.1) N_bonds=N_bonds+1.0 + ! else ! N_bonds=N_bonds+1.0 + ! endif ! current_bond=>current_bond%next_bond ! enddo endif @@ -7343,8 +7366,6 @@ subroutine evolve_icebergs_mts(bergs) skip_first_outer_mts_step=.false. endif - - ! PART 3: MTS fast sub-steps (bonded interactions only): ! Position and velocity are updated by: ! X2 = X1+dt*V1+((dt^2)/2)*a_k +((dt^2)/2)*b_k @@ -7435,10 +7456,6 @@ subroutine evolve_icebergs_mts(bergs) i=berg%ine ; j=berg%jne xi=berg%xi ; yj=berg%yj - ! if (latn>50.e3 .or. latn<10.e3 .or. lonn > 45.e3) then - ! print *,'id,lat,lon',berg%id,latn,lonn - ! endif - if (new_mts) then axn = axn+bxn; ayn=ayn+byn endif @@ -7637,7 +7654,9 @@ subroutine evolve_icebergs_mts(bergs) enddo ! loop over all bergs enddo; enddo ! update 'old' velocities - if (break_bonds_on_sub_steps .and. bergs%bond_break_detected) call break_bonds_dem(bergs) + if (.not. bergs%use_broken_bonds_for_substep_contact) then + if (break_bonds_on_sub_steps .and. bergs%bond_break_detected) call break_bonds_dem(bergs) + endif enddo ! loop over sub-steps !reset interactive_forces diff --git a/src/icebergs_framework.F90 b/src/icebergs_framework.F90 index d4ce54a..3b3618a 100644 --- a/src/icebergs_framework.F90 +++ b/src/icebergs_framework.F90 @@ -422,6 +422,7 @@ module ice_bergs_framework real, allocatable :: tangd2 !< Accumulated tangential displacement, y-component (m) real, allocatable :: nstress !< Normal stress (Pa) real, allocatable :: sstress !< Shear stress (Pa) + integer, allocatable :: broken !< Is a bond broken (1) or not (0)? type(bond), pointer :: other_bond=>null() !< Matching bond ! If save_bond_forces real, allocatable :: F_x !< Zonal force @@ -463,6 +464,7 @@ module ice_bergs_framework real, allocatable :: tangd2 !< Accumulated tangential displacement, y-component (m) real, allocatable :: nstress real, allocatable :: sstress + integer, allocatable :: broken !< Is bond broken (1) or not(0)? end type bond_xyt ! A dynamic buffer, used for communication, that packs types into rectangular memory @@ -648,6 +650,7 @@ module ice_bergs_framework real :: poisson=0.3 ! Poisson's ratio real :: dem_spring_coef=0. real :: dem_damping_coef=0.1 + logical :: use_broken_bonds_for_substep_contact=.false. !only evaluate sub-step contact between particles with a broken bond logical :: bond_break_detected=.false. logical :: dem_shear_for_frac_only=.false. !< If true, DEM shear is calculated for fracture, but zeroed for berg interactions integer :: dem_beam_test=0 !1=Simply supported beam,2=cantilever beam,3=angular vel tes @@ -768,7 +771,6 @@ subroutine ice_bergs_framework_init(bergs, & logical :: uniaxial_test=.false. !adds a tensile stress to the east-most berg element real :: cdrag_grounding=0.0 ! Drag coefficient against ocean bottom real :: h_to_init_grounding=100.0 -!character(len=11) :: fracture_criterion ! 'energy','stress','strain_rate','strain',or 'none' real :: frac_thres_n=0.0 !normal fracture strain threshold real :: frac_thres_t=0.0 !tangential fracture strain threshold logical :: damage_test_1=.false. ! Sets initial damage for bonds that overlap y=17000.0 @@ -876,6 +878,7 @@ subroutine ice_bergs_framework_init(bergs, & real :: poisson=0.3 ! Poisson's ratio real :: dem_spring_coef=0. real :: dem_damping_coef=0.1 +logical :: use_broken_bonds_for_substep_contact=.false. ! only evaluate sub-step contact between particles with a broken bond logical :: dem_shear_for_frac_only=.false. ! If true, DEM shear is calculated for fracture, but zeroed for berg interactions integer :: dem_beam_test=0 !1=Simply supported beam,2=cantilever beam,3=angular vel test ! Element Interactions @@ -933,7 +936,7 @@ subroutine ice_bergs_framework_init(bergs, & fl_l_scale_erosion_only, fl_youngs, fl_strength, save_all_traj_year, save_nonfl_traj_by_class,& save_traj_by_class_start_mass_thres_n, save_traj_by_class_start_mass_thres_s,traj_area_thres_sntbc,& traj_area_thres_fl,tau_is_velocity, ocean_drag_scale, A68_test, & - A68_xdisp,A68_ydisp,& + A68_xdisp,A68_ydisp,use_broken_bonds_for_substep_contact,& orig_dem_moment_of_inertia, break_bonds_on_sub_steps, skip_first_outer_mts_step, rev_mind, & no_frac_first_ts, use_grounding_torque, power_ground, short_step_mts_grounding, radius_based_drag, save_bond_forces @@ -1388,9 +1391,9 @@ subroutine ice_bergs_framework_init(bergs, & endif if (dem) then - buffer_width=buffer_width+3+(max_bonds*5) + buffer_width=buffer_width+3+(max_bonds*6) buffer_width_traj=buffer_width_traj+3 - buffer_width_bond_traj=buffer_width_bond_traj+4+1 + buffer_width_bond_traj=buffer_width_bond_traj+6 elseif (fracture_criterion .ne. 'none') then buffer_width=buffer_width+1+(max_bonds*3) buffer_width_traj=buffer_width_traj+1 @@ -1561,6 +1564,16 @@ subroutine ice_bergs_framework_init(bergs, & bergs%dem_spring_coef=dem_spring_coef bergs%dem_K_damp=2.*bergs%dem_spring_coef/(3.*(1.-bergs%poisson**2)) !damping factor bergs%dem_damping_coef=dem_damping_coef + if (break_bonds_on_sub_steps) then + if (use_broken_bonds_for_substep_contact .and. .not. (bergs%dem .and. bergs%iceberg_bonds_on)) then + call error_mesg('KID, framework', & + 'use_broken_bonds_for_substep_contact requires dem and iceberg_bonds_on', FATAL) + else + bergs%use_broken_bonds_for_substep_contact=use_broken_bonds_for_substep_contact + endif + else + bergs%use_broken_bonds_for_substep_contact=.false. + end if bergs%dem_beam_test=dem_beam_test bergs%constant_interaction_LW=constant_interaction_LW bergs%constant_length=constant_length @@ -2772,6 +2785,18 @@ subroutine set_conglom_ids(bergs) enddo enddo;enddo + if (bergs%use_broken_bonds_for_substep_contact) then + do grdj = grd%jsd,grd%jed; do grdi = grd%isd,grd%ied + this=>bergs%list(grdi,grdj)%first + do while (associated(this)) + if (this%n_bondsthis%next + enddo + enddo;enddo + endif + end subroutine set_conglom_ids !> Assign initial velocities and energies for energy tracking tests @@ -2819,19 +2844,84 @@ recursive subroutine label_conglomerates(berg, new_conglom_id) type(iceberg), pointer :: other_berg type(bond) , pointer :: current_bond + current_bond=>berg%first_bond + + if (dem) then + do while (associated(current_bond)) + if (current_bond%broken.ne.1) then + if (associated(current_bond%other_berg)) then + other_berg=>current_bond%other_berg + if (other_berg%conglom_id.ne.new_conglom_id) then + other_berg%conglom_id=new_conglom_id + if (other_berg%halo_berg==10) other_berg%halo_berg=2 + call label_conglomerates(other_berg,new_conglom_id) + endif + endif + endif + current_bond=>current_bond%next_bond + enddo + else + !same as above with no broken bond check + do while (associated(current_bond)) + if (associated(current_bond%other_berg)) then + other_berg=>current_bond%other_berg + if (other_berg%conglom_id.ne.new_conglom_id) then + other_berg%conglom_id=new_conglom_id + if (other_berg%halo_berg==10) other_berg%halo_berg=2 + call label_conglomerates(other_berg,new_conglom_id) + endif + endif + current_bond=>current_bond%next_bond + enddo + endif +end subroutine label_conglomerates + +!> Remove broken bonds between two different conglomerates +subroutine remove_broken_bonds_between_congloms(berg) + type(iceberg), pointer :: berg !< Berg to start search from + ! Local variables + type(iceberg), pointer :: other_berg + type(bond) , pointer :: current_bond, other_bond, kick_the_bucket + current_bond=>berg%first_bond do while (associated(current_bond)) - if (associated(current_bond%other_berg)) then - other_berg=>current_bond%other_berg - if (other_berg%conglom_id.ne.new_conglom_id) then - other_berg%conglom_id=new_conglom_id - if (other_berg%halo_berg==10) other_berg%halo_berg=2 - call label_conglomerates(other_berg,new_conglom_id) + + if (current_bond%broken.eq.1) then + if (associated(current_bond%other_berg)) then + other_berg=>current_bond%other_berg + if (other_berg%conglom_id.ne.berg%conglom_id) then + !broken bond exists between two different conglomerates + !remove matching bond from other berg + if (save_bond_forces .and. associated(current_bond%other_bond)) then + other_bond=>current_bond%other_bond + call delete_bond_from_list(other_berg,other_bond) + else + other_bond=>other_berg%first_bond + do while (associated(other_bond)) ! loop over all bonds + if (other_bond%other_id.eq.berg%id) then + kick_the_bucket=>other_bond + other_bond=>null() + call delete_bond_from_list(other_berg,kick_the_bucket) + else + other_bond=>other_bond%next_bond + endif + enddo + endif + !remove bond from current berg + kick_the_bucket=>current_bond + current_bond=>current_bond%next_bond + call delete_bond_from_list(berg,kick_the_bucket) + else + current_bond=>current_bond%next_bond + endif + else + current_bond=>current_bond%next_bond endif + else + current_bond=>current_bond%next_bond endif - current_bond=>current_bond%next_bond enddo -end subroutine label_conglomerates +end subroutine remove_broken_bonds_between_congloms !> Keep only the bergs which are part of, or within contact distance of, a conglomerate that overlaps !! the computational domain, or which are halo_berg=1 @@ -3480,6 +3570,7 @@ subroutine pack_berg_into_buffer2(berg, buff, n, max_bonds_in) call push_buffer_value(buff%data(:,n), counter, current_bond%nstress) call push_buffer_value(buff%data(:,n), counter, current_bond%sstress) call push_buffer_value(buff%data(:,n), counter, current_bond%rel_rotation) + call push_buffer_value(buff%data(:,n), counter, current_bond%broken) elseif (fracture_criterion .ne. 'none') then call push_buffer_value(buff%data(:,n), counter, current_bond%rotation) call push_buffer_value(buff%data(:,n), counter, current_bond%rel_rotation) @@ -3519,6 +3610,7 @@ subroutine pack_berg_into_buffer2(berg, buff, n, max_bonds_in) call push_buffer_value(buff%data(:,n), counter, 0) call push_buffer_value(buff%data(:,n), counter, 0) call push_buffer_value(buff%data(:,n), counter, 0) + call push_buffer_value(buff%data(:,n), counter, 0) elseif (fracture_criterion .ne. 'none') then call push_buffer_value(buff%data(:,n), counter, 0) call push_buffer_value(buff%data(:,n), counter, 0) @@ -3901,6 +3993,7 @@ subroutine unpack_berg_from_buffer2(bergs, buff, n, grd, force_append, max_bonds call pull_buffer_value(buff%data(:,n), counter, current_bond%nstress) call pull_buffer_value(buff%data(:,n), counter, current_bond%sstress) call pull_buffer_value(buff%data(:,n), counter, current_bond%rel_rotation) + call pull_buffer_value(buff%data(:,n), counter, current_bond%broken) elseif (fracture_criterion .ne. 'none') then call pull_buffer_value(buff%data(:,n), counter, current_bond%rotation) call pull_buffer_value(buff%data(:,n), counter, current_bond%rel_rotation) @@ -4232,6 +4325,7 @@ subroutine pack_bond_traj_into_buffer2(bond_traj, buff, n) call push_buffer_value(buff%data(:,n),counter,bond_traj%nstress) call push_buffer_value(buff%data(:,n),counter,bond_traj%sstress) call push_buffer_value(buff%data(:,n),counter,bond_traj%rel_rotation) + call push_buffer_value(buff%data(:,n),counter,bond_traj%broken) elseif (fracture_criterion.ne.'none') then call push_buffer_value(buff%data(:,n),counter,bond_traj%rotation) call push_buffer_value(buff%data(:,n),counter,bond_traj%rel_rotation) @@ -4307,6 +4401,7 @@ subroutine unpack_bond_traj_from_buffer2(first, buff, n) call pull_buffer_value(buff%data(:,n),counter,bond_traj%nstress) call pull_buffer_value(buff%data(:,n),counter,bond_traj%sstress) call pull_buffer_value(buff%data(:,n),counter,bond_traj%rel_rotation) + call pull_buffer_value(buff%data(:,n),counter,bond_traj%broken) elseif (fracture_criterion.ne.'none') then call pull_buffer_value(buff%data(:,n),counter,bond_traj%rotation) call pull_buffer_value(buff%data(:,n),counter,bond_traj%rel_rotation) @@ -5413,29 +5508,13 @@ subroutine break_bonds_dem(bergs) if (bergs%fracture_criterion=='stress') then if (current_bond%nstress>frac_thres_n .or. current_bond%sstress>frac_thres_t) then - if (current_bond%nstress>frac_thres_n .and. current_bond%sstress>frac_thres_t) then - print *,'T and S break',this%lat,current_bond%nstress,current_bond%sstress - elseif (current_bond%nstress>frac_thres_n) then - print *,'TENSILE break',this%lat,current_bond%nstress,current_bond%sstress - else - print *,'SHEAR break',this%lat,current_bond%nstress,current_bond%sstress - endif - current_bond%other_id=-1 - endif - elseif (bergs%fracture_criterion=='computed') then - Rsig=current_bond%nstress/frac_thres_n - Rtau=current_bond%sstress/frac_thres_t - frac=0.5*Rsig+sqrt((0.5*Rsig)**2. + Rtau**2.) - if (frac > 1.) then - print *,'BREAK',this%lat,frac,Rsig,Rtau - current_bond%other_id=-1 - endif - elseif (bergs%fracture_criterion=='computed_s') then - Rsig=current_bond%nstress/bergs%frac_thres_n - Rtau=current_bond%sstress/bergs%frac_thres_t - frac=sqrt(Rsig**2. + Rtau**2.) - if (frac > 1.) then - print *,'BREAK',this%lat,frac,Rsig,Rtau + ! if (current_bond%nstress>frac_thres_n .and. current_bond%sstress>frac_thres_t) then + ! print *,'T and S break',this%lat,current_bond%nstress,current_bond%sstress + ! elseif (current_bond%nstress>frac_thres_n) then + ! print *,'TENSILE break',this%lat,current_bond%nstress,current_bond%sstress + ! else + ! print *,'SHEAR break',this%lat,current_bond%nstress,current_bond%sstress + ! endif current_bond%other_id=-1 endif else @@ -5454,21 +5533,7 @@ subroutine break_bonds_dem(bergs) if (other_bond%other_id.eq.this%id) then other_bond%other_id=-1 if (debug) then - if (bergs%fracture_criterion=='computed') then - Rsig=other_bond%nstress/frac_thres_n - Rtau=other_bond%sstress/frac_thres_t - frac=0.5*Rsig+sqrt((0.5*Rsig)**2. + Rtau**2.) - if (frac <= 1.) then - print *,'Other bond did not fracture!' - endif - elseif(bergs%fracture_criterion=='computed_s') then - Rsig=current_bond%nstress/bergs%frac_thres_n - Rtau=current_bond%sstress/bergs%frac_thres_t - frac=sqrt(Rsig**2. + Rtau**2.) - if (frac <= 1.) then - print *,'Other bond did not fracture!' - endif - elseif (other_bond%nstress<=frac_thres_n .and. other_bond%sstress<=frac_thres_t) then + if (other_bond%nstress<=frac_thres_n .and. other_bond%sstress<=frac_thres_t) then print *,'Other bond did not fracture!' print *,'current and other bond nstress,',current_bond%nstress, other_bond%nstress print *,'current and other bond sstress,',current_bond%sstress, other_bond%sstress @@ -5565,6 +5630,8 @@ subroutine form_a_bond(berg, other_id, other_berg_ine, other_berg_jne, other_ber new_bond%nstress=0.; new_bond%sstress=0. allocate(new_bond%rel_rotation) new_bond%rel_rotation=0. + allocate(new_bond%broken) + new_bond%broken=0 if (save_bond_forces) then allocate(new_bond%F_x, new_bond%F_y, new_bond%Fd_x, new_bond%Fd_y, new_bond%T, new_bond%T_d) new_bond%F_x=0.; new_bond%F_y=0.; new_bond%Fd_x=0.; new_bond%Fd_y=0.; new_bond%T=0.; new_bond%T_d=0. diff --git a/src/icebergs_io.F90 b/src/icebergs_io.F90 index d7761af..0e48343 100644 --- a/src/icebergs_io.F90 +++ b/src/icebergs_io.F90 @@ -185,7 +185,8 @@ subroutine write_restart(bergs, time_stamp) first_berg_jne, & first_berg_ine, & other_berg_jne, & - other_berg_ine + other_berg_ine, & + broken integer :: grdi, grdj @@ -481,6 +482,7 @@ subroutine write_restart(bergs, time_stamp) allocate(nstress(nbonds)) allocate(sstress(nbonds)) allocate(rel_rotation(nbonds)) + allocate(broken(nbonds)) elseif (fracture_criterion.ne.'none') then allocate(rotation(nbonds)) allocate(rel_rotation(nbonds)) @@ -525,6 +527,8 @@ subroutine write_restart(bergs, time_stamp) longname='shear stress',units='Pa') id = register_restart_field(bergs_bond_restart,filename_bonds,'rel_rotation',rel_rotation,& longname='relative rotation',units='rad') + id = register_restart_field(bergs_bond_restart,filename_bonds,'broken',broken,& + longname='broken status',units='none') elseif (fracture_criterion.ne.'none') then id = register_restart_field(bergs_bond_restart,filename_bonds,'rotation',rotation,& longname='rotation',units='rad') @@ -564,6 +568,7 @@ subroutine write_restart(bergs, time_stamp) nstress(i)=current_bond%nstress sstress(i)=current_bond%sstress rel_rotation(i)=current_bond%rel_rotation + broken(i)=current_bond%broken elseif (fracture_criterion.ne.'none') then rotation(i)=current_bond%rotation rel_rotation(i)=current_bond%rel_rotation @@ -599,7 +604,8 @@ subroutine write_restart(bergs, time_stamp) tangd2, & nstress, & sstress, & - rel_rotation) + rel_rotation, & + broken) elseif (fracture_criterion.ne.'none') then deallocate( & rotation, & @@ -1278,7 +1284,8 @@ subroutine read_restart_bonds(bergs,Time) first_berg_jne, & first_berg_ine, & other_berg_jne, & - other_berg_ine + other_berg_ine, & + broken real, allocatable, dimension(:) :: tangd1, & tangd2, & nstress, & @@ -1332,6 +1339,7 @@ subroutine read_restart_bonds(bergs,Time) allocate(nstress(nbonds_in_file)) allocate(sstress(nbonds_in_file)) allocate(rel_rotation(nbonds_in_file)) + allocate(broken(nbonds_in_file)) elseif (fracture_criterion.ne.'none') then allocate(rotation(nbonds_in_file)) allocate(rel_rotation(nbonds_in_file)) @@ -1365,6 +1373,7 @@ subroutine read_restart_bonds(bergs,Time) call read_real_vector(filename,'nstress',nstress,grd%domain,value_if_not_in_file=0.) call read_real_vector(filename,'sstress',sstress,grd%domain,value_if_not_in_file=0.) call read_real_vector(filename,'rel_rotation',rel_rotation,grd%domain,value_if_not_in_file=0.) + call read_int_vector(filename,'broken',broken,grd%domain,value_if_not_in_file=0) elseif (fracture_criterion.ne.'none') then call read_real_vector(filename,'rotation',rotation,grd%domain,value_if_not_in_file=0.) call read_real_vector(filename,'rel_rotation',rel_rotation,grd%domain,value_if_not_in_file=0.) @@ -1477,6 +1486,7 @@ subroutine read_restart_bonds(bergs,Time) current_bond%nstress=nstress(k) current_bond%sstress=sstress(k) current_bond%rel_rotation=rel_rotation(k) + current_bond%broken=broken(k) elseif (fracture_criterion.ne.'none') then current_bond%rotation=rotation(k) current_bond%rel_rotation=rel_rotation(k) @@ -1550,7 +1560,8 @@ subroutine read_restart_bonds(bergs,Time) tangd2, & nstress, & sstress, & - rel_rotation) + rel_rotation, & + broken) elseif (fracture_criterion.ne.'none') then deallocate( & rotation, & @@ -2273,7 +2284,7 @@ subroutine write_bond_trajectory(trajectory, bond_traj_name) ! Local variables integer :: iret, ncid, i_dim, i integer :: lonid, latid, yearid, dayid, lenid,n1id, n2id, peid -integer :: rotid,rrotid,nsid,nsrid, damid +integer :: rotid,rrotid,nsid,nsrid, damid, brid integer :: idcnt1_id, idcnt2_id, idij1_id, idij2_id, eeid, edid integer :: axid,ayid,bxid,byid integer :: td1id,td2id,dnsid,dssid @@ -2405,7 +2416,7 @@ subroutine write_bond_trajectory(trajectory, bond_traj_name) if (dem) then td1id = inq_varid(ncid, 'tangd1'); td2id = inq_varid(ncid, 'tangd2') dnsid = inq_varid(ncid, 'nstress'); dssid = inq_varid(ncid, 'sstress') - rrotid = inq_varid(ncid, 'rel_rotation') + rrotid = inq_varid(ncid, 'rel_rotation'); brid = inq_varid(ncid, 'broken') elseif (fracture_criterion.ne.'none') then rotid = inq_varid(ncid, 'rotation'); rrotid = inq_varid(ncid, 'rel_rotation') nsid = inq_varid(ncid, 'n_frac_var') @@ -2444,6 +2455,7 @@ subroutine write_bond_trajectory(trajectory, bond_traj_name) td1id = def_var(ncid, 'tangd1', NF_DOUBLE, i_dim); td2id = def_var(ncid, 'tangd2', NF_DOUBLE, i_dim) dnsid = def_var(ncid, 'nstress', NF_DOUBLE, i_dim); dssid = def_var(ncid, 'sstress', NF_DOUBLE, i_dim) rrotid = def_var(ncid, 'rel_rotation', NF_DOUBLE, i_dim) + brid = def_var(ncid, 'broken', NF_INT, i_dim) elseif (fracture_criterion.ne.'none') then rotid = def_var(ncid, 'rotation', NF_DOUBLE, i_dim); rrotid = def_var(ncid, 'rel_rotation', NF_DOUBLE, i_dim) nsid = def_var(ncid, 'n_frac_var', NF_DOUBLE, i_dim) @@ -2500,6 +2512,7 @@ subroutine write_bond_trajectory(trajectory, bond_traj_name) call put_att(ncid, dnsid, 'long_name', 'normal stress'); call put_att(ncid, dnsid, 'units', 'Pa') call put_att(ncid, dssid, 'long_name', 'shear stress'); call put_att(ncid, dssid, 'units', 'Pa') call put_att(ncid, rrotid, 'long_name', 'relative rotation'); call put_att(ncid, rrotid, 'units', 'rad') + call put_att(ncid, brid, 'long_name', 'broken status'); call put_att(ncid, brid, 'units', 'none') elseif (fracture_criterion.ne.'none') then call put_att(ncid, rotid, 'long_name', 'rotation'); call put_att(ncid, rotid, 'units', 'rad') call put_att(ncid, rrotid, 'long_name', 'relative rotation'); call put_att(ncid, rrotid, 'units', 'rad') @@ -2548,7 +2561,7 @@ subroutine write_bond_trajectory(trajectory, bond_traj_name) if (dem) then call put_double(ncid,td1id,i,this%tangd1); call put_double(ncid,td2id,i,this%tangd2) call put_double(ncid,dnsid,i,this%nstress); call put_double(ncid,dssid,i,this%sstress) - call put_double(ncid,rrotid,i,this%rel_rotation) + call put_double(ncid,rrotid,i,this%rel_rotation); call put_int(ncid,brid,i,this%broken) elseif (fracture_criterion.ne.'none') then call put_double(ncid,rotid,i,this%rotation); call put_double(ncid,rrotid,i,this%rel_rotation) call put_double(ncid,nsid,i,this%n_frac_var) diff --git a/tests/collision_tests/.gitignore b/tests/collision_tests/.gitignore index 68bcbc9..75188aa 100644 --- a/tests/collision_tests/.gitignore +++ b/tests/collision_tests/.gitignore @@ -1 +1,2 @@ -results/ \ No newline at end of file +results/ +input.nml \ No newline at end of file diff --git a/tests/collision_tests/input.nml b/tests/collision_tests/input.nml new file mode 100644 index 0000000..328c2cd --- /dev/null +++ b/tests/collision_tests/input.nml @@ -0,0 +1,161 @@ +!!KID, bergs_chksum: write_restart berg chksum = -1211104969 chksum2= -1211104969 chksum3= 1964715299 chksum4= 1964715299 chksum5= 0#= 16 + +&icebergs_driver_nml + ni=20 !number of elements, x direction + nj=20 !number of elements, y direction + ibdt=60.0 !seconds + ibuo=0.2 !ocean x velocity + ibvo=0.2 !ocean y velocity + ibui=0.0 !ice x velocity + ibvi=0.0 !ice y velocity + ibhrs=48 !total simulation hours + nmax=99999999 !max number of timesteps + saverestart=.true. + halo=3 + debug=.false. + collision_test=.true. !specific to the current test, must be true here +/ + + +&icebergs_nml + + traj_name=KID.nc + bond_traj_name=KID_bonds.nc + + !defaults given in parentheses: + ! explicit_inner_mts=.true. + ! new_mts=.true. + ! mts=.true. !(F) for using the Multiple Time Stepping Velocity Verlet scheme + ! mts_sub_steps=60 !1 !30 !(-1) # of mts sub-steps (-1=automatically determine # from spring const) + ! contact_distance=1.75e3 !2.e3 !(0) unbonded bergs contact at max(contact_dist,sum of their radii) + ! contact_spring_coef=1.e-7!(default=spring_coeff) Berg spring coeff for collisions + ! debug_write=.false. !(F)sets traj_sample_hours=traj_write_hours & includes halo bergs in write + ! remove_unused_bergs=.false. !(T) remove unneeded bergs after PEs transfers + + halo=3 !(4) halo width + Lx=20000. ! (360.) Length of X-domain. For periodicity (use a huge number for non-periodic) + grid_is_latlon=.false. !(T) 'T'=lat/lon grid, uses earth radius to convert to distance + grid_is_regular=.True. !(T) Point in cell can be found assuming regular Cartesian grid + hexagonal_icebergs=.true. !(F)'T'=bergs are rectangles,'F'=bergs are hex (for mass spread) + Static_icebergs=.False. ! (T) True=icebergs do no move + + rho_bergs=850. ! (850) iceberg density + spring_coef=1.e-5 !(1.e-8) Spring const for berg interactions (dflt~highest stable value) + bond_coef=1.e-8 ! (1.e-8) Spring const for berg bonds - not being used right now + radial_damping_coef=1.e-4 ! (1.e-4) Coeff for rel berg motion damping (radial component) + tangental_damping_coef=2.e-5 !(2.e-5) Coeff for rel berg motion damping (tangential comp.) + critical_interaction_damping_on=.true. !(T) Sets damping on rel berg vel to crit val + LoW_ratio=1.5 ! (1.5) Initial ratio L/W for newly calved icebergs + bergy_bit_erosion_fraction=0. ! (0.) Fraction of erosion melt flux to divert to bergy bits + sicn_shift=0. ! (0) Shift of sea-ice concent. in erosion flux modulation (0