Skip to content

Commit adc26e4

Browse files
MathieuMorlighemMathieuMorlighem
committed
CHG: completing first version of nudging code
1 parent d4bfc34 commit adc26e4

9 files changed

Lines changed: 188 additions & 125 deletions

File tree

‎src/c/cores/controlnudging_core.cpp‎

Lines changed: 57 additions & 35 deletions
Original file line numberDiff line numberDiff line change
@@ -12,41 +12,49 @@ void controlnudging_core(FemModel* femmodel){
1212

1313
/*Intermediaries*/
1414
IssmDouble dCdt1, dCdt2, dCdt3, dHdt;
15-
int numvertices=femmodel->vertices->NumberOfVertices();
15+
IssmDouble time;
16+
int maxiter;
1617

1718
/*Hard-coded nudging parameters*/
18-
IssmDouble timeadjust = 1/6; // yr per nudging step (2 months)
19-
IssmDouble H0 = 100; // m — thickness error scale (smaller = more sensitive), 100 m used in van den Akker et al (2025)
20-
IssmDouble r = 0.4; // relaxation strength toward C_inv (0 = none, 1 = strong) , 0.5 used in van den Akker et al (2025)
21-
int nsteps = 12000; // number of nudging steps
19+
IssmDouble minstep = -0.05; // max log10(C) decrease per step (was -0.05, too large), -0.03 or 0.03
20+
IssmDouble maxstep = 0.05; // max log10(C) increase per step
21+
IssmDouble C_log_min = -2; // log10(C) floor = C >= 1, 1 is the typical limit in a budd law
22+
IssmDouble C_log_max = 4.0; // log10(C) ceiling = C <= 10000, 1000 would be the typical limit but I want it to give it a bit more leeway
2223

2324
/*User-defined nudging parameters*/
25+
femmodel->parameters->FindParam(&maxiter, InversionMaxiterEnum);
2426
IssmDouble tau_C = femmodel->parameters->FindParam(InversionTauCEnum);
27+
IssmDouble H0 = femmodel->parameters->FindParam(InversionH0Enum);
28+
IssmDouble r = femmodel->parameters->FindParam(InversionRelaxationEnum);
2529
IssmDouble yts = femmodel->parameters->FindParam(ConstantsYtsEnum);
2630
IssmDouble tmax = femmodel->parameters->FindParam(TimesteppingFinalTimeEnum);
2731
IssmDouble tmin = femmodel->parameters->FindParam(TimesteppingStartTimeEnum);
28-
IssmDouble deltat = (tmax - tmin)*yts; // convert to years
32+
IssmDouble deltat = (tmax - tmin);
2933

3034
/*Fields before/after*/
31-
IssmDouble* C0 = NULL;
32-
IssmDouble* C = NULL;
33-
IssmDouble* Hobs = NULL;
34-
IssmDouble* H = NULL;
35-
IssmDouble* V = NULL;
36-
IssmDouble* O_ls = NULL;
35+
IssmDouble *C0 = NULL;
36+
IssmDouble *C = NULL;
37+
IssmDouble *H_obs = NULL;
38+
IssmDouble *H_old = NULL;
39+
IssmDouble *H = NULL;
40+
IssmDouble *V = NULL;
41+
IssmDouble *O_ls = NULL;
3742

3843
/*Get Fields once and for all*/
44+
int numvertices = femmodel->vertices->NumberOfVertices();
3945
GetVectorFromInputsx(&C0, femmodel, FrictionCoefficientEnum, VertexSIdEnum);
40-
GetVectorFromInputsx(&Hobs, femmodel, ThicknessEnum, VertexSIdEnum);
46+
GetVectorFromInputsx(&H_obs, femmodel, ThicknessEnum, VertexSIdEnum);
4147

4248
femmodel->parameters->SetParam(false,SaveResultsEnum);
43-
for(int m=0;m<nsteps;m++){
44-
_printf0_("\n=== NUDGING STEP "<< m <<"/"<< nsteps << " ===\n");
49+
for(int m=0;m<maxiter;m++){
50+
_printf0_("\n=== NUDGING STEP "<< m+1 <<"/"<< maxiter << " ===\n");
4551

46-
/*we need to make sure we do not modify femmodel at each iteration, make a copy*/
47-
FemModel* femmodel_temp = femmodel->copy();
52+
/*Get ice thickness before we run a transient step*/
53+
GetVectorFromInputsx(&H_old, femmodel, ThicknessEnum, VertexSIdEnum);
4854

4955
/*Solve another transient simulation*/
56+
femmodel->parameters->FindParam(&time,TimeEnum);
57+
femmodel->parameters->SetParam(time+deltat,TimesteppingFinalTimeEnum);
5058
transient_core(femmodel);
5159

5260
/*Extract results*/
@@ -56,55 +64,69 @@ void controlnudging_core(FemModel* femmodel){
5664
GetVectorFromInputsx(&O_ls, femmodel, MaskOceanLevelsetEnum, VertexSIdEnum);
5765

5866
/*Update friction coefficient accordingly*/
67+
IssmDouble RMSE_H = 0.;
68+
IssmDouble RMSE_dHdt = 0.;
5969
for(int i=0;i<numvertices;i++){
6070

71+
/*Compute thickness change for this vertex*/
72+
IssmDouble dH_now = H[i] - H_obs[i];
73+
IssmDouble dHdt_now = (H[i] - H_old[i])/(deltat/yts);
74+
RMSE_H += dH_now*dH_now;
75+
RMSE_dHdt += dHdt_now*dHdt_now;
76+
6177
/*1. : thickness error — push C to reduce H deviation
6278
* Sign: if H > H_obs (too thick), decrease C (less friction → faster ice → larger flux, lower H)*/
63-
dCdt1 = -1*((H[i] - Hobs[i])/ H0) / tau_C;
79+
dCdt1 = -1*(dH_now/H0) / tau_C;
6480

6581
/*2. : tendency — damp ongoing thinning/thickening
6682
* Sign: if dH/dt > 0 (thickening), if it is thickening and already too thick, decrease c faster.
6783
* If it is thinning and the thickness is already to large, the ice is already moving in the correct
6884
* direction, so make C decrease a little bit less */
69-
dHdt = (H[i] - Hobs[i])/deltat;
70-
dCdt2 = -1*(dHdt/H0) /tau_C;
85+
dCdt2 = -1*(dHdt_now/H0);
7186

7287
/*3. Relaxation term*/
73-
IssmDouble C_loginv = log10(max(C[i], 1.));
74-
IssmDouble C_loginv0 = log10(max(C0[i], 1.));
75-
dCdt3 = 1*(r / tau_C)*(C_loginv-C_loginv0);
88+
IssmDouble C_log = log10(max(C[i], 1.));
89+
IssmDouble C_log0 = log10(max(C0[i], 1.));
90+
dCdt3 = 1*(r / tau_C)*(C_log0 - C_log);
7691

7792
/* Do not nudge C where velocity ratio is very low AND ice is too thin
7893
* These cells need upstream replenishment, not local friction changes*/
79-
if(V[i]<15./yts && dHdt<-50/yts && O_ls[i]>0.) {
80-
dCdt2 = 0.;
81-
dCdt3 = 0.;
82-
}
94+
//if(V[i]<15./yts && dHdt<-50/yts && O_ls[i]>0.) {
95+
// dCdt2 = 0.;
96+
// dCdt3 = 0.;
97+
//}
8398

8499
/*Compute total dC in log space*/
85-
IssmDouble dC_log = deltat*(dCdt1 + dCdt2 + dCdt3);
100+
IssmDouble dC_log = deltat/yts*(dCdt1 + dCdt2 + dCdt3);
101+
102+
/*Clip dC_log to not change too much*/
103+
if(dC_log>maxstep) dC_log = maxstep;
104+
if(dC_log<minstep) dC_log = minstep;
86105

87106
/*Correction: prevent C from decreasing where ice is too thin and
88107
* still thinning — the model needs more friction here, not less*/
89-
if( ((H[i] - Hobs[i])< -20.) && (dHdt<-1.) ){
90-
dC_log = max(dC_log, 0.);
91-
}
108+
//if( ((H[i] - H_obs[i])< -20.) && (dHdt<-1.) ){
109+
// dC_log = max(dC_log, 0.);
110+
//}
92111

93112
/*Update friction coefficient now*/
94-
C[i] = pow(10., C_loginv + dC_log);
113+
C[i] = pow(10., C_log + dC_log);
95114
}
96115

97116
InputUpdateFromVectorx(femmodel, C, FrictionCoefficientEnum, VertexSIdEnum);
98-
_error_("not implemented yet");
99117

118+
/*Print statistics*/
119+
_printf0_(" → RMSE H : " << sqrt(RMSE_H/numvertices) << " m\n");
120+
_printf0_(" → RMSE dHdt: " << sqrt(RMSE_dHdt/numvertices) << " m/yr\n");
121+
122+
xDelete<IssmDouble>(H_old);
100123
xDelete<IssmDouble>(H);
101124
xDelete<IssmDouble>(C);
102125
xDelete<IssmDouble>(O_ls);
103126
xDelete<IssmDouble>(V);
104-
delete femmodel_temp;
105127
}
106128

107129
/*Clean up and return*/
108130
xDelete<IssmDouble>(C0);
109-
xDelete<IssmDouble>(Hobs);
131+
xDelete<IssmDouble>(H_obs);
110132
}

‎src/c/modules/ModelProcessorx/Control/CreateParametersControl.cpp‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -153,6 +153,9 @@ void CreateParametersControl(Parameters* parameters,IoModel* iomodel,int solutio
153153
parameters->AddObject(iomodel->CopyConstantObject("md.inversion.maxiter",InversionMaxiterEnum));
154154
break;
155155
case 5: /*Nudging*/
156+
parameters->AddObject(iomodel->CopyConstantObject("md.inversion.maxiter",InversionMaxiterEnum));
157+
parameters->AddObject(iomodel->CopyConstantObject("md.inversion.H0",InversionH0Enum));
158+
parameters->AddObject(iomodel->CopyConstantObject("md.inversion.relaxation",InversionRelaxationEnum));
156159
parameters->AddObject(iomodel->CopyConstantObject("md.inversion.tau_C",InversionTauCEnum));
157160
break;
158161
default:

‎src/c/shared/Enum/Enum.vim‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -345,12 +345,14 @@ syn keyword cConstant InversionGrtolEnum
345345
syn keyword cConstant InversionGttolEnum
346346
syn keyword cConstant InversionIncompleteAdjointEnum
347347
syn keyword cConstant InversionIscontrolEnum
348+
syn keyword cConstant InversionH0Enum
348349
syn keyword cConstant InversionMaxiterEnum
349350
syn keyword cConstant InversionMaxiterPerStepEnum
350351
syn keyword cConstant InversionMaxstepsEnum
351352
syn keyword cConstant InversionNstepsEnum
352353
syn keyword cConstant InversionNumControlParametersEnum
353354
syn keyword cConstant InversionNumCostFunctionsEnum
355+
syn keyword cConstant InversionRelaxationEnum
354356
syn keyword cConstant InversionStepThresholdEnum
355357
syn keyword cConstant InversionStopFlagEnum
356358
syn keyword cConstant InversionTauCEnum

‎src/c/shared/Enum/EnumDefinitions.h‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -339,12 +339,14 @@ enum definitions{
339339
InversionGttolEnum,
340340
InversionIncompleteAdjointEnum,
341341
InversionIscontrolEnum,
342+
InversionH0Enum,
342343
InversionMaxiterEnum,
343344
InversionMaxiterPerStepEnum,
344345
InversionMaxstepsEnum,
345346
InversionNstepsEnum,
346347
InversionNumControlParametersEnum,
347348
InversionNumCostFunctionsEnum,
349+
InversionRelaxationEnum,
348350
InversionStepThresholdEnum,
349351
InversionStopFlagEnum,
350352
InversionTauCEnum,

‎src/c/shared/Enum/EnumToStringx.cpp‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -347,12 +347,14 @@ const char* EnumToStringx(int en){
347347
case InversionGttolEnum : return "InversionGttol";
348348
case InversionIncompleteAdjointEnum : return "InversionIncompleteAdjoint";
349349
case InversionIscontrolEnum : return "InversionIscontrol";
350+
case InversionH0Enum : return "InversionH0";
350351
case InversionMaxiterEnum : return "InversionMaxiter";
351352
case InversionMaxiterPerStepEnum : return "InversionMaxiterPerStep";
352353
case InversionMaxstepsEnum : return "InversionMaxsteps";
353354
case InversionNstepsEnum : return "InversionNsteps";
354355
case InversionNumControlParametersEnum : return "InversionNumControlParameters";
355356
case InversionNumCostFunctionsEnum : return "InversionNumCostFunctions";
357+
case InversionRelaxationEnum : return "InversionRelaxation";
356358
case InversionStepThresholdEnum : return "InversionStepThreshold";
357359
case InversionStopFlagEnum : return "InversionStopFlag";
358360
case InversionTauCEnum : return "InversionTauC";

‎src/c/shared/Enum/Enumjl.vim‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -338,12 +338,14 @@ syn keyword juliaConstC InversionGrtolEnum
338338
syn keyword juliaConstC InversionGttolEnum
339339
syn keyword juliaConstC InversionIncompleteAdjointEnum
340340
syn keyword juliaConstC InversionIscontrolEnum
341+
syn keyword juliaConstC InversionH0Enum
341342
syn keyword juliaConstC InversionMaxiterEnum
342343
syn keyword juliaConstC InversionMaxiterPerStepEnum
343344
syn keyword juliaConstC InversionMaxstepsEnum
344345
syn keyword juliaConstC InversionNstepsEnum
345346
syn keyword juliaConstC InversionNumControlParametersEnum
346347
syn keyword juliaConstC InversionNumCostFunctionsEnum
348+
syn keyword juliaConstC InversionRelaxationEnum
347349
syn keyword juliaConstC InversionStepThresholdEnum
348350
syn keyword juliaConstC InversionStopFlagEnum
349351
syn keyword juliaConstC InversionTauCEnum

0 commit comments

Comments
 (0)