Ice Sheet System Model  4.18
Code documentation
•All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
controlvalidation_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 
11 #ifdef _HAVE_CODIPACK_
12 extern CoDi_global codi_global;
13 #include <sstream> // for output of the CoDiPack tape
14 #endif
15 
16 #ifdef _HAVE_AD_
17 void simul_starttrace2(FemModel* femmodel){/*{{{*/
18 
19  #if defined(_HAVE_ADOLC_)
20  /*Retrive ADOLC parameters*/
21  IssmDouble gcTriggerRatio;
22  IssmDouble gcTriggerMaxSize;
23  IssmDouble obufsize;
24  IssmDouble lbufsize;
25  IssmDouble cbufsize;
26  IssmDouble tbufsize;
27  femmodel->parameters->FindParam(&gcTriggerRatio,AutodiffGcTriggerRatioEnum);
28  femmodel->parameters->FindParam(&gcTriggerMaxSize,AutodiffGcTriggerMaxSizeEnum);
29  femmodel->parameters->FindParam(&obufsize,AutodiffObufsizeEnum);
30  femmodel->parameters->FindParam(&lbufsize,AutodiffLbufsizeEnum);
31  femmodel->parameters->FindParam(&cbufsize,AutodiffCbufsizeEnum);
32  femmodel->parameters->FindParam(&tbufsize,AutodiffTbufsizeEnum);
33 
34  /*Set garbage collection parameters: */
35  setStoreManagerControl(reCast<IssmPDouble>(gcTriggerRatio),reCast<size_t>(gcTriggerMaxSize));
36 
37  /*Start trace: */
38  int skipFileDeletion=1;
39  int keepTaylors=1;
40  int my_rank=IssmComm::GetRank();
41  trace_on(my_rank,keepTaylors,reCast<size_t>(obufsize),reCast<size_t>(lbufsize),reCast<size_t>(cbufsize),reCast<size_t>(tbufsize),skipFileDeletion);
42 
43  #elif defined(_HAVE_CODIPACK_)
44 
45  //fprintf(stderr, "*** Codipack IoModel::StartTrace\n");
46  /*
47  * FIXME codi
48  * - ADOL-C variant uses fine grained tracing with various arguments
49  * - ADOL-C variant sets a garbage collection parameter for its tape
50  * -> These parameters are not read for the CoDiPack ISSM version!
51  */
52  auto& tape_codi = IssmDouble::getGlobalTape();
53  tape_codi.setActive();
54  #if _AD_TAPE_ALLOC_
55  //alloc_profiler.Tag(StartInit, true);
56  IssmDouble x_t(1.0), y_t(1.0);
57  tape_codi.registerInput(y_t);
58  int codi_allocn = 0;
59  femmodel->parameters->FindParam(&codi_allocn,AutodiffTapeAllocEnum);
60  for(int i = 0;i < codi_allocn;++i) {
61  x_t = y_t * y_t;
62  }
63  /*
64  std::stringstream out_s;
65  IssmDouble::getGlobalTape().printStatistics(out_s);
66  _printf0_("StartTrace::Tape Statistics : TapeAlloc count=[" << codi_allocn << "]\n" << out_s.str());
67  */
68  tape_codi.reset();
69  //alloc_profiler.Tag(FinishInit, true);
70  #else
71  tape_codi.reset();
72  #endif
73 
74  #else
75  _error_("not implemented");
76  #endif
77 }/*}}}*/
78 void simul_stoptrace2(){/*{{{*/
79 
80  #if defined(_HAVE_ADOLC_)
81  trace_off();
82  if(VerboseAutodiff()){ /*{{{*/
83 
84  #ifdef _HAVE_ADOLC_
85  int my_rank=IssmComm::GetRank();
86  size_t tape_stats[15];
87  tapestats(my_rank,tape_stats); //reading of tape statistics
88  int commSize=IssmComm::GetSize();
89  int *sstats=new int[7];
90  sstats[0]=tape_stats[NUM_OPERATIONS];
91  sstats[1]=tape_stats[OP_FILE_ACCESS];
92  sstats[2]=tape_stats[NUM_LOCATIONS];
93  sstats[3]=tape_stats[LOC_FILE_ACCESS];
94  sstats[4]=tape_stats[NUM_VALUES];
95  sstats[5]=tape_stats[VAL_FILE_ACCESS];
96  sstats[6]=tape_stats[TAY_STACK_SIZE];
97  int *rstats=NULL;
98  if (my_rank==0) rstats=new int[commSize*7];
99  ISSM_MPI_Gather(sstats,7,ISSM_MPI_INT,rstats,7,ISSM_MPI_INT,0,IssmComm::GetComm());
100  if (my_rank==0) {
101  int offset=50;
102  int rOffset=(commSize/10)+1;
103  _printf_(" ADOLC statistics: \n");
104  _printf_(" "<<setw(offset)<<left<<"#independents: " <<setw(12)<<right<<tape_stats[NUM_INDEPENDENTS] << "\n");
105  _printf_(" "<<setw(offset)<<left<<"#dependents: " <<setw(12)<<right<<tape_stats[NUM_DEPENDENTS] << "\n");
106  _printf_(" "<<setw(offset)<<left<<"max #live active variables: " <<setw(12)<<right<<tape_stats[NUM_MAX_LIVES] << "\n");
107  _printf_(" operations: entry size "<< sizeof(unsigned char) << " Bytes \n");
108  _printf_(" "<<setw(offset)<<left<<" #entries in buffer (AutodiffObufsizeEnum) " <<setw(12)<<right<<tape_stats[OP_BUFFER_SIZE] << "\n");
109  for (int r=0;r<commSize;++r)
110  _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");
111  _printf_(" locations: entry size " << sizeof(locint) << " Bytes\n");
112  _printf_(" "<<setw(offset)<<left<<" #entries in buffer (AutodiffLbufsizeEnum) " <<setw(12)<<right<<tape_stats[LOC_BUFFER_SIZE] << "\n");
113  for (int r=0;r<commSize;++r)
114  _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");
115  _printf_(" constant values: entry size " << sizeof(double) << " Bytes\n");
116  _printf_(" "<<setw(offset)<<left<<" #entries in buffer (AutodiffCbufsizeEnum) " <<setw(12)<<right<<tape_stats[VAL_BUFFER_SIZE] << "\n");
117  for (int r=0;r<commSize;++r)
118  _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");
119  _printf_(" Taylor stack: entry size " << sizeof(revreal) << " Bytes\n");
120  _printf_(" "<<setw(offset)<<left<<" #entries in buffer (AutodiffTbufsizeEnum) " <<setw(12)<<right<<tape_stats[TAY_BUFFER_SIZE] << "\n");
121  for (int r=0;r<commSize;++r)
122  _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");
123  delete []rstats;
124  }
125  delete [] sstats;
126  #endif
127 
128  #ifdef _HAVE_CODIPACK_
129  #ifdef _AD_TAPE_ALLOC_
130  //_printf_("Allocation time P(" << my_rank << "): " << alloc_profiler.DeltaTime(StartInit, FinishInit) << "\n");
131  #endif
132  std::stringstream out_s;
133  IssmDouble::getGlobalTape().printStatistics(out_s);
134  _printf0_("CoDiPack Profiling::Tape Statistics :\n" << out_s.str());
135  #endif
136  } /*}}}*/
137 
138  #elif defined(_HAVE_CODIPACK_)
139  auto& tape_codi = IssmDouble::getGlobalTape();
140  tape_codi.setPassive();
141  if(VerboseAutodiff()){
142  int my_rank=IssmComm::GetRank();
143  if(my_rank == 0) {
144  // FIXME codi "just because" for now
145  tape_codi.printStatistics(std::cout);
146  codi_global.print(std::cout);
147  }
148  }
149  #else
150  _error_("not implemented");
151  #endif
152 }/*}}}*/
153 #endif
154 
155 void controlvalidation_core(FemModel* femmodel){
156 
157  int solution_type,n;
158  int num_responses;
159  IssmDouble j0,j;
160  IssmDouble Ialpha,exponent,alpha;
161  IssmDouble* scaling_factors = NULL;
162  IssmDouble* jlist = NULL;
163  int my_rank=IssmComm::GetRank();
164 
165  /*Recover parameters used throughout the solution*/
166  femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
167  femmodel->parameters->SetParam(false,SaveResultsEnum);
168  femmodel->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum);
169  femmodel->parameters->FindParam(&scaling_factors,NULL,InversionControlScalingFactorsEnum);
170 
171  /*Get initial guess*/
172  IssmPDouble* X0 = NULL;
173  GetPassiveVectorFromControlInputsx(&X0,&n,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"value");
174 
175  /*Allocate vectors*/
176  IssmDouble* X = xNew<IssmDouble>(n);
177  IssmPDouble* G = xNew<IssmPDouble>(n);
178 
179  /*out of solution_type, figure out solution core and adjoint function pointer*/
180  void (*solutioncore)(FemModel*)=NULL;
181  CorePointerFromSolutionEnum(&solutioncore,femmodel->parameters,solution_type);
182 
183  #if defined(_HAVE_ADOLC_)
184  /*{{{*/
185  IssmDouble* aX=xNew<IssmDouble>(n);
186  if(my_rank==0){
187  for(int i=0;i<n;i++){
188  aX[i]<<=X0[i];
189  }
190  }
191  _error_("not implemented yet...");
192  /*}}}*/
193  #elif defined(_HAVE_CODIPACK_)
194  simul_starttrace2(femmodel);
195  IssmDouble* aX=xNew<IssmDouble>(n);
196  auto& tape_codi = IssmDouble::getGlobalTape();
197  codi_global.input_indices.clear();
198  if(my_rank==0){
199  for (int i=0;i<n;i++) {
200  aX[i]=X0[i];
201  tape_codi.registerInput(aX[i]);
202  codi_global.input_indices.push_back(aX[i].getGradientData());
203  }
204  }
205  SetControlInputsFromVectorx(femmodel,aX);
206  xDelete(aX);
207 
208  if(VerboseControl()) _printf0_(" Compute Initial cost function\n");
209  solutioncore(femmodel);
210 
211  /*Get Dependents*/
212  IssmDouble output_value;
213  int num_dependents;
214  IssmPDouble *dependents;
215  DataSet* dependent_objects=NULL;
216  IssmDouble J=0.;
217  femmodel->parameters->FindParam(&num_dependents,AutodiffNumDependentsEnum);
218  femmodel->parameters->FindParam(&dependent_objects,AutodiffDependentObjectsEnum);
219 
220  /*Go through our dependent variables, and compute the response:*/
221  dependents=xNew<IssmPDouble>(num_dependents);
222  codi_global.output_indices.clear();
223  for(int i=0;i<dependent_objects->Size();i++){
224  DependentObject* dep=(DependentObject*)dependent_objects->GetObjectByOffset(i);
225  if(solution_type==TransientSolutionEnum){
226  output_value = dep->GetValue();
227  }
228  else{
229  dep->Responsex(&output_value,femmodel);
230  }
231  _printf0_("=== output ="<<output_value<<" \n");
232  if(my_rank==0) {
233  tape_codi.registerOutput(output_value);
234  dependents[i] = output_value.getValue();
235  codi_global.output_indices.push_back(output_value.getGradientData());
236  J+=output_value;
237  }
238  }
239  j0 = J;
240  _printf0_("Initial cost function J(x) = "<<setw(12)<<setprecision(7)<<j0<<"\n");
241  _assert_(j0>0.);
242  simul_stoptrace2();
243  /*initialize direction index in the weights vector: */
244  if(my_rank==0){
245  tape_codi.setGradient(codi_global.output_indices[0],1.0);
246  }
247  tape_codi.evaluate();
248 
249  /*Get gradient for this dependent */
250  auto in_size = codi_global.input_indices.size();
251  for(size_t i = 0; i < in_size; ++i) {
252  G[i] = tape_codi.getGradient(codi_global.input_indices[i]);
253  }
254  #else
255  /*{{{*/
256  void (*adjointcore)(FemModel*) = NULL;
257  AdjointCorePointerFromSolutionEnum(&adjointcore,solution_type);
258 
259  if(VerboseControl()) _printf0_(" Compute Initial solution\n");
260  solutioncore(femmodel);
261  if(VerboseControl()) _printf0_(" Compute Adjoint\n");
262  adjointcore(femmodel);
263 
264  if(VerboseControl()) _printf0_(" Compute Initial cost function\n");
265  femmodel->CostFunctionx(&j0,&jlist,NULL);
266  _printf0_("Initial cost function J(x) = "<<setw(12)<<setprecision(7)<<j0<<"\n");
267  xDelete<IssmDouble>(jlist);
268 
269  if(VerboseControl()) _printf0_(" Compute Gradient\n");
270  Gradjx(&G,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
271  for(int i=0;i<n;i++) G[i] = -G[i];
272  /*}}}*/
273  #endif
274 
275  /*Allocate output*/
276  int num = 26;
277  IssmDouble* output = xNew<IssmDouble>(2*num);
278 
279  /*Start loop*/
280  _printf0_(" alpha Ialpha \n");
281  _printf0_("_________________________\n");
282  for(int m=0;m<num;m++){
283 
284  /*Create new vector*/
285  alpha = pow(2.,-m);
286  for(int i=0;i<n;i++) X[i] = X0[i] + alpha*scaling_factors[0];
287 
288  /*Calculate j(k+alpha delta k) */
289  SetControlInputsFromVectorx(femmodel,X);
290  solutioncore(femmodel);
291 
292  #if defined(_HAVE_CODIPACK_)
293  j=0.;
294  for(int i=0;i<dependent_objects->Size();i++){
295  DependentObject* dep=(DependentObject*)dependent_objects->GetObjectByOffset(i);
296  if(solution_type==TransientSolutionEnum){
297  output_value = dep->GetValue();
298  }
299  else{
300  dep->Responsex(&output_value,femmodel);
301  }
302  j+=output_value;
303  }
304  #else
305  femmodel->CostFunctionx(&j,NULL,NULL);
306  #endif
307 
308  IssmDouble Den = 0.;
309  for(int i=0;i<n;i++) Den += alpha* G[i] * scaling_factors[0];
310  Ialpha = fabs((j - j0)/Den - 1.);
311  _assert_(fabs(Den)>0.);
312 
313  _printf0_(" " << setw(11) << setprecision (5)<<alpha<<" " << setw(11) << setprecision (5)<<Ialpha<<"\n");
314  output[m*2+0] = alpha;
315  output[m*2+1] = Ialpha;
316  }
317 
318  /*output*/
319  #ifdef _HAVE_AD_
320  IssmPDouble* J_passive=xNew<IssmPDouble>(2*num);
321  for(int i=0;i<2*num;i++) J_passive[i]=reCast<IssmPDouble>(output[i]);
322  femmodel->results->AddObject(new GenericExternalResult<IssmPDouble*>(femmodel->results->Size()+1,JEnum,J_passive,num,2,0,0));
323  xDelete<IssmPDouble>(J_passive);
324  IssmDouble* aG=xNew<IssmDouble>(n);
325  for(int i=0;i<n;i++) aG[i] = G[i];
326  ControlInputSetGradientx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,aG);
327  xDelete<IssmDouble>(aG);
328  #else
329  femmodel->results->AddObject(new GenericExternalResult<IssmPDouble*>(femmodel->results->Size()+1,JEnum,output,num,2,0,0));
330  ControlInputSetGradientx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,G);
331  #endif
332  femmodel->OutputControlsx(&femmodel->results);
333 
334  /*Clean up and return*/
335  xDelete<IssmDouble>(output);
336  xDelete<IssmPDouble>(G);
337  xDelete<IssmDouble>(X);
338  xDelete<double>(X0);
339  xDelete<IssmDouble>(scaling_factors);
340 }
DataSet::Size
int Size()
Definition: DataSet.cpp:399
controlvalidation_core
void controlvalidation_core(FemModel *femmodel)
Definition: controlvalidation_core.cpp:155
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
_printf0_
#define _printf0_(StreamArgs)
Definition: Print.h:29
DataSet::AddObject
int AddObject(Object *object)
Definition: DataSet.cpp:252
_printf_
#define _printf_(StreamArgs)
Definition: Print.h:22
AdjointCorePointerFromSolutionEnum
void AdjointCorePointerFromSolutionEnum(void(**padjointcore)(FemModel *), int solutiontype)
Definition: AdjointCorePointerFromSolutionEnum.cpp:18
cores.h
AutodiffObufsizeEnum
@ AutodiffObufsizeEnum
Definition: EnumDefinitions.h:54
TransientSolutionEnum
@ TransientSolutionEnum
Definition: EnumDefinitions.h:1317
FemModel::parameters
Parameters * parameters
Definition: FemModel.h:46
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
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
FemModel::OutputControlsx
void OutputControlsx(Results **presults)
Definition: FemModel.cpp:2171
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
alpha
IssmDouble alpha(IssmDouble x, IssmDouble y, IssmDouble z, int testid)
Definition: fsanalyticals.cpp:221
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
SolutionTypeEnum
@ SolutionTypeEnum
Definition: EnumDefinitions.h:398
AutodiffCbufsizeEnum
@ AutodiffCbufsizeEnum
Definition: EnumDefinitions.h:42
FemModel::loads
Loads * loads
Definition: FemModel.h:54
FemModel::elements
Elements * elements
Definition: FemModel.h:44
DependentObject::Responsex
void Responsex(IssmDouble *poutput_value, FemModel *femmodel)
Definition: DependentObject.cpp:93
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
AutodiffNumDependentsEnum
@ AutodiffNumDependentsEnum
Definition: EnumDefinitions.h:52
Parameters::FindParam
void FindParam(bool *pinteger, int enum_type)
Definition: Parameters.cpp:262
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
AutodiffTapeAllocEnum
@ AutodiffTapeAllocEnum
Definition: EnumDefinitions.h:55
IssmPDouble
IssmDouble IssmPDouble
Definition: types.h:38
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