Ice Sheet System Model  4.18
Code documentation
UpdateElementsAndMaterialsControl.cpp
Go to the documentation of this file.
1 /*
2  * UpdateElementsAndMaterialsControl:
3  */
4 
5 #include "../../../toolkits/toolkits.h"
6 #include "../../../classes/classes.h"
7 #include "../../../shared/shared.h"
8 #include "../../MeshPartitionx/MeshPartitionx.h"
9 #include "../ModelProcessorx.h"
10 
11 void UpdateElementsAndMaterialsControl(Elements* elements,Parameters* parameters,Inputs2* inputs2,Materials* materials, IoModel* iomodel){
12  /*Intermediary*/
13  bool control_analysis;
14  int M,N;
15  int control,cost_function,domaintype;
16  int num_controls,num_cost_functions;
17  IssmDouble yts,scale;
18  Element *element = NULL;
19  Material *material = NULL;
20  int *control_enums = NULL;
21  char **controls = NULL;
22  char **cost_functions = NULL;
23  IssmDouble *independent = NULL;
24  IssmDouble *independents_min = NULL;
25  IssmDouble *independents_max = NULL;
26  IssmDouble *weights = NULL;
27 
28  /*Fetch parameters: */
29  iomodel->FindConstant(&control_analysis,"md.inversion.iscontrol");
30  if(!control_analysis) return;
31 
32  /*Fetch parameters: */
33  bool isautodiff;
34  iomodel->FindConstant(&isautodiff,"md.autodiff.isautodiff");
35  if(isautodiff){
36  UpdateElementsAndMaterialsControlAD(elements,parameters,inputs2,materials,iomodel);
37  return;
38  }
39 
40  /*Process controls and convert from string to enums*/
41  iomodel->FindConstant(&num_controls,"md.inversion.num_control_parameters");
42  iomodel->FindConstant(&controls,&num_controls,"md.inversion.control_parameters");
43  if(num_controls<1) _error_("no controls found");
44  control_enums=xNew<int>(num_controls);
45  for(int i=0;i<num_controls;i++){
46  control_enums[i]=StringToEnumx(controls[i]);
47  }
48 
49  /*Process cost functions and convert from string to enums*/
50  iomodel->FindConstant(&num_cost_functions,"md.inversion.num_cost_functions");
51  iomodel->FindConstant(&cost_functions,&num_cost_functions,"md.inversion.cost_functions");
52  if(num_cost_functions<1) _error_("No cost functions found");
53  int* cost_function_enums=xNew<int>(num_cost_functions);
54  for(int i=0;i<num_cost_functions;++i){
55  cost_function_enums[i]=StringToEnumx(cost_functions[i]);
56  }
57 
58  /*Fetch Observations and add to inputs*/
59  iomodel->FindConstant(&domaintype,"md.mesh.domain_type");
60  iomodel->FindConstant(&yts,"md.constants.yts");
61  iomodel->FetchData(&weights,&M,&N,"md.inversion.cost_functions_coefficients");
62 
63  /*Transpose weights for simplicity!*/
64  if(M*N && N>1){
65  IssmDouble* weights_transp = xNew<IssmDouble>(M*N);
66  for(int i=0;i<M;i++) for(int j=0;j<N;j++) weights_transp[j*M+i] = weights[i*N+j];
67  xDelete<IssmDouble>(weights);
68  weights = weights_transp;
69  }
70 
71  if(M!=iomodel->numberofvertices && N!=num_cost_functions) _error_("not supported");
72  for(int i=0;i<num_cost_functions;i++){
73  cost_function=cost_function_enums[i];
74  if( cost_function==ThicknessAbsMisfitEnum) iomodel->FetchDataToInput(inputs2,elements,"md.inversion.thickness_obs",InversionThicknessObsEnum);
75  else if(cost_function==SurfaceAbsMisfitEnum) iomodel->FetchDataToInput(inputs2,elements,"md.inversion.surface_obs",InversionSurfaceObsEnum);
76  else if(cost_function==RheologyBInitialguessMisfitEnum) iomodel->FetchDataToInput(inputs2,elements,"md.materials.rheology_B",RheologyBInitialguessEnum);
77  else if(cost_function==SurfaceAbsVelMisfitEnum
78  || cost_function==SurfaceRelVelMisfitEnum
79  || cost_function==SurfaceLogVelMisfitEnum
80  || cost_function==SurfaceLogVxVyMisfitEnum
81  || cost_function==SurfaceAverageVelMisfitEnum){
82  iomodel->FetchDataToInput(inputs2,elements,"md.inversion.vx_obs",InversionVxObsEnum);
83  if(domaintype!=Domain2DverticalEnum) iomodel->FetchDataToInput(inputs2,elements,"md.inversion.vy_obs",InversionVyObsEnum);
84  }
85  for(int j=0;j<elements->Size();j++){
86  Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(j));
87  element->DatasetInputAdd(InversionCostFunctionsCoefficientsEnum,&weights[i*iomodel->numberofvertices],inputs2,iomodel,M,1,1,cost_function,7,cost_function);
88  }
89  }
90  parameters->AddObject(new IntParam(ControlInputSizeMEnum,iomodel->numberofvertices));
91  xDelete<IssmDouble>(weights);
92 
93  /*Get controls*/
94  iomodel->FetchData(&independents_min,&M,&N,"md.inversion.min_parameters");
95  if(M!=iomodel->numberofvertices && N!=num_controls) _error_("not supported");
96  iomodel->FetchData(&independents_max,&M,&N,"md.inversion.max_parameters");
97  if(M!=iomodel->numberofvertices && N!=num_controls) _error_("not supported");
98 
99  /*Transpose weights for simplicity!*/
100  if(M*N && N>1){
101  IssmDouble* independents_min_transp = xNew<IssmDouble>(M*N);
102  for(int i=0;i<M;i++) for(int j=0;j<N;j++) independents_min_transp[j*M+i] = independents_min[i*N+j];
103  xDelete<IssmDouble>(independents_min);
104  independents_min = independents_min_transp;
105 
106  IssmDouble* independents_max_transp = xNew<IssmDouble>(M*N);
107  for(int i=0;i<M;i++) for(int j=0;j<N;j++) independents_max_transp[j*M+i] = independents_max[i*N+j];
108  xDelete<IssmDouble>(independents_max);
109  independents_max = independents_max_transp;
110  }
111 
112  for(int i=0;i<num_controls;i++){
113  control = control_enums[i];
114  scale = 1.;
115 
116  switch(control){
117  /*List of supported controls*/
118  case BalancethicknessThickeningRateEnum: iomodel->FetchData(&independent,&M,&N,"md.balancethickness.thickening_rate");scale = 1./yts; break;
119  case BalancethicknessSpcthicknessEnum: iomodel->FetchData(&independent,&M,&N,"md.balancethickness.spcthickness"); break;
120  case VxEnum: iomodel->FetchData(&independent,&M,&N,"md.initialization.vx");scale = 1./yts; break;
121  case VyEnum: iomodel->FetchData(&independent,&M,&N,"md.initialization.vy");scale = 1./yts; break;
122  case ThicknessEnum: iomodel->FetchData(&independent,&M,&N,"md.geometry.thickness"); break;
123  case FrictionCoefficientEnum: iomodel->FetchData(&independent,&M,&N,"md.friction.coefficient"); break;
124  case FrictionAsEnum: iomodel->FetchData(&independent,&M,&N,"md.friction.As"); break;
125  case BalancethicknessApparentMassbalanceEnum: iomodel->FetchData(&independent,&M,&N,"md.balancethickness.apparent_massbalance"); break;
126  case BalancethicknessOmegaEnum: iomodel->FetchData(&independent,&M,&N,"md.balancethickness.omega"); break;
127  case MaterialsRheologyBEnum: iomodel->FetchData(&independent,&M,&N,"md.materials.rheology_B"); break;
128  /*Special cases*/
129  case MaterialsRheologyBbarEnum: iomodel->FetchData(&independent,&M,&N,"md.materials.rheology_B"); break;
130  case DamageDbarEnum: iomodel->FetchData(&independent,&M,&N,"md.damage.D"); break;
131  default:
132  _error_("Control " << EnumToStringx(control) << " not implemented yet");
133  }
134  if(M!=iomodel->numberofvertices && N!=1) _error_("not supported yet");
135 
136  /*Special case if 3d*/
137  if(iomodel->domaintype==Domain3DEnum){
138  if(control==MaterialsRheologyBbarEnum) control=MaterialsRheologyBEnum;
139  if(control==DamageDbarEnum) control=DamageDEnum;
140  }
141 
142  for(int j=0;j<elements->Size();j++){
143  Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(j));
144  element->ControlInputCreate(independent,&independents_min[i*iomodel->numberofvertices],&independents_max[i*iomodel->numberofvertices],inputs2,iomodel,M,N,scale,control,i+1);
145  }
146  xDelete<IssmDouble>(independent);
147  }
148  xDelete<IssmDouble>(independents_min);
149  xDelete<IssmDouble>(independents_max);
150 
151 
152  /*Free data: */
153  for(int i=0;i<num_controls;i++){
154  switch(control_enums[i]){
155  /*List of supported controls*/
156  case BalancethicknessThickeningRateEnum: iomodel->DeleteData(1,"md.balancethickness.thickening_rate"); break;
157  case BalancethicknessSpcthicknessEnum: iomodel->DeleteData(1,"md.balancethickness.spcthickness"); break;
158  case VxEnum: iomodel->DeleteData(1,"md.initialization.vx"); break;
159  case VyEnum: iomodel->DeleteData(1,"md.initialization.vy"); break;
160  case ThicknessEnum: iomodel->DeleteData(1,"md.geometry.thickness"); break;
161  case FrictionCoefficientEnum: iomodel->DeleteData(1,"md.friction.coefficient"); break;
162  case FrictionAsEnum: iomodel->DeleteData(1,"md.friction.As"); break;
163  case BalancethicknessApparentMassbalanceEnum: iomodel->DeleteData(1,"md.balancethickness.apparent_massbalance"); break;
164  case BalancethicknessOmegaEnum: iomodel->DeleteData(1,"md.balancethickness.omega"); break;
165  case MaterialsRheologyBEnum: iomodel->DeleteData(1,"md.materials.rheology_B"); break;
166  /*Special cases*/
167  case MaterialsRheologyBbarEnum: iomodel->DeleteData(1,"md.materials.rheology_B"); break;
168  case DamageDbarEnum: iomodel->DeleteData(1,"md.damage.D"); break;
169  default:
170  _error_("Control " << EnumToStringx(control_enums[i]) << " not implemented yet");
171  }
172  }
173 
174  xDelete<int>(control_enums);
175  xDelete<int>(cost_function_enums);
176  for(int i=0;i<num_cost_functions;i++) xDelete<char>(cost_functions[i]);
177  xDelete<char*>(cost_functions);
178  for(int i=0;i<num_controls;i++) xDelete<char>(controls[i]);
179  xDelete<char*>(controls);
180 }
181 void UpdateElementsAndMaterialsControlAD(Elements* elements,Parameters* parameters,Inputs2* inputs2,Materials* materials, IoModel* iomodel){
182 
183  #if defined(_HAVE_AD_)
184  /*Intermediaries*/
185  int num_independent_objects,M,N,M_par,N_par;
186  char** names = NULL;
187  int* types = NULL;
188  int* control_sizes = NULL;
189  int* M_all = NULL;
190  IssmDouble* independent = NULL;
191  IssmDouble* independents_fullmin = NULL;
192  IssmDouble* independents_fullmax = NULL;
193  bool control_analysis =false;
194 
195  iomodel->FindConstant(&control_analysis,"md.inversion.iscontrol");
196 
197  /*Now, return if no control*/
198  if(!control_analysis) return;
199 
200  /*Step1: create controls (independents)*/
201  iomodel->FetchData(&num_independent_objects,"md.autodiff.num_independent_objects");
202  _assert_(num_independent_objects>0);
203  iomodel->FetchData(&names,&M,"md.autodiff.independent_object_names");
204  _assert_(M==num_independent_objects);
205  iomodel->FetchData(&types,NULL,NULL,"md.autodiff.independent_object_types");
206 
207  M_all = xNew<int>(num_independent_objects);
208 
209  /*create independent objects, and at the same time, fetch the corresponding independent variables,
210  *and declare them as such in ADOLC: */
211  iomodel->FetchData(&independents_fullmin,&M_par,&N_par,"md.autodiff.independent_min_parameters");
212  iomodel->FetchData(&independents_fullmax,&M_par,&N_par,"md.autodiff.independent_max_parameters");
213  iomodel->FetchData(&control_sizes,NULL,NULL,"md.autodiff.independent_control_sizes");
214 
215  int* start_point = NULL;
216  start_point = xNew<int>(num_independent_objects);
217  int counter = 0;
218  for(int i=0;i<num_independent_objects;i++){
219  start_point[i]=counter;
220  counter+=control_sizes[i];
221  }
222 
223  for(int i=0;i<num_independent_objects;i++){
224 
225  if(types[i]==1){ /* vector:*/
226 
227  /*Get field name and input Enum from independent name*/
228  char* iofieldname = NULL;
229  int input_enum;
230  IssmDouble* independents_min = NULL;
231  IssmDouble* independents_max = NULL;
232 
233  FieldAndEnumFromCode(&input_enum,&iofieldname,names[i]);
234 
235  /*Fetch required data*/
236  iomodel->FetchData(&independent,&M,&N,iofieldname);
237  _assert_(independent);
238  _assert_(N==control_sizes[i]);
239 
240  independents_min = NULL; independents_min = xNew<IssmDouble>(M*N);
241  independents_max = NULL; independents_max = xNew<IssmDouble>(M*N);
242  for(int m=0;m<M;m++){
243  for(int n=0;n<N;n++){
244  independents_min[N*m+n]=independents_fullmin[N_par*m+start_point[i]+n];
245  independents_max[N*m+n]=independents_fullmax[N_par*m+start_point[i]+n];
246  }
247  }
248  if(N!=1) M_all[i]=M-1;
249  else M_all[i]=M;
250 
251  for(int j=0;j<elements->Size();j++){
252  Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(j));
253  element->ControlInputCreate(independent,independents_min,independents_max,inputs2,iomodel,M,N,1.,input_enum,i+1);
254  }
255  xDelete<IssmDouble>(independent);
256  xDelete<IssmDouble>(independents_min);
257  xDelete<IssmDouble>(independents_max);
258 
259  }
260  else{
261  _error_("Independent cannot be of size " << types[i]);
262  }
263  }
264  parameters->AddObject(new IntVecParam(ControlInputSizeNEnum,control_sizes,num_independent_objects));
265  parameters->AddObject(new IntVecParam(ControlInputSizeMEnum,M_all,num_independent_objects));
266 
267  /*cleanup*/
268  for(int i=0;i<num_independent_objects;i++){
269  xDelete<char>(names[i]);
270  }
271  xDelete<char*>(names);
272  xDelete<int>(types);
273  xDelete<int>(M_all);
274  xDelete<IssmDouble>(independents_fullmin);
275  xDelete<IssmDouble>(independents_fullmax);
276  xDelete<int>(start_point);
277  xDelete<int>(control_sizes);
278  /*Step2: create cost functions (dependents)*/
279 
280  return;
281 #else
282  _error_("AD not compiled");
283 #endif
284 }
DataSet::Size
int Size()
Definition: DataSet.cpp:399
InversionThicknessObsEnum
@ InversionThicknessObsEnum
Definition: EnumDefinitions.h:631
FrictionCoefficientEnum
@ FrictionCoefficientEnum
Definition: EnumDefinitions.h:574
_assert_
#define _assert_(ignore)
Definition: exceptions.h:37
IssmDouble
double IssmDouble
Definition: types.h:37
InversionVyObsEnum
@ InversionVyObsEnum
Definition: EnumDefinitions.h:634
BalancethicknessApparentMassbalanceEnum
@ BalancethicknessApparentMassbalanceEnum
Definition: EnumDefinitions.h:982
BalancethicknessThickeningRateEnum
@ BalancethicknessThickeningRateEnum
Definition: EnumDefinitions.h:474
InversionCostFunctionsCoefficientsEnum
@ InversionCostFunctionsCoefficientsEnum
Definition: EnumDefinitions.h:629
Parameters
Declaration of Parameters class.
Definition: Parameters.h:18
Elements
Declaration of Elements class.
Definition: Elements.h:17
SurfaceRelVelMisfitEnum
@ SurfaceRelVelMisfitEnum
Definition: EnumDefinitions.h:828
Material
Definition: Material.h:21
ControlInputSizeMEnum
@ ControlInputSizeMEnum
Definition: EnumDefinitions.h:105
UpdateElementsAndMaterialsControl
void UpdateElementsAndMaterialsControl(Elements *elements, Parameters *parameters, Inputs2 *inputs2, Materials *materials, IoModel *iomodel)
Definition: UpdateElementsAndMaterialsControl.cpp:11
IntVecParam
Definition: IntVecParam.h:20
VyEnum
@ VyEnum
Definition: EnumDefinitions.h:850
Parameters::AddObject
void AddObject(Param *newparam)
Definition: Parameters.cpp:67
UpdateElementsAndMaterialsControlAD
void UpdateElementsAndMaterialsControlAD(Elements *elements, Parameters *parameters, Inputs2 *inputs2, Materials *materials, IoModel *iomodel)
Definition: UpdateElementsAndMaterialsControl.cpp:181
FrictionAsEnum
@ FrictionAsEnum
Definition: EnumDefinitions.h:571
SurfaceLogVxVyMisfitEnum
@ SurfaceLogVxVyMisfitEnum
Definition: EnumDefinitions.h:826
FieldAndEnumFromCode
void FieldAndEnumFromCode(int *out_enum, char **pfield, const char *string_in)
Definition: IoCodeConversions.cpp:9
MaterialsRheologyBbarEnum
@ MaterialsRheologyBbarEnum
Definition: EnumDefinitions.h:644
SurfaceAbsVelMisfitEnum
@ SurfaceAbsVelMisfitEnum
Definition: EnumDefinitions.h:818
Element
Definition: Element.h:41
IoModel::numberofvertices
int numberofvertices
Definition: IoModel.h:99
RheologyBInitialguessEnum
@ RheologyBInitialguessEnum
Definition: EnumDefinitions.h:672
IoModel::DeleteData
void DeleteData(int num,...)
Definition: IoModel.cpp:500
BalancethicknessOmegaEnum
@ BalancethicknessOmegaEnum
Definition: EnumDefinitions.h:473
SurfaceAbsMisfitEnum
@ SurfaceAbsMisfitEnum
Definition: EnumDefinitions.h:817
Materials
Declaration of Materials class.
Definition: Materials.h:16
EnumToStringx
const char * EnumToStringx(int enum_in)
Definition: EnumToStringx.cpp:15
IoModel::FindConstant
void FindConstant(bool *pvalue, const char *constant_name)
Definition: IoModel.cpp:2362
DamageDEnum
@ DamageDEnum
Definition: EnumDefinitions.h:516
IntParam
Definition: IntParam.h:20
ThicknessAbsMisfitEnum
@ ThicknessAbsMisfitEnum
Definition: EnumDefinitions.h:837
InversionSurfaceObsEnum
@ InversionSurfaceObsEnum
Definition: EnumDefinitions.h:630
IoModel::FetchData
void FetchData(bool *pboolean, const char *data_name)
Definition: IoModel.cpp:933
Inputs2
Declaration of Inputs class.
Definition: Inputs2.h:23
Domain3DEnum
@ Domain3DEnum
Definition: EnumDefinitions.h:536
DamageDbarEnum
@ DamageDbarEnum
Definition: EnumDefinitions.h:518
StringToEnumx
int StringToEnumx(const char *string_in, bool notfounderror=true)
Definition: StringToEnumx.cpp:14
_error_
#define _error_(StreamArgs)
Definition: exceptions.h:49
MaterialsRheologyBEnum
@ MaterialsRheologyBEnum
Definition: EnumDefinitions.h:643
DataSet::GetObjectByOffset
Object * GetObjectByOffset(int offset)
Definition: DataSet.cpp:334
VxEnum
@ VxEnum
Definition: EnumDefinitions.h:846
SurfaceLogVelMisfitEnum
@ SurfaceLogVelMisfitEnum
Definition: EnumDefinitions.h:825
Element::ControlInputCreate
void ControlInputCreate(IssmDouble *doublearray, IssmDouble *independents_min, IssmDouble *independents_max, Inputs2 *inputs2, IoModel *iomodel, int M, int N, IssmDouble scale, int input_enum, int id)
Definition: Element.cpp:1753
ThicknessEnum
@ ThicknessEnum
Definition: EnumDefinitions.h:840
IoModel::FetchDataToInput
void FetchDataToInput(Inputs2 *inputs2, Elements *elements, const char *vector_name, int input_enum)
Definition: IoModel.cpp:1651
InversionVxObsEnum
@ InversionVxObsEnum
Definition: EnumDefinitions.h:633
SurfaceAverageVelMisfitEnum
@ SurfaceAverageVelMisfitEnum
Definition: EnumDefinitions.h:821
IoModel
Definition: IoModel.h:48
Element::DatasetInputAdd
void DatasetInputAdd(int enum_type, IssmDouble *vector, Inputs2 *inputs2, IoModel *iomodel, int M, int N, int vector_type, int vector_enum, int code, int input_enum)
Definition: Element.cpp:1831
ControlInputSizeNEnum
@ ControlInputSizeNEnum
Definition: EnumDefinitions.h:106
IoModel::domaintype
int domaintype
Definition: IoModel.h:78
Domain2DverticalEnum
@ Domain2DverticalEnum
Definition: EnumDefinitions.h:535
RheologyBInitialguessMisfitEnum
@ RheologyBInitialguessMisfitEnum
Definition: EnumDefinitions.h:673
BalancethicknessSpcthicknessEnum
@ BalancethicknessSpcthicknessEnum
Definition: EnumDefinitions.h:986
Generated on Thu Jul 2 2020 08:09:19 for Ice Sheet System Model by   doxygen 1.8.19