[MITgcm-devel] seaice dynamics in ocean 2-D
Jean-Michel Campin
jmc at mit.edu
Wed Jun 29 14:10:59 EDT 2016
Hi Martin,
Thanks for your suggestions.
I've just tried to comment out:
SEAICE_OLx=1
SEAICE_OLy=1
and the big differences between c65j and new code are gone.
I will check more carefully what the modifs from around checkpoint65r
are doing (specially in a 2-D set-up, i.e. 1-D seaice dyn).
Cheers,
Jean-Michel
On Wed, Jun 29, 2016 at 09:54:50AM +0200, Martin Losch wrote:
> Hi Jean-Michel,
>
> that???s annoying. I have changed the dynamics solver(s) of pkg/seaice a lot since then, see attached file changes_to_seaice.txt, where I compiled all (documented) changes to pkg/seaice since ckpt65j. Below are the ones that I think may have affected your experiement, with comments; the most likely candidate is this: checkpoint65r (2015/12/21), seaice_lsr.F r1.119. I hope this helps.
>
> Martin
>
> *** this will definitely change your results at the truncation level ***
> checkpoint65m (2015/06/15)
> o pkg/seaice:
> - seaice_calc_viscosity.F (r1.30): change computation of deltaC**2 to ensure
> positiveness, modify a few comments and improve variable names,
> - the changed computation of deltaC affects most seaice-related
> experiments at the truncation level.
> seaice_itd.thermo is affected by far the most, with
> only 7 digits of agreement remaining; update experiments.
>
> *** this is unlikely a problem, but you never know, especially with the new parameter that may lead to a different default ***
> checkpoint65n (2015/07/29)
> o pkg/seaice:
> - Rename SEAICE_VECTORIZE_LSR_ZEBRA to SEAICE_LSR_ZEBRA and define it
> in lab_sea forward experiment (changes results -> updated)
> o pkg/seaice:
> add new parameter SEAICE_deltaMin that is used ***only***
> for regularizing Delta (and nothing else, like the parameter
> which is also used for all sorts of things SEAICE_EPS).
> Defaults to SEAICE_EPS for backward compatibility only.
>
> *** this should not affect your results, but it???s quite a big change to many subroutines (seaice_calc_lhs.F, seaice_dynsolver.F, seaice_evp.F, seaice_lsr.F, seaice_preconditioner.F) ****
> checkpoint65o (2015/09/14)
> o pkg/seaice: introduce new parameter SEAICEscaleSurfStress (default = .FALSE.)
> - if TRUE scale ice-ocean and ice-atmosphere stress acting on sea ice
> by AREA according to Connelly et al. (2004), JPO.
> - For EVP and the LSR solver, the implicit part of the drag term is
> the only non-zero term in the denominator when there is no ice (which
> prevented running the model with zero ice-ocean drag). If
> SEAICEscaleSurfStress = .TRUE., this results in a division by zero
> (or zero main diagonals BU/BV) which need to be caught. A practical
> consequence is that for open water, the momentum equation reduce
> to trivially 0 = 0 + 0 (for EVP). For LSR, BU/BV are reset to 1, if
> they happen to be zero, often leading to a non-zero solution over
> open water.
>
> *** since you are using SEAICE_OLx/y=1, this is actually a very likely candidate, before that all LSR solutions with SEAICE_OLx/y>0 where clearly wrong (stripes in solutions along tile edges), JFNK solutions also improved, but not significantly, because LSR is only used as a preconditioner and the Newton step fixes this ???automatically"***
> checkpoint65r (2015/12/21)
> o pkg/seaice: (seaice_lsr.F r1.119)
> fix bug in tridiagonal solvers for SEAICE_OLx/y>0, this affects the JFNK
> solver (offline.dyn_jfnk) and global_ocean.cs32x15.seaice -> updated
>
> *** another candidate, but unlikely to have caused the large differences ***
> checkpoint65t (2016/02/10)
> o pkg/seaice/seaice_lsr.F: replace 1./SEAICE_deltaTdyn with recip_deltaT,
> affects some experiments at truncation level:
> - global_ocean.cs32x15.icedyn (11 digits),
> - global_ocean.cs32x15.seaice (11 digits),
> - seaice_itd (12 digits), seaice_itd.lipscomb07 (12 digits) -> update
> - also update global_ocean.cs32x15 adjoint and tangent-linear results
>
>
> checkpoint65k (2015/04/01)
> o pkg/seaice: add CPP brackets (cost function related).
>
> checkpoint65l (2015/05/04)
> o pkg/seaice/seaice_growth.F (in case of SEAICE_ITD)
> - replace tmpscal1**1.36 by faster exp(1.36*log(tmpscal1))
>
> checkpoint65m (2015/06/15)
> o pkg/seaice:
> - seaice_calc_viscosity.F: change computation of deltaC**2 to ensure
> positiveness, modify a few comments and improve variable names,
> - the changed computation of deltaC affects most seaice-related
> experiments at the truncation level.
> seaice_itd.thermo is affected by far the most, with
> only 7 digits of agreement remaining; update experiments.
>
> checkpoint65n (2015/07/29)
> o pkg/seaice:
> - initialise deltaC in seaice_lsr, cosmetic changes in seaice_lsr
> - add cpp-flag SEAICE_DELTA_SMOOTHREG for option of regularising
> delta with a smooth function in s/r seaice_calc_viscosity
> (no effect on EVP)
> - rename local suffix "sqr" into "sq" for more consistent variable names
> o pkg/seaice:
> add new parameter SEAICE_deltaMin that is used ***only***
> for regularizing Delta (and nothing else, like the parameter
> which is also used for all sorts of things SEAICE_EPS).
> Defaults to SEAICE_EPS for backward compatibility only.
> o pkg/seaice: replace erroneously storing e12Csqr on the wrong tape with
> the wrong key by inialising it before use -> fixes a recomputation
> problem but does not fix the gradient
> o pkg/seaice:
> - Rename SEAICE_VECTORIZE_LSR_ZEBRA to SEAICE_LSR_ZEBRA and define it
> in lab_sea forward experiment (changes results -> updated)
> - update adjoint experiment lab_sea after cleaning up seaice_lsr_tridiagu/v
> o pkg/seaice: modify seaice_lsr.F in order to improve the gradient
> computations; for SEAICE_VECTORIZE_LSR
> - move the loops over which the tridiagonal solvers (seaice_lsr_tridiagu/v)
> sweep to the outside,
> - remove store directives and add "CADJ loop sequential" directives
> in analogy to model/src/solve_tridiagonal.F
> - replace many "#ifdef SEAICE_VECTORIZE_LSR_ZEBRA" by variable loop
> boundaries to yield more readable code. This has also the charming side
> effect that your can use SEAICE_VECTORIZE_LSR_ZEBRA without
> SEAICE_VECTORIZE_LSR (i.e. adjoint without recomputations in
> seaice_lsr_tridiagu/v now requires either of these flags, vectorization
> still requires SEAICE_VECTORIZE_LSR with SEAICE_VECTORIZE_LSR_ZEBRA as
> an additional option)
> The tridiagonal solvers are now completely analoguous to solve_tridiagonal.F
> but the adjoint gradients (which are affected by this change) still explode.
>
> checkpoint65o (2015/09/14)
> o pkg/seaice: introduce fast ice parameterisation following Itkin et al. (2015)
> and Koenig-Beatty+Holland (2012)
> - two new parameters SEAICE_tensilFac and SEAICE_tensilDepth
> - global 2D field tensileStrength, computed in seaice_calc_ice_strength
> - requires an extra input argument for S/R seaice_calc_viscosities
> - tensileStrength's depth dependence is different from Itkin et al. (2015)
> (to be changed and played with, for now exponential)
> o pkg/seaice: introduce new parameter SEAICEscaleSurfStress (default = .FALSE.)
> - if TRUE scale ice-ocean and ice-atmosphere stress acting on sea ice
> by AREA according to Connelly et al. (2004), JPO.
> - For EVP and the LSR solver, the implicit part of the drag term is
> the only non-zero term in the denominator when there is no ice (which
> prevented running the model with zero ice-ocean drag). If
> SEAICEscaleSurfStress = .TRUE., this results in a division by zero
> (or zero main diagonals BU/BV) which need to be caught. A practical
> consequence is that for open water, the momentum equation reduce
> to trivially 0 = 0 + 0 (for EVP). For LSR, BU/BV are reset to 1, if
> they happen to be zero, often leading to a non-zero solution over
> open water.
> o pkg/seaice: revise EVP code (seaice_evp.F) in preparation for a more
> efficient EVP method.
> - introduce area-weighted averages for e12 (analogously to
> seaice_calc_viscosities.F); this totally changes lab_sea.hb87 (2 digits
> remain, but also changing the order of terms in the previous average
> changed the results as much) -> update experiment
> - add SEAICE_DELTA_SMOOTHREG code
> - disentangle computation of zetaC/deltaC and zetaZ/deltaZ, also for
> TEM-code
> - adjust loop bounds so that only required fields are computed
> - add new averaging code for zetaZ/deltaZ (again similar to
> seaice_calc_viscosities.F), inactive by default for now (will change
> results again), but will replace old code eventually
>
> checkpoint65p (2015/10/23)
> o pkg/seaice/seaice_evp.F: fix store directives and key computations, since
> there is no adjoint test for EVP, we don't know what it does to the adjoint
> solution, but at least there are no more recomputation warnings left
> o pkg/seaice/seaice_evp.F: slight modification of averaging so that EVP stands
> a chance of giving the same results as a fully converged VP (JFNK) solver
> o pkg/seaice:
> - replace SEAICE_EVP_USE_ADAPTATION by run time parameters, requires
> 5 additionals 2D fields and 6 new store directives
> - change logic in seaice_readparms.F: EVP code can now be turned on in
> various ways, setting SEAICE_deltaTev is no longer required,
> alternatively set SEAICE_evpAlpha/Beta, SEAICEuseEVPstar, SEAICEuseEVPrev,
> or SEAICEaEVPcoeff
> o pkg/seaice/seaice_evp.F: introduce "adaptive EVP" by Madlen Kimmritz
> - this is an intermediate state of a truely converging EVP that has
> the potential of outperforming implicit VP methods, I do not recommend
> using it yet.
> - for now, "adaptive EVP" is turned on by setting SEAICE_EVP_USE_ADAPTATION
> and all parameters are hardwired
> - for convergence it still requires very smooth regularisations, that means,
> if SEAICE_DELTA_SMOOTHREG is defined, deltaCreg=deltaC+deltaMin, which is
> different from the VP case with deltaCreg=sqrt(deltaC**2+deltaMin**2); also
> specific averaging appears to be required, still subject to tuning
> - still needs cleaning and I need to turn parameter into runtime parameters
> - add code to compute and print residuals of iteration, compile with
> defining ALLOW_SEAICE_EVP_RESIDUAL to enable it
> o pkg/seaice/seaice_evp.F: fix a bug (factor of four in front of e12Csq)
> that entered with r1.54 on Sep04, 2015, -> update results
>
> checkpoint65q (2015/11/18)
> o pkg/seaice/seaice_evp:
> - fix residual computations so that they no longer affect the actual
> simulation code
> - fix seaice_readparams to be able to actually set new parameters
> SEAICEaEVPcStar and SEAICEaEVPalphaMin, also fix logic a little
>
> checkpoint65r (2015/12/21)
> o pkg/seaice:
> fix bug in tridiagonal solvers for SEAICE_OLx/y>0, this affects the JFNK
> solver (offline.dyn_jfnk) and global_ocean.cs32x15.seaice -> updated
> o pkg/seaice:
> - fix a bug in computing areaS in (J-J -> J-1) that affects runs
> when SEAICEscaleSurfStress=.TRUE.,
> - apply areaW and areaS in S/R seaice_calc_lhs
> o pkg/seaice/seaice_growth:
> - catch potential division by zero in ITD code, does not change
> verification results
> o pkg/seaice/seaice_evp:
> - refine residual computation
>
> checkpoint65s (2016/01/13)
> o pkg/seaice: add new parameter SEAICEpressReplFac to choose between
> pressure replacement method (=1., default) or original Hibler (1979) (=0.)
> value can be from domain [0,1] to average between the two methods
>
> checkpoint65t (2016/02/10)
> o pkg/seaice:
> - change the termination criterion so that maxits (SEAICElinearIterMax)
> has an effect
> - unify iteration parameters for implicit solvers (JFNK and Picard)
> SEAICEnonLinIterMax replaces SEAICEnewtonIterMax/NPSEUDOTIMESTEPS
> SEAICElinearIterMax replaces SEAICEkrylovIterMax/SOLV_MAX_ITER
> SEAICEpreLinIterMax replaces SOLV_MAX_ITER in preconditioner
> SEAICEpreNL_IterMax replaces NPSEUDOTIMESTEPS in preconditioner
> SEAICEnonLinTol replaces JFNKgamma_nonlin
> - remove S/R SEAICE_FGMRES_DRIVER and put content into S/R SEAICE_JFKN
> - add new Picard-Krylov solver, compile with SEAICE_ALLOW_KRYLOV and
> use with SEAICEuseKrylov
> o pkg/seaice/seaice_lsr.F
> fix small bug: replace zetaZ by zetaZloc in S/R seaice_lsr_rhsu/v; only
> relevant if SEAICEuseStrImpCpl=T, and even then it does not seem to have
> an effect.
> o pkg/seaice: add two new parameters
> - SEAICEuseLSR to simplify the logic in seaice_dynsolver (this
> parameter is not in any namelist)
> - SEAICEusrPicardAsPrecon to allow s/r seaice_lsr to be used as
> a preconditioner for the non-linear Newton iteration of the JFNK
> solver
> o pkg/seaice/seaice_lsr.F: replace 1./SEAICE_deltaTdyn with recip_deltaT,
> affects some experiments at truncation level:
> - global_ocean.cs32x15.icedyn (11 digits),
> - global_ocean.cs32x15.seaice (11 digits),
> - seaice_itd (12 digits), seaice_itd.lipscomb07 (12 digits) -> update
> - also update global_ocean.cs32x15 adjoint and tangent-linear results
>
> checkpoint65v (2016/04/08)
> o pkg/seaice/seaice_evp.F
> - simplify computation of deltaZ and zetaZ to simple averaging following
> Bouillon et al. (2013). This removes a lot of code, but also changes
> verification experiment lab_sea.hb87 -> updated
>
> checkpoint65w (2016/05/12)
> o pkg/seaice: ITD code
> - fix picking up from a single category model by initialising
> heff/area/hsnowITD = 0. in seaice_itd_pickup.F
> - split a long warning message in seaice_check_pickup.F into two lines,
> because NEC SX-ACE cannot deal with the long line.
> o pkg/seaice:
> - add missing checks for SEAICE_ALLOW_KRYLOV and SEAICE_ALLOW_BOTTOMDRAG
> o pkg/seaice:
> - add basal stress as parameterisation for grounding ice
> following Lemieux et al. (2015)
> - turn on by setting new parameter SEAICEbasalDragK2 to a value
> larger than 0. Lemieux et al. (2015) recommend 15.
> - The adjoint complains with extra recomputations so there is also a
> new cpp-flag SEAICE_ALLOW_BOTTOMDRAG that is undefined by default in
> order to postpone the problem
> - compile the new code in lab_sea
>
> checkpoint65x (2016/06/11)
> o pkg/seaice: fix tensile stength computation
> - replace global field tensileStrength by tensileStrFac
> - tensileStrFac can be computed once in S/R seaice_init_varia
> - fortunately that does not affect any verification results
>
> o pkg/dic & seaice:
> - fix compilation of ocean component for coupled atm+ocn model with
> seaice dynamics & dic.
>
> > On 29 Jun 2016, at 01:34, Jean-Michel Campin <jmc at mit.edu> wrote:
> >
> > Hi Martin,
> >
> > I am trying to reproduce results from ~1.5 year ago (close to checkpoint65j),
> > with a 2-D (y-z) southern ocean channel, and am am getting very different seaice
> > velocity.
> > I am using pkg/thsice thermodyn with pkg/seaice dynamics, deltaT=1.h
> > and my data.seaice is attached.
> >
> > I can reproduce the old results when I use:
> > - full code from checkpoint65j (Feb 25, 2015)
> > - or just pkg/seaice from checkpoint65j into latest code
> >
> > I made a short (10.iter) comparison between c65j (d12a) and latest code (d12c),
> > starting from a pickup of the old run. The largest difference are on vIce as
> > seem on this plot (x-axis = latitude, magenta and red are d12a results
> > at iter 1 and 9 respectively; green, yellow, ceilan and blue are d12c results
> > at iter 1 3 5 and 9).
> >
> > Do you have any suggestion of what could cause these large difference ?
> >
> > Cheers,
> > Jean-Michel
> > <data.seaice><fig_uIce_d12.jpg>_______________________________________________
> > MITgcm-devel mailing list
> > MITgcm-devel at mitgcm.org
> > http://mitgcm.org/mailman/listinfo/mitgcm-devel
>
> _______________________________________________
> MITgcm-devel mailing list
> MITgcm-devel at mitgcm.org
> http://mitgcm.org/mailman/listinfo/mitgcm-devel
More information about the MITgcm-devel
mailing list