diff --git a/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 b/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 index 5e14d9686..23251b2d1 100755 --- a/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 @@ -63,8 +63,11 @@ module ice_dyn_shared real (kind=dbl_kind), public :: & revp , & ! 0 for classic EVP, 1 for revised EVP - e_ratio , & ! e = EVP ellipse aspect ratio - ecci , & ! 1/e^2 + e_yieldcurve, & ! VP aspect ratio of elliptical yield curve + e_plasticpot, & ! VP aspect ratio of elliptical plastic potential + epp2i , & ! 1/(e_plasticpot)^2 + e_factor , & ! (e_yieldcurve)^2/(e_plasticpot)^4 + ecci , & ! temporary for 1d evp dtei , & ! 1/dte, where dte is subcycling timestep (1/s) ! dte2T , & ! dte/2T denom1 ! constants for stress equation @@ -220,7 +223,6 @@ subroutine set_evp_parameters (dt) !real (kind=dbl_kind) :: & !dte , & ! subcycling timestep for EVP dynamics, s - !ecc , & ! (ratio of major to minor ellipse axes)^2 !tdamp2 ! 2*(wave damping time scale T) character(len=*), parameter :: subname = '(set_evp_parameters)' @@ -230,10 +232,10 @@ subroutine set_evp_parameters (dt) !dtei = c1/dte ! 1/s dtei = real(ndte,kind=dbl_kind)/dt - ! major/minor axis length ratio, squared - !ecc = e_ratio**2 - !ecci = c1/ecc ! 1/ecc - ecci = c1/e_ratio**2 ! 1/ecc + ! variables for elliptical yield curve and plastic potential + epp2i = c1/e_plasticpot**2 + e_factor = e_yieldcurve**2 / e_plasticpot**4 + ecci = c1/e_yieldcurve**2 ! temporary for 1d evp ! constants for stress equation !tdamp2 = c2*eyc*dt ! s @@ -1352,10 +1354,10 @@ subroutine strain_rates (nx_block, ny_block, & - cxp(i,j)*uvel(i ,j-1) + dxt(i,j)*uvel(i ,j ) ! Delta (in the denominator of zeta, eta) - Deltane = sqrt(divune**2 + ecci*(tensionne**2 + shearne**2)) - Deltanw = sqrt(divunw**2 + ecci*(tensionnw**2 + shearnw**2)) - Deltasw = sqrt(divusw**2 + ecci*(tensionsw**2 + shearsw**2)) - Deltase = sqrt(divuse**2 + ecci*(tensionse**2 + shearse**2)) + Deltane = sqrt(divune**2 + e_factor*(tensionne**2 + shearne**2)) + Deltanw = sqrt(divunw**2 + e_factor*(tensionnw**2 + shearnw**2)) + Deltasw = sqrt(divusw**2 + e_factor*(tensionsw**2 + shearsw**2)) + Deltase = sqrt(divuse**2 + e_factor*(tensionse**2 + shearse**2)) end subroutine strain_rates @@ -1419,19 +1421,19 @@ subroutine viscous_coeffs_and_rep_pressure (strength, tinyarea, & zetax2ne = (c1+Ktens)*tmpcalcne ! northeast rep_prsne = (c1-Ktens)*tmpcalcne*Deltane - etax2ne = ecci*zetax2ne ! CHANGE FOR e_plasticpot + etax2ne = epp2i*zetax2ne zetax2nw = (c1+Ktens)*tmpcalcnw ! northwest rep_prsnw = (c1-Ktens)*tmpcalcnw*Deltanw - etax2nw = ecci*zetax2nw ! CHANGE FOR e_plasticpot + etax2nw = epp2i*zetax2nw zetax2sw = (c1+Ktens)*tmpcalcsw ! southwest rep_prssw = (c1-Ktens)*tmpcalcsw*Deltasw - etax2sw = ecci*zetax2sw ! CHANGE FOR e_plasticpot + etax2sw = epp2i*zetax2sw zetax2se = (c1+Ktens)*tmpcalcse ! southeast rep_prsse = (c1-Ktens)*tmpcalcse*Deltase - etax2se = ecci*zetax2se ! CHANGE FOR e_plasticpot + etax2se = epp2i*zetax2se ! else diff --git a/cicecore/cicedynB/dynamics/ice_dyn_vp.F90 b/cicecore/cicedynB/dynamics/ice_dyn_vp.F90 index 1a6c68548..2f1285084 100644 --- a/cicecore/cicedynB/dynamics/ice_dyn_vp.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_vp.F90 @@ -46,7 +46,7 @@ module ice_dyn_vp use ice_domain, only: nblocks, distrb_info use ice_domain_size, only: max_blocks use ice_dyn_shared, only: dyn_prep1, dyn_prep2, dyn_finish, & - ecci, cosw, sinw, fcor_blk, uvel_init, vvel_init, & + cosw, sinw, fcor_blk, uvel_init, vvel_init, & seabed_stress_factor_LKD, seabed_stress_factor_prob, seabed_stress_method, & seabed_stress, Ktens, stack_velocity_field, unstack_velocity_field use ice_fileunits, only: nu_diag @@ -1320,6 +1320,9 @@ end subroutine calc_zeta_dPr ! Computes the VP stresses (as diagnostic) +! Lemieux, J.-F., and Dupont, F. (2020), On the calculation of normalized +! viscous-plastic sea ice stresses, Geosci. Model Dev., 13, 1763–1769, + subroutine stress_vp (nx_block , ny_block , & icellt , & indxti , indxtj , & diff --git a/cicecore/cicedynB/general/ice_init.F90 b/cicecore/cicedynB/general/ice_init.F90 index 2e67af51c..b299ef77f 100644 --- a/cicecore/cicedynB/general/ice_init.F90 +++ b/cicecore/cicedynB/general/ice_init.F90 @@ -105,9 +105,10 @@ subroutine input_data use ice_dyn_shared, only: ndte, kdyn, revised_evp, yield_curve, & evp_algorithm, & seabed_stress, seabed_stress_method, & - k1, k2, alphab, threshold_hw, & - Ktens, e_ratio, coriolis, ssh_stress, & - kridge, brlx, arlx + k1, k2, alphab, threshold_hw, Ktens, & + e_yieldcurve, e_plasticpot, coriolis, & + ssh_stress, kridge, brlx, arlx + use ice_dyn_vp, only: maxits_nonlin, precond, dim_fgmres, dim_pgmres, maxits_fgmres, & maxits_pgmres, monitor_nonlin, monitor_fgmres, & monitor_pgmres, reltol_nonlin, reltol_fgmres, reltol_pgmres, & @@ -208,20 +209,20 @@ subroutine input_data namelist /dynamics_nml/ & kdyn, ndte, revised_evp, yield_curve, & - evp_algorithm, & + evp_algorithm, & brlx, arlx, ssh_stress, & advection, coriolis, kridge, ktransport, & kstrength, krdg_partic, krdg_redist, mu_rdg, & - e_ratio, Ktens, Cf, seabed_stress, & - k1, maxits_nonlin, precond, dim_fgmres, & + e_yieldcurve, e_plasticpot, Ktens, & + maxits_nonlin, precond, dim_fgmres, & dim_pgmres, maxits_fgmres, maxits_pgmres, monitor_nonlin, & monitor_fgmres, monitor_pgmres, reltol_nonlin, reltol_fgmres, & reltol_pgmres, algo_nonlin, dim_andacc, reltol_andacc, & damping_andacc, start_andacc, fpfunc_andacc, use_mean_vrel, & - ortho_type, & - k2, alphab, threshold_hw, & - seabed_stress_method, Pstar, Cstar - + ortho_type, seabed_stress, seabed_stress_method, & + k1, k2, alphab, threshold_hw, & + Cf, Pstar, Cstar + namelist /shortwave_nml/ & shortwave, albedo_type, & albicev, albicei, albsnowv, albsnowi, & @@ -367,7 +368,8 @@ subroutine input_data alphab = 20.0_dbl_kind ! alphab=Cb factor in Lemieux et al 2015 threshold_hw = 30.0_dbl_kind ! max water depth for grounding Ktens = 0.0_dbl_kind ! T=Ktens*P (tensile strength: see Konig and Holland, 2010) - e_ratio = 2.0_dbl_kind ! VP ellipse aspect ratio + e_yieldcurve = 2.0_dbl_kind ! VP aspect ratio of elliptical yield curve + e_plasticpot = 2.0_dbl_kind ! VP aspect ratio of elliptical plastic potential maxits_nonlin = 4 ! max nb of iteration for nonlinear solver precond = 'pgmres' ! preconditioner for fgmres: 'ident' (identity), 'diag' (diagonal), 'pgmres' (Jacobi-preconditioned GMRES) dim_fgmres = 50 ! size of fgmres Krylov subspace @@ -729,7 +731,8 @@ subroutine input_data call broadcast_scalar(alphab, master_task) call broadcast_scalar(threshold_hw, master_task) call broadcast_scalar(Ktens, master_task) - call broadcast_scalar(e_ratio, master_task) + call broadcast_scalar(e_yieldcurve, master_task) + call broadcast_scalar(e_plasticpot, master_task) call broadcast_scalar(advection, master_task) call broadcast_scalar(conserv_check, master_task) call broadcast_scalar(shortwave, master_task) @@ -1440,9 +1443,10 @@ subroutine input_data if (kdyn == 1 .or. kdyn == 3) then write(nu_diag,1030) ' yield_curve = ', trim(yield_curve) if (trim(yield_curve) == 'ellipse') & - write(nu_diag,1002) ' e_ratio = ', e_ratio, ' : aspect ratio of ellipse' + write(nu_diag,1002) ' e_yieldcurve = ', e_yieldcurve, ' : aspect ratio of yield curve' + write(nu_diag,1002) ' e_plasticpot = ', e_plasticpot, ' : aspect ratio of plastic potential' endif - + if (trim(coriolis) == 'latitude') then tmpstr2 = ' : latitude-dependent Coriolis parameter' elseif (trim(coriolis) == 'contant') then diff --git a/configuration/scripts/ice_in b/configuration/scripts/ice_in index 443ff1cbb..bb44663eb 100644 --- a/configuration/scripts/ice_in +++ b/configuration/scripts/ice_in @@ -141,7 +141,8 @@ Cstar = 20 Cf = 17. Ktens = 0. - e_ratio = 2. + e_yieldcurve = 2. + e_plasticpot = 2. seabed_stress = .false. seabed_stress_method = 'LKD' k1 = 7.5 diff --git a/configuration/scripts/options/set_nml.alt03 b/configuration/scripts/options/set_nml.alt03 index 255d77261..98e794735 100644 --- a/configuration/scripts/options/set_nml.alt03 +++ b/configuration/scripts/options/set_nml.alt03 @@ -19,7 +19,7 @@ sw_dtemp = 0.02d0 tfrz_option = 'linear_salt' revised_evp = .false. Ktens = 0. -e_ratio = 2. +e_yieldcurve = 2. seabed_stress = .true. use_bathymetry = .true. l_mpond_fresh = .true. diff --git a/doc/source/cice_index.rst b/doc/source/cice_index.rst index 12bc8d32e..85acbece3 100644 --- a/doc/source/cice_index.rst +++ b/doc/source/cice_index.rst @@ -205,7 +205,9 @@ either Celsius or Kelvin units). "etax2", "2 x eta (shear viscous coefficient)", "kg/s" "evap", "evaporative water flux", "kg/m\ :math:`^2`/s" "ew_boundary_type", "type of east-west boundary condition", "" - "eyc", "coefficient for calculating the parameter E, 0\ :math:`<` eyc :math:`<`\ 1", "0.36" + "eyc", "coefficient for calculating the parameter E, 0\ :math:`<` eyc :math:`<`\ 1", "0.36" + "e_yieldcurve", "yield curve minor/major axis ratio", "2" + "e_plasticpot", "plastic potential minor/major axis ratio", "2" "**F**", "", "" "faero_atm", "aerosol deposition rate", "kg/m\ :math:`^2`/s" "faero_ocn", "aerosol flux to the ocean", "kg/m\ :math:`^2`/s" diff --git a/doc/source/master_list.bib b/doc/source/master_list.bib index 7c2a45a35..e4afa0c46 100644 --- a/doc/source/master_list.bib +++ b/doc/source/master_list.bib @@ -60,7 +60,7 @@ @string{CRST @string{IJHPCA={Int. J High Perform. Comput. Appl}} @string{PTRSA={Philos. Trans. Royal Soc. A}} @string{SIAMJCP={SIAM J. Sci. Comput.}} - +@string{TC={The Cryosphere}} % ********************************************** @@ -979,6 +979,17 @@ @article{Roach18 volume = {123}, year = {2018} } + +@Article{Ringeisen21 + author = "D. Ringeisen and L.B. Tremblay and M. Losch", + title = "{Non-normal flow rules affect fracture angles in sea ice viscous-plastic rheologies}", + journal = TC, + year = {2021}, + volume = {15}, + pages = {2873-2888}, + url = {https://doi.org/10.5194/tc-15-2873-2021} +} + @Article{Tsujino18, author = "H. Tsujino and S. Urakawa and R.J. Small and W.M. Kim and S.G. Yeager and et al.", title = "{JRA‐55 based surface dataset for driving ocean–sea‐ice models (JRA55‐do)}", diff --git a/doc/source/science_guide/sg_dynamics.rst b/doc/source/science_guide/sg_dynamics.rst index 6d7d32976..9c529b8ec 100644 --- a/doc/source/science_guide/sg_dynamics.rst +++ b/doc/source/science_guide/sg_dynamics.rst @@ -350,9 +350,9 @@ Following again the LKD method, the seabed stress coefficients are finally expre .. _internal-stress: -*************** +******************************** Internal stress -*************** +******************************** For convenience we formulate the stress tensor :math:`\bf \sigma` in terms of :math:`\sigma_1=\sigma_{11}+\sigma_{22}`, @@ -378,15 +378,6 @@ CICE can output the internal ice pressure which is an important field to support The internal ice pressure (``sigP``) is the average of the normal stresses multiplied by :math:`-1` and is therefore simply equal to :math:`-\sigma_1/2`. -Following the approach of :cite:`Konig10` (see also :cite:`Lemieux16`), the -elliptical yield curve can be modified such that the ice has isotropic tensile strength. -The tensile strength :math:`T_p` is expressed as a fraction of the ice strength :math:`P`, that is :math:`T_p=k_t P` -where :math:`k_t` should be set to a value between 0 and 1 (this can -be changed at runtime with the namelist parameter ``Ktens``). The ice -strength :math:`P` is a function of the ice thickness distribution as -described in the `Icepack -Documentation `_. - .. _stress-vp: Viscous-Plastic @@ -395,28 +386,45 @@ Viscous-Plastic The VP constitutive law is given by .. math:: - \sigma_{ij} = 2 \eta \dot{\epsilon}_{ij} + (\zeta - \eta) D_D - P_R(1 - k_t)\frac{\delta_{ij}}{2} + \sigma_{ij} = 2 \eta \dot{\epsilon}_{ij} + (\zeta - \eta) D_D - P_R\frac{\delta_{ij}}{2} :label: vp-const -where :math:`\eta` and :math:`\zeta` are the bulk and shear viscosities. +where :math:`\eta` and :math:`\zeta` are the bulk and shear viscosities and +:math:`P_R` is a “replacement pressure” (see :cite:`Geiger98`, for example), +which serves to prevent residual ice motion due to spatial +variations of the ice strength :math:`P` when the strain rates are exactly zero. + An elliptical yield curve is used, with the viscosities given by .. math:: \zeta = {P(1+k_t)\over 2\Delta}, .. math:: - \eta = e^{-2} \zeta, + \eta = e_g^{-2} \zeta, where .. math:: - \Delta = \left[D_D^2 + {1\over e^2}\left(D_T^2 + D_S^2\right)\right]^{1/2} + \Delta = \left[D_D^2 + {e_f^2\over e_g^4}\left(D_T^2 + D_S^2\right)\right]^{1/2}. + +The ice strength :math:`P` is a function of the ice thickness distribution as +described in the `Icepack Documentation `_. + +Two modifications to the standard VP rheology of :cite:`Hibler79` are available. +First, following the approach of :cite:`Konig10` (see also :cite:`Lemieux16`), the +elliptical yield curve can be modified such that the ice has isotropic tensile strength. +The tensile strength is expressed as a fraction of :math:`P`, that is :math:`k_t P` +where :math:`k_t` should be set to a value between 0 and 1 (this can +be changed at runtime with the namelist parameter ``Ktens``). + +Second, while :math:`e_f` is the ratio of the major and minor axes of the elliptical yield curve, the parameter +:math:`e_g` characterizes the plastic potential, i.e. another ellipse that decouples the flow rule from the +yield curve (:cite:`Ringeisen21`). :math:`e_f` and :math:`e_g` are respectively called ``e_yieldcurve`` and ``e_plasticpot`` in the code and +can be set in the namelist. The plastic potential can lead to more realistic fracture angles between linear kinematic features. :cite:`Ringeisen21` suggest to set :math:`e_f` to a value larger than 1 and to have :math:`e_g < e_f`. -and :math:`P_R` is a “replacement pressure” (see :cite:`Geiger98`, for -example), which serves to prevent residual ice motion due to spatial -variations of :math:`P` when the rates of strain are exactly zero. +By default, the namelist parameters are set to :math:`e_f=e_g=2` and :math:`k_t=0` which correspond to the standard VP rheology. -The parameter :math:`e` is the ratio of the major and minor axes of the elliptical yield curve, also called the ellipse aspect ratio. It can be changed using the namelist parameter ``e_ratio``. +There are four options in the code for solving the sea ice momentum equation with a VP formulation: the standard EVP approach, a 1d EVP solver, the revised EVP approach and an implicit Picard solver. The modifications to the yield curve and to the flow rule described above are available for these four different solution methods. .. _stress-evp: @@ -428,7 +436,7 @@ regularized version of the VP constitutive law :eq:`vp-const`. The constitutive .. math:: {1\over E}{\partial\sigma_1\over\partial t} + {\sigma_1\over 2\zeta} - + {P_R(1-k_t)\over 2\zeta} = D_D, \\ + + {P_R\over 2\zeta} = D_D, \\ :label: sig1 .. math:: @@ -455,7 +463,7 @@ parameter less than one. Including the modification proposed by :cite:`Bouillon1 .. math:: \begin{aligned} {\partial\sigma_1\over\partial t} + {\sigma_1\over 2T} - + {P_R(1-k_t)\over 2T} &=& {\zeta \over T} D_D, \\ + + {P_R\over 2T} &=& {\zeta \over T} D_D, \\ {\partial\sigma_2\over\partial t} + {\sigma_2\over 2T} &=& {\eta \over T} D_T,\\ {\partial\sigma_{12}\over\partial t} + {\sigma_{12}\over 2T} &=& @@ -466,7 +474,7 @@ Once discretized in time, these last three equations are written as .. math:: \begin{aligned} {(\sigma_1^{k+1}-\sigma_1^{k})\over\Delta t_e} + {\sigma_1^{k+1}\over 2T} - + {P_R^k(1-k_t)\over 2T} &=& {\zeta^k\over T} D_D^k, \\ + + {P_R^k\over 2T} &=& {\zeta^k\over T} D_D^k, \\ {(\sigma_2^{k+1}-\sigma_2^{k})\over\Delta t_e} + {\sigma_2^{k+1}\over 2T} &=& {\eta^k \over T} D_T^k,\\ {(\sigma_{12}^{k+1}-\sigma_{12}^{k})\over\Delta t_e} + {\sigma_{12}^{k+1}\over 2T} &=& @@ -480,6 +488,78 @@ the viscosity terms in the subcycling. Choices of the parameters used to define :math:`T` and :math:`\Delta t_e` are discussed in Sections :ref:`revp` and :ref:`parameters`. +.. _evp1d: + +1d EVP solver +~~~~~~~~~~~~~ + +The standard EVP solver iterates hundreds of times, where each iteration includes a communication through MPI and a limited number of calculations. This limits how much the solver can be optimized as the speed is primarily determined by the communication. The 1d EVP solver avoids the communication by utilizing shared memory, which removes the requirement for calls to the MPI communicator. As a consequence of this the potential scalability of the code is improved. The performance is best on shared memory but the solver is also functional on MPI and hybrid MPI/OpenMP setups as it will run on the master processor alone. + +The scalability of geophysical models is in general terms limited by the memory usage. In order to optimize this the 1d EVP solver solves the same equations that are outlined in the section :ref:`stress-evp` but it transforms all matrices to vectors (1d matrices) as this compiles better with the computer hardware. The vectorization and the contiguous placement of arrays in the memory makes it easier for the compiler to optimize the code and pass pointers instead of copying the vectors. The 1d solver is not supported for tripole grids and the code will abort if this combination is attempted. + +.. _revp: + +Revised EVP approach +~~~~~~~~~~~~~~~~~~~~ + +The revised EVP approach is based on a pseudo-time iterative scheme :cite:`Lemieux12`, :cite:`Bouillon13`, :cite:`Kimmritz15`. By construction, the revised EVP approach should lead to the VP solution +(given the right numerical parameters and a sufficiently large number of iterations). To do so, the inertial term is formulated such that it matches the backward Euler approach of +implicit solvers and there is an additional term for the pseudo-time iteration. Hence, with the revised approach, the discretized momentum equations :eq:`umom` and :eq:`vmom` become + +.. math:: + {\beta^*(u^{k+1}-u^k)\over\Delta t_e} + {m(u^{k+1}-u^n)\over\Delta t} + {\left({\tt vrel} \cos\theta + C_b \right)} u^{k+1} + - {\left(mf+{\tt vrel}\sin\theta\right)} v^{k+1} + = {{\partial\sigma_{1j}^{k+1}\over\partial x_j}} + + {\tau_{ax} - mg{\partial H_\circ\over\partial x} } + + {\tt vrel} {\left(U_w\cos\theta-V_w\sin\theta\right)}, + :label: umomr + +.. math:: + {\beta^*(v^{k+1}-v^k)\over\Delta t_e} + {m(v^{k+1}-v^n)\over\Delta t} + {\left({\tt vrel} \cos\theta + C_b \right)}v^{k+1} + + {\left(mf+{\tt vrel}\sin\theta\right)} u^{k+1} + = {{\partial\sigma_{2j}^{k+1}\over\partial x_j}} + + {\tau_{ay} - mg{\partial H_\circ\over\partial y} } + + {\tt vrel}{\left(U_w\sin\theta+V_w\cos\theta\right)}, + :label: vmomr + +where :math:`\beta^*` is a numerical parameter and :math:`u^n, v^n` are the components of the previous time level solution. +With :math:`\beta=\beta^* \Delta t \left( m \Delta t_e \right)^{-1}` :cite:`Bouillon13`, these equations can be written as + +.. math:: + \underbrace{\left((\beta+1){m\over\Delta t}+{\tt vrel} \cos\theta\ + C_b \right)}_{\tt cca} u^{k+1} + - \underbrace{\left(mf+{\tt vrel}\sin\theta\right)}_{\tt ccb}v^{k+1} + = \underbrace{{\partial\sigma_{1j}^{k+1}\over\partial x_j}}_{\tt strintx} + + \underbrace{\tau_{ax} - mg{\partial H_\circ\over\partial x} }_{\tt forcex} + + {\tt vrel}\underbrace{\left(U_w\cos\theta-V_w\sin\theta\right)}_{\tt waterx} + {m\over\Delta t}(\beta u^k + u^n), + :label: umomr2 + +.. math:: + \underbrace{\left(mf+{\tt vrel}\sin\theta\right)}_{\tt ccb} u^{k+1} + + \underbrace{\left((\beta+1){m\over\Delta t}+{\tt vrel} \cos\theta + C_b \right)}_{\tt cca}v^{k+1} + = \underbrace{{\partial\sigma_{2j}^{k+1}\over\partial x_j}}_{\tt strinty} + + \underbrace{\tau_{ay} - mg{\partial H_\circ\over\partial y} }_{\tt forcey} + + {\tt vrel}\underbrace{\left(U_w\sin\theta+V_w\cos\theta\right)}_{\tt watery} + {m\over\Delta t}(\beta v^k + v^n), + :label: vmomr2 + +At this point, the solutions :math:`u^{k+1}` and :math:`v^{k+1}` are obtained in the same manner as for the standard EVP approach (see equations :eq:`cevpuhat` to :eq:`cevpb`). + +Introducing another numerical parameter :math:`\alpha=2T \Delta t_e ^{-1}` :cite:`Bouillon13`, the stress equations in :eq:`sigdisc` become + +.. math:: + \begin{aligned} + {\alpha (\sigma_1^{k+1}-\sigma_1^{k})} + {\sigma_1^{k}} + + {P_R^k} &=& 2 \zeta^k D_D^k, \\ + {\alpha (\sigma_2^{k+1}-\sigma_2^{k})} + {\sigma_2^{k}} &=& 2 \eta^k D_T^k,\\ + {\alpha (\sigma_{12}^{k+1}-\sigma_{12}^{k})} + {\sigma_{12}^{k}} &=& + \eta^k D_S^k,\end{aligned} + +where as opposed to the classic EVP, the second term in each equation is at iteration :math:`k` :cite:`Bouillon13`. Also, as opposed to the classic EVP, +:math:`\Delta t_e` times the number of subcycles (or iterations) does not need to be equal to the advective time step :math:`\Delta t`. +Finally, as with the classic EVP approach, the stresses are initialized using the previous time level values. +The revised EVP is activated by setting the namelist parameter ``revised_evp = true``. +In the code :math:`\alpha` is ``arlx`` and :math:`\beta` is ``brlx``. The values of ``arlx`` and ``brlx`` can be set in the namelist. +It is recommended to use large values of these parameters and to set :math:`\alpha=\beta` :cite:`Kimmritz15`. + .. _stress-eap: Elastic-Anisotropic-Plastic @@ -667,77 +747,3 @@ rheology we compute the area loss rate due to ridging as Both ridging rate and sea ice strength are computed in the outer loop of the dynamics. - -.. _revp: - -**************** -Revised approach -**************** - -The revised EVP approach is based on a pseudo-time iterative scheme :cite:`Lemieux12`, :cite:`Bouillon13`, :cite:`Kimmritz15`. By construction, the revised EVP approach should lead to the VP solution -(given the right numerical parameters and a sufficiently large number of iterations). To do so, the inertial term is formulated such that it matches the backward Euler approach of -implicit solvers and there is an additional term for the pseudo-time iteration. Hence, with the revised approach, the discretized momentum equations :eq:`umom` and :eq:`vmom` become - -.. math:: - {\beta^*(u^{k+1}-u^k)\over\Delta t_e} + {m(u^{k+1}-u^n)\over\Delta t} + {\left({\tt vrel} \cos\theta + C_b \right)} u^{k+1} - - {\left(mf+{\tt vrel}\sin\theta\right)} v^{k+1} - = {{\partial\sigma_{1j}^{k+1}\over\partial x_j}} - + {\tau_{ax} - mg{\partial H_\circ\over\partial x} } - + {\tt vrel} {\left(U_w\cos\theta-V_w\sin\theta\right)}, - :label: umomr - -.. math:: - {\beta^*(v^{k+1}-v^k)\over\Delta t_e} + {m(v^{k+1}-v^n)\over\Delta t} + {\left({\tt vrel} \cos\theta + C_b \right)}v^{k+1} - + {\left(mf+{\tt vrel}\sin\theta\right)} u^{k+1} - = {{\partial\sigma_{2j}^{k+1}\over\partial x_j}} - + {\tau_{ay} - mg{\partial H_\circ\over\partial y} } - + {\tt vrel}{\left(U_w\sin\theta+V_w\cos\theta\right)}, - :label: vmomr - -where :math:`\beta^*` is a numerical parameter and :math:`u^n, v^n` are the components of the previous time level solution. -With :math:`\beta=\beta^* \Delta t \left( m \Delta t_e \right)^{-1}` :cite:`Bouillon13`, these equations can be written as - -.. math:: - \underbrace{\left((\beta+1){m\over\Delta t}+{\tt vrel} \cos\theta\ + C_b \right)}_{\tt cca} u^{k+1} - - \underbrace{\left(mf+{\tt vrel}\sin\theta\right)}_{\tt ccb}v^{k+1} - = \underbrace{{\partial\sigma_{1j}^{k+1}\over\partial x_j}}_{\tt strintx} - + \underbrace{\tau_{ax} - mg{\partial H_\circ\over\partial x} }_{\tt forcex} - + {\tt vrel}\underbrace{\left(U_w\cos\theta-V_w\sin\theta\right)}_{\tt waterx} + {m\over\Delta t}(\beta u^k + u^n), - :label: umomr2 - -.. math:: - \underbrace{\left(mf+{\tt vrel}\sin\theta\right)}_{\tt ccb} u^{k+1} - + \underbrace{\left((\beta+1){m\over\Delta t}+{\tt vrel} \cos\theta + C_b \right)}_{\tt cca}v^{k+1} - = \underbrace{{\partial\sigma_{2j}^{k+1}\over\partial x_j}}_{\tt strinty} - + \underbrace{\tau_{ay} - mg{\partial H_\circ\over\partial y} }_{\tt forcey} - + {\tt vrel}\underbrace{\left(U_w\sin\theta+V_w\cos\theta\right)}_{\tt watery} + {m\over\Delta t}(\beta v^k + v^n), - :label: vmomr2 - -At this point, the solutions :math:`u^{k+1}` and :math:`v^{k+1}` are obtained in the same manner as for the standard EVP approach (see equations :eq:`cevpuhat` to :eq:`cevpb`). - -Introducing another numerical parameter :math:`\alpha=2T \Delta t_e ^{-1}` :cite:`Bouillon13`, the stress equations in :eq:`sigdisc` become - -.. math:: - \begin{aligned} - {\alpha (\sigma_1^{k+1}-\sigma_1^{k})} + {\sigma_1^{k}} - + {P_R^k(1-k_t)} &=& 2 \zeta^k D_D^k, \\ - {\alpha (\sigma_2^{k+1}-\sigma_2^{k})} + {\sigma_2^{k}} &=& 2 \eta^k D_T^k,\\ - {\alpha (\sigma_{12}^{k+1}-\sigma_{12}^{k})} + {\sigma_{12}^{k}} &=& - \eta^k D_S^k,\end{aligned} - -where as opposed to the classic EVP, the second term in each equation is at iteration :math:`k` :cite:`Bouillon13`. Also, as opposed to the classic EVP, -:math:`\Delta t_e` times the number of subcycles (or iterations) does not need to be equal to the advective time step :math:`\Delta t`. -Finally, as with the classic EVP approach, the stresses are initialized using the previous time level values. -The revised EVP is activated by setting the namelist parameter ``revised_evp = true``. -In the code :math:`\alpha` is ``arlx`` and :math:`\beta` is ``brlx``. The values of ``arlx`` and ``brlx`` can be set in the namelist. -It is recommended to use large values of these parameters and to set :math:`\alpha=\beta` :cite:`Kimmritz15`. - -.. _evp1d: - -**************** -1d EVP solver -**************** - -The standard EVP solver iterates hundreds of times, where each iteration includes a communication through MPI and a limited number of calculations. This limits how much the solver can be optimized as the speed is primarily determined by the communication. The 1d EVP solver avoids the communication by utilizing shared memory, which removes the requirement for calls to the MPI communicator. As a consequence of this the potential scalability of the code is improved. The performance is best on shared memory but the solver is also functional on MPI and hybrid MPI/OpenMP setups as it will run on the master processor alone. - -The scalability of geophysical models is in general terms limited by the memory usage. In order to optimize this the 1d EVP solver solves the same equations that are outlined in the section :ref:`stress-evp` but it transforms all matrices to vectors (1d matrices) as this compiles better with the computer hardware. The vectorization and the contiguous placement of arrays in the memory makes it easier for the compiler to optimize the code and pass pointers instead of copying the vectors. The 1d solver is not supported for tripole grids and the code will abort if this combination is attempted.