Skip to content

Commit 3fc9d34

Browse files
MathieuMorlighemMathieuMorlighem
committed
CHG: added more options for nudging
1 parent fc36666 commit 3fc9d34

11 files changed

Lines changed: 172 additions & 125 deletions

File tree

‎src/c/cores/controlnudging_core.cpp‎

Lines changed: 23 additions & 20 deletions
Original file line numberDiff line numberDiff line change
@@ -15,24 +15,21 @@ void controlnudging_core(FemModel* femmodel){
1515
IssmDouble time;
1616
int maxiter;
1717

18-
/*Hard-coded nudging parameters*/
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
23-
2418
/*User-defined nudging parameters*/
2519
femmodel->parameters->FindParam(&maxiter, InversionMaxiterEnum);
26-
IssmDouble tau_C = femmodel->parameters->FindParam(InversionTauCEnum);
27-
IssmDouble H0 = femmodel->parameters->FindParam(InversionH0Enum);
28-
IssmDouble r = femmodel->parameters->FindParam(InversionRelaxationEnum);
29-
IssmDouble yts = femmodel->parameters->FindParam(ConstantsYtsEnum);
30-
IssmDouble tmax = femmodel->parameters->FindParam(TimesteppingFinalTimeEnum);
31-
IssmDouble tmin = femmodel->parameters->FindParam(TimesteppingStartTimeEnum);
32-
IssmDouble deltat = (tmax - tmin);
20+
IssmDouble tau_C = femmodel->parameters->FindParam(InversionTauCEnum);
21+
IssmDouble max_inc = femmodel->parameters->FindParam(InversionMaxIncrementEnum);
22+
IssmDouble H0 = femmodel->parameters->FindParam(InversionH0Enum);
23+
IssmDouble r = femmodel->parameters->FindParam(InversionRelaxationEnum);
24+
IssmDouble yts = femmodel->parameters->FindParam(ConstantsYtsEnum);
25+
IssmDouble tmax = femmodel->parameters->FindParam(TimesteppingFinalTimeEnum);
26+
IssmDouble tmin = femmodel->parameters->FindParam(TimesteppingStartTimeEnum);
27+
IssmDouble deltat = (tmax - tmin);
3328

3429
/*Fields before/after*/
3530
IssmDouble *C0 = NULL;
31+
IssmDouble *Cmin = NULL;
32+
IssmDouble *Cmax = NULL;
3633
IssmDouble *C = NULL;
3734
IssmDouble *H_obs = NULL;
3835
IssmDouble *H_old = NULL;
@@ -43,6 +40,8 @@ void controlnudging_core(FemModel* femmodel){
4340
/*Get Fields once and for all*/
4441
int numvertices = femmodel->vertices->NumberOfVertices();
4542
GetVectorFromInputsx(&C0, femmodel, FrictionCoefficientEnum, VertexSIdEnum);
43+
GetVectorFromInputsx(&Cmin,femmodel, InversionMinParameterEnum, VertexSIdEnum);
44+
GetVectorFromInputsx(&Cmax,femmodel, InversionMaxParameterEnum, VertexSIdEnum);
4645
GetVectorFromInputsx(&H_obs, femmodel, ThicknessEnum, VertexSIdEnum);
4746

4847
femmodel->parameters->SetParam(false,SaveResultsEnum);
@@ -59,10 +58,10 @@ void controlnudging_core(FemModel* femmodel){
5958

6059
/*Extract results*/
6160
xDelete<IssmDouble>(C);
62-
GetVectorFromInputsx(&H, femmodel, ThicknessEnum, VertexSIdEnum);
63-
GetVectorFromInputsx(&C, femmodel, FrictionCoefficientEnum, VertexSIdEnum);
64-
GetVectorFromInputsx(&V, femmodel, VelEnum, VertexSIdEnum);
65-
GetVectorFromInputsx(&O_ls, femmodel, MaskOceanLevelsetEnum, VertexSIdEnum);
61+
GetVectorFromInputsx(&H, femmodel, ThicknessEnum, VertexSIdEnum);
62+
GetVectorFromInputsx(&C, femmodel, FrictionCoefficientEnum, VertexSIdEnum);
63+
GetVectorFromInputsx(&V, femmodel, VelEnum, VertexSIdEnum);
64+
GetVectorFromInputsx(&O_ls,femmodel, MaskOceanLevelsetEnum, VertexSIdEnum);
6665

6766
/*Update friction coefficient accordingly*/
6867
IssmDouble RMSE_H = 0.;
@@ -101,8 +100,8 @@ void controlnudging_core(FemModel* femmodel){
101100
IssmDouble dC_log = deltat*(dCdt1 + dCdt2 + dCdt3);
102101

103102
/*Clip dC_log to not change too much*/
104-
if(dC_log>maxstep) dC_log = maxstep;
105-
if(dC_log<minstep) dC_log = minstep;
103+
if(dC_log> max_inc) dC_log = max_inc;
104+
if(dC_log<-max_inc) dC_log = -max_inc;
106105

107106
/*Correction: prevent C from decreasing where ice is too thin and
108107
* still thinning — the model needs more friction here, not less*/
@@ -112,12 +111,14 @@ void controlnudging_core(FemModel* femmodel){
112111

113112
/*Update friction coefficient now*/
114113
C[i] = pow(10., C_log + dC_log);
114+
if(C[i] > Cmax[i]) C[i] = Cmax[i];
115+
if(C[i] < Cmin[i]) C[i] = Cmin[i];
115116
}
116117

117118
InputUpdateFromVectorx(femmodel, C, FrictionCoefficientEnum, VertexSIdEnum);
118119

119120
/*Print statistics*/
120-
_printf0_(" → RMSE H : " << sqrt(RMSE_H/numvertices) << " m\n");
121+
_printf0_(" → RMSE H : " << sqrt(RMSE_H/numvertices) << " m\n");
121122
_printf0_(" → RMSE dHdt: " << sqrt(RMSE_dHdt/numvertices) << " m/yr\n");
122123

123124
xDelete<IssmDouble>(H_old);
@@ -132,5 +133,7 @@ void controlnudging_core(FemModel* femmodel){
132133
/*Clean up and return*/
133134
xDelete<IssmDouble>(C);
134135
xDelete<IssmDouble>(C0);
136+
xDelete<IssmDouble>(Cmin);
137+
xDelete<IssmDouble>(Cmax);
135138
xDelete<IssmDouble>(H_obs);
136139
}

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

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -154,6 +154,7 @@ void CreateParametersControl(Parameters* parameters,IoModel* iomodel,int solutio
154154
break;
155155
case 5: /*Nudging*/
156156
parameters->AddObject(iomodel->CopyConstantObject("md.inversion.maxiter",InversionMaxiterEnum));
157+
parameters->AddObject(iomodel->CopyConstantObject("md.inversion.max_increment",InversionMaxIncrementEnum));
157158
parameters->AddObject(iomodel->CopyConstantObject("md.inversion.H0",InversionH0Enum));
158159
parameters->AddObject(iomodel->CopyConstantObject("md.inversion.relaxation",InversionRelaxationEnum));
159160
parameters->AddObject(iomodel->CopyConstantObject("md.inversion.tau_C",InversionTauCEnum));

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

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -27,6 +27,8 @@ void UpdateElementsAndMaterialsControl(Elements* elements,Parameters* parameters
2727
iomodel->FindConstant(&inversiontype,"md.inversion.type");
2828
if(inversiontype==5){
2929
iomodel->FetchDataToInput(inputs,elements,"md.inversion.dhdt_obs",BalancethicknessThickeningRateEnum);
30+
iomodel->FetchDataToInput(inputs,elements,"md.inversion.max_parameters",InversionMaxParameterEnum);
31+
iomodel->FetchDataToInput(inputs,elements,"md.inversion.min_parameters",InversionMinParameterEnum);
3032
return;
3133
}
3234

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

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -346,6 +346,7 @@ syn keyword cConstant InversionGttolEnum
346346
syn keyword cConstant InversionIncompleteAdjointEnum
347347
syn keyword cConstant InversionIscontrolEnum
348348
syn keyword cConstant InversionH0Enum
349+
syn keyword cConstant InversionMaxIncrementEnum
349350
syn keyword cConstant InversionMaxiterEnum
350351
syn keyword cConstant InversionMaxiterPerStepEnum
351352
syn keyword cConstant InversionMaxstepsEnum
@@ -1008,6 +1009,8 @@ syn keyword cConstant IceEnum
10081009
syn keyword cConstant IceMaskNodeActivationEnum
10091010
syn keyword cConstant InputEnum
10101011
syn keyword cConstant InversionCostFunctionsCoefficientsEnum
1012+
syn keyword cConstant InversionMaxParameterEnum
1013+
syn keyword cConstant InversionMinParameterEnum
10111014
syn keyword cConstant InversionSurfaceObsEnum
10121015
syn keyword cConstant InversionThicknessObsEnum
10131016
syn keyword cConstant InversionVelObsEnum

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

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -340,6 +340,7 @@ enum definitions{
340340
InversionIncompleteAdjointEnum,
341341
InversionIscontrolEnum,
342342
InversionH0Enum,
343+
InversionMaxIncrementEnum,
343344
InversionMaxiterEnum,
344345
InversionMaxiterPerStepEnum,
345346
InversionMaxstepsEnum,
@@ -1004,6 +1005,8 @@ enum definitions{
10041005
IceMaskNodeActivationEnum,
10051006
InputEnum,
10061007
InversionCostFunctionsCoefficientsEnum,
1008+
InversionMaxParameterEnum,
1009+
InversionMinParameterEnum,
10071010
InversionSurfaceObsEnum,
10081011
InversionThicknessObsEnum,
10091012
InversionVelObsEnum,

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

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -348,6 +348,7 @@ const char* EnumToStringx(int en){
348348
case InversionIncompleteAdjointEnum : return "InversionIncompleteAdjoint";
349349
case InversionIscontrolEnum : return "InversionIscontrol";
350350
case InversionH0Enum : return "InversionH0";
351+
case InversionMaxIncrementEnum : return "InversionMaxIncrement";
351352
case InversionMaxiterEnum : return "InversionMaxiter";
352353
case InversionMaxiterPerStepEnum : return "InversionMaxiterPerStep";
353354
case InversionMaxstepsEnum : return "InversionMaxsteps";
@@ -1010,6 +1011,8 @@ const char* EnumToStringx(int en){
10101011
case IceMaskNodeActivationEnum : return "IceMaskNodeActivation";
10111012
case InputEnum : return "Input";
10121013
case InversionCostFunctionsCoefficientsEnum : return "InversionCostFunctionsCoefficients";
1014+
case InversionMaxParameterEnum : return "InversionMaxParameter";
1015+
case InversionMinParameterEnum : return "InversionMinParameter";
10131016
case InversionSurfaceObsEnum : return "InversionSurfaceObs";
10141017
case InversionThicknessObsEnum : return "InversionThicknessObs";
10151018
case InversionVelObsEnum : return "InversionVelObs";

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

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -339,6 +339,7 @@ syn keyword juliaConstC InversionGttolEnum
339339
syn keyword juliaConstC InversionIncompleteAdjointEnum
340340
syn keyword juliaConstC InversionIscontrolEnum
341341
syn keyword juliaConstC InversionH0Enum
342+
syn keyword juliaConstC InversionMaxIncrementEnum
342343
syn keyword juliaConstC InversionMaxiterEnum
343344
syn keyword juliaConstC InversionMaxiterPerStepEnum
344345
syn keyword juliaConstC InversionMaxstepsEnum
@@ -1001,6 +1002,8 @@ syn keyword juliaConstC IceEnum
10011002
syn keyword juliaConstC IceMaskNodeActivationEnum
10021003
syn keyword juliaConstC InputEnum
10031004
syn keyword juliaConstC InversionCostFunctionsCoefficientsEnum
1005+
syn keyword juliaConstC InversionMaxParameterEnum
1006+
syn keyword juliaConstC InversionMinParameterEnum
10041007
syn keyword juliaConstC InversionSurfaceObsEnum
10051008
syn keyword juliaConstC InversionThicknessObsEnum
10061009
syn keyword juliaConstC InversionVelObsEnum

0 commit comments

Comments
 (0)