Ice Sheet System Model  4.18
Code documentation
controladm1qn3_core.cpp
Go to the documentation of this file.
1 
4 #include <ctime>
5 #include <config.h>
6 #include "./cores.h"
7 #include "../toolkits/toolkits.h"
8 #include "../classes/classes.h"
9 #include "../shared/shared.h"
10 #include "../modules/modules.h"
11 #include "../solutionsequences/solutionsequences.h"
12 
13 #ifdef _HAVE_CODIPACK_
14 extern CoDi_global codi_global;
15 #include <sstream> // for output of the CoDiPack tape
16 #endif
17 
18 #if defined (_HAVE_M1QN3_) && defined(_HAVE_AD_)
19 /*m1qn3 prototypes {{{*/
20 extern "C" void *ctonbe_; // DIS mode : Conversion
21 extern "C" void *ctcabe_; // DIS mode : Conversion
22 extern "C" void *euclid_; // Scalar product
23 typedef void (*SimulFunc) (long* indic,long* n, double* x,double* pf,double* g,long [],float [],void* dzs);
24 extern "C" void m1qn3_ (void f(long* indic,long* n, double* x,double* pf,double* g,long [],float [],void* dzs),
25  void **, void **, void **,
26  long *, double [],double *, double[], double*, double *,
27  double *, char [], long *, long *, long *, long *, long *, long *, long [],double [], long *,
28  long *, long *, long [], float [],void* );
29 
30 /*Use struct to provide arguments*/
31 typedef struct{
32  FemModel * femmodel;
33  IssmPDouble* Jlist;
34  int M;
35  int N;
36  int* i;
37 } m1qn3_struct;
38 /*}}}*/
39 
40 /*m1qm3 functions*/
41 void simul_starttrace(FemModel* femmodel){/*{{{*/
42 
43  #if defined(_HAVE_ADOLC_)
44  /*Retrive ADOLC parameters*/
45  IssmDouble gcTriggerRatio;
46  IssmDouble gcTriggerMaxSize;
47  IssmDouble obufsize;
48  IssmDouble lbufsize;
49  IssmDouble cbufsize;
50  IssmDouble tbufsize;
51  femmodel->parameters->FindParam(&gcTriggerRatio,AutodiffGcTriggerRatioEnum);
52  femmodel->parameters->FindParam(&gcTriggerMaxSize,AutodiffGcTriggerMaxSizeEnum);
53  femmodel->parameters->FindParam(&obufsize,AutodiffObufsizeEnum);
54  femmodel->parameters->FindParam(&lbufsize,AutodiffLbufsizeEnum);
55  femmodel->parameters->FindParam(&cbufsize,AutodiffCbufsizeEnum);
56  femmodel->parameters->FindParam(&tbufsize,AutodiffTbufsizeEnum);
57 
58  /*Set garbage collection parameters: */
59  setStoreManagerControl(reCast<IssmPDouble>(gcTriggerRatio),reCast<size_t>(gcTriggerMaxSize));
60 
61  /*Start trace: */
62  int skipFileDeletion=1;
63  int keepTaylors=1;
64  int my_rank=IssmComm::GetRank();
65  trace_on(my_rank,keepTaylors,reCast<size_t>(obufsize),reCast<size_t>(lbufsize),reCast<size_t>(cbufsize),reCast<size_t>(tbufsize),skipFileDeletion);
66 
67  #elif defined(_HAVE_CODIPACK_)
68 
69  //fprintf(stderr, "*** Codipack IoModel::StartTrace\n");
70  /*
71  * FIXME codi
72  * - ADOL-C variant uses fine grained tracing with various arguments
73  * - ADOL-C variant sets a garbage collection parameter for its tape
74  * -> These parameters are not read for the CoDiPack ISSM version!
75  */
76  auto& tape_codi = IssmDouble::getGlobalTape();
77  tape_codi.setActive();
78  #if _AD_TAPE_ALLOC_
79  //alloc_profiler.Tag(StartInit, true);
80  IssmDouble x_t(1.0), y_t(1.0);
81  tape_codi.registerInput(y_t);
82  int codi_allocn = 0;
83  femmodel->parameters->FindParam(&codi_allocn,AutodiffTapeAllocEnum);
84  for(int i = 0;i < codi_allocn;++i) {
85  x_t = y_t * y_t;
86  }
87  /*
88  std::stringstream out_s;
89  IssmDouble::getGlobalTape().printStatistics(out_s);
90  _printf0_("StartTrace::Tape Statistics : TapeAlloc count=[" << codi_allocn << "]\n" << out_s.str());
91  */
92  tape_codi.reset();
93  //alloc_profiler.Tag(FinishInit, true);
94  #else
95  tape_codi.reset();
96  #endif
97 
98  #else
99  _error_("not implemented");
100  #endif
101 }/*}}}*/
102 void simul_stoptrace(){/*{{{*/
103 
104  #if defined(_HAVE_ADOLC_)
105  trace_off();
106  if(VerboseAutodiff()){ /*{{{*/
107 
108  #ifdef _HAVE_ADOLC_
109  int my_rank=IssmComm::GetRank();
110  size_t tape_stats[15];
111  tapestats(my_rank,tape_stats); //reading of tape statistics
112  int commSize=IssmComm::GetSize();
113  int *sstats=new int[7];
114  sstats[0]=tape_stats[NUM_OPERATIONS];
115  sstats[1]=tape_stats[OP_FILE_ACCESS];
116  sstats[2]=tape_stats[NUM_LOCATIONS];
117  sstats[3]=tape_stats[LOC_FILE_ACCESS];
118  sstats[4]=tape_stats[NUM_VALUES];
119  sstats[5]=tape_stats[VAL_FILE_ACCESS];
120  sstats[6]=tape_stats[TAY_STACK_SIZE];
121  int *rstats=NULL;
122  if (my_rank==0) rstats=new int[commSize*7];
123  ISSM_MPI_Gather(sstats,7,ISSM_MPI_INT,rstats,7,ISSM_MPI_INT,0,IssmComm::GetComm());
124  if (my_rank==0) {
125  int offset=50;
126  int rOffset=(commSize/10)+1;
127  _printf_(" ADOLC statistics: \n");
128  _printf_(" "<<setw(offset)<<left<<"#independents: " <<setw(12)<<right<<tape_stats[NUM_INDEPENDENTS] << "\n");
129  _printf_(" "<<setw(offset)<<left<<"#dependents: " <<setw(12)<<right<<tape_stats[NUM_DEPENDENTS] << "\n");
130  _printf_(" "<<setw(offset)<<left<<"max #live active variables: " <<setw(12)<<right<<tape_stats[NUM_MAX_LIVES] << "\n");
131  _printf_(" operations: entry size "<< sizeof(unsigned char) << " Bytes \n");
132  _printf_(" "<<setw(offset)<<left<<" #entries in buffer (AutodiffObufsizeEnum) " <<setw(12)<<right<<tape_stats[OP_BUFFER_SIZE] << "\n");
133  for (int r=0;r<commSize;++r)
134  _printf_(" ["<<setw(rOffset)<<right<<r<<"]"<<setw(offset-rOffset-4)<<left<<" #entries total" <<setw(12)<<right<<rstats[r*7+0] << (rstats[r*7+1]?" ->file":"") << "\n");
135  _printf_(" locations: entry size " << sizeof(locint) << " Bytes\n");
136  _printf_(" "<<setw(offset)<<left<<" #entries in buffer (AutodiffLbufsizeEnum) " <<setw(12)<<right<<tape_stats[LOC_BUFFER_SIZE] << "\n");
137  for (int r=0;r<commSize;++r)
138  _printf_(" ["<<setw(rOffset)<<right<<r<<"]"<<setw(offset-rOffset-4)<<left<<" #entries total" <<setw(12)<<right<<rstats[r*7+2] << (rstats[r*7+3]?" ->file":"") << "\n");
139  _printf_(" constant values: entry size " << sizeof(double) << " Bytes\n");
140  _printf_(" "<<setw(offset)<<left<<" #entries in buffer (AutodiffCbufsizeEnum) " <<setw(12)<<right<<tape_stats[VAL_BUFFER_SIZE] << "\n");
141  for (int r=0;r<commSize;++r)
142  _printf_(" ["<<setw(rOffset)<<right<<r<<"]"<<setw(offset-rOffset-4)<<left<<" #entries total" <<setw(12)<<right<<rstats[r*7+4] << (rstats[r*7+5]?" ->file":"") << "\n");
143  _printf_(" Taylor stack: entry size " << sizeof(revreal) << " Bytes\n");
144  _printf_(" "<<setw(offset)<<left<<" #entries in buffer (AutodiffTbufsizeEnum) " <<setw(12)<<right<<tape_stats[TAY_BUFFER_SIZE] << "\n");
145  for (int r=0;r<commSize;++r)
146  _printf_(" ["<<setw(rOffset)<<right<<r<<"]"<<setw(offset-rOffset-4)<<left<<" #entries total" <<setw(12)<<right<<rstats[r*7+6] << (rstats[r*7+6]>tape_stats[TAY_BUFFER_SIZE]?" ->file":"") << "\n");
147  delete []rstats;
148  }
149  delete [] sstats;
150  #endif
151 
152  #ifdef _HAVE_CODIPACK_
153  #ifdef _AD_TAPE_ALLOC_
154  //_printf_("Allocation time P(" << my_rank << "): " << alloc_profiler.DeltaTime(StartInit, FinishInit) << "\n");
155  #endif
156  std::stringstream out_s;
157  IssmDouble::getGlobalTape().printStatistics(out_s);
158  _printf0_("CoDiPack Profiling::Tape Statistics :\n" << out_s.str());
159  #endif
160  } /*}}}*/
161 
162  #elif defined(_HAVE_CODIPACK_)
163  auto& tape_codi = IssmDouble::getGlobalTape();
164  tape_codi.setPassive();
165  if(VerboseAutodiff()){
166  int my_rank=IssmComm::GetRank();
167  if(my_rank == 0) {
168  // FIXME codi "just because" for now
169  tape_codi.printStatistics(std::cout);
170  codi_global.print(std::cout);
171  }
172  }
173  #else
174  _error_("not implemented");
175  #endif
176 }/*}}}*/
177 void simul_ad(long* indic,long* n,double* X,double* pf,double* G,long izs[1],float rzs[1],void* dzs){/*{{{*/
178 
179  /*Get rank*/
180  int my_rank=IssmComm::GetRank();
181 
182  /*Recover Arguments*/
183  m1qn3_struct *input_struct = (m1qn3_struct*)dzs;
184 
185  FemModel* femmodel = input_struct->femmodel;
186  int num_responses,num_controls,solution_type;
187  femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
188 
189  /*In transient, we need to make sure we do not modify femmodel at each iteration, make a copy*/
190  if(solution_type == TransientSolutionEnum) femmodel = input_struct->femmodel->copy();
191 
192  IssmPDouble* Jlist = input_struct->Jlist;
193  int JlistM = input_struct->M;
194  int JlistN = input_struct->N;
195  int* Jlisti = input_struct->i;
196  int intn = (int)*n;
197 
198  /*Recover some parameters*/
199  IssmDouble *scaling_factors = NULL;
200  int *N = NULL;
201  int *control_enum = NULL;
202  femmodel->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum);
203  femmodel->parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
204  femmodel->parameters->FindParam(&scaling_factors,NULL,InversionControlScalingFactorsEnum);
205  femmodel->parameters->FindParam(&N,NULL,ControlInputSizeNEnum);
206  femmodel->parameters->FindParam(&control_enum,NULL,InversionControlParametersEnum);
207  int numberofvertices=femmodel->vertices->NumberOfVertices();
208 
209  /*Constrain input vector and update controls*/
210  double *XL = NULL;
211  double *XU = NULL;
212  GetPassiveVectorFromControlInputsx(&XL,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"lowerbound");
213  GetPassiveVectorFromControlInputsx(&XU,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"upperbound");
214 
215  int N_add = 0;
216  for (int c=0;c<num_controls;c++){
217  for(int i=0;i<numberofvertices*N[c];i++){
218  int index = N_add*numberofvertices+i;
219  X[index] = X[index]*reCast<double>(scaling_factors[c]);
220  if(X[index]>XU[index]) X[index]=XU[index];
221  if(X[index]<XL[index]) X[index]=XL[index];
222  }
223  N_add+=N[c];
224  }
225 
226  /*Start Tracing*/
227  simul_starttrace(femmodel);
228  /*Set X as our new control input and as INDEPENDENT*/
229 #ifdef _HAVE_AD_
230  IssmDouble* aX=xNew<IssmDouble>(intn,"t");
231 #else
232  IssmDouble* aX=xNew<IssmDouble>(intn);
233 #endif
234 
235  #if defined(_HAVE_ADOLC_)
236  if(my_rank==0){
237  for(int i=0;i<intn;i++){
238  aX[i]<<=X[i];
239  }
240  }
241  #elif defined(_HAVE_CODIPACK_)
242  auto& tape_codi = IssmDouble::getGlobalTape();
243  codi_global.input_indices.clear();
244  if(my_rank==0){
245  for (int i=0;i<intn;i++) {
246  aX[i]=X[i];
247  tape_codi.registerInput(aX[i]);
248  codi_global.input_indices.push_back(aX[i].getGradientData());
249  }
250  }
251  #else
252  _error_("not suppoted");
253  #endif
254 
255  ISSM_MPI_Bcast(aX,intn,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
256  SetControlInputsFromVectorx(femmodel,aX);
257  xDelete<IssmDouble>(aX);
258 
259  /*Compute solution (forward)*/
260  void (*solutioncore)(FemModel*)=NULL;
261  CorePointerFromSolutionEnum(&solutioncore,femmodel->parameters,solution_type);
262  solutioncore(femmodel);
263 
264  /*Reset the time to zero for next optimization*/
265  if(solution_type==TransientSolutionEnum){
266  IssmDouble restart_time;
267  femmodel->parameters->FindParam(&restart_time,TimesteppingStartTimeEnum);
268  femmodel->parameters->SetParam(restart_time,TimeEnum);
269  }
270 
271  /*Get Dependents*/
272  IssmDouble output_value;
273  int num_dependents;
274  IssmPDouble *dependents;
275  DataSet *dependent_objects = NULL;
276  IssmDouble J = 0.;
277  femmodel->parameters->FindParam(&num_dependents,AutodiffNumDependentsEnum);
278  femmodel->parameters->FindParam(&dependent_objects,AutodiffDependentObjectsEnum);
279 
280  /*Go through our dependent variables, and compute the response:*/
281  dependents=xNew<IssmPDouble>(num_dependents);
282  #if defined(_HAVE_CODIPACK_)
283  codi_global.output_indices.clear();
284  #endif
285  for(int i=0;i<dependent_objects->Size();i++){
286  DependentObject* dep=(DependentObject*)dependent_objects->GetObjectByOffset(i);
287  if(solution_type==TransientSolutionEnum) output_value = dep->GetValue();
288  if(solution_type!=TransientSolutionEnum) dep->Responsex(&output_value,femmodel);
289  if(my_rank==0) {
290 
291  #if defined(_HAVE_CODIPACK_)
292  tape_codi.registerOutput(output_value);
293  dependents[i] = output_value.getValue();
294  codi_global.output_indices.push_back(output_value.getGradientData());
295 
296  #elif defined(_HAVE_ADOLC_)
297  output_value>>=dependents[i];
298 
299  #else
300  _error_("not suppoted");
301  #endif
302  J+=output_value;
303  }
304  }
305 
306  /*Turning off trace tape*/
307  simul_stoptrace();
308  //time_t now = time(NULL);
309  //if(my_rank==0) _printf_("\nTIME: "<<now<<"\n");
310 
311  /*intermediary: */
312  int num_independents=intn;
313  IssmPDouble *aWeightVector=NULL;
314  IssmPDouble *weightVectorTimesJac=NULL;
315  IssmPDouble *totalgradient=xNewZeroInit<IssmPDouble>(num_independents);
316 
317  /*if no dependents, no point in running a driver: */
318  if(!(num_dependents*num_independents)) _error_("this is not allowed");
319 
320  /*for adolc to run in parallel, we 0 out on rank~=0. But we still keep track of num_dependents:*/
321  int num_dependents_old = num_dependents;
322  int num_independents_old = num_independents;
323 
324  #if defined(_HAVE_ADOLC_)
325  /*Get gradient for ADOLC {{{*/
326  if(my_rank!=0){
327  num_dependents = 0;
328  num_independents = 0;
329  }
330 
331  /*get the EDF pointer:*/
332  ext_diff_fct *anEDF_for_solverx_p=xDynamicCast<GenericParam<Adolc_edf> * >(femmodel->parameters->FindParamObject(AdolcParamEnum))->GetParameterValue().myEDF_for_solverx_p;
333 
334  /* these are always needed regardless of the interpreter */
335  anEDF_for_solverx_p->dp_x=xNew<double>(anEDF_for_solverx_p->max_n);
336  anEDF_for_solverx_p->dp_y=xNew<double>(anEDF_for_solverx_p->max_m);
337 
338  /* Ok, now we are going to call the fos_reverse in a loop on the index, from 0 to num_dependents, so
339  * as to generate num_dependents gradients: */
340  for(int aDepIndex=0;aDepIndex<num_dependents_old;aDepIndex++){
341 
342  /*initialize direction index in the weights vector: */
343  aWeightVector=xNewZeroInit<IssmPDouble>(num_dependents);
344  if (my_rank==0) aWeightVector[aDepIndex]=1.;
345 
346  /*initialize output gradient: */
347  weightVectorTimesJac=xNew<IssmPDouble>(num_independents);
348 
349  /*set the forward method function pointer: */
350  #ifdef _HAVE_GSL_
351  anEDF_for_solverx_p->fos_reverse=EDF_fos_reverse_for_solverx;
352  #endif
353  #ifdef _HAVE_MUMPS_
354  anEDF_for_solverx_p->fos_reverse_iArr=fos_reverse_mumpsSolveEDF;
355  #endif
356 
357  anEDF_for_solverx_p->dp_U=xNew<IssmPDouble>(anEDF_for_solverx_p->max_m);
358  anEDF_for_solverx_p->dp_Z=xNew<IssmPDouble>(anEDF_for_solverx_p->max_n);
359 
360  /*call driver: */
361  fos_reverse(my_rank,num_dependents,num_independents, aWeightVector, weightVectorTimesJac );
362 
363  /*Add to totalgradient: */
364  if(my_rank==0) for(int i=0;i<num_independents;i++) {
365  totalgradient[i]+=weightVectorTimesJac[i];
366  }
367 
368  /*free resources :*/
369  xDelete(weightVectorTimesJac);
370  xDelete(aWeightVector);
371  }
372  /*}}}*/
373  #elif defined(_HAVE_CODIPACK_)
374  /*Get gradient for CoDiPack{{{*/
375  if(VerboseAutodiff())_printf0_(" CoDiPack fos_reverse\n");
376 
377  /* call the fos_reverse in a loop on the index, from 0 to num_dependents, so
378  * as to generate num_dependents gradients: */
379  for(int dep_index=0;dep_index<num_dependents_old;dep_index++){
380 
381  /*initialize direction index in the weights vector: */
382  if(my_rank==0){
383  if(dep_index<0 || dep_index>=num_dependents || codi_global.output_indices.size() <= dep_index){
384  _error_("index value for dependent index should be in [0,num_dependents-1]");
385  }
386  tape_codi.setGradient(codi_global.output_indices[dep_index],1.0);
387  }
388  tape_codi.evaluate();
389 
390  /*Get gradient for this dependent */
391  weightVectorTimesJac=xNew<IssmPDouble>(num_independents);
392  auto in_size = codi_global.input_indices.size();
393  for(size_t i = 0; i < in_size; ++i){
394  _assert_(i<num_independents);
395  weightVectorTimesJac[i] = tape_codi.getGradient(codi_global.input_indices[i]);
396  }
397  if(my_rank==0) for(int i=0;i<num_independents;i++){
398  totalgradient[i]+=weightVectorTimesJac[i];
399  }
400  xDelete(weightVectorTimesJac);
401  }
402  /*}}}*/
403  #else
404  _error_("not suppoted");
405  #endif
406 
407  /*Broadcast gradient to other ranks*/
408  ISSM_MPI_Bcast(totalgradient,num_independents_old,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm());
409  /*Check size of Jlist to avoid crashes*/
410  _assert_((*Jlisti)<JlistM);
411  _assert_(JlistN==num_responses+1);
412 
413  /*Compute objective function and broadcast it to other cpus*/
414  *pf = reCast<double>(J);
415  ISSM_MPI_Bcast(pf,1,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm());
416  _printf0_("f(x) = "<<setw(12)<<setprecision(7)<<*pf<<" | ");
417 
418  /*Record cost function values and delete Jtemp*/
419  for(int i=0;i<num_responses;i++) Jlist[(*Jlisti)*JlistN+i] = reCast<IssmPDouble>(dependents[i]);
420  Jlist[(*Jlisti)*JlistN+num_responses] = reCast<IssmPDouble>(J);
421 
422  if(*indic==0){
423  /*dry run, no gradient required*/
424 
425  /*Retrieve objective functions independently*/
426  _printf0_(" N/A |\n");
427  for(int i=0;i<num_responses;i++) _printf0_(" "<<setw(12)<<setprecision(7)<<Jlist[(*Jlisti)*JlistN+i]);
428  _printf0_("\n");
429 
430  *Jlisti = (*Jlisti) +1;
431  xDelete<double>(XU);
432  xDelete<double>(XL);
433  return;
434  }
435 
436  /*Compute gradient*/
437  for(long i=0;i<num_independents_old;i++) G[i] = totalgradient[i];
438 
439  /*Constrain Gradient*/
440  IssmDouble Gnorm = 0.;
441  N_add = 0;
442  for (int c=0;c<num_controls;c++){
443  for(int i=0;i<numberofvertices*N[c];i++){
444  int index = N_add*numberofvertices+i;
445  if(X[index]>=XU[index]) G[index]=0.;
446  if(X[index]<=XL[index]) G[index]=0.;
447  G[index] = G[index]*reCast<double>(scaling_factors[c]);
448  X[index] = X[index]/reCast<double>(scaling_factors[c]);
449  Gnorm += G[index]*G[index];
450  }
451  N_add+=N[c];
452  }
453  Gnorm = sqrt(Gnorm);
454 
455  /*Print info*/
456  _printf0_(" "<<setw(12)<<setprecision(7)<<Gnorm<<" |");
457  for(int i=0;i<num_responses;i++) _printf0_(" "<<setw(12)<<setprecision(7)<<Jlist[(*Jlisti)*JlistN+i]);
458  _printf0_("\n");
459 
460  /*Clean-up and return*/
461  *Jlisti = (*Jlisti) +1;
462  xDelete<double>(XU);
463  xDelete<double>(XL);
464  xDelete<int>(control_enum);
465  xDelete<int>(N);
466  xDelete<IssmDouble>(scaling_factors);
467  xDelete<IssmPDouble>(totalgradient);
468 }/*}}}*/
469 void controladm1qn3_core(FemModel* femmodel){/*{{{*/
470 
471  /*Intermediaries*/
472  long omode;
473  double f;
474  double dxmin,gttol;
475  int maxsteps,maxiter;
476  int intn,numberofvertices,num_controls,num_cost_functions,solution_type;
477  IssmDouble *scaling_factors = NULL;
478  double *X = NULL;
479  double *G = NULL;
480  int* N = NULL;
481  int offset = 0;
482  int N_add;
483  int* control_enum;
484 
485  /*Recover some parameters*/
486  femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
487  femmodel->parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
488  femmodel->parameters->FindParam(&num_cost_functions,InversionNumCostFunctionsEnum);
489  femmodel->parameters->FindParam(&maxsteps,InversionMaxstepsEnum);
490  femmodel->parameters->FindParam(&maxiter,InversionMaxiterEnum);
491  femmodel->parameters->FindParamAndMakePassive(&dxmin,InversionDxminEnum);
492  femmodel->parameters->FindParamAndMakePassive(&gttol,InversionGttolEnum);
493  femmodel->parameters->FindParam(&scaling_factors,NULL,InversionControlScalingFactorsEnum);
494  femmodel->parameters->FindParam(&control_enum,NULL,InversionControlParametersEnum);
495  femmodel->parameters->SetParam(false,SaveResultsEnum);
496  femmodel->parameters->FindParam(&N,NULL,ControlInputSizeNEnum);
497  numberofvertices=femmodel->vertices->NumberOfVertices();
498 
499  /*Initialize M1QN3 parameters*/
500  if(VerboseControl())_printf0_(" Initialize M1QN3 parameters\n");
501  SimulFunc simul_ptr = &simul_ad; /*Cost function address*/
502  void** prosca = &euclid_; /*Dot product function (euclid is the default)*/
503  char normtype[] = "dfn"; /*Norm type: dfn = scalar product defined by prosca*/
504  long izs[5]; /*Arrays used by m1qn3 subroutines*/
505  long iz[5]; /*Integer m1qn3 working array of size 5*/
506  float rzs[1]; /*Arrays used by m1qn3 subroutines*/
507  long impres = 0; /*verbosity level*/
508  long imode[3] = {0}; /*scaling and starting mode, 0 by default*/
509  long indic = 4; /*compute f and g*/
510  long reverse = 0; /*reverse or direct mode*/
511  long io = 6; /*Channel number for the output*/
512 
513  /*Optimization criterions*/
514  long niter = long(maxsteps); /*Maximum number of iterations*/
515  long nsim = long(maxiter);/*Maximum number of function calls*/
516 
517  /*Get initial guess*/
518  GetPassiveVectorFromControlInputsx(&X,&intn,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"value");
519  //_assert_(intn==numberofvertices*num_controls);
520 
521  /*Get problem dimension and initialize gradient and initial guess*/
522  long n = long(intn);
523  G = xNew<double>(n);
524 
525  /*Scale control for M1QN3*/
526  N_add = 0;
527  for (int c=0;c<num_controls;c++){
528  for(int i=0;i<numberofvertices*N[c];i++){
529  int index = N_add*numberofvertices+i;
530  X[index] = X[index]/reCast<IssmPDouble>(scaling_factors[c]);
531  }
532  N_add+=N[c];
533  }
534 
535  /*Allocate m1qn3 working arrays (see documentation)*/
536  long m = 100;
537  long ndz = 4*n+m*(2*n+1);
538  double* dz = xNew<double>(ndz);
539  if(VerboseControl())_printf0_(" Computing initial solution\n");
540  _printf0_("\n");
541  _printf0_("Cost function f(x) | Gradient norm |g(x)| | List of contributions\n");
542  _printf0_("____________________________________________________________________\n");
543 
544  /*Prepare structure for m1qn3*/
545  m1qn3_struct mystruct;
546  mystruct.femmodel = femmodel;
547  mystruct.M = maxiter;
548  mystruct.N = num_cost_functions+1;
549  mystruct.Jlist = xNewZeroInit<IssmPDouble>(mystruct.M*mystruct.N);
550  mystruct.i = xNewZeroInit<int>(1);
551  /*Initialize Gradient and cost function of M1QN3*/
552  indic = 4; /*gradient required*/
553  simul_ad(&indic,&n,X,&f,G,izs,rzs,(void*)&mystruct);
554  /*Estimation of the expected decrease in f during the first iteration*/
555  double df1=f;
556 
557  /*Call M1QN3 solver*/
558  m1qn3_(simul_ptr,prosca,&ctonbe_,&ctcabe_,
559  &n,X,&f,G,&dxmin,&df1,
560  &gttol,normtype,&impres,&io,imode,&omode,&niter,&nsim,iz,dz,&ndz,
561  &reverse,&indic,izs,rzs,(void*)&mystruct);
562  switch(int(omode)){
563  case 0: _printf0_(" Stop requested (indic = 0)\n"); break;
564  case 1: _printf0_(" Convergence reached (gradient satisfies stopping criterion)\n"); break;
565  case 2: _printf0_(" Bad initialization\n"); break;
566  case 3: _printf0_(" Line search failure\n"); break;
567  case 4: _printf0_(" Maximum number of iterations exceeded\n");break;
568  case 5: _printf0_(" Maximum number of function calls exceeded\n"); break;
569  case 6: _printf0_(" stopped on dxmin during line search\n"); break;
570  case 7: _printf0_(" <g,d> > 0 or <y,s> <0\n"); break;
571  default: _printf0_(" Unknown end condition\n");
572  }
573 
574  /*Constrain solution vector*/
575  double *XL = NULL;
576  double *XU = NULL;
577  GetPassiveVectorFromControlInputsx(&XL,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"lowerbound");
578  GetPassiveVectorFromControlInputsx(&XU,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"upperbound");
579 
580  N_add = 0;
581  for (int c=0;c<num_controls;c++){
582  for(int i=0;i<numberofvertices*N[c];i++){
583  int index = N_add*numberofvertices+i;
584  X[index] = X[index]*reCast<double>(scaling_factors[c]);
585  if(X[index]>XU[index]) X[index]=XU[index];
586  if(X[index]<XL[index]) X[index]=XL[index];
587  }
588  N_add +=N[c];
589  }
590 
591  /*Set X as our new control*/
592  IssmDouble* aX=xNew<IssmDouble>(intn);
593  IssmDouble* aG=xNew<IssmDouble>(intn);
594 
595  for(int i=0;i<intn;i++) {
596  aX[i] = reCast<IssmDouble>(X[i]);
597  aG[i] = reCast<IssmDouble>(G[i]);
598  }
599 
600  ControlInputSetGradientx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,aG);
601  SetControlInputsFromVectorx(femmodel,aX);
602 
603  xDelete(aX);
604 
605  if (solution_type == TransientSolutionEnum){
606  int step = 1;
607  femmodel->parameters->SetParam(step,StepEnum);
608  femmodel->results->AddObject(new GenericExternalResult<IssmPDouble*>(femmodel->results->Size()+1,JEnum,mystruct.Jlist,(*mystruct.i),mystruct.N,1,0));
609 
610  int offset = 0;
611  for(int i=0;i<num_controls;i++){
612 
613  /*Disect results*/
614  GenericExternalResult<IssmPDouble*>* G_output = new GenericExternalResult<IssmPDouble*>(femmodel->results->Size()+1,Gradient1Enum+i,&G[offset],N[i],numberofvertices,1,0.);
615  GenericExternalResult<IssmPDouble*>* X_output = new GenericExternalResult<IssmPDouble*>(femmodel->results->Size()+1,control_enum[i],&X[offset],N[i],numberofvertices,1,0.);
616 
617  /*transpose for consistency with MATLAB's formating*/
618  G_output->Transpose();
619  X_output->Transpose();
620 
621  /*Add to results*/
622  femmodel->results->AddObject(G_output);
623  femmodel->results->AddObject(X_output);
624 
625  offset += N[i]*numberofvertices;
626  }
627  }
628  else{
629  femmodel->results->AddObject(new GenericExternalResult<IssmPDouble*>(femmodel->results->Size()+1,JEnum,mystruct.Jlist,(*mystruct.i),mystruct.N,0,0));
630 
631  femmodel->OutputControlsx(&femmodel->results);
632  }
633 
634  xDelete(aG);
635 
636  /*Add last cost function to results*/
637 
638  /*Finalize*/
639  if(VerboseControl()) _printf0_(" preparing final solution\n");
640  femmodel->parameters->SetParam(true,SaveResultsEnum);
641  void (*solutioncore)(FemModel*)=NULL;
642  CorePointerFromSolutionEnum(&solutioncore,femmodel->parameters,solution_type);
643  solutioncore(femmodel);
644 
645  /*Clean-up and return*/
646  xDelete<double>(G);
647  xDelete<double>(X);
648  xDelete<double>(dz);
649  xDelete<double>(XU);
650  xDelete<double>(XL);
651  xDelete<IssmDouble>(scaling_factors);
652  xDelete<IssmPDouble>(mystruct.Jlist);
653  xDelete<int>(mystruct.i);
654 }/*}}}*/
655 
656 #else
657 void controladm1qn3_core(FemModel* femmodel){_error_("M1QN3 or ADOLC/CoDiPack not installed");}
658 #endif //_HAVE_M1QN3_
DataSet::Size
int Size()
Definition: DataSet.cpp:399
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
_assert_
#define _assert_(ignore)
Definition: exceptions.h:37
FemModel::vertices
Vertices * vertices
Definition: FemModel.h:49
GetPassiveVectorFromControlInputsx
void GetPassiveVectorFromControlInputsx(IssmPDouble **pvector, int *pN, Elements *elements, Nodes *nodes, Vertices *vertices, Loads *loads, Materials *materials, Parameters *parameters, const char *data)
Definition: GetVectorFromControlInputsx.cpp:127
IssmDouble
double IssmDouble
Definition: types.h:37
AutodiffDependentObjectsEnum
@ AutodiffDependentObjectsEnum
Definition: EnumDefinitions.h:43
InversionNumControlParametersEnum
@ InversionNumControlParametersEnum
Definition: EnumDefinitions.h:223
_printf0_
#define _printf0_(StreamArgs)
Definition: Print.h:29
StepEnum
@ StepEnum
Definition: EnumDefinitions.h:403
DataSet::AddObject
int AddObject(Object *object)
Definition: DataSet.cpp:252
_printf_
#define _printf_(StreamArgs)
Definition: Print.h:22
InversionControlParametersEnum
@ InversionControlParametersEnum
Definition: EnumDefinitions.h:209
cores.h
AutodiffObufsizeEnum
@ AutodiffObufsizeEnum
Definition: EnumDefinitions.h:54
TransientSolutionEnum
@ TransientSolutionEnum
Definition: EnumDefinitions.h:1317
FemModel::parameters
Parameters * parameters
Definition: FemModel.h:46
TimeEnum
@ TimeEnum
Definition: EnumDefinitions.h:427
IssmComm::GetComm
static ISSM_MPI_Comm GetComm(void)
Definition: IssmComm.cpp:30
CorePointerFromSolutionEnum
void CorePointerFromSolutionEnum(void(**psolutioncore)(FemModel *), Parameters *parameters, int solutiontype)
Definition: CorePointerFromSolutionEnum.cpp:18
JEnum
@ JEnum
Definition: EnumDefinitions.h:1134
FemModel::results
Results * results
Definition: FemModel.h:48
Vertices::NumberOfVertices
int NumberOfVertices(void)
Definition: Vertices.cpp:255
InversionMaxstepsEnum
@ InversionMaxstepsEnum
Definition: EnumDefinitions.h:221
Gradient1Enum
@ Gradient1Enum
Definition: EnumDefinitions.h:1086
FemModel::OutputControlsx
void OutputControlsx(Results **presults)
Definition: FemModel.cpp:2171
ISSM_MPI_DOUBLE
#define ISSM_MPI_DOUBLE
Definition: issmmpi.h:125
Parameters::SetParam
void SetParam(bool boolean, int enum_type)
Definition: Parameters.cpp:441
ISSM_MPI_INT
#define ISSM_MPI_INT
Definition: issmmpi.h:127
VerboseAutodiff
bool VerboseAutodiff(void)
Definition: Verbosity.cpp:29
DependentObject
Definition: DependentObject.h:15
AutodiffLbufsizeEnum
@ AutodiffLbufsizeEnum
Definition: EnumDefinitions.h:51
FemModel::nodes
Nodes * nodes
Definition: FemModel.h:56
xDelete
void xDelete(T **&aT_pp)
Definition: MemOps.h:97
FemModel::materials
Materials * materials
Definition: FemModel.h:45
IssmComm::GetSize
static int GetSize(void)
Definition: IssmComm.cpp:46
controladm1qn3_core
void controladm1qn3_core(FemModel *femmodel)
Definition: controladm1qn3_core.cpp:657
GenericExternalResult::Transpose
void Transpose(void)
Definition: GenericExternalResult.h:202
SolutionTypeEnum
@ SolutionTypeEnum
Definition: EnumDefinitions.h:398
AutodiffCbufsizeEnum
@ AutodiffCbufsizeEnum
Definition: EnumDefinitions.h:42
FemModel::loads
Loads * loads
Definition: FemModel.h:54
Parameters::FindParamAndMakePassive
void FindParamAndMakePassive(IssmPDouble *pscalar, int enum_type)
Definition: Parameters.cpp:379
ISSM_MPI_PDOUBLE
#define ISSM_MPI_PDOUBLE
Definition: issmmpi.h:126
FemModel::elements
Elements * elements
Definition: FemModel.h:44
ISSM_MPI_Bcast
int ISSM_MPI_Bcast(void *buffer, int count, ISSM_MPI_Datatype datatype, int root, ISSM_MPI_Comm comm)
Definition: issmmpi.cpp:162
DependentObject::Responsex
void Responsex(IssmDouble *poutput_value, FemModel *femmodel)
Definition: DependentObject.cpp:93
TimesteppingStartTimeEnum
@ TimesteppingStartTimeEnum
Definition: EnumDefinitions.h:432
FemModel
Definition: FemModel.h:31
GenericExternalResult
Definition: GenericExternalResult.h:21
IssmComm::GetRank
static int GetRank(void)
Definition: IssmComm.cpp:34
InversionControlScalingFactorsEnum
@ InversionControlScalingFactorsEnum
Definition: EnumDefinitions.h:210
_error_
#define _error_(StreamArgs)
Definition: exceptions.h:49
DataSet::GetObjectByOffset
Object * GetObjectByOffset(int offset)
Definition: DataSet.cpp:334
InversionNumCostFunctionsEnum
@ InversionNumCostFunctionsEnum
Definition: EnumDefinitions.h:224
AdolcParamEnum
@ AdolcParamEnum
Definition: EnumDefinitions.h:12
InversionMaxiterEnum
@ InversionMaxiterEnum
Definition: EnumDefinitions.h:219
AutodiffNumDependentsEnum
@ AutodiffNumDependentsEnum
Definition: EnumDefinitions.h:52
Parameters::FindParam
void FindParam(bool *pinteger, int enum_type)
Definition: Parameters.cpp:262
FemModel::copy
FemModel * copy()
Definition: FemModel.cpp:320
InversionDxminEnum
@ InversionDxminEnum
Definition: EnumDefinitions.h:212
SetControlInputsFromVectorx
void SetControlInputsFromVectorx(FemModel *femmodel, IssmDouble *vector)
Definition: SetControlInputsFromVectorx.cpp:9
AutodiffGcTriggerRatioEnum
@ AutodiffGcTriggerRatioEnum
Definition: EnumDefinitions.h:49
DependentObject::GetValue
IssmDouble GetValue(void)
Definition: DependentObject.cpp:102
InversionGttolEnum
@ InversionGttolEnum
Definition: EnumDefinitions.h:216
ControlInputSizeNEnum
@ ControlInputSizeNEnum
Definition: EnumDefinitions.h:106
AutodiffTapeAllocEnum
@ AutodiffTapeAllocEnum
Definition: EnumDefinitions.h:55
IssmPDouble
IssmDouble IssmPDouble
Definition: types.h:38
Parameters::FindParamObject
Param * FindParamObject(int enum_type)
Definition: Parameters.cpp:588
DataSet
Declaration of DataSet class.
Definition: DataSet.h:14
AutodiffGcTriggerMaxSizeEnum
@ AutodiffGcTriggerMaxSizeEnum
Definition: EnumDefinitions.h:48
VerboseControl
bool VerboseControl(void)
Definition: Verbosity.cpp:27
AutodiffTbufsizeEnum
@ AutodiffTbufsizeEnum
Definition: EnumDefinitions.h:56
ISSM_MPI_Gather
int ISSM_MPI_Gather(void *sendbuf, int sendcnt, ISSM_MPI_Datatype sendtype, void *recvbuf, int recvcnt, ISSM_MPI_Datatype recvtype, int root, ISSM_MPI_Comm comm)
Definition: issmmpi.cpp:242
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