diff --git a/driver/icebergs_driver.F90 b/driver/icebergs_driver.F90 index 42815b4..e112545 100644 --- a/driver/icebergs_driver.F90 +++ b/driver/icebergs_driver.F90 @@ -18,6 +18,7 @@ program icebergs_driver use mpp_mod, only : mpp_pe use mpp_mod, only : mpp_root_pe use mpp_mod, only: mpp_sync +use mpp_mod, only: mpp_exit use time_manager_mod, only : time_type use time_manager_mod, only : set_date use time_manager_mod, only : set_calendar_type @@ -60,6 +61,7 @@ program icebergs_driver logical :: chaotic_test=.false. logical :: grounding_test=.false. logical :: big_grounding_test=.false. +logical :: dem_test_contact=.false. logical :: a68_test=.false. logical :: vert_int_ocean_vel=.false. logical :: tau_is_velocity=.true. @@ -87,7 +89,7 @@ program icebergs_driver saverestart,ibui,ibvi,collision_test,chaotic_test,grounding_test,& gridres,write_time_inc,bump_depth, sst, a68_test, data_dir, vert_int_ocean_vel,& reverse_a68_forcings,tau_is_velocity,transient_a68_data_start_ind,& - Old_wrong_Rearth_and_SSH,no_wind,REarth,big_grounding_test !mom_Rearth + Old_wrong_Rearth_and_SSH,no_wind,REarth,big_grounding_test,dem_test_contact !mom_Rearth ! For loops integer :: isc !< Start of i-index for computational domain (used for loops) integer :: iec !< End of i-index for computational domain (used for loops) @@ -334,6 +336,32 @@ program icebergs_driver enddo; enddo print *,'min, max depth',minval(depth(:,:)),maxval(depth(:,:)) depth = 1000.0-depth + endif + + if (dem_test_contact) then + lat(:,:)=lat(:,:)-0.45 !-25000.1 + lon(:,:)=lon(:,:)-0.45 + !assign land cells at the N and S portions of the domain. + !input.nml: set coastal drift parameter to prevent grounding + do j=jsd,jed; do i=isd,ied + if (lat(i,j)<=(-5.e3) .or. lat(i,j)>=(220.e3)) wet(i,j)=0.0 + enddo;enddo + + mid=50.e3 + + do j=jsd,jed; do i=isd,ied + if (lon(i,j)>65 .or. lat(i,j)==mid) then + vo(i,j)=0.0 + else + if (lat(i,j)>mid) then + vo(i,j)=-ibvo + else + vo(i,j)=ibvo + endif + endif + enddo;enddo + + endif if (grounding_test) then @@ -513,13 +541,15 @@ program icebergs_driver if (saverestart) call icebergs_save_restart(bergs) -!Deallocate all memory and disassociated pointer -call icebergs_end(bergs) - call cpu_time(time_finish) if (mpp_pe()==0) then write(*,*) 'Simulation finished in ', (time_finish - time_begin)/60., 'minutes' end if +call mpp_exit() + +!Deallocate all memory and disassociated pointer +call icebergs_end(bergs) + end program icebergs_driver diff --git a/src/icebergs.F90 b/src/icebergs.F90 index fedc1b2..f13a2eb 100644 --- a/src/icebergs.F90 +++ b/src/icebergs.F90 @@ -1176,6 +1176,7 @@ subroutine calculate_force_dem(bergs, berg, other_berg, current_bond, & 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)) Fd_x = Fd_x-damping_coef*ur; Fd_y = Fd_y-damping_coef*vr !linear damping force F_x = F_x + Fn_x; F_y = F_y + Fn_y !linear forces endif @@ -2089,7 +2090,7 @@ subroutine accel_explicit_inner_mts(bergs, berg, i, j, xi, yj, lat, uvel, vvel, do grdi = max(berg%ine-1,bergs%grd%isd+1),min(berg%ine+1,bergs%grd%ied) other_berg=>bergs%list(grdi,grdj)%first do while (associated(other_berg)) - if (other_berg%id>0 .and. other_berg%conglom_id.eq.berg%conglom_id) then + if (other_berg%id>0 .and. other_berg%conglom_id.eq.berg%conglom_id .and. other_berg%n_bondsthis%first_bond - do while (associated(current_bond)) ! loop over all bonds - N_bonds=N_bonds+1.0 - current_bond=>current_bond%next_bond - enddo + ! current_bond=>this%first_bond + ! do while (associated(current_bond)) ! loop over all bonds + ! N_bonds=N_bonds+1.0 + ! current_bond=>current_bond%next_bond + ! enddo endif if (this%static_berg .eq. 1) N_bonds=N_max !Static icebergs melt like ice shelves endif @@ -7303,35 +7305,40 @@ subroutine evolve_icebergs_mts(bergs) do while (associated(berg)) ! loop over all bergs !only evolve non-static bergs that overlap, or are part of a conglom that overlaps, the computational domain: if (berg%static_berg .lt. 0.5 .and. (berg%conglom_id.ne.0 .or. bergs%force_convergence)) then - latn = berg%lat ; lonn = berg%lon - uvel1=berg%uvel ; vvel1=berg%vvel !V_n (previous cycle) - axn = 0.0 ; ayn = 0.0 !note these are redefined at the start of subroutine accel_mts - bxn = 0.0 ; byn = 0.0 - i=berg%ine ; j=berg%jne - xi=berg%xi ; yj=berg%yj - - call accel_mts(bergs, berg, i, j, xi, yj, latn, uvel1, vvel1, uvel1, vvel1, dt, rx, ry, & - ax1, ay1, axn, ayn, bxn, byn, save_bond_energy,Fec_x, Fec_y, Fdc_x, Fdc_y) - - if (Fdc_x .ne. 0. .or. Fdc_y .ne. 0.) had_collision=.true. - if (monitor_energy .and. last_iter) & - call mts_energy_part_1(bergs,berg,ax1,ay1,axn,ayn,bxn,byn,Fec_x,Fec_y,Fdc_x,Fdc_y) - - ! Saving all the iceberg variables - berg%axn=axn; berg%ayn=ayn - berg%bxn=bxn; berg%byn=byn - - if (bergs%force_convergence) then - berg%uvel_prev=berg%uvel+(dt*ax1); berg%vvel_prev=berg%vvel+(dt*ay1) !the new velocity - if (ii==1) usum=usum+berg%uvel_old**2 + berg%vvel_old**2 - usum1=usum1+berg%uvel_prev**2+berg%vvel_prev**2 - usum2=usum2+(berg%uvel_prev-berg%uvel_old)**2+(berg%vvel_prev-berg%vvel_old)**2 - else - berg%uvel=berg%uvel+(dt*ax1); berg%vvel=berg%vvel+(dt*ay1) + if (ii==1 .or. berg%static_berg==0.1) then + latn = berg%lat ; lonn = berg%lon + uvel1=berg%uvel ; vvel1=berg%vvel !V_n (previous cycle) + axn = 0.0 ; ayn = 0.0 !note these are redefined at the start of subroutine accel_mts + bxn = 0.0 ; byn = 0.0 + i=berg%ine ; j=berg%jne + xi=berg%xi ; yj=berg%yj + + call accel_mts(bergs, berg, i, j, xi, yj, latn, uvel1, vvel1, uvel1, vvel1, dt, rx, ry, & + ax1, ay1, axn, ayn, bxn, byn, save_bond_energy,Fec_x, Fec_y, Fdc_x, Fdc_y) + + if (Fdc_x .ne. 0. .or. Fdc_y .ne. 0. .and. bergs%force_convergence) then + had_collision=.true. + berg%static_berg=0.1 + endif + if (monitor_energy .and. last_iter) & + call mts_energy_part_1(bergs,berg,ax1,ay1,axn,ayn,bxn,byn,Fec_x,Fec_y,Fdc_x,Fdc_y) + + ! Saving all the iceberg variables + berg%axn=axn; berg%ayn=ayn + berg%bxn=bxn; berg%byn=byn + + if (bergs%force_convergence) then + berg%uvel_prev=berg%uvel+(dt*ax1); berg%vvel_prev=berg%vvel+(dt*ay1) !the new velocity + if (ii==1) usum=usum+berg%uvel_old**2 + berg%vvel_old**2 + usum1=usum1+berg%uvel_prev**2+berg%vvel_prev**2 + usum2=usum2+(berg%uvel_prev-berg%uvel_old)**2+(berg%vvel_prev-berg%vvel_old)**2 + else + berg%uvel=berg%uvel+(dt*ax1); berg%vvel=berg%vvel+(dt*ay1) - !The final velocity from the previous cycle, which can be used - !during post-processing to calculate kinetic energy and momentum - berg%uvel_prev=berg%uvel ; berg%vvel_prev=berg%vvel + !The final velocity from the previous cycle, which can be used + !during post-processing to calculate kinetic energy and momentum + berg%uvel_prev=berg%uvel ; berg%vvel_prev=berg%vvel + endif endif endif berg=>berg%next @@ -7345,6 +7352,7 @@ subroutine evolve_icebergs_mts(bergs) do while (associated(berg)) ! loop over all bergs if (berg%static_berg .lt. 0.5) then berg%uvel_old=berg%uvel_prev; berg%vvel_old=berg%vvel_prev + if (last_iter .and. berg%static_berg==0.1) berg%static_berg=0. endif berg=>berg%next enddo diff --git a/tests/a68/long_run.nml b/tests/a68/long_run.nml index 507561d..e85f80a 100644 --- a/tests/a68/long_run.nml +++ b/tests/a68/long_run.nml @@ -26,7 +26,7 @@ &icebergs_nml Rearth= - !6363827 !6.378e6 !6.371229e6 !6.36e6 !6.378e6 + A68_test=.true. A68_xdisp= A68_ydisp= @@ -40,9 +40,9 @@ constant_interaction_LW=.true. - traj_name=_1p.nc - bond_traj_name=_bonds_1p.nc - + traj_name=.nc + bond_traj_name=_bonds.nc + break_bonds_on_sub_steps=.true. power_ground=.false. @@ -121,7 +121,7 @@ coastal_drift=0. !(0.)Vel added to currents to cause bergs to drift away from land cells tidal_drift=0. !(0.)Amp of stoch. tidal vel added to currents to cause random berg drift Runge_not_Verlet=.false. !(T) True=Runge Kutta, False=Verlet. - use_mixed_melting=.false.!(F) 'T' partly determ. melt frm 3eq model(uses BergParams&Bond #) + use_mixed_melting=.true. !.false.!(F) 'T' partly determ. melt frm 3eq model(uses BergParams&Bond #) apply_thickness_cutoff_to_gridded_melt=.true.!(F)Prevnt melt for Ho < melt_cuttoff on grid apply_thickness_cutoff_to_bergs_melt=.true. !(F)Prevnt melt for Ho < melt_cuttoff on bergs melt_cutoff=10.0 ! (-1.0) Min ocean thickness for melting to occur (no apply for vals < 0) @@ -153,7 +153,7 @@ verbose_hrs=1!24000 ! (24) Period between verbose messages traj_sample_hrs=0!1 !(24) sampling period for trajectory storage traj_write_hrs=0!6 !(480) Period for writing sampled trajectories to disk - save_bond_traj=.true. !.false. !(F) If T, saves bond trajectories in a separate file + save_bond_traj=.false. !(F) If T, saves bond trajectories in a separate file budget=.false. !(T) Calculate budgets debug=.false. !(F) Turn on debugging really_debug=.false. !(F) Turn on debugging @@ -211,7 +211,7 @@ / &fms_nml - clock_grain='NONE' + clock_grain='MODULE' domains_stack_size = 2000000 clock_flags='SYNC' / \ No newline at end of file diff --git a/tests/a68/input.nml b/tests/a68/long_run_no_inc_contact_dist.nml similarity index 89% rename from tests/a68/input.nml rename to tests/a68/long_run_no_inc_contact_dist.nml index f216136..12de105 100644 --- a/tests/a68/input.nml +++ b/tests/a68/long_run_no_inc_contact_dist.nml @@ -4,11 +4,11 @@ data_dir='/lustre/f2/dev/gfdl/Alexander.Huth/MOM6-examples/src/icebergs/tests/a68/data/' Old_wrong_Rearth_and_SSH=.false. !mom_Rearth=.true. - Rearth=6363827 + Rearth= !6363827 !6.378e6 no_wind=.false. tau_is_velocity=.true.!.false. - ibdt=1800.0 !seconds + ibdt= !3600.0 !1800.0 !seconds ibuo=0.0 !ocean x velocity ibvo=0.0 !ocean y velocity ibui=0.0 !ice x velocity @@ -25,24 +25,24 @@ &icebergs_nml - Rearth=6363827 - !6363827 !6.378e6 !6.371229e6 !6.36e6 !6.378e6 + Rearth= + A68_test=.true. - A68_xdisp=-37.51 - A68_ydisp=-55.2166 + A68_xdisp= + A68_ydisp= short_step_mts_grounding=.true. - tau_is_velocity=.true.!.false. - ocean_drag_scale=17.8 !20. !1 + tau_is_velocity=.true. + ocean_drag_scale= skip_first_outer_mts_step=.true. add_curl_to_torque=.false. constant_interaction_LW=.true. - traj_name=longrun_od_17.8_tn_ng_ydn.1_ss90_ns18_rx-37.51_ry-55.2166_1p.nc - bond_traj_name=longrun_od_17.8_tn_ng_ydn.1_ss90_ns18_rx-37.51_ry-55.2166_bonds_1p.nc - + traj_name=.nc + bond_traj_name=_bonds.nc + break_bonds_on_sub_steps=.true. power_ground=.false. @@ -68,27 +68,23 @@ !defaults given in parentheses: new_mts=.true. mts=.true. !(F) for using the Multiple Time Stepping Velocity Verlet scheme - mts_sub_steps=90 !200 !75 !75 !150 !50 !(-1) # of mts sub-steps + mts_sub_steps= !200 !75 !75 !150 !50 !(-1) # of mts sub-steps !mts_sub_steps=10 force_convergence=.true. !(F) experimental MTS convergence scheme to better conserve momentum during collisions convergence_tolerance=1e-4 !8 !(1.e-8) tolerance for the MTS force_convergence scheme - contact_distance=4.e3 !(0) unbonded bergs contact at max(contact_dist,sum of their radii) + contact_distance=0.!4.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 - cdrag_grounding=1.e4 !0.!1.e6!3.16e6!1.e6!1.e3!3.16e6 !for power-law-type friction + cdrag_grounding= !0.!1.e6!3.16e6!1.e6!1.e3!3.16e6 !for power-law-type friction h_to_init_grounding=0!30!15!50. !15. ewsame=.false. !(F) set T if periodic and 2 PEs along the x (zonal) direction remove_unused_bergs=.true. !(T) remove unneeded bergs after PEs transfers fracture_criterion='stress' !('none') 'strain_rate','strain',or 'none' - frac_thres_scaling = 1 !0000000 - frac_thres_n = 18e3 !25000.0 !5000.0 !4250. !12000.0 !8000.0 !5575. !3000. !2.1e4 !12. - frac_thres_t = 100e3 !1000. !2000.0 !2500. !6000.0 !4000.0 !5575. !1750. !35000. !12. - - !!fracture_criterion='strain' - !!frac_thres_n = 0.0162 - !!frac_thres_t = 0.0457 + frac_thres_scaling = 1 + frac_thres_n = e3 + frac_thres_t = 1000e3 halo=3 !(4) halo width Lx=360 !28.875!45000. ! (360.) Length of X-domain. For periodicity (use a huge number for non-periodic) @@ -125,7 +121,7 @@ coastal_drift=0. !(0.)Vel added to currents to cause bergs to drift away from land cells tidal_drift=0. !(0.)Amp of stoch. tidal vel added to currents to cause random berg drift Runge_not_Verlet=.false. !(T) True=Runge Kutta, False=Verlet. - use_mixed_melting=.false.!(F) 'T' partly determ. melt frm 3eq model(uses BergParams&Bond #) + use_mixed_melting=.true. !.false.!(F) 'T' partly determ. melt frm 3eq model(uses BergParams&Bond #) apply_thickness_cutoff_to_gridded_melt=.true.!(F)Prevnt melt for Ho < melt_cuttoff on grid apply_thickness_cutoff_to_bergs_melt=.true. !(F)Prevnt melt for Ho < melt_cuttoff on bergs melt_cutoff=10.0 ! (-1.0) Min ocean thickness for melting to occur (no apply for vals < 0) @@ -155,9 +151,9 @@ verbose=.false. !(F) Be verbose to stderr verbose_hrs=1!24000 ! (24) Period between verbose messages - traj_sample_hrs=0!1 !(24) sampling period for trajectory storage - traj_write_hrs=0!6 !(480) Period for writing sampled trajectories to disk - save_bond_traj=.true. !.false. !(F) If T, saves bond trajectories in a separate file + traj_sample_hrs=1 !0!1 !(24) sampling period for trajectory storage + traj_write_hrs=1 !0!6 !(480) Period for writing sampled trajectories to disk + save_bond_traj=.false. !(F) If T, saves bond trajectories in a separate file budget=.false. !(T) Calculate budgets debug=.false. !(F) Turn on debugging really_debug=.false. !(F) Turn on debugging @@ -215,7 +211,7 @@ / &fms_nml - clock_grain='NONE' + clock_grain='MODULE' domains_stack_size = 2000000 clock_flags='SYNC' / \ No newline at end of file diff --git a/tests/a68/long_run_no_inc_contact_dist.sh b/tests/a68/long_run_no_inc_contact_dist.sh new file mode 100755 index 0000000..9d02863 --- /dev/null +++ b/tests/a68/long_run_no_inc_contact_dist.sh @@ -0,0 +1,40 @@ +#!/bin/bash + +Rearth=6363827 +odrag=(17.8) +ss=90 +name=(lrun_nicd_od_17.8_tn_ng_ydn.1_ss90_ns18_rx-37.51_ry-55.2166) +ibdt=1800.0 + +gcoef=(1.e4) +nstress=(18) +A68dx=-37.51 +A68dy=-55.2166 + +odrag=(10.35) +ss=180 +name=(lrun_nicd_hrtest) +ibdt=3600.0 +nstress=(24.) #(22.5) +A68dx=-37.5 +A68dy=-55.1 + +for ((i = 0; i < ${#odrag[@]}; ++i)); do + drag=(${odrag[$i]}) + fname=(${name[$i]}) + gc=(${gcoef[$i]}) + ns=(${nstress[$i]}) + + echo $drag + echo $fname + echo $Rearth + echo $ss + echo $A68dx + echo $A68dy + echo $ns + + sed "s//$drag/; s//$Rearth/; s//$ss/; s//$ns/; s//$gc/; s//$A68dx/; s//$A68dy/; s//$ibdt/; s//$fname/g" \ + long_run_no_inc_contact_dist.nml > input.nml + ./RUN1p + # &>/dev/null & +done diff --git a/tests/chaotic_test/.gitignore b/tests/chaotic_test/.gitignore index 68bcbc9..f395a63 100644 --- a/tests/chaotic_test/.gitignore +++ b/tests/chaotic_test/.gitignore @@ -1 +1,3 @@ -results/ \ No newline at end of file +results/ +csvresults/ +csvresults2/ \ No newline at end of file diff --git a/tests/dem_test_contact/.gitignore b/tests/dem_test_contact/.gitignore new file mode 100644 index 0000000..68bcbc9 --- /dev/null +++ b/tests/dem_test_contact/.gitignore @@ -0,0 +1 @@ +results/ \ No newline at end of file diff --git a/tests/dem_test_contact/README b/tests/dem_test_contact/README new file mode 100644 index 0000000..ec575a8 --- /dev/null +++ b/tests/dem_test_contact/README @@ -0,0 +1,7 @@ +1. Make the iceberg input files: cd makeberg; ./RUN; cd .. +2.(optional) edit input.nml +3. run the model: ./RUN +4, animate the results: animate_trajectories.py + +This is a good test for visualizing fracture dependence on particle size. Try running makeberg.py with the full radius and then half the radius, and with regular vs dem nmls. Make sure to edit the nmls so that contact distance is halved when radius is halved + diff --git a/tests/dem_test_contact/RUN b/tests/dem_test_contact/RUN new file mode 100755 index 0000000..3dfc97d --- /dev/null +++ b/tests/dem_test_contact/RUN @@ -0,0 +1,7 @@ +rm iceberg_trajectories.nc +rm iceberg_trajectories.nc.* +rm bond_trajectories.nc +rm bond_trajectories.nc.* +srun -n1 ./../../build/bergs.x + + diff --git a/tests/dem_test_contact/animate_trajectories.py b/tests/dem_test_contact/animate_trajectories.py new file mode 100755 index 0000000..d163b57 --- /dev/null +++ b/tests/dem_test_contact/animate_trajectories.py @@ -0,0 +1,161 @@ +#!/usr/bin/env python + +from netCDF4 import Dataset +import pdb +import netCDF4 as nc +import numpy as np +import matplotlib +import math +import os +matplotlib.use("tkagg") +from pylab import * +import matplotlib.pyplot as plt +from matplotlib.animation import FuncAnimation +import argparse + +def parseCommandLine(): + parser = argparse.ArgumentParser(description= + '''Generate animation of iceberg trajectories.''', + epilog='Written by Alex Huth, 2020') + parser.add_argument('-fname', type=str, default='iceberg_trajectories.nc', + help=''' provide filename to plot''') + parser.add_argument('-s', type=int, default='14000', + help='''plotted particle size''') + optCmdLineArgs = parser.parse_args() + return optCmdLineArgs + + + +def main(args): + + print('reading file') + + filename=args.fname + psize=args.s + #filename = 'iceberg_trajectories.nc' + field = 'lon' + field2 = 'lat' + + with nc.Dataset(filename) as file: + x = file.variables['lon'][:]/1.e3 + y = file.variables['lat'][:]/1.e3 + day = file.variables['day'][:] + length = file.variables['length'][:] + width = file.variables['width'][:] + thick = file.variables['thickness'][:] + od = file.variables['od'][:] + + ud = np.unique(day) + t = ud[0] + + radius = length*width*(1./(2*np.sqrt(3))) + + #for determining if grounded: + rho_bergs=850. + rho_seawater=1025. + h_to_ground=20 #200. + draught=(rho_bergs/rho_seawater)*thick + if (h_to_ground==0): + groundfrac=draught-od + else: + groundfrac=1.-(od-draught)/h_to_ground + groundfrac[groundfrac<0.]=0. + groundfrac[groundfrac>1.]=1. + + #groundfrac[thick<200.]=0. + #groundfrac[thick>=200.]=1. + + + # frame info + num_frames = len(ud) + movie_len = 5 #seconds + #frame_len = 1000.0*movie_len/num_frames/10000000 + frame_len = movie_len/num_frames/1000 + + xmin = 40 + xmax = 100 + ymin = 20 #xmin + ymax = 80 #xmax + + anim_running = True + + def animate(i): + + x1 = x[day == ud[i]] + y1 = y[day == ud[i]] + data = np.hstack((x1[:,np.newaxis],y1[:,np.newaxis])) + + scat.set_offsets(data) + + t = ud[i] + time_text.set_text('time = %.1f days' % t ) + + r1 = radius[day == ud[i]] + scat.set_sizes(r1/psize) + + hstat = groundfrac[day == ud[i]] + cstring=[] + for j in range(len(hstat)): + if (hstat[j]==0): + cstring.append("none") + else: + cstring.append("b") + scat.set_color(cstring) + scat.set_edgecolor('r') + + return scat,time_text + + def init(): + scat.set_offsets([]) + return scat, + + def onClick(event): + global anim_running + if anim_running: + ani.event_source.stop() + anim_running = False + else: + ani.event_source.start() + anim_running = True + + # Now we can do the plotting! + f = plt.figure(figsize=(5,5)) + f.tight_layout() + ax1 = plt.subplot(111,xlim=(xmin, xmax), ylim=(ymin, ymax)) + scat = ax1.scatter([],[],marker='o')#,edgecolor='red') + time_text = ax1.text(0.02, 0.95, '', transform=ax1.transAxes) + + # Change major ticks to show every 20. + ax1.xaxis.set_major_locator(MultipleLocator(5)) + ax1.yaxis.set_major_locator(MultipleLocator(5)) + + ax1.xaxis.set_minor_locator(MultipleLocator(1)) + ax1.yaxis.set_minor_locator(MultipleLocator(1)) + + # Turn grid on for both major and minor ticks and style minor slightly + # differently. + ax1.grid(which='major', color='#CCCCCC', linestyle=':') + ax1.grid(which='minor', color='#CCCCCC', linestyle=':') + + ax1.set_xlabel('x (km)') + ax1.set_ylabel('y (km)') + ax1.set_title('Iceberg trajectory') + + f.canvas.mpl_connect('button_press_event', onClick) + + ani = FuncAnimation( + f,animate,init_func=init,frames=num_frames, + interval=frame_len,blit=True,repeat=True) + + print('time',t) + + plt.show() + #animation.save("iceberg_traj_animation.mp4") + + +print('Script complete') + +if __name__ == '__main__': + optCmdLineArgs= parseCommandLine() + anim_running = True + main(optCmdLineArgs) diff --git a/tests/dem_test_contact/diag_table b/tests/dem_test_contact/diag_table new file mode 100644 index 0000000..c1c8819 --- /dev/null +++ b/tests/dem_test_contact/diag_table @@ -0,0 +1,3 @@ +"icebergs test" +1 1 1 0 0 0 + diff --git a/tests/dem_test_contact/input.nml b/tests/dem_test_contact/input.nml new file mode 100644 index 0000000..e3e1133 --- /dev/null +++ b/tests/dem_test_contact/input.nml @@ -0,0 +1,200 @@ + +&icebergs_driver_nml + ni=45 !30 !number of elements, x direction + nj=45 !30 !number of elements, y direction + gridres=5000.0 + ibdt=1800.0 !seconds + ibuo=0.2 !ocean x velocity + ibvo=0.2 !-0.1 !ocean y velocity + ibui=0.0 !ice x velocity + ibvi=0.0 !ice y velocity + ibhrs=212!480 !20 days !1440 !960 !total simulation hours (60 !40 days) + write_time_inc = 48 !increment at which time is written to screen + nmax=999999999 !max number of timesteps + saverestart=.true. + halo=3 + debug=.true. + collision_test=.false. + !big_grounding_test=.true. + dem_test_contact=.true. + bump_depth=50. !200.0 !(m) depth of top of gaussian bump below sea level +/ + +&icebergs_nml + + constant_interaction_LW=.true. + !traj_name=test_dem.nc + + dem=.true. + poisson=0.3 + dem_damping_coef=1.0 + explicit_inner_mts=.true. + break_bonds_on_sub_steps=.true. + + rho_bergs=850. ! (850) iceberg density + !dem_spring_coef=E=youngs. For hexagonal packing and a constant radius and density for all elements, + !a dem_spring_coef equivalent to the spring_coef (k) used for non-DEM would be: + !E=k*rho_bergs*2*sqrt(3)*(radius**2). If rho_bergs=850, radius=389.71 m, k=1.e-5, then E=4471.94 + !if square packing, use A=4*(radius**2) instead of A=2*sqrt(3)*(radius**2) + !E=34000 for square packing with 1000 m radius, rho_bergs=850, and k=1.e-5. A=L*W always works? + dem_spring_coef=5e6 !5e6 !1e6 !1e7!34000.0 !4471.938678723104 !1000.0 + !spring_coef=0.00015094124684696098 !dem_spring_coef=1.e6, r=1500, and hex pack + !spring_coef=0.0030188249369392196 !dem_spring_coef=5.e6, r=750, and hex pack! + !spring_coef=0.0006037649873878439 !dem_spring_coef=1.e6, r=750, and hex pack + spring_coef=0.0007547062342348049 !if dem_spring_coef=5e6, r=1500, and hex pack? + !spring_coef=0.0015094124684696098 !dem_spring_coef=1.e7, r=1500, and hex pack + + monitor_energy=.false. + !defaults given in parentheses: + new_mts=.true. + mts=.true. !(F) for using the Multiple Time Stepping Velocity Verlet scheme + mts_sub_steps=500 !90 # of mts sub-steps #use 90 for 1.5km radii, 500 works for 500 m radii + explicit_inner_mts=.true. + force_convergence=.true. !(F) experimental MTS convergence scheme to better conserve momentum during collisions + convergence_tolerance=1e-8 !(1.e-8) tolerance for the MTS force_convergence scheme + + contact_distance=4.e3 !(a68) !(0) unbonded bergs contact at max(contact_dist,sum of their radii) + contact_spring_coef=1.e-7!1.e-7!(a68)!6!8 !(default=spring_coeff) Berg spring coeff for collisions + + debug_write=.false. !(F)sets traj_sample_hours=traj_write_hours & includes halo bergs in write + cdrag_grounding=1.e5 !3.16e6 !for power-law-type friction + h_to_init_grounding=20 !200.0 + ewsame=.false. !(F) set T if periodic and 2 PEs along the x (zonal) direction + remove_unused_bergs=.false.!.true. !(T) remove unneeded bergs after PEs transfers + + fracture_criterion='stress' !none' !'stress' !('none') 'strain_rate','strain',or 'none' + frac_thres_scaling = 99999999 !1.05 !0.75 !0.5 + frac_thres_n = 2000 !1850 !!!1000 !2250 !4100 + frac_thres_t = 1000 !1200 !1450 !500 !!!1500 !800 !10000000 !800 !1400 + + + halo=3 !(4) halo width + Lx=225000. ! (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 + + bond_coef=1.e-8 ! (1.e-8) Spring const for berg bonds - not being used right now + radial_damping_coef=0.!1.e-4 ! (1.e-4) Coeff for rel berg motion damping (radial component) + tangental_damping_coef=0.0!2.e-5 !(2.e-5) Coeff for rel berg motion damping (tangential comp.) + scale_damping_by_pmag=.false. !(T) Scales damping by mag of (proj matrix \cdot relvel) + critical_interaction_damping_on=.true. !(T) Sets damping on rel berg vel to crit val + tang_crit_int_damp_on=.false. !(F) critical interaction tamping for tangential component? + 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 (0grdxmax]=grdxmax +CBymin[CBxmingrdymax]=grdymax +CBxc=0.5*(CBxmin+CBxmax); CByc=0.5*(CBymin+CBymax) + +# #if just want one particle per CB +# if just2particles: +# CBxmin=CBxc; CBxmax=CBxc +# CBymin=CByc; CBymax=CByc + +#--- min and max thicknesses for each CB --- +h1=thickness1 +h2=thickness2 +CBhmax=np.array([h1,h2]); CBhmin=np.array([h1,h2]) + +#radii +CBrad=np.array([radius,radius]) +#print('radii',CBrad) + +berg_x=[]; berg_y=[] +berg_id=[]; berg_static=[] +berg_width=[]; berg_bonds=[] +berg_h=[]; berg_mass_scaling=[] +berg_mass=[] +berg_uvel=[]; berg_vvel=[] +berg_uvel_in=np.array([0.2,0.2]); berg_vvel_in=np.array([0.2,-0.2]) + +berg_count=0 + +for i in range(nbergs): + # if (i>0): + # print('num bergs in 1st conglom',berg_count) + + x_start=CBxmin[i]+(CBrad[i]*2./np.sqrt(3)) + + if (i>0): + x_start=x_start+1.5e3 + + # if flip2ndconglom and i>0: + # x_start=CBxmax[i]-(CBrad[i]*2./np.sqrt(3)) + + #pdb.set_trace() + if x_start>CBxmax[i]: + x_start=CBxmax[i] + y_start0=CBymin[i]+CBrad[i] + if y_start0>CBymax[i]: + y_start0=CBymax[i] + element_area=(3.*np.sqrt(3.)/2.)*((4./3.)*(CBrad[i])**2) + #H (thickness) is a linear function of a berg element's position from the + #center of the CB. H=Hmax at the center of the CB and H=Hmin at a corner of the + #CB, where the dist of a corner from the center is: + cdistb=np.sqrt((CBxmin[i]-CBxc[i])**2+(CBymin[i]-CByc[i])**2) + j=0 + x_val=x_start + offset=0.0 + # if (i==0): + uvel=berg_uvel_in[i] + #0.1#1 + # else: + # uvel=-0.15#-0.05 + # if (offset2ndconglom): + # offset=250.0 + vvel=berg_vvel_in[i] + + #print('uvel vvel',uvel,vvel) + #0.0#25 + #uvel=0;vvel=0 + #berg_count_start=berg_count + while x_val<=CBxmax[i] and x_val>=CBxmin[i]: + y_start=y_start0+((j%2)*CBrad[i])+offset + k=0 + y_val=y_start + while y_val<=(CBymax[i]+offset): + berg_count=berg_count+1 + berg_id.append(berg_count) + berg_x.append(x_val) + berg_y.append(y_val) + berg_width.append(sqrt(element_area)) + #dist of berg elem from center of CB + bdistc=np.sqrt((x_val-CBxc[i])**2+(y_val-CByc[i])**2) + bh=CBhmin[i]*bdistc/cdistb + CBhmax[i]*(1-bdistc/cdistb) + # if (just2particles): + # bh=h1 + berg_h.append(bh) #thickness + berg_mass_scaling.append(1) + berg_mass.append(bh*rho_ice*element_area) + berg_static.append(0) + berg_uvel.append(uvel) + berg_vvel.append(vvel) + # if (x_val<22000): + # berg_vvel.append(vvel) + # else: + # berg_vvel.append(0.4) + #berg_CBid.append(i) + k=k+1 + y_val=y_start+(2*k*CBrad[i]) + j=j+1 + # if flip2ndconglom and i>0: + # x_val=x_start-(np.sqrt(3)*CBrad[i]*j) + # else: + x_val=x_start+(np.sqrt(3)*CBrad[i]*j) + +#print('Number of bergs',berg_count) +print 'Number of bergs',berg_count + +#pdb.set_trace() + +#Create iceberg restart file +Ice_geometry_source='Generic' +Create_iceberg_restart_file(berg_count,berg_x,berg_y,berg_h,berg_width,berg_mass,berg_mass_scaling,berg_id,Ice_geometry_source,berg_static,berg_uvel,berg_vvel) diff --git a/tests/nc2csv.py b/tests/nc2csv.py index d5d390c..05a9602 100755 --- a/tests/nc2csv.py +++ b/tests/nc2csv.py @@ -150,9 +150,9 @@ def str2bool(string): elif string.lower() in ("no", "false", "f", "0"): Value=False else: - print '**********************************************************************' - print 'The input variable ' ,str(string) , ' is not suitable for boolean conversion, using default' - print '**********************************************************************' + print( '**********************************************************************') + print( 'The input variable ' ,str(string) , ' is not suitable for boolean conversion, using default') + print( '**********************************************************************') Value=None return