From 8c49091491fad174bf04e576d25964021b828b3c Mon Sep 17 00:00:00 2001 From: Gillian Smith Date: Thu, 9 Apr 2026 11:52:51 +0100 Subject: [PATCH 1/3] fix: apply ave4() to state.flowdirx and state.flowdiry to make dims match --- igm/processes/data_assimilation/cost_terms/regu_thk.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/igm/processes/data_assimilation/cost_terms/regu_thk.py b/igm/processes/data_assimilation/cost_terms/regu_thk.py index 2eb72583..58219269 100755 --- a/igm/processes/data_assimilation/cost_terms/regu_thk.py +++ b/igm/processes/data_assimilation/cost_terms/regu_thk.py @@ -5,6 +5,7 @@ import numpy as np import tensorflow as tf +from igm.processes.data_assimilation.utils import ave4 def regu_thk(cfg,state): if cfg.processes.data_assimilation.regularization.thk_version == 1: @@ -49,6 +50,9 @@ def regu_thk_v1(cfg,state): - gamma * tf.math.reduce_sum(state.thk) ) else: + state.flowdirx = ave4(state.flowdirx) + state.flowdiry = ave4(state.flowdiry) + dbdx = (field[:, 1:] - field[:, :-1])/state.dx dbdx = (dbdx[1:, :] + dbdx[:-1, :]) / 2.0 dbdy = (field[1:, :] - field[:-1, :])/state.dx From ff93aba20fc76b56f4cdbad78e21ed71e450a15d Mon Sep 17 00:00:00 2001 From: Gillian Smith Date: Thu, 9 Apr 2026 12:54:10 +0100 Subject: [PATCH 2/3] additional term in 2nd derivative regu; choice of new or old discretization for bx and by --- igm/conf/processes/data_assimilation.yaml | 1 + .../data_assimilation/cost_terms/regu_thk.py | 90 +++++++++++++------ 2 files changed, 66 insertions(+), 25 deletions(-) diff --git a/igm/conf/processes/data_assimilation.yaml b/igm/conf/processes/data_assimilation.yaml index a2adc773..4c36a2eb 100644 --- a/igm/conf/processes/data_assimilation.yaml +++ b/igm/conf/processes/data_assimilation.yaml @@ -32,6 +32,7 @@ data_assimilation: thk_version: 1 thk_2nd_der: 1.0e+6 thk_1st_der: 1.0e+2 + thk_1st_der_disc_version: new abl_acc_balance: 1.0 slidingco: 1.0 arrhenius: 10.0 diff --git a/igm/processes/data_assimilation/cost_terms/regu_thk.py b/igm/processes/data_assimilation/cost_terms/regu_thk.py index 58219269..39ec82e2 100755 --- a/igm/processes/data_assimilation/cost_terms/regu_thk.py +++ b/igm/processes/data_assimilation/cost_terms/regu_thk.py @@ -102,42 +102,82 @@ def regu_thk_v2(cfg,state): w_acc = 0.5 * (1.0 + tf.math.tanh((state.usurf - ELA) / 100.0)) rect = (r_acc * w_acc + r_abl * (1.0 - w_acc)) - # Compute derivatives directly on 2D tensors - kx, ky, kxx, kyy, kxy = _kernels(state.dx) # Derivative stencils - bx = _conv2(field, kx); by = _conv2(field, ky) # Derivatives of field - bxx = _conv2(field, kxx); byy = _conv2(field, kyy); bxy = _conv2(field, kxy) + anis_factor = cfg.processes.data_assimilation.regularization.smooth_anisotropy_factor + # Derivative stencils + kx, ky, kxx, kyy, kxy = _kernels(state.dx) + + # 1st derivatives (gradient) + if cfg.processes.data_assimilation.regularization.thk_1st_der_disc_version == "new": + bx = _conv2(field, kx); by = _conv2(field, ky) # Derivatives of field + if cfg.processes.data_assimilation.optimization.sole_mask: + bx = tf.where( state.icemaskobs > 0.0, bx, 0.0) + by = tf.where( state.icemaskobs > 0.0, by, 0.0) + else: # do as in regu_thk_v1() + if anis_factor == 1: + bx = (field[:, 1:] - field[:, :-1])/state.dx + by = (field[1:, :] - field[:-1, :])/state.dx + if cfg.processes.data_assimilation.optimization.sole_mask: + bx = tf.where( (state.icemaskobs[:, 1:] > 0.5) & (state.icemaskobs[:, :-1] > 0.5) , bx, 0.0) + by = tf.where( (state.icemaskobs[1:, :] > 0.5) & (state.icemaskobs[:-1, :] > 0.5) , by, 0.0) + else: + bx = (field[:, 1:] - field[:, :-1])/state.dx + bx = (bx[1:, :] + bx[:-1, :]) / 2.0 + by = (field[1:, :] - field[:-1, :])/state.dx + by = (by[:, 1:] + by[:, :-1]) / 2.0 + if cfg.processes.data_assimilation.optimization.sole_mask: + MASK = (state.icemaskobs[1:, 1:] > 0.5) & (state.icemaskobs[1:, :-1] > 0.5) & (state.icemaskobs[:-1, 1:] > 0.5) & (state.icemaskobs[:-1, :-1] > 0.5) + bx = tf.where( MASK, bx, 0.0) + by = tf.where( MASK, by, 0.0) + + # 2nd derivatives (Hessian) + bxx = _conv2(field, kxx); byy = _conv2(field, kyy); bxy = _conv2(field, kxy) if cfg.processes.data_assimilation.optimization.sole_mask: - bx = tf.where( state.icemaskobs > 0.0, bx, 0.0) - by = tf.where( state.icemaskobs > 0.0, by, 0.0) bxx = tf.where( state.icemaskobs > 0.0, bxx, 0.0) byy = tf.where( state.icemaskobs > 0.0, byy, 0.0) bxy = tf.where( state.icemaskobs > 0.0, bxy, 0.0) - if cfg.processes.data_assimilation.regularization.smooth_anisotropy_factor == 1: - Dnn = byy - Dtautau = bxx - else: - tx, ty = state.flowdirx, state.flowdiry # along-flow - nx, ny = -ty, tx - Dtautau = (tx*tx)*bxx + 2.0*(tx*ty)*bxy + (ty*ty)*byy # along-flow curvature - Dnn = (nx*nx)*bxx + 2.0*(nx*ny)*bxy + (ny*ny)*byy # cross-flow curvature - - alpha = cfg.processes.data_assimilation.regularization.thk_2nd_der - beta = cfg.processes.data_assimilation.regularization.thk_1st_der - anis_factor = cfg.processes.data_assimilation.regularization.smooth_anisotropy_factor + alpha_2 = cfg.processes.data_assimilation.regularization.thk_2nd_der + alpha_1 = cfg.processes.data_assimilation.regularization.thk_1st_der gamma = cfg.processes.data_assimilation.regularization.convexity_weight if cfg.processes.data_assimilation.optimization.fix_opti_normalization_issue: - J = alpha * rect * (tf.square(Dtautau) + anis_factor * tf.square(Dnn) ) \ - + beta * (tf.square(bx) + tf.square(by)) \ - - gamma * state.thk - return tf.reduce_mean( J) + if anis_factor == 1: + R_1 = 0.5 * tf.reduce_mean(tf.square(bx)) + tf.reduce_mean(tf.square(by)) + R_2 = 0.5 * rect * tf.reduce_mean(tf.square(bxx) + 2 * tf.square(bxy) + tf.square(byy)) + else: + ux, uy = state.flowdirx, state.flowdiry + if cfg.processes.data_assimilation.regularization.thk_1st_der_disc_version == "new": + ux1, uy1 = ux, uy + else: + ux1, uy1 = ave4(ux), ave4(uy) # so that dims match bx and by + + R_1 = 0.5 * tf.reduce_mean(tf.square(ux1*bx + uy1*by) + anis_factor * tf.square(uy1*bx - ux1*by)) + R_2 = 0.5 * tf.reduce_mean( + tf.square(ux*ux*bxx + 2*ux*uy*bxy + uy*uy*byy) \ + + 2 * anis_factor * tf.square(ux*uy*bxx + (uy*uy-ux*ux)*bxy + ux*uy*byy) \ + + anis_factor * anis_factor * tf.square(uy*uy*bxx - 2*ux*uy*bxy + ux*ux*byy) + ) + convexity_term = tf.reduce_mean(state.thk) else: - return alpha * (tf.nn.l2_loss( tf.square(Dtautau) + anis_factor * tf.square(Dnn) )) \ - + beta * tf.nn.l2_loss( tf.square(bx) + tf.square(by) ) \ - - gamma * tf.reduce_sum(state.thk) + if anis_factor == 1: + R_1 = tf.nn.l2_loss(bx) + tf.nn.l2_loss(by) + R_2 = tf.nn.l2_loss(bxx) + 2 * tf.nn.l2_loss(bxy) + tf.nn.l2_loss(byy) + else: + ux, uy = state.flowdirx, state.flowdiry + if cfg.processes.data_assimilation.regularization.thk_1st_der_disc_version == "new": + ux1, uy1 = ux, uy + else: + ux1, uy1 = ave4(ux), ave4(uy) # so that dims match bx and by + R_1 = tf.nn.l2_loss(ux1*bx + uy1*by) + anis_factor * tf.nn.l2_loss(uy1*bx - ux1*by) + R_2 = tf.nn.l2_loss(ux*ux*bxx + 2*ux*uy*bxy + uy*uy*byy) \ + + 2 * anis_factor * tf.nn.l2_loss(ux*uy*bxx + (uy*uy-ux*ux)*bxy + ux*uy*byy) \ + + anis_factor * anis_factor * tf.nn.l2_loss(uy*uy*bxx - 2*ux*uy*bxy + ux*ux*byy) + convexity_term = tf.reduce_sum(state.thk) + + R = alpha_1 * R_1 + alpha_2 * R_2 - gamma * convexity_term + return R def map_range(x,source_min, source_max, target_min, target_max): return target_min + (x - source_min) / (source_max - source_min) * (target_max - target_min) From 6fef41803b7fcbc6ae1c7f75b7ec11af7ddc4d71 Mon Sep 17 00:00:00 2001 From: Gillian Smith Date: Wed, 22 Apr 2026 13:04:19 +0100 Subject: [PATCH 3/3] fix: missing factor of 0.5 from by term --- igm/processes/data_assimilation/cost_terms/regu_thk.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/igm/processes/data_assimilation/cost_terms/regu_thk.py b/igm/processes/data_assimilation/cost_terms/regu_thk.py index 39ec82e2..62782787 100755 --- a/igm/processes/data_assimilation/cost_terms/regu_thk.py +++ b/igm/processes/data_assimilation/cost_terms/regu_thk.py @@ -143,7 +143,7 @@ def regu_thk_v2(cfg,state): if cfg.processes.data_assimilation.optimization.fix_opti_normalization_issue: if anis_factor == 1: - R_1 = 0.5 * tf.reduce_mean(tf.square(bx)) + tf.reduce_mean(tf.square(by)) + R_1 = 0.5 * (tf.reduce_mean(tf.square(bx)) + tf.reduce_mean(tf.square(by))) R_2 = 0.5 * rect * tf.reduce_mean(tf.square(bxx) + 2 * tf.square(bxy) + tf.square(byy)) else: ux, uy = state.flowdirx, state.flowdiry