Skip to content

Fixes to thickness regularization#45

Merged
jouvetg merged 3 commits into
instructed-glacier-model:mainfrom
gillian-smith:fix/regu_thk
May 8, 2026
Merged

Fixes to thickness regularization#45
jouvetg merged 3 commits into
instructed-glacier-model:mainfrom
gillian-smith:fix/regu_thk

Conversation

@gillian-smith
Copy link
Copy Markdown
Contributor

Here I propose a few fixes/changes/additions to the thickness regularization in the data_assimilation module.

regu_thk_v1()

  • when setting processes.data_assimilation.regularization.smooth_anisotropy_factor less than 1.0, ave4() needs to be applied to make dims of state.flowdirx, state.flowdiry match dims of dbdx, dbdy

regu_thk_v2()

  • The 2nd derivative term has a term missing, which is required to make it rotationally invariant/independent of choice of $x$ and $y$ axes: see https://hackmd.io/@gillian-smith/r1luUbVh-g for derivation (many thanks to James Maddison for spotting this)

  • The 1st derivative term is isotropic even if processes.data_assimilation.regularization.smooth_anisotropy_factor < 1.0 - I have changed this to match behaviour of regu_thk_v1()

  • For the 1st derivatives, the finite difference kernels kx and ky decouple odd- and even-indexed cells. This induces chessboard modes when the 2nd derivative term is inactive (processes.data_assimilation.regularization.thk_2nd_der=0). I propose reverting to the discretization used for regu_thk_v1(), or at least give users the option, with new parameter processes.data_assimilation.regularization.thk_1st_der_disc_version which can be "new" or "old".

  • Currently the kernel kxy also decouples odd- and even-numbered cells, but since bxy is always combined with other terms, this appears to not produce any visible numerical artefacts. I have left this as-is, but if we would prefer to change it, perhaps kxy = tf.constant([[ 1., -1., 0.],[ -1., 2., -1.],[0., -1., 1.]], tf.float32) / (2.0*dx*dx) is an option? (see https://en.wikipedia.org/wiki/Finite_difference#Multivariate_finite_differences)

  • suggest naming the coefficient variables alpha_2 and alpha_1 instead of alpha and beta, just because $\beta$ represents the anisotropy factor in https://doi.org/10.1017/jog.2022.41 and https://doi.org/10.5194/egusphere-2026-788, and it would be less confusing to keep this notation

After these fixes,
processes.data_assimilation.regularization.thk_version=2 processes.data_assimilation.regularization.abl_acc_balance=1 processes.data_assimilation.regularization.thk_2nd_der=0 processes.data_assimilation.regularization.thk_1st_der_disc_version=old

has effectively the same behaviour as
processes.data_assimilation.regularization.thk_version=1,

and the parameters processes.data_assimilation.regularization.thk_1st_der and processes.data_assimilation.regularization.thk are equivalent (up to a constant factor).

I think the regu_thk_v1() and regu_thk_v2() can be unified into one function, but perhaps later in another pull request.

@jouvetg jouvetg merged commit 7e10ca7 into instructed-glacier-model:main May 8, 2026
3 checks passed
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants