From 165fb8a98de9ecb750124d6829a654e806cc5c59 Mon Sep 17 00:00:00 2001 From: William Lipscomb Date: Tue, 20 Apr 2021 17:51:26 -0600 Subject: [PATCH 1/5] Made the Glen flow-law exponent 'n' a config variable Until now, the exponent n in the Glen flow law has been set in glimmer_physcon.F90 as an integer parameter, gn = 3. With this commit, n is renamed n_glen (a more greppable name) for use in Glissade. It is declared in glimmer_physcon.F90 as real(dp) with default value 3.0d0. Since n_glen is no longer a parameter, it can now be set in the config file like other physical 'constants' (e.g., rhoi and rhoo) that are not truly constant, but can take different values in different models or benchmarking experiments. To avoid changing answers in old Glide code, I retained the integer parameter gn = 3 in glimmer_physcon.F90. This parameter is now used only in the Glide dycore (e.g., glide_velo.F90). In Glissade, I replaced gn with n_glen in several places: (1) In subroutine glissade_interior_dissipation_sia, the dissipation factor includes n_glen. Note: n_glen does not appear explicitly in the 1st-order dissipation, which is proportional to tau_eff^2/efvs. (2) In glissade_velo_sia, n_glen appears in the vertical integral for the velocity. (3) In glissade_velo_higher, flwafact = flwa**(-1./n_glen) where flwa = A. (4) In glissade_velo_higher, the exponent p_effstr = (1.d0 - n_glen)/n_glen is used to compute effective_viscosity for BP, DIVA, or L1L2. (5) In glissade_velo_higher, subroutine compute_3d_velocity_l1l2 has a factor proportional to tau**((n_glen-1.)/2.) in the vertical integral. For (1) and (2), n_glen was previously assumed to be an integer. Switching it to real(dp) is answer-changing at the machine roundoff level for runs using the local SIA solver in glissade_velo_sia.F90. For (3), (4) and (5), n_glen replaces the equivalent expression real(gn,dp). For this reason, answers are BFB when using the BP, DIVA or L1L2 solver. In glissade_basal_traction, n_glen appeared in the expression for beta in the Coulomb sliding option, HO_BABC_COULOMB_FRICTION. Here, I replaced n_glen with powerlaw_m (also = 3.0d0 by default) to be consistent with the expressions for beta in the Schoof and Tsai laws. In glimmer_scales, Glen's n appears in the expressions for the scaling parameters vis0 and vis_scale, Here, I defined a local integer parameter gn = 3 so that the scales are unchanged. It is now possible for the user to specify arbitrary values of n_glen in tests such as the slab test. Another minor change: I added some logic to the subroutine that computes L1L2 velocities. For which_ho_efvs = 2 = HO_EFVS_NONLINEAR, the effective viscosity (efvs) is computed from the effective strain rate using an equation from Perego et al. (2012). But for option 0 (efvs = constant) and option 1 (efvs = multiple of flow factor), this strain rate equation in the code does not apply. For these options, I added an alternative that computes velocity in terms of the strain-rate-independent efvs. This allows us to use L1L2 for problems with constant efvs (e.g., the slab test). --- libglide/felix_dycore_interface.F90 | 2 +- libglide/glide_setup.F90 | 10 ++- libglimmer/glimmer_paramets.F90 | 3 +- libglimmer/glimmer_physcon.F90 | 8 ++- libglimmer/glimmer_scales.F90 | 2 +- libglissade/glissade_basal_traction.F90 | 9 ++- libglissade/glissade_therm.F90 | 14 ++-- libglissade/glissade_velo.F90 | 3 - libglissade/glissade_velo_higher.F90 | 85 ++++++++++++++++--------- libglissade/glissade_velo_sia.F90 | 13 ++-- 10 files changed, 93 insertions(+), 56 deletions(-) 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..909244fb 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 @@ -1996,7 +1996,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 +2021,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) @@ -2206,7 +2207,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 +2372,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/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_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..f6364650 100644 --- a/libglissade/glissade_therm.F90 +++ b/libglissade/glissade_therm.F90 @@ -1640,10 +1640,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 +2251,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 +2267,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 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..e9845437 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 @@ -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 From 44f771c90c2df1df1c19632d7d05d86b20abc726 Mon Sep 17 00:00:00 2001 From: William Lipscomb Date: Wed, 21 Apr 2021 10:37:52 -0600 Subject: [PATCH 2/5] Do not call init_isostasy unless isostasy = 1 The code was calling subroutine init_isostasy with isostasy = 0 = ISOSTASY_NONE. This subroutine is now called only if isostasy = 1 = ISOSTASY_COMPUTE. --- libglissade/glissade.F90 | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/libglissade/glissade.F90 b/libglissade/glissade.F90 index ef0dc5dd..4093a137 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) From 8614ae3dede10932afd039c5e912423930ba1a93 Mon Sep 17 00:00:00 2001 From: William Lipscomb Date: Sat, 31 Jul 2021 12:17:28 -0600 Subject: [PATCH 3/5] Minor code changes to support the slab test I modified glissade.F90 to abort cleanly with a call to glide_finalise after an advective CFL error. This is done only when the user does *not* specify adaptive subcycling. The clean abort allows the new slabStability script to keep going, launching a new run with a shorter timestep. In subroutine glissade_flow_factor of glissade_therm.F90, I removed the FLWA_INPUT option (option 3 of whichflwa). This option is redundant with option 0, FLWA_CONST_FLWA, in which the user can set default_flwa in the parameters section of the config file, and otherwise CISM uses default_flwa = 1.0e-16 Pa^-n yr^-1. We probably should rename default_flwa to constant_flwa, but leaving it for now to avoid breaking config files in test cases. Note: This option was added by Matt Hoffman in 2014. I am unaware of test cases that use this option (most have flow_law = 0 or 2), but if there are any, we will need to fix them by switching to whichflwa = 0. In subroutine glissade_therm_driver of glissade_therm.F90, I increased the threshold for column energy conservation errors from 1.0d-8 to 1.0d-7 W/m^2. For very small timesteps of ~1.0e-6 yr, the smaller threshold can be exceeded because of roundoff errors. In subroutine glissade_check_cfl of glissade_transport.F90, I modified the fatal abort for large CFL violations (advective CFL > 10). Now, CISM aborts for CFL > 10 only when adaptive_cfl_threshold > 0, i.e. transport subcycling is enabled. This prevents excessive subcycling for egregious CFL violations. If adaptive_cfl_threshold = 0. (the default), i.e. transport subcycling is not enabled, then the code exits cleanly (with a call to glide_finalise) in glissade.F90 when advective CFL > 1. I added a TODO note to set the default value of geot (the geothermal heat flux) to 0 instead of 0.05 W/m^2. With the default value, simple prognostic tests like the dome are not mass-conserving. Not changing just yet because answers will change for the dome test. --- libglide/glide_setup.F90 | 9 +++---- libglide/glide_types.F90 | 3 +-- libglissade/glissade.F90 | 40 +++++++++++++++++++++------- libglissade/glissade_therm.F90 | 37 +++++++++++-------------- libglissade/glissade_transport.F90 | 20 +++++++++----- libglissade/glissade_velo_higher.F90 | 4 +-- 6 files changed, 66 insertions(+), 47 deletions(-) diff --git a/libglide/glide_setup.F90 b/libglide/glide_setup.F90 index 909244fb..c5f4cc3d 100644 --- a/libglide/glide_setup.F90 +++ b/libglide/glide_setup.F90 @@ -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 = (/ & @@ -2034,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) 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/libglissade/glissade.F90 b/libglissade/glissade.F90 index 4093a137..efa2d483 100644 --- a/libglissade/glissade.F90 +++ b/libglissade/glissade.F90 @@ -1899,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 @@ -1977,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 @@ -2165,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 @@ -2192,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_therm.F90 b/libglissade/glissade_therm.F90 index f6364650..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) @@ -2416,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) @@ -2490,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) @@ -2560,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_higher.F90 b/libglissade/glissade_velo_higher.F90 index e9845437..77b2f9f5 100644 --- a/libglissade/glissade_velo_higher.F90 +++ b/libglissade/glissade_velo_higher.F90 @@ -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. From 840c406071abc6df939de3655eab35abe63a6413 Mon Sep 17 00:00:00 2001 From: William Lipscomb Date: Sat, 31 Jul 2021 13:25:12 -0600 Subject: [PATCH 4/5] Debugged and extended the Dukowicz slab test case This commit modifies the run and plot scripts for the Dukowicz slab test case, as described in Section 5 of this paper: 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. 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 exponent 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). These solutions can be compared to those simulated by CISM. The original scripts, runSlab.py and plotSlab.py, were written by Matt Hoffman with support for 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 value of n >= 1 (not necessarily an integer). The tests assume which_ho_efvs = 2 (nonlinear viscosity) with flow_law = 0 (constant). * 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. 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. (Note: 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.) I modified the plotting script, plotSlab.py, to look in the config and output files for physics parameters that are no longer hardwired. I slightly modified the analytic formulas to give the correct solution for non-integer n. This script creates two plots. The first plot shows 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. I also added a new script, stabilitySlab.py, to carry out stability tests as described in: Robinson, A., D. Goldberg, and W. H. Lipscomb, A comparison of the performance of depth-integrated ice-dynamics solvers, to be submitted to The Cryosphere. The idea is that 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). I have run the stability script for several solvers (DIVA, L1L2, SIA, SSA) for each of two physical cases: one with thick shearing ice and one with thin sliding ice. Each suite of experiments runs in a few minutes on 4 Macbook cores for solvers other than BP. As expected, DIVA and SSA are much more stable than L1L2 and SIA. This and the previous commit correspond to the CISM code and scripts used for the initial submission by Robinson et al. (2021). --- tests/slab/plotSlab.py | 167 ++++++++++++---- tests/slab/runSlab.py | 274 +++++++++++++++++++------ tests/slab/slab.config | 22 +- tests/slab/stabilitySlab.py | 387 ++++++++++++++++++++++++++++++++++++ 4 files changed, 742 insertions(+), 108 deletions(-) create mode 100644 tests/slab/stabilitySlab.py 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()) From bd9b1f4794476d8fab4881d81a84d1c2f2028729 Mon Sep 17 00:00:00 2001 From: William Lipscomb Date: Thu, 5 Aug 2021 15:16:08 -0600 Subject: [PATCH 5/5] Updated the slab README file Rewrote the slab README file to describe the new command line options for runSlab.py, and the new script stabilitySlab.py. --- tests/slab/README.md | 90 +++++++++++++++++++++++++++++++++++++++----- 1 file changed, 80 insertions(+), 10 deletions(-) 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).