Ice Sheet System Model  4.18
Code documentation
control_core.cpp
Go to the documentation of this file.
1 
5 #include "./cores.h"
6 #include "../toolkits/toolkits.h"
7 #include "../classes/classes.h"
8 #include "../shared/shared.h"
9 #include "../modules/modules.h"
10 #include "../solutionsequences/solutionsequences.h"
11 
12 /*Local prototypes*/
13 /*{{{*/
14 IssmDouble FormFunction(IssmDouble* X,void* usr);
15 IssmDouble FormFunctionGradient(IssmDouble** pG,IssmDouble* X,void* usr);
16 typedef struct {
17  FemModel* femmodel;
18  int nsize;
19 } AppCtx;
20 /*}}}*/
21 
22 void control_core(FemModel* femmodel){/*{{{*/
23 
24  /*parameters: */
25  int num_controls,nsize,nsteps;
26  int solution_type;
27  bool isFS,dakota_analysis;
28  int *control_type = NULL;
29  int *maxiter = NULL;
30  IssmDouble *cm_jump = NULL;
31  IssmDouble *J = NULL;
32 
33  /*Solution and Adjoint core pointer*/
34  void (*solutioncore)(FemModel*) = NULL;
35  void (*adjointcore)(FemModel*) = NULL;
36 
37  /*Recover parameters used throughout the solution*/
38  femmodel->parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
39  femmodel->parameters->FindParam(&control_type,NULL,InversionControlParametersEnum);
40  femmodel->parameters->FindParam(&nsteps,InversionNstepsEnum);
41  femmodel->parameters->FindParam(&maxiter,NULL,InversionMaxiterPerStepEnum);
42  femmodel->parameters->FindParam(&cm_jump,NULL,InversionStepThresholdEnum);
43  femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
44  femmodel->parameters->FindParam(&isFS,FlowequationIsFSEnum);
45  femmodel->parameters->FindParam(&dakota_analysis,QmuIsdakotaEnum);
46  femmodel->parameters->SetParam(false,SaveResultsEnum);
47 
48  /*out of solution_type, figure out solution core and adjoint function pointer*/
49  CorePointerFromSolutionEnum(&solutioncore,femmodel->parameters,solution_type);
50  AdjointCorePointerFromSolutionEnum(&adjointcore,solution_type);
51 
52  /*Launch once a complete solution to set up all inputs*/
53  if(VerboseControl()) _printf0_(" preparing initial solution\n");
54  if(isFS) solutioncore(femmodel);
55 
56  /*Get initial guess*/
57  Vector<IssmDouble> *Xpetsc = NULL;
58  GetVectorFromControlInputsx(&Xpetsc,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"value");
59  IssmDouble* X0 = Xpetsc->ToMPISerial();
60  Xpetsc->GetSize(&nsize);
61  delete Xpetsc;
62 
63  /*Initialize some of the BrentSearch arguments: */
64  OptPars optpars;
65  optpars.xmin = 0;
66  optpars.xmax = 1;
67  optpars.nsteps = nsteps;
68  optpars.nsize = nsize;
69  optpars.maxiter = maxiter;
70  optpars.cm_jump = cm_jump;
71 
72  /*Initialize function argument*/
73  AppCtx usr;
74  usr.femmodel = femmodel;
75  usr.nsize = nsize;
76 
77  /*Call Brent optimization*/
78  BrentSearch(&J,optpars,X0,&FormFunction,&FormFunctionGradient,(void*)&usr);
79 
80  if(VerboseControl()) _printf0_(" preparing final solution\n");
81  IssmDouble *XL = NULL;
82  IssmDouble *XU = NULL;
83  GetVectorFromControlInputsx(&XL,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"lowerbound");
84  GetVectorFromControlInputsx(&XU,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"upperbound");
85  for(long i=0;i<nsize;i++){
86  if(X0[i]>XU[i]) X0[i]=XU[i];
87  if(X0[i]<XL[i]) X0[i]=XL[i];
88  }
89  xDelete<IssmDouble>(XU);
90  xDelete<IssmDouble>(XL);
91  SetControlInputsFromVectorx(femmodel,X0);
92  femmodel->parameters->SetParam(true,SaveResultsEnum);
93  solutioncore(femmodel);
94 
95  /*some results not computed by steadystate_core or stressbalance_core: */
96  if(!dakota_analysis){ //do not save this if we are running the control core from a qmu run!
97  femmodel->OutputControlsx(&femmodel->results);
98 
99  #ifdef _HAVE_AD_
100  IssmPDouble* J_passive=xNew<IssmPDouble>(nsteps);
101  for(int i=0;i<nsteps;i++) J_passive[i]=reCast<IssmPDouble>(J[i]);
102  femmodel->results->AddObject(new GenericExternalResult<IssmPDouble*>(femmodel->results->Size()+1,JEnum,J_passive,nsteps,1,0,0));
103  xDelete<IssmPDouble>(J_passive);
104  #else
105  femmodel->results->AddObject(new GenericExternalResult<IssmPDouble*>(femmodel->results->Size()+1,JEnum,J,nsteps,1,0,0));
106  #endif
107  }
108 
109  /*Free ressources: */
110  xDelete<int>(control_type);
111  xDelete<int>(maxiter);
112  xDelete<IssmDouble>(cm_jump);
113  xDelete<IssmDouble>(J);
114  xDelete<IssmDouble>(X0);
115 }/*}}}*/
116 IssmDouble FormFunction(IssmDouble* X,void* usrvoid){/*{{{*/
117 
118  /*output: */
119  IssmDouble J;
120 
121  /*parameters: */
122  int solution_type,analysis_type,num_responses;
123  bool conserve_loads = true;
124  AppCtx* usr = (AppCtx*)usrvoid;
125  FemModel *femmodel = usr->femmodel;
126  int nsize = usr->nsize;
127 
128  /*Recover parameters: */
129  femmodel->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
130  femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
131  femmodel->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum);
132 
133  /*Constrain input vector*/
134  IssmDouble *XL = NULL;
135  IssmDouble *XU = NULL;
136  GetVectorFromControlInputsx(&XL,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"lowerbound");
137  GetVectorFromControlInputsx(&XU,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"upperbound");
138  for(long i=0;i<nsize;i++){
139  if(X[i]>XU[i]) X[i]=XU[i];
140  if(X[i]<XL[i]) X[i]=XL[i];
141  }
142 
143  /*Update control input*/
144  SetControlInputsFromVectorx(femmodel,X);
145 
146  /*solve forward: */
147  switch(solution_type){
148  case SteadystateSolutionEnum:
149  femmodel->SetCurrentConfiguration(StressbalanceAnalysisEnum);
150  stressbalance_core(femmodel); //We need a 3D velocity!! (vz is required for the next thermal run)
151  break;
152  case StressbalanceSolutionEnum:{
153  femmodel->SetCurrentConfiguration(StressbalanceAnalysisEnum);
154 
155  bool is_schur_cg_solver = false;
156  #ifdef _HAVE_PETSC_
157  int solver_type;
158  PetscOptionsDetermineSolverType(&solver_type);
159  if(solver_type==FSSolverEnum) is_schur_cg_solver = true;
160  #endif
161 
162  if(is_schur_cg_solver){
163  solutionsequence_schurcg(femmodel);
164  }else{
165  solutionsequence_nonlinear(femmodel,conserve_loads);
166  }
167  }
168  break;
169  case BalancethicknessSolutionEnum:
170  femmodel->SetCurrentConfiguration(BalancethicknessAnalysisEnum);
171  solutionsequence_linear(femmodel);
172  break;
173  case BalancethicknessSoftSolutionEnum:
174  /*NOTHING*/
175  break;
176  case Balancethickness2SolutionEnum:
177  femmodel->SetCurrentConfiguration(Balancethickness2AnalysisEnum);
178  solutionsequence_linear(femmodel);
179  break;
180  default:
181  _error_("Solution " << EnumToStringx(solution_type) << " not implemented yet");
182  }
183 
184  /*Compute misfit for this velocity field.*/
185  IssmDouble* Jlist = NULL;
186  femmodel->CostFunctionx(&J,&Jlist,NULL);
187  _printf0_("f(x) = "<<setw(12)<<setprecision(7)<<J<<" | ");
188  for(int i=0;i<num_responses;i++) _printf0_(" "<<setw(12)<<setprecision(7)<<Jlist[i]);
189  _printf0_("\n");
190 
191  /*Free ressources:*/
192  xDelete<IssmDouble>(XU);
193  xDelete<IssmDouble>(XL);
194  xDelete<IssmDouble>(Jlist);
195  return J;
196 }/*}}}*/
197 IssmDouble FormFunctionGradient(IssmDouble** pG,IssmDouble* X,void* usrvoid){/*{{{*/
198 
199  /*output: */
200  IssmDouble J;
201  int temp;
202 
203  /*parameters: */
204  void (*adjointcore)(FemModel*)=NULL;
205  int solution_type,analysis_type,num_responses,num_controls,numvertices;
206  bool conserve_loads = true;
207  IssmDouble *scalar_list = NULL;
208  IssmDouble *Jlist = NULL;
209  IssmDouble *G = NULL;
210  IssmDouble *norm_list = NULL;
211  AppCtx *usr = (AppCtx*)usrvoid;
212  FemModel *femmodel = usr->femmodel;
213  int nsize = usr->nsize;
214 
215  /*Recover parameters: */
216  femmodel->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
217  femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
218  femmodel->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum);
219  femmodel->parameters->FindParam(&scalar_list,&temp,&temp,InversionGradientScalingEnum);
220  femmodel->parameters->FindParam(&num_controls,InversionNumControlParametersEnum); _assert_(num_controls);
221  numvertices=femmodel->vertices->NumberOfVertices();
222 
223  /*Constrain input vector*/
224  IssmDouble *XL = NULL;
225  IssmDouble *XU = NULL;
226  GetVectorFromControlInputsx(&XL,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"lowerbound");
227  GetVectorFromControlInputsx(&XU,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"upperbound");
228  for(long i=0;i<nsize;i++){
229  if(X[i]>XU[i]) X[i]=XU[i];
230  if(X[i]<XL[i]) X[i]=XL[i];
231  }
232 
233  /*Update control input*/
234  SetControlInputsFromVectorx(femmodel,X);
235 
236  /*Compute new temperature at this point*/
237  if(solution_type==SteadystateSolutionEnum) steadystate_core(femmodel);
238 
239  /*Compute Adjoint*/
240  AdjointCorePointerFromSolutionEnum(&adjointcore,solution_type);
241  adjointcore(femmodel);
242 
243  /*Compute gradient*/
244  Gradjx(&G,&norm_list,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
245 
246  /*Compute scaling factor*/
247  IssmDouble scalar = scalar_list[0]/norm_list[0];
248  for(int i=1;i<num_controls;i++) scalar=min(scalar,scalar_list[i]/norm_list[i]);
249 
250  /*Constrain Gradient*/
251  for(int i=0;i<num_controls;i++){
252  for(int j=0;j<numvertices;j++){
253  G[i*numvertices+j] = scalar*G[i*numvertices+j];
254  }
255  }
256 
257  for(long i=0;i<nsize;i++){
258  if(X[i]>=XU[i]) G[i]=0.;
259  if(X[i]<=XL[i]) G[i]=0.;
260  }
261 
262  /*Needed for output results*/
263  ControlInputSetGradientx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,G);
264 
265  /*Compute misfit for this velocity field.*/
266  femmodel->CostFunctionx(&J,&Jlist,NULL);
267  _printf0_("f(x) = "<<setw(12)<<setprecision(7)<<J<<" | ");
268  for(int i=0;i<num_responses;i++) _printf0_(" "<<setw(12)<<setprecision(7)<<Jlist[i]);
269  _printf0_("\n");
270 
271  /*Clean-up and return*/
272  xDelete<IssmDouble>(XU);
273  xDelete<IssmDouble>(XL);
274  xDelete<IssmDouble>(norm_list);
275  xDelete<IssmDouble>(scalar_list);
276  xDelete<IssmDouble>(Jlist);
277  *pG = G;
278  return J;
279 }/*}}}*/
DataSet::Size
int Size()
Definition: DataSet.cpp:399
GetVectorFromControlInputsx
void GetVectorFromControlInputsx(Vector< IssmDouble > **pvector, Elements *elements, Nodes *nodes, Vertices *vertices, Loads *loads, Materials *materials, Parameters *parameters, const char *data)
Definition: GetVectorFromControlInputsx.cpp:9
BalancethicknessAnalysisEnum
@ BalancethicknessAnalysisEnum
Definition: EnumDefinitions.h:981
SaveResultsEnum
@ SaveResultsEnum
Definition: EnumDefinitions.h:302
ControlInputSetGradientx
void ControlInputSetGradientx(Elements *elements, Nodes *nodes, Vertices *vertices, Loads *loads, Materials *materials, Parameters *parameters, IssmDouble *gradient)
Definition: ControlInputSetGradientx.cpp:9
OptPars::nsize
int nsize
Definition: OptPars.h:17
_assert_
#define _assert_(ignore)
Definition: exceptions.h:37
FemModel::vertices
Vertices * vertices
Definition: FemModel.h:49
Balancethickness2SolutionEnum
@ Balancethickness2SolutionEnum
Definition: EnumDefinitions.h:980
IssmDouble
double IssmDouble
Definition: types.h:37
InversionNumControlParametersEnum
@ InversionNumControlParametersEnum
Definition: EnumDefinitions.h:223
_printf0_
#define _printf0_(StreamArgs)
Definition: Print.h:29
PetscOptionsDetermineSolverType
void PetscOptionsDetermineSolverType(int *psolver_type)
Definition: PetscOptionsDetermineSolverType.cpp:20
StressbalanceAnalysisEnum
@ StressbalanceAnalysisEnum
Definition: EnumDefinitions.h:1285
AppCtx
Definition: control_core.cpp:16
InversionNstepsEnum
@ InversionNstepsEnum
Definition: EnumDefinitions.h:222
FormFunction
IssmDouble FormFunction(IssmDouble *X, void *usr)
Definition: control_core.cpp:116
DataSet::AddObject
int AddObject(Object *object)
Definition: DataSet.cpp:252
InversionControlParametersEnum
@ InversionControlParametersEnum
Definition: EnumDefinitions.h:209
AdjointCorePointerFromSolutionEnum
void AdjointCorePointerFromSolutionEnum(void(**padjointcore)(FemModel *), int solutiontype)
Definition: AdjointCorePointerFromSolutionEnum.cpp:18
cores.h
Balancethickness2AnalysisEnum
@ Balancethickness2AnalysisEnum
Definition: EnumDefinitions.h:979
BalancethicknessSoftSolutionEnum
@ BalancethicknessSoftSolutionEnum
Definition: EnumDefinitions.h:984
FemModel::parameters
Parameters * parameters
Definition: FemModel.h:46
CorePointerFromSolutionEnum
void CorePointerFromSolutionEnum(void(**psolutioncore)(FemModel *), Parameters *parameters, int solutiontype)
Definition: CorePointerFromSolutionEnum.cpp:18
solutionsequence_schurcg
void solutionsequence_schurcg(FemModel *femmodel)
Definition: solutionsequence_schurcg.cpp:834
FemModel::CostFunctionx
void CostFunctionx(IssmDouble *pJ, IssmDouble **pJlist, int *pn)
Definition: FemModel.cpp:1062
JEnum
@ JEnum
Definition: EnumDefinitions.h:1134
FemModel::results
Results * results
Definition: FemModel.h:48
OptPars::maxiter
int * maxiter
Definition: OptPars.h:15
Vertices::NumberOfVertices
int NumberOfVertices(void)
Definition: Vertices.cpp:255
BalancethicknessSolutionEnum
@ BalancethicknessSolutionEnum
Definition: EnumDefinitions.h:985
FSSolverEnum
@ FSSolverEnum
Definition: EnumDefinitions.h:1061
FemModel::OutputControlsx
void OutputControlsx(Results **presults)
Definition: FemModel.cpp:2171
solutionsequence_linear
void solutionsequence_linear(FemModel *femmodel)
Definition: solutionsequence_linear.cpp:10
steadystate_core
void steadystate_core(FemModel *femmodel)
Definition: steadystate_core.cpp:20
InversionMaxiterPerStepEnum
@ InversionMaxiterPerStepEnum
Definition: EnumDefinitions.h:220
solutionsequence_nonlinear
void solutionsequence_nonlinear(FemModel *femmodel, bool conserve_loads)
Definition: solutionsequence_nonlinear.cpp:11
AppCtx::nsize
int nsize
Definition: control_core.cpp:18
BrentSearch
void BrentSearch(IssmDouble **pJ, OptPars optpars, IssmDouble *X0, IssmDouble(*f)(IssmDouble *, void *), IssmDouble(*g)(IssmDouble **, IssmDouble *, void *), void *usr)
Definition: BrentSearch.cpp:23
Vector::GetSize
void GetSize(int *pM)
Definition: Vector.h:185
control_core
void control_core(FemModel *femmodel)
Definition: control_core.cpp:22
OptPars::xmax
IssmDouble xmax
Definition: OptPars.h:13
Parameters::SetParam
void SetParam(bool boolean, int enum_type)
Definition: Parameters.cpp:441
FemModel::nodes
Nodes * nodes
Definition: FemModel.h:56
FemModel::materials
Materials * materials
Definition: FemModel.h:45
EnumToStringx
const char * EnumToStringx(int enum_in)
Definition: EnumToStringx.cpp:15
StressbalanceSolutionEnum
@ StressbalanceSolutionEnum
Definition: EnumDefinitions.h:1288
Gradjx
void Gradjx(Vector< IssmDouble > **pgradient, IssmDouble **pnorm_list, Elements *elements, Nodes *nodes, Vertices *vertices, Loads *loads, Materials *materials, Parameters *parameters)
Definition: Gradjx.cpp:9
QmuIsdakotaEnum
@ QmuIsdakotaEnum
Definition: EnumDefinitions.h:288
SolutionTypeEnum
@ SolutionTypeEnum
Definition: EnumDefinitions.h:398
FemModel::loads
Loads * loads
Definition: FemModel.h:54
OptPars
Definition: OptPars.h:10
FemModel::elements
Elements * elements
Definition: FemModel.h:44
FemModel
Definition: FemModel.h:31
InversionStepThresholdEnum
@ InversionStepThresholdEnum
Definition: EnumDefinitions.h:225
InversionGradientScalingEnum
@ InversionGradientScalingEnum
Definition: EnumDefinitions.h:214
GenericExternalResult
Definition: GenericExternalResult.h:21
OptPars::nsteps
int nsteps
Definition: OptPars.h:16
OptPars::xmin
IssmDouble xmin
Definition: OptPars.h:12
_error_
#define _error_(StreamArgs)
Definition: exceptions.h:49
FlowequationIsFSEnum
@ FlowequationIsFSEnum
Definition: EnumDefinitions.h:138
min
IssmDouble min(IssmDouble a, IssmDouble b)
Definition: extrema.cpp:14
InversionNumCostFunctionsEnum
@ InversionNumCostFunctionsEnum
Definition: EnumDefinitions.h:224
AppCtx::femmodel
FemModel * femmodel
Definition: control_core.cpp:17
FemModel::SetCurrentConfiguration
void SetCurrentConfiguration(int configuration_type)
Definition: FemModel.cpp:634
Parameters::FindParam
void FindParam(bool *pinteger, int enum_type)
Definition: Parameters.cpp:262
OptPars::cm_jump
IssmDouble * cm_jump
Definition: OptPars.h:14
SetControlInputsFromVectorx
void SetControlInputsFromVectorx(FemModel *femmodel, IssmDouble *vector)
Definition: SetControlInputsFromVectorx.cpp:9
AnalysisTypeEnum
@ AnalysisTypeEnum
Definition: EnumDefinitions.h:36
FormFunctionGradient
IssmDouble FormFunctionGradient(IssmDouble **pG, IssmDouble *X, void *usr)
Definition: control_core.cpp:197
SteadystateSolutionEnum
@ SteadystateSolutionEnum
Definition: EnumDefinitions.h:1283
Vector::ToMPISerial
doubletype * ToMPISerial(void)
Definition: Vector.h:277
IssmPDouble
IssmDouble IssmPDouble
Definition: types.h:38
stressbalance_core
void stressbalance_core(FemModel *femmodel)
Definition: stressbalance_core.cpp:13
VerboseControl
bool VerboseControl(void)
Definition: Verbosity.cpp:27
Vector< IssmDouble >
femmodel
FemModel * femmodel
Definition: esmfbinders.cpp:16
Generated on Thu Jul 2 2020 08:09:18 for Ice Sheet System Model by   doxygen 1.8.19