diff --git a/libglide/felix_dycore_interface.F90 b/libglide/felix_dycore_interface.F90 index f702d9a1..ce1ed0dc 100644 --- a/libglide/felix_dycore_interface.F90 +++ b/libglide/felix_dycore_interface.F90 @@ -146,7 +146,7 @@ end subroutine felix_velo_init subroutine felix_velo_driver(model) use glimmer_global, only : dp - use glimmer_physcon, only: gn, scyr + use glimmer_physcon, only: scyr use glimmer_paramets, only: thk0, len0, vel0, vis0 use glimmer_log use glide_types diff --git a/libglide/glide_setup.F90 b/libglide/glide_setup.F90 index 3b967d05..c5f4cc3d 100644 --- a/libglide/glide_setup.F90 +++ b/libglide/glide_setup.F90 @@ -160,7 +160,7 @@ end subroutine glide_printconfig subroutine glide_scale_params(model) !> scale parameters use glide_types - use glimmer_physcon, only: scyr, gn + use glimmer_physcon, only: scyr use glimmer_paramets, only: thk0, tim0, len0, vel0, vis0, acc0, tau0 implicit none @@ -883,11 +883,10 @@ subroutine print_options(model) 'advective-diffusive balance ',& 'temp from external file ' /) - character(len=*), dimension(0:3), parameter :: flow_law = (/ & - 'const 1e-16 Pa^-n a^-1 ', & + character(len=*), dimension(0:2), parameter :: flow_law = (/ & + 'uniform factor flwa ', & 'Paterson and Budd (T = -5 C)', & - 'Paterson and Budd ', & - 'read flwa/flwastag from file' /) + 'Paterson and Budd ' /) !TODO - Rename slip_coeff to which_btrc? character(len=*), dimension(0:5), parameter :: slip_coeff = (/ & @@ -1996,7 +1995,7 @@ subroutine handle_parameters(section, model) use glimmer_config use glide_types use glimmer_log - use glimmer_physcon, only: rhoi, rhoo, grav, shci, lhci, trpt + use glimmer_physcon, only: rhoi, rhoo, grav, shci, lhci, trpt, n_glen implicit none type(ConfigSection), pointer :: section @@ -2021,6 +2020,7 @@ subroutine handle_parameters(section, model) call GetValue(section,'lhci', lhci) call GetValue(section,'trpt', trpt) #endif + call GetValue(section,'n_glen', n_glen) loglevel = GM_levels-GM_ERROR call GetValue(section,'log_level',loglevel) @@ -2033,9 +2033,9 @@ subroutine handle_parameters(section, model) call GetValue(section,'pmp_offset', model%temper%pmp_offset) call GetValue(section,'pmp_threshold', model%temper%pmp_threshold) call GetValue(section,'geothermal', model%paramets%geot) - !TODO - Change default_flwa to flwa_constant? Would have to change config files. call GetValue(section,'flow_factor', model%paramets%flow_enhancement_factor) call GetValue(section,'flow_factor_float', model%paramets%flow_enhancement_factor_float) + !TODO - Change default_flwa to flwa_constant? Would have to change config files. call GetValue(section,'default_flwa', model%paramets%default_flwa) call GetValue(section,'efvs_constant', model%paramets%efvs_constant) call GetValue(section,'effstrain_min', model%paramets%effstrain_min) @@ -2206,7 +2206,7 @@ end subroutine handle_parameters subroutine print_parameters(model) - use glimmer_physcon, only: rhoi, rhoo, lhci, shci, trpt, grav + use glimmer_physcon, only: rhoi, rhoo, lhci, shci, trpt, grav, n_glen use glide_types use glimmer_log implicit none @@ -2371,6 +2371,9 @@ subroutine print_parameters(model) write(message,*) 'triple point of water (K) : ', trpt call write_log(message) + write(message,*) 'Glen flow law exponent : ', n_glen + call write_log(message) + write(message,*) 'geothermal flux (W/m^2) : ', model%paramets%geot call write_log(message) diff --git a/libglide/glide_types.F90 b/libglide/glide_types.F90 index 147608e9..d55a727e 100644 --- a/libglide/glide_types.F90 +++ b/libglide/glide_types.F90 @@ -104,7 +104,6 @@ module glide_types integer, parameter :: FLWA_CONST_FLWA = 0 integer, parameter :: FLWA_PATERSON_BUDD_CONST_TEMP = 1 integer, parameter :: FLWA_PATERSON_BUDD = 2 - integer, parameter :: FLWA_INPUT = 3 integer, parameter :: BTRC_ZERO = 0 integer, parameter :: BTRC_CONSTANT = 1 @@ -470,7 +469,6 @@ module glide_types !> \item[1] \emph{Paterson and Budd} relationship, !> with temperature set to $-5^{\circ}\mathrm{C}$ !> \item[2] \emph{Paterson and Budd} relationship - !> \item[3] Read flwa/flwastag from file !> \end{description} integer :: whichbtrc = 0 @@ -2148,6 +2146,7 @@ module glide_types !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !TODO - Move these parameters to types associated with a certain kind of physics + !TODO - Set default geot = 0, so that idealized tests by default have no mass loss type glide_paramets real(dp),dimension(5) :: bpar = (/ 0.2d0, 0.5d0, 0.0d0 ,1.0d-2, 1.0d0/) real(dp) :: btrac_const = 0.d0 ! m yr^{-1} Pa^{-1} (gets scaled during init) diff --git a/libglimmer/glimmer_paramets.F90 b/libglimmer/glimmer_paramets.F90 index 019f44e6..aa8b595d 100644 --- a/libglimmer/glimmer_paramets.F90 +++ b/libglimmer/glimmer_paramets.F90 @@ -33,7 +33,7 @@ module glimmer_paramets use glimmer_global, only : dp - use glimmer_physcon, only : scyr, gn + use glimmer_physcon, only : scyr implicit none save @@ -118,6 +118,7 @@ module glimmer_paramets real(dp), parameter :: grav_glam = 9.81d0 ! m s^{-2} ! GLAM scaling parameters; units are correct if thk0 has units of meters + integer, parameter :: gn = 3 ! Glen flow exponent; fixed at 3 for purposes of setting vis0 real(dp), parameter :: tau0 = rhoi_glam*grav_glam*thk0 ! stress scale in GLAM ( Pa ) real(dp), parameter :: evs0 = tau0 / (vel0/len0) ! eff. visc. scale in GLAM ( Pa s ) real(dp), parameter :: vis0 = tau0**(-gn) * (vel0/len0) ! rate factor scale in GLAM ( Pa^-3 s^-1 ) diff --git a/libglimmer/glimmer_physcon.F90 b/libglimmer/glimmer_physcon.F90 index 0c990d29..f697bf3e 100644 --- a/libglimmer/glimmer_physcon.F90 +++ b/libglimmer/glimmer_physcon.F90 @@ -79,11 +79,17 @@ module glimmer_physcon ! real(dp) :: trpt = 273.16d0 !< Triple point of water (K) #endif + ! The Glen flow-law exponent, n_glen, is used in Glissade. + ! It is not a parameter, because the default can be overridden in the config file. + ! TODO: Allow setting n_glen independently for each ice sheet instance? + ! Note: Earlier code used an integer parameter, gn = 3, for all flow-law calculations. + ! For backward compatiblity, gn = 3 is retained for Glide. + real(dp) :: n_glen = 3.0d0 !< Exponent in Glen's flow law; user-configurable real(dp) in Glissade + integer, parameter :: gn = 3 !< Exponent in Glen's flow law; fixed integer parameter in Glide real(dp),parameter :: celsius_to_kelvin = 273.15d0 !< Note: Not quite equal to trpt real(dp),parameter :: scyr = 31536000.d0 !< Number of seconds in a year of exactly 365 days real(dp),parameter :: rhom = 3300.0d0 !< The density of magma(?) (kg m-3) real(dp),parameter :: rhos = 2600.0d0 !< The density of solid till (kg m$^{-3}$) - integer, parameter :: gn = 3 !< The power dependency of Glen's flow law. real(dp),parameter :: actenh = 139.0d3 !< Activation energy in Glen's flow law for \f$T^{*}\geq263\f$K. (J mol-1) real(dp),parameter :: actenl = 60.0d3 !< Activation energy in Glen's flow law for \f$T^{*}<263\f$K. (J mol-1) real(dp),parameter :: arrmlh = 1.733d3 !< Constant of proportionality in Arrhenius relation diff --git a/libglimmer/glimmer_scales.F90 b/libglimmer/glimmer_scales.F90 index f95c6a86..380ff2b8 100644 --- a/libglimmer/glimmer_scales.F90 +++ b/libglimmer/glimmer_scales.F90 @@ -45,7 +45,7 @@ subroutine glimmer_init_scales ! set scale factors for I/O (can't have non-integer powers) - use glimmer_physcon, only : scyr, gn + use glimmer_physcon, only : scyr use glimmer_paramets, only : thk0, tim0, vel0, vis0, len0, acc0, tau0, evs0 implicit none diff --git a/libglissade/glissade.F90 b/libglissade/glissade.F90 index ef0dc5dd..efa2d483 100644 --- a/libglissade/glissade.F90 +++ b/libglissade/glissade.F90 @@ -457,7 +457,11 @@ subroutine glissade_initialise(model, evolve_ice) ! handle relaxed/equilibrium topo ! Initialise isostasy first - call init_isostasy(model) + if (model%options%isostasy == ISOSTASY_COMPUTE) then + + call init_isostasy(model) + + endif select case(model%isostasy%whichrelaxed) @@ -1895,7 +1899,7 @@ subroutine glissade_thermal_solve(model, dt) model%temper%btemp_ground, & ! deg C model%temper%btemp_float, & ! deg C bmlt_ground_unscaled) ! m/s - + ! Update basal hydrology, if needed ! Note: glissade_calcbwat uses SI units @@ -1973,6 +1977,7 @@ subroutine glissade_thickness_tracer_solve(model) use glissade_bmlt_float, only: verbose_bmlt_float use glissade_calving, only: verbose_calving use glissade_grid_operators, only: glissade_vertical_interpolate + use glide_stop, only: glide_finalise implicit none @@ -2161,21 +2166,25 @@ subroutine glissade_thickness_tracer_solve(model) model%geomderv%dusrfdew*thk0/len0, model%geomderv%dusrfdns*thk0/len0, & model%velocity%uvel * scyr * vel0, model%velocity%vvel * scyr * vel0, & model%numerics%dt_transport * tim0 / scyr, & + model%numerics%adaptive_cfl_threshold, & model%numerics%adv_cfl_dt, model%numerics%diff_cfl_dt) ! Set the transport timestep. ! The timestep is model%numerics%dt by default, but optionally can be reduced for subcycling + !WHL - debug +! if (main_task) then +! print*, 'Checked advective CFL threshold' +! print*, 'model dt (yr) =', model%numerics%dt * tim0/scyr +! print*, 'adv_cfl_dt =', model%numerics%adv_cfl_dt +! endif + + advective_cfl = model%numerics%dt*(tim0/scyr) / model%numerics%adv_cfl_dt + if (model%numerics%adaptive_cfl_threshold > 0.0d0) then - !WHL - debug -! if (main_task) then -! print*, 'Check advective CFL threshold' -! print*, 'model dt (yr) =', model%numerics%dt * tim0/scyr -! print*, 'adv_cfl_dt =', model%numerics%adv_cfl_dt -! endif + ! subcycle the transport when advective_cfl exceeds the threshold - advective_cfl = model%numerics%dt*(tim0/scyr) / model%numerics%adv_cfl_dt if (advective_cfl > model%numerics%adaptive_cfl_threshold) then ! compute the number of subcycles @@ -2188,14 +2197,29 @@ subroutine glissade_thickness_tracer_solve(model) print*, 'Ratio =', advective_cfl / model%numerics%adaptive_cfl_threshold print*, 'nsubcyc =', nsubcyc endif + else nsubcyc = 1 endif dt_transport = model%numerics%dt * tim0 / real(nsubcyc,dp) ! convert to s else ! no adaptive subcycling - nsubcyc = model%numerics%subcyc - dt_transport = model%numerics%dt_transport * tim0 ! convert to s + + advective_cfl = model%numerics%dt*(tim0/scyr) / model%numerics%adv_cfl_dt + + ! If advective_cfl exceeds 1.0, then abort cleanly. Otherwise, set dt_transport and proceed. + ! Note: Usually, it would be enough to write a fatal abort message. + ! The call to glide_finalise was added to allow CISM to finish cleanly when running + ! a suite of automated stability tests, e.g. with the stabilitySlab.py script. + if (advective_cfl > 1.0d0) then + if (main_task) print*, 'advective CFL violation; call glide_finalise and exit cleanly' + call glide_finalise(model, crash=.true.) + stop + else + nsubcyc = model%numerics%subcyc + dt_transport = model%numerics%dt_transport * tim0 ! convert to s + endif + endif !------------------------------------------------------------------------- diff --git a/libglissade/glissade_basal_traction.F90 b/libglissade/glissade_basal_traction.F90 index 69f1fd75..5dc1e2cb 100644 --- a/libglissade/glissade_basal_traction.F90 +++ b/libglissade/glissade_basal_traction.F90 @@ -440,9 +440,12 @@ subroutine calcbeta (whichbabc, & ! If this factor is not present in the input file, it is set to 1 everywhere. ! Compute beta - ! gn = Glen's n from physcon module - beta(:,:) = coulomb_c * basal_physics%effecpress_stag(:,:) * speed(:,:)**(1.0d0/gn - 1.0d0) * & - (speed(:,:) + basal_physics%effecpress_stag(:,:)**gn * big_lambda)**(-1.0d0/gn) + ! Note: Where this equation has powerlaw_m, we used to have Glen's flow exponent n, + ! following the notation of Leguy et al. (2014). + ! Changed to powerlaw_m to be consistent with the Schoof and Tsai laws. + m = basal_physics%powerlaw_m + beta(:,:) = coulomb_c * basal_physics%effecpress_stag(:,:) * speed(:,:)**(1.0d0/m - 1.0d0) * & + (speed(:,:) + basal_physics%effecpress_stag(:,:)**m * big_lambda)**(-1.0d0/m) ! If c_space_factor /= 1.0 everywhere, then multiply beta by c_space_factor if (maxval(abs(basal_physics%c_space_factor_stag(:,:) - 1.0d0)) > tiny(0.0d0)) then diff --git a/libglissade/glissade_therm.F90 b/libglissade/glissade_therm.F90 index 10b1fcca..f155521b 100644 --- a/libglissade/glissade_therm.F90 +++ b/libglissade/glissade_therm.F90 @@ -1115,11 +1115,14 @@ subroutine glissade_therm_driver(whichtemp, & dissipcol(ew,ns) = dissipcol(ew,ns) * thck(ew,ns)*rhoi*shci ! Verify that the net input of energy into the column is equal to the change in - ! internal energy. + ! internal energy. delta_e = (ucondflx(ew,ns) - lcondflx(ew,ns) + dissipcol(ew,ns)) * dttem - if (abs((efinal-einit-delta_e)/dttem) > 1.0d-8) then + ! Note: For very small dttem (e.g., 1.0d-6 year or less), this error can be triggered + ! by roundoff error. In that case, the user may need to increase the threshold. + ! July 2021: Increased from 1.0d-8 to 1.0d-7 to allow smaller dttem. + if (abs((efinal-einit-delta_e)/dttem) > 1.0d-7) then if (verbose_column) then print*, 'Ice thickness:', thck(ew,ns) @@ -1640,10 +1643,11 @@ subroutine glissade_enthalpy_matrix_elements(dttem, & ! At each temperature point, compute the temperature part of the enthalpy. ! enth_T = enth for cold ice, enth_T < enth for temperate ice - enth_T(0) = rhoi*shci*temp(0) !WHL - not sure enth_T(0) is needed - do up = 1, upn + do up = 1, upn-1 enth_T(up) = (1.d0 - waterfrac(up)) * rhoi*shci*temp(up) enddo + enth_T(0) = rhoi*shci*temp(0) + enth_T(up) = rhoi*shci*temp(up) !WHL - debug if (verbose_column) then @@ -2250,8 +2254,7 @@ subroutine glissade_interior_dissipation_sia(ewn, nsn, & ! Compute the dissipation source term associated with strain heating, ! based on the shallow-ice approximation. - use glimmer_physcon, only : gn ! Glen's n - use glimmer_physcon, only: rhoi, shci, grav + use glimmer_physcon, only: rhoi, shci, grav, n_glen integer, intent(in) :: ewn, nsn, upn ! grid dimensions @@ -2267,12 +2270,14 @@ subroutine glissade_interior_dissipation_sia(ewn, nsn, & real(dp), dimension(:,:,:), intent(out) :: & dissip ! interior heat dissipation (deg/s) - integer, parameter :: p1 = gn + 1 - integer :: ew, ns real(dp), dimension(upn-1) :: sia_dissip_fact ! factor in SIA dissipation calculation real(dp) :: geom_fact ! geometric factor + real(dp) :: p1 ! exponent = n_glen + 1 + + p1 = n_glen + 1.0d0 + ! Two methods of doing this calculation: ! 1. find dissipation at u-pts and then average ! 2. find dissipation at H-pts by averaging quantities from u-pts @@ -2414,7 +2419,7 @@ subroutine glissade_flow_factor(whichflwa, whichtemp, & integer, intent(in) :: whichflwa !> which method of calculating A integer, intent(in) :: whichtemp !> which method of calculating temperature; - !> include waterfrac in calculation if using enthalpy method + !> include waterfrac in calculation if using enthalpy method real(dp),dimension(:), intent(in) :: stagsigma !> vertical coordinate at layer midpoints real(dp),dimension(:,:), intent(in) :: thck !> ice thickness (m) real(dp),dimension(:,:,:), intent(in) :: temp !> 3D temperature field (deg C) @@ -2488,17 +2493,16 @@ subroutine glissade_flow_factor(whichflwa, whichtemp, & endif ! Multiply the default rate factor by the enhancement factor if applicable - ! Note: Here, default_flwa is assumed to have units of Pa^{-n} s^{-1}, + ! Note: Here, the input default_flwa is assumed to have units of Pa^{-n} s^{-1}, ! whereas model%paramets%default_flwa has units of Pa^{-n} yr^{-1}. ! initialize - if (whichflwa /= FLWA_INPUT) then - do ns = 1, nsn - do ew = 1, ewn - flwa(:,ew,ns) = enhancement_factor(ew,ns) * default_flwa - enddo + !TODO - Move the next few lines inside the select case construct. + do ns = 1, nsn + do ew = 1, ewn + flwa(:,ew,ns) = enhancement_factor(ew,ns) * default_flwa enddo - endif + enddo select case(whichflwa) @@ -2558,21 +2562,12 @@ subroutine glissade_flow_factor(whichflwa, whichtemp, & end do case(FLWA_CONST_FLWA) - - ! do nothing (flwa is initialized to default_flwa above) - - case(FLWA_INPUT) - ! do nothing - use flwa from input or forcing file - print *, 'FLWA', minval(flwa), maxval(flwa) + ! do nothing (flwa is set above, with units Pa^{-n} s^{-1}) end select - ! This logic assumes that the input flwa is already in dimensionless model units. - ! TODO: Make a different assumption about input units? - if (whichflwa /= FLWA_INPUT) then - ! Change flwa to model units (glissade_flow_factor assumes SI units of Pa{-n} s^{-1}) - flwa(:,:,:) = flwa(:,:,:) / vis0 - endif + ! Change flwa to model units (glissade_flow_factor assumes SI units of Pa{-n} s^{-1}) + flwa(:,:,:) = flwa(:,:,:) / vis0 deallocate(enhancement_factor) diff --git a/libglissade/glissade_transport.F90 b/libglissade/glissade_transport.F90 index a1c57219..e0974b96 100644 --- a/libglissade/glissade_transport.F90 +++ b/libglissade/glissade_transport.F90 @@ -979,6 +979,7 @@ subroutine glissade_check_cfl(ewn, nsn, nlyr, & parallel, & stagthk, dusrfdew, dusrfdns, & uvel, vvel, deltat, & + adaptive_cfl_threshold, & allowable_dt_adv, allowable_dt_diff) ! Calculate maximum allowable time step based on both @@ -1015,6 +1016,10 @@ subroutine glissade_check_cfl(ewn, nsn, nlyr, & real(dp), intent(in) :: & deltat ! model deltat (yrs) + real(dp), intent(in) :: & + adaptive_cfl_threshold ! threshold for adaptive subcycling + ! if = 0, there is no adaptive subcycling; code aborts when CFL > 1 + real(dp), intent(out) :: & allowable_dt_adv ! maximum allowable dt (yrs) based on advective CFL @@ -1159,13 +1164,16 @@ subroutine glissade_check_cfl(ewn, nsn, nlyr, & write(message,*) 'Advective CFL violation! Maximum allowable time step for advective CFL condition is ' & // trim(adjustl(dt_string)) // ' yr, limited by global position i=' & // trim(adjustl(xpos_string)) // ' j=' //trim(adjustl(ypos_string)) + call write_log(trim(message),GM_WARNING) - ! If the violation is egregious (defined as deltat > 10 * allowable_dt_adv), then abort. - ! Otherwise, write a warning and proceed. - if (deltat > 10.d0 * allowable_dt_adv) then - call write_log(trim(message),GM_FATAL) - else - call write_log(trim(message),GM_WARNING) + ! If adaptive subcyling is allowed, then make the code abort for egregious CFL violations, + ! (defined as deltat > 10 * allowable_dt_adv), to prevent excessive subcycling. + + if (main_task .and. adaptive_cfl_threshold > 0.0d0) then + if (deltat > 10.d0 * allowable_dt_adv) then + print*, 'deltat, allowable_dt_adv, ratio =', deltat, allowable_dt_adv, deltat/allowable_dt_adv + call write_log('Aborting with CFL violation', GM_FATAL) + endif endif endif diff --git a/libglissade/glissade_velo.F90 b/libglissade/glissade_velo.F90 index bb6e6e38..a9dadd17 100644 --- a/libglissade/glissade_velo.F90 +++ b/libglissade/glissade_velo.F90 @@ -43,9 +43,6 @@ subroutine glissade_velo_driver(model) ! Glissade higher-order velocity driver - use glimmer_global, only : dp - use glimmer_physcon, only: gn, scyr - use glimmer_paramets, only: thk0, len0, vel0, vis0, tau0, evs0 use glimmer_log use glide_types use glissade_velo_higher, only: glissade_velo_higher_solve diff --git a/libglissade/glissade_velo_higher.F90 b/libglissade/glissade_velo_higher.F90 index 033ca254..77b2f9f5 100644 --- a/libglissade/glissade_velo_higher.F90 +++ b/libglissade/glissade_velo_higher.F90 @@ -57,7 +57,7 @@ module glissade_velo_higher use glimmer_global, only: dp - use glimmer_physcon, only: gn, rhoi, rhoo, grav, scyr, pi + use glimmer_physcon, only: n_glen, rhoi, rhoo, grav, scyr, pi use glimmer_paramets, only: eps08, eps10, thk0, len0, tim0, tau0, vel0, vis0, evs0 use glimmer_paramets, only: vel_scale, len_scale ! used for whichefvs = HO_EFVS_FLOWFACT use glimmer_log @@ -200,8 +200,8 @@ module glissade_velo_higher ! logical :: verbose = .true. logical :: verbose_init = .false. ! logical :: verbose_init = .true. -! logical :: verbose_solver = .false. - logical :: verbose_solver = .true. + logical :: verbose_solver = .false. +! logical :: verbose_solver = .true. logical :: verbose_Jac = .false. ! logical :: verbose_Jac = .true. logical :: verbose_residual = .false. @@ -2082,7 +2082,7 @@ subroutine glissade_velo_higher_solve(model, & ! gn = exponent in Glen's flow law (= 3 by default) do k = 1, nz-1 if (flwa(k,i,j) > 0.0d0) then - flwafact(k,i,j) = 0.5d0 * flwa(k,i,j)**(-1.d0/real(gn,dp)) + flwafact(k,i,j) = 0.5d0 * flwa(k,i,j)**(-1.d0/n_glen) endif enddo enddo @@ -4222,6 +4222,7 @@ subroutine glissade_velo_higher_solve(model, & usrf, & dusrf_dx, dusrf_dy, & flwa, efvs, & + whichefvs, efvs_constant, & whichgradient_margin, & max_slope, & uvel, vvel) @@ -6426,6 +6427,7 @@ subroutine compute_3d_velocity_L1L2(nx, ny, & usrf, & dusrf_dx, dusrf_dy, & flwa, efvs, & + whichefvs, efvs_constant, & whichgradient_margin, & max_slope, & uvel, vvel) @@ -6486,6 +6488,12 @@ subroutine compute_3d_velocity_L1L2(nx, ny, & flwa, & ! temperature-based flow factor A, Pa^{-n} yr^{-1} efvs ! effective viscosity, Pa yr + integer, intent(in) :: & + whichefvs ! option for effective viscosity calculation + + real(dp), intent(in) :: & + efvs_constant ! constant value of effective viscosity (Pa yr) + integer, intent(in) :: & whichgradient_margin ! option for computing gradient at ice margin ! 0 = include all neighbor cells in gradient calculation @@ -6840,7 +6848,7 @@ subroutine compute_3d_velocity_L1L2(nx, ny, & ! Compute vertical integration factor at each active vertex ! This is int_b_to_z{-2 * A * tau^2 * rho*g*(s-z) * dz}, - ! similar to the factor computed in Glide and glissade_velo_sia.. + ! similar to the factor computed in Glide and glissade_velo_sia. ! Note: tau_xz ~ rho*g*(s-z)*ds_dx; ds_dx term is computed on edges below do j = 1, ny-1 @@ -6921,9 +6929,27 @@ subroutine compute_3d_velocity_L1L2(nx, ny, & tau_eff_sq = stagtau_parallel_sq(i,j) & + tau_xz(k,i,j)**2 + tau_yz(k,i,j)**2 - ! Note: This formula is correct for any value of Glen's n, but currently efvs is computed - ! only for gn = 3 (in which case (n-1)/2 = 1). - fact = 2.d0 * stagflwa(i,j) * tau_eff_sq**((gn-1.d0)/2.d0) * (sigma(k+1) - sigma(k))*stagthck(i,j) + ! Note: The first formula below is correct for whichefvs = 2 (efvs computed from effective strain rate), + ! but not for whichefvs = 0 (constant efvs) or whichefvs = 1 (multiple of flow factor). + ! For these options we need a modified formula. + ! + ! Recall: efvs = 1/2 * A^(-1/n) * eps_e^[(1-n)/n] + ! = 1/2 * A^(-1/n) * [A tau_e^n]^[(1-n)/n] + ! = 1/2 * A^(-1) * tau_e^(1-n) + ! => 1/efvs = 2 * A * tau_e(n-1) + ! + ! Thus, for options 0 and 1, we can replace 2 * A * tau_e^(n-1) below with 1/efvs. + + if (whichefvs == HO_EFVS_NONLINEAR) then + fact = 2.d0 * stagflwa(i,j) * tau_eff_sq**((n_glen-1.d0)/2.d0) & + * (sigma(k+1) - sigma(k))*stagthck(i,j) + else ! HO_EFVS_CONSTANT, HO_EFVS_FLOWFACT + if (efvs(k,i,j) > 0.0d0) then + fact = (sigma(k+1) - sigma(k))*stagthck(i,j) / efvs(k,i,j) + else + fact = 0.0d0 + endif + endif ! reset velocity to prescribed basal value if Dirichlet condition applies ! else compute velocity at this level @@ -7875,15 +7901,6 @@ subroutine compute_effective_viscosity (whichefvs, whichapprox, & integer, intent(in) :: i, j, k, p - !---------------------------------------------------------------- - ! Local parameters - !---------------------------------------------------------------- - - real(dp), parameter :: & - p_effstr = (1.d0 - real(gn,dp))/real(gn,dp), &! exponent (1-n)/n in effective viscosity relation - p2_effstr = p_effstr/2 ! exponent (1-n)/(2n) in effective viscosity relation - - !---------------------------------------------------------------- ! Local variables !---------------------------------------------------------------- @@ -7896,8 +7913,14 @@ subroutine compute_effective_viscosity (whichefvs, whichapprox, & integer :: n + real(dp) :: & + p_effstr ! exponent (1-n)/n in effective viscosity relation + real(dp), parameter :: p2 = -1.d0/3.d0 + ! Set exponent that depends on Glen's exponent + p_effstr = (1.d0 - n_glen)/n_glen + select case(whichefvs) case(HO_EFVS_CONSTANT) @@ -7988,11 +8011,11 @@ subroutine compute_effective_viscosity (whichefvs, whichapprox, & ! Compute effective viscosity (PGB 2012, eq. 4) ! Units: flwafact has units Pa yr^{1/n} ! effstrain has units yr^{-1} - ! p2_effstr = (1-n)/(2n) - ! = -1/3 for n=3 + ! p_effstr = (1-n)/n + ! = -2/3 for n=3 ! Thus efvs has units Pa yr - efvs = flwafact * effstrainsq**p2_effstr + efvs = flwafact * effstrainsq**(p_effstr/2.d0) if (verbose_efvs .and. this_rank==rtest .and. i==itest .and. j==jtest .and. k==ktest .and. p==ptest) then print*, ' ' @@ -8081,8 +8104,8 @@ subroutine compute_effective_viscosity_L1L2(whichefvs, ! Local parameters !---------------------------------------------------------------- - real(dp), parameter :: & - p_effstr = (1.d0 - real(gn,dp)) / real(gn,dp) ! exponent (1-n)/n in effective viscosity relation + real(dp) :: & + p_effstr ! exponent (1-n)/n in effective viscosity relation !---------------------------------------------------------------- ! Local variables @@ -8107,6 +8130,9 @@ subroutine compute_effective_viscosity_L1L2(whichefvs, integer :: n, k + ! Set exponent that depends on Glen's exponent + p_effstr = (1.d0 - n_glen) / n_glen + select case(whichefvs) case(HO_EFVS_CONSTANT) @@ -8125,7 +8151,7 @@ subroutine compute_effective_viscosity_L1L2(whichefvs, ! ! Units: flwafact has units Pa yr^{1/n} ! effstrain has units yr^{-1} - ! p_effstr = (1-n)/n + ! p_effstr = (1-n)/n ! = -2/3 for n=3 ! Thus efvs has units Pa yr @@ -8320,14 +8346,6 @@ subroutine compute_effective_viscosity_diva(whichefvs, integer, intent(in) :: i, j, p - !---------------------------------------------------------------- - ! Local parameters - !---------------------------------------------------------------- - - real(dp), parameter :: & - p_effstr = (1.d0 - real(gn,dp))/real(gn,dp), &! exponent (1-n)/n in effective viscosity relation - p2_effstr = p_effstr/2 ! exponent (1-n)/(2n) in effective viscosity relation - !---------------------------------------------------------------- ! Local variables !---------------------------------------------------------------- @@ -8346,11 +8364,17 @@ subroutine compute_effective_viscosity_diva(whichefvs, integer :: n, k real(dp) :: du_dz, dv_dz + real(dp) :: & + p_effstr ! exponent (1-n)/n in effective viscosity relation + !WHL - For ISMIP-HOM, the cubic solve is not robust. It leads to oscillations ! in successive iterations between uvel_2d/vvel_2d and btractx/btracty !TODO - Remove the cubic solve for efvs, unless we find a way to make it robust? logical, parameter :: cubic = .false. + ! Set exponent that depends on Glen's exponent + p_effstr = (1.d0 - n_glen)/n_glen + select case(whichefvs) case(HO_EFVS_CONSTANT) @@ -8493,7 +8517,8 @@ subroutine compute_effective_viscosity_diva(whichefvs, effstrainsq = effstrain_min**2 & + du_dx**2 + dv_dy**2 + du_dx*dv_dy + 0.25d0*(dv_dx + du_dy)**2 & + 0.25d0 * (du_dz**2 + dv_dz**2) - efvs(k) = flwafact(k) * effstrainsq**p2_effstr + efvs(k) = flwafact(k) * effstrainsq**(p_effstr/2.d0) + enddo endif ! cubic diff --git a/libglissade/glissade_velo_sia.F90 b/libglissade/glissade_velo_sia.F90 index 66884aaf..1efcb554 100644 --- a/libglissade/glissade_velo_sia.F90 +++ b/libglissade/glissade_velo_sia.F90 @@ -55,7 +55,7 @@ module glissade_velo_sia use glimmer_global, only: dp - use glimmer_physcon, only: gn, rhoi, grav, scyr + use glimmer_physcon, only: n_glen, rhoi, grav, scyr use glimmer_paramets, only: thk0, len0, vel0, vis0, tau0 ! use glimmer_log, only: write_log @@ -881,16 +881,15 @@ subroutine glissade_velo_sia_interior(nx, ny, nz, & if (stagthck(i,j) > thklim) then - siafact = 2.d0 * (rhoi*grav)**gn * stagthck(i,j)**(gn+1) & - * (dusrf_dx(i,j)**2 + dusrf_dy(i,j)**2) ** ((gn-1)/2) - + siafact = 2.d0 * (rhoi*grav)**n_glen * stagthck(i,j)**(n_glen+1) & + * (dusrf_dx(i,j)**2 + dusrf_dy(i,j)**2) ** ((n_glen-1)/2) vintfact(nz,i,j) = 0.d0 do k = nz-1, 1, -1 - vintfact(k,i,j) = vintfact(k+1,i,j) - & - siafact * stagflwa(k,i,j) & - * ((sigma(k) + sigma(k+1))/2.d0) ** gn & + vintfact(k,i,j) = vintfact(k+1,i,j) - & + siafact * stagflwa(k,i,j) & + * ((sigma(k) + sigma(k+1))/2.d0) ** n_glen & * (sigma(k+1) - sigma(k)) enddo ! k diff --git a/tests/slab/README.md b/tests/slab/README.md index c3767f86..71b95feb 100644 --- a/tests/slab/README.md +++ b/tests/slab/README.md @@ -1,18 +1,88 @@ Slab test case ============== -WARNING: THIS TEST CASE IS IN DEVELOPMENT AND HAS NOT BEEN SCIENTIFICALLY VALIDATED. -USE AT YOUR OWN RISK! +This directory contains python scripts for running an experiment involving a +uniform, infinite ice sheet ("slab") on an inclined plane. +The test case is described in sections 5.1-5.2 of: + Dukowicz, J. K., 2012, Reformulating the full-Stokes ice sheet model for a + more efficient computational solution. The Cryosphere, 6, 21-34, + doi:10.5194/tc-6-21-2012. -This directory contains python scripts for running an experiment involving a -uniform and infinite ice sheet ("slab") on an inclined plane. +Some results from this test case are described in Sect. 3.4 of: + Robinson, A., D. Goldberg, and W. H. Lipscomb, A comparison of the performance + of depth-integrated ice-dynamics solvers. Submitted to The Cryosphere, Aug. 2021. + +The test case consists of an ice slab of uniform thickness moving down an +inclined plane by a combination of sliding and shearing. +Analytic Stokes and first-order velocity solutions exist for all values of Glen's n >= 1. +The solutions for n = 1 are derived in Dukowicz (2012), and solutions for n > 1 +are derived in an unpublished manuscript by Dukowicz (2013). + +The original scripts, runSlab.py and plotSlab.py, were written by Matt Hoffman +with support for Glens' n = 1. They came with warnings that the test is not supported. +The test is now supported, and the scripts include some new features: + +* The user may specify any n >= 1 (not necessarily an integer). + The tests assume which_ho_efvs = 2 (nonlinear viscosity) with flow_law = 0 (constant A). +* Physics parameters are no longer hard-coded. The user can enter the ice thickness, + beta, viscosity coefficient (mu_n), and slope angle (theta) on the command line. +* The user can specify time parameters dt (the dynamic time step) and nt (number of steps). + The previous version did not support transient runs. +* The user can specify a small thickness perturbation dh, which is added to the initial + uniform thickness via random sampling from a Gaussian distribution. + The perturbation will grow or decay, depending on the solver stability for given dx and dt. + +The run script is executed by a command like the following: + +> python runSlab.py -n 4 -a DIVA -theta 0.0375 -thk 1000. -mu 1.e5 -beta 1000. + +In this case, the user runs on 4 processors with the DIVA solver, a slope angle of 0.0375 degrees, +Glen's n = 1 (the default), slab thickness H = 1000 m, sliding coefficient beta = 1000 Pa (m/yr)^{-1}, +and viscosity coefficient 1.e5 Pa yr. +These parameters correspond to the thick shearing test case described by Robinson et al. (2021). + +To see the full set of command-line options, type 'python runSlab.py -h'. + +Notes on effective viscosity: + * For n = 1, the viscosity coefficient mu_1 has a default value of 1.e6 Pa yr in the relation + mu = mu_1 * eps((1-n)/n), where eps is the effective strain rate. + * For n > 1, the user can specify a coefficient mu_n; otherwise the run script computes mu_n + such that the basal and surface speeds are nearly the same as for an n = 1 case with the + mu_1 = 1.e6 Pa yr and the same values of thickness, beta, and theta. + * There is a subtle difference between the Dukowicz and CISM definitions of the + effective strain rate; the Dukowicz value is twice as large. Later, it might be helpful + to make the Dukowicz convention consistent with CISM.) + +The plotting script, plotSlab.py, is run by typing 'python plotSlab.py'. It creates two plots. +The first plot shows the vertical velocity profile in nondimensional units and in units of m/yr. +There is excellent agreement between higher-order CISM solutions and the analytic solution +for small values of the slope angle theta. For steep slopes, the answers diverge as expected. + +For the second plot, the extent of the y-axis is wrong. This remains to be fixed. + +This directory also includes a new script, stabilitySlab.py, to carry out the stability tests +described in Robinson et al. (2021). + +For a given set of physics parameters and stress-balance approximation (DIVA, L1L2, etc.), +the script launches multiple CISM runs at a range of grid resolutions. +At each grid resolution, the script determines the maximum stable time step. +A run is deemed stable when the standard deviation of an initial small thickness perturbation +is reduced over the course of 100 time steps. A run is unstable if the standard deviation +increases or if the model aborts (usually with a CFL violation). + +To run the stability script, type a command like the following: + +> python stabilitySlab.py -n 4 -a DIVA -theta 0.0375 -thk 1000. -mu 1.e5 -beta 1000. \ + -dh 0.1 -nt 100 -nr 12 -rmin 10. -rmax 40000. + +Here, the first few commands correspond to the thick shearing test case and are passed repeatedly +to the run script. The remaining commands specify that each run will be initialized +with a Gaussian perturbation of amplitude 0.1 m and run for 100 timesteps. +The maximum stable timestep will be determined at 12 resolutions ranging from 10m to 40 km. +This test takes several minutes to complete on a Macbook Pro with 4 cores. -The test case is described in sections 5.1-2 of: - J.K. Dukoqicz, 2012. Reformulating the full-Stokes ice sheet model for a - more efficient computational solution. The Cryosphere, 6, 21-34. - www.the-cryosphere.net/6/21/2012/ +To see the full set of commmand line options, type 'python stabilitySlab.py -h'. -Blatter-Pattyn First-order solution is described in J.K. Dukowicz, manuscript -in preparation. +For questions, please contact Willian Lipscomb (lipscomb@ucar.edu) or Gunter Leguy (gunterl@ucar.edu). diff --git a/tests/slab/plotSlab.py b/tests/slab/plotSlab.py index 214c6531..6bfa7663 100755 --- a/tests/slab/plotSlab.py +++ b/tests/slab/plotSlab.py @@ -1,10 +1,9 @@ #!/usr/bin/env python2 - """ This script plots the results of an experiment with an ice "slab" on an inclined plane. Test case is described in sections 5.1-2 of: - J.K. Dukoqicz, 2012. Reformulating the full-Stokes ice sheet model for a + J.K. Dukowicz, 2012. Reformulating the full-Stokes ice sheet model for a more efficient computational solution. The Cryosphere, 6, 21-34. www.the-cryosphere.net/6/21/2012/ @@ -12,13 +11,12 @@ in preparation. """ #FIXME: Manuscript likely not in prep anymore -- JHK, 08/07/2015 +# Not published as of July 2021 -- WHL # Written by Matt Hoffman, Dec. 16, 2013 # Reconfigured by Joseph H Kennedy at ORNL on August 7, 2015 to work with the regression testing # NOTE: Did not adjust inner workings except where needed. - - -#NOTE: this script is assuming n=3, but more general solutions are available. +# Revised by William Lipscomb in 2021 to support more options, including general values of Glen's n. import os import sys @@ -28,8 +26,12 @@ import matplotlib.pyplot as plt from netCDF import * -from math import tan, pi, sin, cos -from runSlab import n, rhoi, grav, theta, beta, efvs, thickness # Get the values used to run the experiment +from math import tan, pi, sin, cos, atan + +# Get hard-coded parameters from the run script. +from runSlab import rhoi, grav + +import ConfigParser import argparse parser = argparse.ArgumentParser(description=__doc__, @@ -46,16 +48,15 @@ help="The tests output file you would like to plot. If a path is" \ +"passed via this option, the -o/--output-dir option will be ignored.") +parser.add_argument('-c','--config-file', + help="The configure file used to set up the test case and run CISM") + # =========================================================== # Define some variables and functions used in the main script # =========================================================== -# Calculate scales from Ducowicz unpub. man. -eta = beta * thickness * efvs**-n * (rhoi * grav * thickness)**(n-1) -velscale = (rhoi * grav * thickness / efvs)**n * thickness -thetar = theta * pi/180.0 # theta in radians - +#WHL args.output-file with a hyphen? def get_in_file(): if args.output_file: out_d, out_f = os.path.split(args.output_file) @@ -76,7 +77,7 @@ def get_in_file(): newest = max(matching, key=os.path.getmtime) print("\nWARNING: MULTIPLE *.out.nc FILES DETECTED!") print( "==========================================") - print( "Ploting the most recently modified file in the output directory:") + print( "Plotting the most recently modified file in the output directory:") print( " "+newest) print( "To plot another file, specify it with the -f/--outfile option.\n") @@ -94,6 +95,25 @@ def get_in_file(): return filein +def split_file_name(file_name): + """ + Get the root name, size, and number of processors from an out.nc filename. + #WHL - Adapted from plotISMIP_HOM.py + """ + root = '' + size = '' + proc = '' + + file_details = file_name.replace('.out.nc','') .split('.') +# print(file_details) +# print('len = ' + str(len(file_details))) + + if len(file_details) > 2: + proc = '.'+file_details[2] + size = '.'+file_details[1] + root = file_details[0] + + return (root, size, proc) # ========================= # Actual script starts here @@ -103,10 +123,7 @@ def main(): Plot the slab test results. """ - print("WARNING: THIS TEST CASE IS IN DEVELOPMENT. USE AT YOUR OWN RISK!") - - - filein = get_in_file() + filein = get_in_file() # Get needed variables from the output file x1 = filein.variables['x1'][:] @@ -120,28 +137,96 @@ def main(): # use integer floor division operator to get an index close to the center xp = len(x0)//2 yp = len(y0)//2 - #yp = 15 - #xp = 15 # ===================================================================== - print 'Using x index of '+str(xp)+'='+str(x0[xp]) - print 'Using y index of '+str(yp)+'='+str(y0[yp]) + + print('Using x index of '+str(xp)+'='+str(x0[xp])) + print('Using y index of '+str(yp)+'='+str(y0[yp])) thk = filein.variables['thk'][:] if netCDF_module == 'Scientific.IO.NetCDF': - thk = thk * filein.variables['thk'].scale_factor + thk = thk * filein.variables['thk'].scale_factor topg = filein.variables['topg'][:] if netCDF_module == 'Scientific.IO.NetCDF': - topg = topg * filein.variables['topg'].scale_factor + topg = topg * filein.variables['topg'].scale_factor uvel = filein.variables['uvel'][:] if netCDF_module == 'Scientific.IO.NetCDF': - uvel = uvel * filein.variables['uvel'].scale_factor - + uvel = uvel * filein.variables['uvel'].scale_factor + beta_2d = filein.variables['beta'][:] + if netCDF_module == 'Scientific.IO.NetCDF': + beta_2d = beta_2d * filein.variables['beta'].scale_factor + + # Get the name of the config file + # If not entered on the command line, then construct from the output filename + + if not args.config_file: + root, size, proc = split_file_name(args.output_file) + args.config_file = root + size + proc + '.config' + + configpath = os.path.join(args.output_dir, args.config_file) + print('configpath = ' + configpath) + + # Get gn and default_flwa from the config file + + try: + config_parser = ConfigParser.SafeConfigParser() + config_parser.read( configpath ) + + gn = float(config_parser.get('parameters','n_glen')) + flwa = float(config_parser.get('parameters', 'default_flwa')) + + except ConfigParser.Error as error: + print("Error parsing " + args.config ) + print(" "), + print(error) + sys.exit(1) + + # Derive the viscosity constant mu_n from flwa + # This expression is derived in the comments on flwa in runSlab.py. + mu_n = 1.0 / (2.0**((1.0+gn)/(2.0*gn)) * flwa**(1.0/gn)) + + # Get the ice thickness from the output file. + # If thickness = constant (i.e., the optional perturbation dh = 0), it does not matter where we sample. + # Note: In general, this thickness will differ from the baseline 'thk' that is used in runSlab.py + # to create the input file. + # This is because the baseline value is measured perpendicular to the sloped bed, + # whereas the CISM value is in the vertical direction, which is not perpendicular to the bed. + thickness = thk[0,yp,xp] + + # Get beta from the output file. + # Since beta = constant, it does not matter where we sample. + beta = beta_2d[0,yp,xp] + + # Derive theta from the output file as atan(slope(topg)) + # Since the slope is constant, it does not matter where we sample. + slope = (topg[0,yp,xp] - topg[0,yp,xp+1]) / (x0[xp+1] - x0[xp]) + thetar = atan(slope) + theta = thetar * 180.0/pi + + # Compute the dimensionless parameter eta and the velocity scale, + # which appear in the scaled velocity solution. + eta = (beta * thickness / mu_n**gn) * (rhoi * grav * thickness)**(gn-1) + velscale = (rhoi * grav * thickness / mu_n)**gn * thickness + + print('gn = ' + str(gn)) + print('rhoi = ' + str(rhoi)) + print('grav = ' + str(grav)) + print('thck = ' + str(thickness)) + print('mu_n = ' + str(mu_n)) + print('flwa = ' + str(flwa)) + print('beta = ' + str(beta)) + print('eta = ' + str(eta)) + print('theta= ' + str(theta)) + print('velscale = ' + str(velscale)) # === Plot the results at the given location === # Note we are not plotting like in Fig 3 of paper. # That figure plotted a profile against zprime. # It seemed more accurate to plot a profile against z to avoid interpolating model results (analytic solution can be calculated anywhere). - # Also, the analytic solution calculates the bed-parallel u velocity, but CISM calculates u as parallel to the geoid, so we need to transform the analytic solution to the CISM coordinate system. + # Also, the analytic solution calculates the bed-parallel u velocity, but CISM calculates u as parallel to the geoid, + # so we need to transform the analytic solution to the CISM coordinate system. + + #WHL - I think the analytic solution is actually for u(z'), which is not bed-parallel. + # The bed-parallel solution would be u'(z'), with w'(z') = 0. fig = plt.figure(1, facecolor='w', figsize=(12, 6)) @@ -151,24 +236,23 @@ def main(): x = (x0-x0[xp]) / thickness # calculate rotated zprime coordinates for this column (we assume the solution truly is spatially uniform) zprime = x[xp] * sin(thetar) + z * cos(thetar) - #print 'zprime', zprime # Calculate analytic solution for x-component of velocity (eq. 39 in paper) for the CISM-column - #uvelStokesAnalyticScaled = sin(theta * pi/180.0) * cos(theta * pi/180.0) * (0.5 * zprime**2 - zprime - 1.0/eta) - uvelStokesAnalyticScaled = (-1)**n * 2**((1.0-n)/2.0) * sin(thetar)**n * cos(thetar) / (n+1) \ - * ( (zprime - 1.0)**(n+1) - (-1.0)**(n+1) ) + sin(thetar) * cos(thetar) / eta + uvelStokesAnalyticScaled = sin(thetar) * cos(thetar) / eta \ + - 2**((1.0-gn)/2.0) * sin(thetar)**gn * cos(thetar) / (gn+1) * ( (1.0 - zprime)**(gn+1) - 1.0 ) - # Calculate the BP FO solution for x-component of velocity (Ducowicz, in prep. paper, Eq.30, n=3) - #uvelFOAnalyticScaled = (tan(theta * pi/180.0))**3 / (8.0 * (1.0 + 3.0 * (sin(theta * pi/180.0)**2))**2) \ - uvelFOAnalyticScaled = (-1)**n * 2**((1.0-n)/2.0) * tan(thetar)**n / \ - ( (n + 1) * (1.0 + 3.0 * sin(thetar)**2)**((n+1.0)/2.0) ) \ - * ( (zprime - 1.0)**(n+1) - (-1.0)**(n+1) ) + tan(thetar) / eta + # Calculate the BP FO solution for x-component of velocity (Dukowicz, in prep. paper, Eq.30, n=3) + uvelFOAnalyticScaled = + tan(thetar) / eta \ + - 2**((1.0-gn)/2.0) * tan(thetar)**gn / \ + ( (gn + 1) * (1.0 + 3.0 * sin(thetar)**2)**((gn+1.0)/2.0) ) \ + * ( (1.0 - zprime)**(gn+1) - 1.0 ) ### 1. Plot as nondimensional variables # Plot analytic solution fig.add_subplot(1,2,1) plt.plot(uvelStokesAnalyticScaled, z, '-kx', label='Analytic Stokes') plt.plot(uvelFOAnalyticScaled, z, '-ko', label='Analytic FO') + # Plot model results plt.plot(uvel[0,:,yp,xp] / velscale, z, '--ro', label='CISM') plt.ylim((-0.05, 1.05)) @@ -191,7 +275,16 @@ def main(): plt.title('Velocity profile at x=' + str(x0[xp]) + ' m, y=' + str(y0[yp]) + ' m\n(Unscaled coordinates)') ################# +# print('y0_min:') +# print(y0.min()) +# print('y0_max:') +# print(y0.max()) + # Now plot maps to show if the velocities vary over the domain (they should not) + # For some reason, the y-axis has a greater extent than the range (y0.min, y0.max). + #TODO - Fix the y-axis extent. Currently, the extent is too large for small values of ny. + #TODO - Plot the thickness relative to the initial thickness. + fig = plt.figure(2, facecolor='w', figsize=(12, 6)) fig.add_subplot(1,2,1) uvelDiff = uvel[0,0,:,:] - uvel[0,0,yp,xp] @@ -224,14 +317,11 @@ def main(): #plt.plot(level, tan(thetar)**3 / (8.0 * (1.0 + 3.0 * sin(thetar)**2)**2) * (1.0 - (level-1.0)**4 ) + tan(thetar)/eta, 'b--' , label='nonlinear fo') #plt.ylim((0.0, 0.04)); plt.xlabel("z'"); plt.ylabel('u'); plt.legend() - plt.draw() plt.show() filein.close() - print("WARNING: THIS TEST CASE IS IN DEVELOPMENT. USE AT YOUR OWN RISK!") - # Run only if this is being run as a script. if __name__=='__main__': @@ -240,4 +330,3 @@ def main(): # run the script sys.exit(main()) - diff --git a/tests/slab/runSlab.py b/tests/slab/runSlab.py index 2fc0217a..b6009ed5 100755 --- a/tests/slab/runSlab.py +++ b/tests/slab/runSlab.py @@ -1,6 +1,5 @@ #!/usr/bin/env python2 -#FIXME: More detailed description of this test case!!! """ Run an experiment with an ice "slab". """ @@ -8,10 +7,12 @@ # Authors # ------- # Modified from dome.py by Matt Hoffman, Dec. 16, 2013 -# Test case described in sections 5.1-2 of: -# J.K. Dukoqicz, 2012. Reformulating the full-Stokes ice sheet model for a -# more efficient computational solution. The Cryosphere, 6, 21-34. www.the-cryosphere.net/6/21/2012/ -# Reconfigured by Joseph H Kennedy at ORNL on April 27, 2015 to work with the regression testing +# Test case described in sections 5.1- 5.2 of: +# J.K. Dukowicz, 2012. Reformulating the full-Stokes ice sheet model for a +# more efficient computational solution. The Cryosphere, 6, 21-34, +# https://doi.org/10.5194/tc-6-21-2012. +# Reconfigured by Joseph H Kennedy at ORNL on April 27, 2015 to work with the regression testing. +# Revised by William Lipscomb in 2021 to support more options. import os import sys @@ -19,10 +20,10 @@ import subprocess import ConfigParser -import numpy +import numpy as np import netCDF -from math import sqrt, tan, pi, cos +from math import sqrt, sin, cos, tan, pi # Parse the command line options # ------------------------------ @@ -56,11 +57,36 @@ def unsigned_int(x): parser.add_argument('-s','--setup-only', action='store_true', help="Set up the test, but don't actually run it.") - -# Additional test specific options: -#parser.add_argument('--scale', type=unsigned_int, default=0, -# help="Scales the problem size by 2**SCALE. SCALE=0 creates a 31x31 grid, SCALE=1 " -# +"creates a 62x62 grid, and SCALE=2 creates a 124x124 grid.") +# Additional options to set grid, solver, physics parameters, etc.: +#Note: The default mu_n = 0.0 is not actually used. +# Rather, mu_n is computed below, unless mu_n > 0 is specified in the command line. +# For n = 1, the default is mu_1 = 1.0e6 Pa yr. +parser.add_argument('-a','--approx', default='BP', + help="Stokes approximation (SIA, SSA, BP, L1L2, DIVA)") +parser.add_argument('-beta','--beta', default=2000.0, + help="Friction parameter beta (Pa (m/yr)^{-1})") +parser.add_argument('-dh','--delta_thck', default=0.0, + help="Thickness perturbation (m)") +parser.add_argument('-dt','--tstep', default=0.01, + help="Time step (yr)") +parser.add_argument('-gn','--glen_exponent', default=1, + help="Exponent in Glen flow law") +parser.add_argument('-l','--levels', default=10, + help="Number of vertical levels") +parser.add_argument('-mu','--mu_n', default=0.0, + help="Viscosity parameter mu_n (Pa yr^{1/n})") +parser.add_argument('-nt','--n_tsteps', default=0, + help="Number of timesteps") +parser.add_argument('-nx','--nx_grid', default=50, + help="Number of grid cells in x direction") +parser.add_argument('-ny','--ny_grid', default=5, + help="Number of grid cells in y direction") +parser.add_argument('-r','--resolution', default=100.0, + help="Grid resolution (m)") +parser.add_argument('-theta','--theta', default=5.0, + help="Slope angle (deg)") +parser.add_argument('-thk','--thickness', default=1000.0, + help="Ice thickness (m)") # Some useful functions @@ -112,28 +138,11 @@ def prep_commands(args, config_name): return commands - -# Hard coded test specific parameters # ----------------------------------- -#FIXME: Some of these could just be options! - -# Physical parameters -n = 1 # flow law parameter - only the n=1 case is currently supported -# (implementing the n=3 case would probably require implementing a new efvs option in CISM) -rhoi = 910.0 # kg/m3 -grav = 9.1801 # m^2/s - -# Test case parameters -theta = 18 # basal inclination angle (degrees) unpub. man. uses example with theta=18 -thickness = 1000.0 # m thickness in the rotated coordinate system, not in CISM coordinates +# Hard-cosed test case parameters +rhoi = 917.0 # kg/m^3 +grav = 9.81 # m^2/s baseElevation = 1000.0 # arbitrary height to keep us well away from sea level - -efvs = 2336041.42829 # hardcoded in CISM for constant viscosity setting (2336041.42829 Pa yr) - -eta = 10.0 # unpub. man. uses example with eta=10.0 -beta = eta / thickness / efvs**-n / (rhoi * grav * thickness)**(n-1) # Pa yr m^-1 -# Note: Fig. 3 in Ducowicz (2013) uses eta=18, where eta=beta*H/efvs - # the main script function # ------------------------ @@ -142,24 +151,24 @@ def main(): Run the slab test. """ - print("WARNING: THIS TEST CASE IS IN DEVELOPMENT. USE AT YOUR OWN RISK!") - # check that file name modifier, if it exists, starts with a '-' if not (args.modifier == '') and not args.modifier.startswith('-') : args.modifier = '-'+args.modifier # get the configuration # --------------------- + + dx = float(args.resolution) + dy = dx + nx = int(args.nx_grid) + ny = int(args.ny_grid) + nz = int(args.levels) + dt = float(args.tstep) + tend = float(args.n_tsteps) * dt + try: config_parser = ConfigParser.SafeConfigParser() config_parser.read( args.config ) - - nz = int(config_parser.get('grid','upn')) - nx = int(config_parser.get('grid','ewn')) - ny = int(config_parser.get('grid','nsn')) - dx = float(config_parser.get('grid','dew')) - dy = float(config_parser.get('grid','dns')) - file_name = config_parser.get('CF input', 'name') root, ext = os.path.splitext(file_name) @@ -169,7 +178,8 @@ def main(): print(error) sys.exit(1) - res = str(nx).zfill(4) + res=str(int(dx)).zfill(5) # 00100 for 100m, 01000 for 1000m, etc. + if args.parallel > 0: mod = args.modifier+'.'+res+'.p'+str(args.parallel).zfill(3) else: @@ -180,32 +190,146 @@ def main(): out_name = root+mod+'.out'+ext - # Setup the domain + # Set up the domain # ---------------- - offset = 1.0 * float(nx)*dx * tan(theta * pi/180.0) - - # create the new config file + # Create the new config file # -------------------------- if not args.quiet: print("\nCreating config file: "+config_name) + # Set grid and time parameters config_parser.set('grid', 'upn', str(nz)) config_parser.set('grid', 'ewn', str(nx)) config_parser.set('grid', 'nsn', str(ny)) config_parser.set('grid', 'dew', str(dx)) config_parser.set('grid', 'dns', str(dy)) + config_parser.set('time', 'dt', str(dt)) + config_parser.set('time', 'tend',str(tend)) + + # Set physics parameters that are needed to create the config file and the input netCDF file. + # Note: rhoi and grav are hardwired above. + + # ice thickness + thickness = float(args.thickness) + + # friction parameter beta (Pa (m/yr)^{-1}) + beta = float(args.beta) + + # basal inclination angle (degrees) + theta = float(args.theta) + theta_rad = theta * pi/180.0 # convert to radians + + # flow law exponent + # Any value n >= 1 is supported. + gn = float(args.glen_exponent) + + # Note: Fig. 3 in Dukowicz (2012) uses eta = 18 and theta = 18 deg. + # This gives u(1) = 10.0 * u(0), where u(1) = usfc and u(0) = ubas. + + # viscosity coefficient mu_n, dependent on n (Pa yr^{1/n}) + # The nominal default is mu_n = 0.0, but this value is never used. + # If a nonzero value is specified on the command line, it is used; + # else, mu_n is computed here. The goal is to choose a value mu_n(n) that + # will result in vertical shear similar to a default case with n = 1 and mu_1, + # provided we have similar values of H and theta. + # In the Dukowicz unpublished manuscript, the viscosity mu is given by + # mu = mu_n * eps_e^[(1-n)/n], where eps_e is the effective strain rate. + # For n = 1, we choose a default value of 1.0e6 Pa yr. + # For n > 1, we choose mu_n (units of Pa yr^{1/n}) to match the surface velocity + # we would have with n = 1 and the same values of H and theta. + # The general velocity solution is + # u(z') = u_b + du(z') + # where u_b = rhoi * grav * sin(theta) * cos(theta) / beta + # and du(z') = 2^{(1-n)/2}/(n+1) * sin^n(theta) * cos(theta) + # * (rhoi*grav*H/mu_n)^n * H * [1 - (1 - z'/H)^{n+1}] + # For z' = H and general n, we have + # du_n(H) = 2^{(1-n)/2}/(n+1) * sin^n(theta) * cos(theta) + # * (rhoi*grav*H/mu_n)^n * H + # For z' = H and n = 1, we have + # du_1(H) = (1/2) * sin(theta) * cos(theta) * (rhoi*grav*H/mu_1) * H + # If we equate du_n(H) with du_1(H), we can solve for mu_n: + # mu_n = [ 2^{(3-n)/(2n)}/(n+1) * sin^{n-1}(theta) * (rhoi*grav*H)^{n-1} * mu_1 ]^{1/n} + # with units Pa yr^{1/n} + # This value should give nearly the same shearing velocity du(H) for exponent n > 1 + # as we would get for n = 1, given mu_1 and the same values of H and theta. + + if float(args.mu_n) > 0.0: + mu_n = float(args.mu_n) + mu_n_pwrn = mu_n**gn + else: + mu_1 = 1.0e6 # default value for mu_1 (Pa yr) + mu_n_pwrn = 2.0**((3.0-gn)/2.0)/(gn+1.0) * sin(theta_rad)**(gn-1.0) \ + * (rhoi*grav*thickness)**(gn-1.0) * mu_1 # (mu_n)^n + mu_n = mu_n_pwrn**(1.0/gn) + + # Given mu_n, compute the temperature-independent flow factor A = 1 / [2^((1+n)/2) * mu_n^n]. + # This is how CISM incorporates a prescribed mu_n (with flow_law = 0, i.e. constant flwa). + # Note: The complicated exponent of 2 in the denominator results from CISM and the Dukowicz papers + # having different conventions for the squared effective strain rate, eps_sq. + # In CISM: mu = 1/2 * A^(-1/n) * eps_sq_c^((1-n)/(2n)) + # where eps_sq_c = 1/2 * eps_ij * eps_ij + # eps_ij = strain rate tensor + # In Dukowicz: mu = mu_n * eps_sq_d^((1-n)/(2n)) + # where eps_sq_d = eps_ij * eps_ij = 2 * eps_sq_c + # Equating the two values of mu, we get mu_n * 2^((1-n)/(2n)) = 1/2 * A^(-1/n), + # which we solve to find A = 1 / [2^((1+n)/2) * mu_n^n] + # Conversely, we have mu_n = 1 / [2^((1+n)/(2n)) * A^(1/n)] + #TODO: Modify the Dukowicz derivations to use the same convention as CISM. + flwa = 1.0 / (2.0**((1.0+gn)/2.0) * mu_n_pwrn) + + # Find the dimensionless parameter eta + # This is diagnostic only; not used directly by CISM + eta = (beta * thickness / mu_n**gn) * (rhoi * grav * thickness)**(gn-1) + + # periodic offset; depends on theta and dx + offset = 1.0 * float(nx)*dx * tan(theta_rad) + + # Print some values + print('nx = ' + str(nx)) + print('ny = ' + str(ny)) + print('nz = ' + str(nz)) + print('dt = ' + str(dt)) + print('tend = ' + str(tend)) + print('rhoi = ' + str(rhoi)) + print('grav = ' + str(grav)) + print('thck = ' + str(thickness)) + print('beta = ' + str(beta)) + print('gn = ' + str(gn)) + print('mu_n = ' + str(mu_n)) + print('flwa = ' + str(flwa)) + print('eta = ' + str(eta)) + print('theta = ' + str(theta)) + print('offset = ' + str(offset)) + + # Write some options and parameters to the config file config_parser.set('parameters', 'periodic_offset_ew', str(offset)) - + config_parser.set('parameters', 'rhoi', str(rhoi)) + config_parser.set('parameters', 'grav', str(grav)) + config_parser.set('parameters', 'n_glen', str(gn)) + config_parser.set('parameters', 'default_flwa', str(flwa)) + + if (args.approx == 'SIA'): + approx = 0 + elif (args.approx == 'SSA'): + approx = 1 + elif (args.approx == 'BP'): + approx = 2 + elif (args.approx == 'L1L2'): + approx = 3 + elif (args.approx == 'DIVA'): + approx = 4 + config_parser.set('ho_options', 'which_ho_approx', str(approx)) + config_parser.set('CF input', 'name', file_name) config_parser.set('CF output', 'name', out_name) config_parser.set('CF output', 'xtype', 'double') - + config_parser.set('CF output', 'frequency', str(tend)) # write output at start and end of run + with open(config_name, 'wb') as config_file: config_parser.write(config_file) - # create the input netCDF file # ---------------------------- if not args.quiet: @@ -222,8 +346,8 @@ def main(): nc_file.createDimension('x0',nx-1) # staggered grid nc_file.createDimension('y0',ny-1) - x = dx*numpy.arange(nx,dtype='float32') - y = dx*numpy.arange(ny,dtype='float32') + x = dx*np.arange(nx,dtype='float32') + y = dx*np.arange(ny,dtype='float32') nc_file.createVariable('time','f',('time',))[:] = [0] nc_file.createVariable('x1','f',('x1',))[:] = x @@ -231,20 +355,49 @@ def main(): nc_file.createVariable('x0','f',('x0',))[:] = dx/2 + x[:-1] # staggered grid nc_file.createVariable('y0','f',('y0',))[:] = dy/2 + y[:-1] - # Calculate values for the required variables. - thk = numpy.zeros([1,ny,nx],dtype='float32') - topg = numpy.zeros([1,ny,nx],dtype='float32') - artm = numpy.zeros([1,ny,nx],dtype='float32') - unstagbeta = numpy.zeros([1,ny,nx],dtype='float32') + thk = np.zeros([1,ny,nx],dtype='float32') + topg = np.zeros([1,ny,nx],dtype='float32') + artm = np.zeros([1,ny,nx],dtype='float32') + unstagbeta = np.zeros([1,ny,nx],dtype='float32') # Calculate the geometry of the slab of ice - thk[:] = thickness / cos(theta * pi/180.0) + # Note: Thickness is divided by cos(theta), since thck in CISM is the vertical distance + # from bed to surface. On a slanted bed, this is greater than the distance measured + # in the direction perpendicular to the bed. + # Why does topg use a tan function? Is the bed slanted? + # Do we need unstagbeta instead of beta? Compare to ISMIP-HOM tests. + + thk[:] = thickness / cos(theta_rad) xmax = x[:].max() for i in range(nx): - topg[0,:,i] = (xmax - x[i]) * tan(theta * pi/180.0) + baseElevation + topg[0,:,i] = (xmax - x[i]) * tan(theta_rad) + baseElevation unstagbeta[:] = beta + # Optionally, add a small perturbation to the thickness field + + if args.delta_thck: + dh = float(args.delta_thck) + for i in range(nx): + + # Apply a Gaussian perturbation, using the Box-Mueller algorithm: + # https://en.wikipedia.org/wiki/Normal_distribution#Generating_values_from_normal_distribution + + mu = 0.0 # mean of normal distribution + sigma = 1.0 # stdev of normal distribution + + rnd1 = np.random.random() # two random numbers between 0 and 1 + rnd2 = np.random.random() + + # Either of the next two lines will sample a number at random from a normal distribution + rnd_normal = mu + sigma * sqrt(-2.0 * np.log(rnd1)) * cos(2.0*pi*rnd2) +# rnd_normal = mu + sigma * sqrt(-2.0 * np.log(rnd2)) * sin(2.0*pi*rnd1) + + dthk = dh * rnd_normal + thk[0,:,i] = thk[0,:,i] + dthk + print(i, dthk, thk[0,ny/2,i]) + thk_in = thk # for comparing later to final thk + # Create the required variables in the netCDF file. nc_file.createVariable('thk', 'f',('time','y1','x1'))[:] = thk nc_file.createVariable('topg','f',('time','y1','x1'))[:] = topg @@ -274,6 +427,8 @@ def main(): print("\nRunning CISM slab test") print( "======================\n") + print('command_list =' + str(command_list)) + process = subprocess.check_call(str.join("; ",command_list), shell=True) try: @@ -289,6 +444,7 @@ def main(): if not args.quiet: print("\nFinished running the CISM slab test") print( "===================================\n") + else: run_script = args.output_dir+os.sep+root+mod+".run" @@ -304,7 +460,6 @@ def main(): print( "======================================") print( " To run the test, use: "+run_script) - print("WARNING: THIS TEST CASE IS IN DEVELOPMENT. USE AT YOUR OWN RISK!") # Run only if this is being run as a script. if __name__=='__main__': @@ -314,4 +469,3 @@ def main(): # run the script sys.exit(main()) - diff --git a/tests/slab/slab.config b/tests/slab/slab.config index d9ffcd61..fbba9139 100644 --- a/tests/slab/slab.config +++ b/tests/slab/slab.config @@ -1,30 +1,34 @@ [grid] -upn = 50 +upn = 20 ewn = 30 -nsn = 20 +nsn = 5 dew = 50 dns = 50 [time] tstart = 0. tend = 0. -dt = 1. +dt = 0.01 +dt_diag = 0.01 +idiag = 15 +jdiag = 5 [options] -dycore = 2 # 1 = glam, 2 = glissade -flow_law = 0 # 0 = constant +dycore = 2 # 2 = glissade +flow_law = 0 # 0 = constant flwa (default = 1.e-16 Pa-n yr-1) evolution = 3 # 3 = remapping -temperature = 1 # 1 = prognostic, 3 = enthalpy +temperature = 1 # 1 = prognostic +basal_mass_balance = 0 # 0 = basal mbal not in continuity eqn [ho_options] which_ho_babc = 5 # 5 = externally-supplied beta(required by test case) -which_ho_efvs = 0 # 0 = constant (required by test case - makes n effectively 1) -which_ho_sparse = 3 # 1 = SLAP GMRES, 3 = glissade parallel PCG, 4 = Trilinos for linear solver +which_ho_sparse = 3 # 1 = SLAP GMRES, 3 = glissade parallel PCG which_ho_nonlinear = 0 # 0 = Picard, 1 = accelerated Picard +which_ho_approx = 4 # 2 = BP, 3 = L1L2, 4 = DIVA [parameters] ice_limit = 1. # min thickness (m) for dynamics -periodic_offset_ew = 487.379544349 +geothermal = 0. [CF default] comment = created with slab.py diff --git a/tests/slab/stabilitySlab.py b/tests/slab/stabilitySlab.py new file mode 100644 index 00000000..5529896a --- /dev/null +++ b/tests/slab/stabilitySlab.py @@ -0,0 +1,387 @@ +#!/usr/bin/env python2 +# -*- coding: utf-8 -*- + +""" +This script runs a series of CISM experiments at different resolutions. +At each resolution, it determines the maximum stable time step. +A run is deemed to be stable if the standard deviation of a small thickness perturbation +decreases during a transient run (100 timesteps by default). + +Used to obtain the CISM stability results described in: +Robinson, A., D. Goldberg, and W. H. Lipscomb, A comparison of the performance +of depth-integrated ice-dynamics solvers, to be submitted. +""" + +# Authors +# ------- +# Created by William Lipscomb, July 2021 + +import os +import sys +import errno +import subprocess +import ConfigParser + +import numpy as np +import netCDF +from math import sqrt, log10 + +# Parse the command line options +# ------------------------------ +import argparse +parser = argparse.ArgumentParser(description=__doc__, + formatter_class=argparse.ArgumentDefaultsHelpFormatter) + +# small helper function so argparse will understand unsigned integers +def unsigned_int(x): + x = int(x) + if x < 1: + raise argparse.ArgumentTypeError("This argument is an unsigned int type! Should be an integer greater than zero.") + return x + +# The following command line arguments determine the set of resolutions to run the slab test. +# At each resolution, we aim to find the maximum stable time step. +# Note: If args.n_resolution > 1, then args.resolution (see below) is ignored. + +parser.add_argument('-nr','--n_resolution', default=1, + help="number of resolutions") +parser.add_argument('-rmin','--min_resolution', default=10.0, + help="minimum resolution (m)") +parser.add_argument('-rmax','--max_resolution', default=40000.0, + help="minimum resolution (m)") + +# The following command line arguments are the same as in runSlab.py. +# Not sure how to avoid code repetition. + +parser.add_argument('-c','--config', default='./slab.config', + help="The configure file used to setup the test case and run CISM") +parser.add_argument('-e','--executable', default='./cism_driver', + help="The CISM driver") +parser.add_argument('-m', '--modifier', metavar='MOD', default='', + help="Add a modifier to file names. FILE.EX will become FILE.MOD.EX") +parser.add_argument('-n','--parallel', metavar='N', type=unsigned_int, default=0, + help="Run in parallel using N processors.") +parser.add_argument('-o', '--output_dir',default='./output', + help="Write all created files here.") +parser.add_argument('-a','--approx', default='BP', + help="Stokes approximation (SIA, SSA, BP, L1L2, DIVA)") +parser.add_argument('-beta','--beta', default=2000.0, + help="Friction parameter beta (Pa (m/yr)^{-1})") +parser.add_argument('-dh','--delta_thck', default=0.0, + help="Thickness perturbation (m)") +parser.add_argument('-dt','--tstep', default=0.01, + help="Time step (yr)") +parser.add_argument('-gn','--glen_exponent', default=1, + help="Exponent in Glen flow law") +parser.add_argument('-l','--levels', default=10, + help="Number of vertical levels") +parser.add_argument('-mu','--mu_n', default=0.0, + help="Viscosity parameter mu_n (Pa yr^{1/n})") +parser.add_argument('-nt','--n_tsteps', default=0, + help="Number of timesteps") +parser.add_argument('-nx','--nx_grid', default=50, + help="Number of grid cells in x direction") +parser.add_argument('-ny','--ny_grid', default=5, + help="Number of grid cells in y direction") +parser.add_argument('-r','--resolution', default=100.0, + help="Grid resolution (m)") +parser.add_argument('-theta','--theta', default=5.0, + help="Slope angle (deg)") +parser.add_argument('-thk','--thickness', default=1000.0, + help="Ice thickness") + + ############ + # Functions # + ############ + +def reading_file(inputFile): + + #Check whether a netCDF file exists, and return a list of times in the file + + ReadVarFile = True + try: + filein = netCDF.NetCDFFile(inputFile,'r') + time = filein.variables['time'][:] + filein.close() + print('Was able to read file ' + inputFile) + print(time) + except: + ReadVarFile = False + time = [0.] + print('Was not able to read file' + inputFile) + + return time, ReadVarFile + + +def check_output_file(outputFile, time_end): + + # Check that the output file exists with the expected final time slice + + # Path to experiment + path_to_slab_output = './output/' + + # File to check + filename = path_to_slab_output + outputFile + + # Read the output file + print('Reading file ' + str(filename)) + time_var, VarRead = reading_file(filename) + +# print(time_var) + + # Checking that the last time entry is the same as we expect from time_end + # Allow for a small roundoff difference. + if abs(time_var[-1] - time_end) < 1.0e-7: + check_time_var = True + else: + check_time_var = False + + print('time_end = ' + str(time_end)) + print('last time in file = ' + str(time_var[-1])) + + # Creating the status of both checks + check_passed = check_time_var and VarRead + + if check_passed: + print('Found output file with expected file time slice') + else: + if (not VarRead): + print('Output file cannot be read') + else: + if not check_time_var: + print('Output file is missing time slices') + + return check_passed + + +def main(): + + print('In main') + + """ + For each of several values of the horizontal grid resolution, determine the maximum + stable time step for a given configuration of the slab test. + """ + + resolution = [] + + # Based on the input arguments, make a list of resolutions at which to run the test. + # The formula and the default values of rmin and rmax give resolutions agreeing with + # those used by Alex Robinson for Yelmo, for the case nres = 12: + # resolution = [10., 21., 45., 96., 204., 434., 922., 1960., 4170., 8850., 18800., 40000.] + + print('Computing resolutions') + print(args.n_resolution) + if int(args.n_resolution) > 1: + nres = int(args.n_resolution) + resolution = [0. for n in range(nres)] + rmin = float(args.min_resolution) + rmax = float(args.max_resolution) + for n in range(nres): + res = 10.0**(log10(rmin) + (log10(rmax) - log10(rmin))*float(n)/float(nres-1)) + # Round to 3 significant figures (works for log10(res) < 5) + if log10(res) > 4.: + resolution[n] = round(res, -2) + elif log10(res) > 3.: + resolution[n] = round(res, -1) + else: + resolution[n] = round(res) + else: + nres = 1 + resolution.append(float(args.resolution)) + + print('nres = ' + str(nres)) + print(resolution) + + # Create an array to store max time step for each resolution + rows, cols = (nres, 2) + res_tstep = [[0. for i in range(cols)] for j in range(rows)] + for n in range(nres): + res_tstep[n][0] = resolution[n] + + for n in range(nres): + + print('output_dir: ' + args.output_dir) + + # Construct the command for calling the main runSlab script + run_command = 'python runSlab.py' + run_command = run_command + ' -c ' + args.config + run_command = run_command + ' -e ' + args.executable + if args.parallel > 0: + run_command = run_command + ' -n ' + str(args.parallel) + run_command = run_command + ' -o ' + args.output_dir + run_command = run_command + ' -a ' + args.approx + run_command = run_command + ' -beta ' + str(args.beta) + run_command = run_command + ' -dh ' + str(args.delta_thck) + run_command = run_command + ' -gn ' + str(args.glen_exponent) + run_command = run_command + ' -l ' + str(args.levels) + run_command = run_command + ' -mu ' + str(args.mu_n) + run_command = run_command + ' -nt ' + str(args.n_tsteps) + run_command = run_command + ' -nx ' + str(args.nx_grid) + run_command = run_command + ' -ny ' + str(args.ny_grid) + run_command = run_command + ' -theta '+ str(args.theta) + run_command = run_command + ' -thk '+ str(args.thickness) + + tend = float(args.n_tsteps) * args.tstep + + res = resolution[n] + run_command = run_command + ' -r ' + str(res) + + # Choose the time step. + # Start by choosing a very small timestep that can be assumed stable + # and a large step that can be assumed unstable. + # Note: SIA-type solvers at 10m resolution can require dt <~ 1.e-6 yr. + + tstep_lo = 1.0e-7 + tstep_hi = 1.0e+5 + tstep_log_precision = 1.0e-4 + print('Initial tstep_lo = ' + str(tstep_lo)) + print('Initial tstep_hi = ' + str(tstep_hi)) + print('Log precision = ' + str(tstep_log_precision)) + + while (log10(tstep_hi) - log10(tstep_lo)) > tstep_log_precision: + + # Compute the time step as the geometric mean of the tstep_lo and tstep_hi. + # tstep_lo is the largest time step known to be stable. + # tstep_hi is the smallest time step known to be unstable. + + tstep = sqrt(tstep_lo*tstep_hi) + + run_command_full = run_command + ' -dt ' + str(tstep) + + print("\nRunning CISM slab test...") + print('resolution = ' + str(res)) + print('tstep = ' + str(tstep)) + print('run_command = ' + run_command_full) + + process = subprocess.check_call(run_command_full, shell=True) + + print("\nFinished running the CISM slab test") + + # Determine the name of the output file. + # Must agree with naming conventions in runSlab.py + + file_name = args.config + root, ext = os.path.splitext(file_name) + + res=str(int(res)).zfill(5) # 00100 for 100m, 01000 for 1000m, etc. + + if args.parallel > 0: + mod = args.modifier + '.' + res + '.p' + str(args.parallel).zfill(3) + else: + mod = args.modifier + '.' + res + + outputFile = root + mod + '.out.nc' + + # Check whether the output file exists with the desired final time slice. + + time_end = float(args.n_tsteps) * tstep + + print('outputFile = ' + str(outputFile)) + print('n_tsteps = ' + str(float(args.n_tsteps))) + print('tstep = ' + str(tstep)) + print('time_end = ' + str(time_end)) + + check_passed = check_output_file(outputFile, time_end) + + if check_passed: + + print('Compute stdev of initial and final thickness for j = ny/2') + nx = int(args.nx_grid) + ny = int(args.ny_grid) + + # Read initial and final thickness from output file + outpath = os.path.join(args.output_dir, outputFile) + print('outpath = ' + outpath) + filein = netCDF.NetCDFFile(outpath,'r') + thk = filein.variables['thk'][:] + + j = ny/2 + thk_in = thk[0,j,:] + thk_out = thk[1,j,:] + + # Compute + Hav_in = 0.0 + Hav_out = 0.0 + for i in range(nx): + Hav_in = Hav_in + thk_in[i] + Hav_out = Hav_out + thk_out[i] + Hav_in = Hav_in / nx + Hav_out = Hav_out / nx + + # Compute + H2av_in = 0.0 + H2av_out = 0.0 + for i in range(nx): + H2av_in = H2av_in + thk_in[i]**2 + H2av_out = H2av_out + thk_out[i]**2 + H2av_in = H2av_in / nx + H2av_out = H2av_out / nx + + print('H2av_out =' + str(H2av_out)) + print('Hav_out^2 =' + str(Hav_out**2)) + + # Compute stdev = sqrt( - ^2) + var_in = H2av_in - Hav_in**2 + var_out = H2av_out - Hav_out**2 + + if var_in > 0.: + stdev_in = sqrt(H2av_in - Hav_in**2) + else: + stdev_in = 0. + + if var_out > 0.: + stdev_out = sqrt(H2av_out - Hav_out**2) + else: + stdev_out = 0. + + if stdev_in > 0.: + ratio = stdev_out/stdev_in + else: + ratio = 0. + + print('stdev_in = ' + str(stdev_in)) + print('stdev_out = ' + str(stdev_out)) + print('ratio = ' + str(ratio)) + + # Determine whether the run was stable. + # A run is defined to be stable if the final standard deviation of thickness + # is less than the initial standard deviation + + if ratio < 1.: + tstep_lo = max(tstep_lo, tstep) + print('Stable, new tstep_lo =' + str(tstep_lo)) + else: + tstep_hi = min(tstep_hi, tstep) + print('Unstable, new tstep_hi =' + str(tstep_hi)) + + else: # check_passed = F; not stable + tstep_hi = min(tstep_hi, tstep) + print('Unstable, new tstep_hi =' + str(tstep_hi)) + + print('Latest tstep_lo = ' + str(tstep_lo)) + print('Latest tstep_hi = ' + str(tstep_hi)) + + # Add to the array containing the max stable timestep at each resolution. + # Take the max stable timestep to be the average of tstep_lo and tstep_hi. + res_tstep[n][1] = 0.5 * (tstep_lo + tstep_hi) + + print('New res_tstep, res #' + str(n)) + print(res_tstep) + + # Print a table containing the max timestep for each resolution + for n in range(nres): + float_res = res_tstep[n][0] + float_dt = res_tstep[n][1] + formatted_float_res = "{:8.1f}".format(float_res) + formatted_float_dt = "{:.3e}".format(float_dt) # exponential notation with 3 decimal places + print(formatted_float_res + ' ' + formatted_float_dt) + +# Run only if this is being run as a script. +if __name__=='__main__': + + # get the command line arguments + args = parser.parse_args() + + # run the script + sys.exit(main())