Ice Sheet System Model  4.18
Code documentation
BalancethicknessAnalysis.cpp
Go to the documentation of this file.
1 #include "./BalancethicknessAnalysis.h"
2 #include "../toolkits/toolkits.h"
3 #include "../classes/classes.h"
4 #include "../shared/shared.h"
5 #include "../modules/modules.h"
6 #include "../classes/Inputs2/DatasetInput2.h"
7 
8 /*Model processing*/
9 void BalancethicknessAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
10 
11  /*Fetch parameters: */
12  int stabilization;
13  iomodel->FindConstant(&stabilization,"md.balancethickness.stabilization");
14 
15  /*Do not add constraints in DG*/
16  if(stabilization!=3){
17  IoModelToConstraintsx(constraints,iomodel,"md.balancethickness.spcthickness",BalancethicknessAnalysisEnum,P1Enum);
18  }
19 
20 }/*}}}*/
21 void BalancethicknessAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
22 
23  /*Intermediary*/
24  int element;
25  int stabilization;
26 
27  /*Fetch parameters: */
28  iomodel->FindConstant(&stabilization,"md.balancethickness.stabilization");
29 
30  /*Loads only in DG*/
31  if (stabilization==3){
32 
33  /*Get faces and elements*/
34  CreateFaces(iomodel);
35  iomodel->FetchData(1,"md.geometry.thickness");
36 
37  /*First load data:*/
38  for(int i=0;i<iomodel->numberoffaces;i++){
39 
40  /*Get left and right elements*/
41  element=iomodel->faces[4*i+2]-1; //faces are [node1 node2 elem1 elem2]
42 
43  /*Now, if this element is not in the partition, pass: */
44  if(!iomodel->my_elements[element]) continue;
45 
46  /* Add load */
47  loads->AddObject(new Numericalflux(i+1,i,i,iomodel));
48  }
49 
50  /*Free data: */
51  iomodel->DeleteData(1,"md.geometry.thickness");
52  }
53 }/*}}}*/
54 void BalancethicknessAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel,bool isamr){/*{{{*/
55 
56  int stabilization;
57  iomodel->FindConstant(&stabilization,"md.balancethickness.stabilization");
58 
59  /*Check in 3d*/
60  if(stabilization==3 && iomodel->domaintype==Domain3DEnum) _error_("DG 3d not implemented yet");
61 
62  /*First fetch data: */
63  if(iomodel->domaintype==Domain3DEnum) iomodel->FetchData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface");
64  if(stabilization!=3){
65  ::CreateNodes(nodes,iomodel,BalancethicknessAnalysisEnum,P1Enum);
66  }
67  else{
68  ::CreateNodes(nodes,iomodel,BalancethicknessAnalysisEnum,P1DGEnum);
69  }
70  iomodel->DeleteData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface");
71 }/*}}}*/
72 int BalancethicknessAnalysis::DofsPerNode(int** doflist,int domaintype,int approximation){/*{{{*/
73  return 1;
74 }/*}}}*/
75 void BalancethicknessAnalysis::UpdateElements(Elements* elements,Inputs2* inputs2,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
76 
77  int stabilization,finiteelement;
78 
79  /*Fetch data needed: */
80  iomodel->FindConstant(&stabilization,"md.balancethickness.stabilization");
81 
82  /*Finite element type*/
83  finiteelement = P1Enum;
84  if(stabilization==3){
85  finiteelement = P1DGEnum;
86  }
87 
88  /*Update elements: */
89  int counter=0;
90  for(int i=0;i<iomodel->numberofelements;i++){
91  if(iomodel->my_elements[i]){
92  Element* element=(Element*)elements->GetObjectByOffset(counter);
93  element->Update(inputs2,i,iomodel,analysis_counter,analysis_type,finiteelement);
94  counter++;
95  }
96  }
97 
98  iomodel->FetchDataToInput(inputs2,elements,"md.geometry.thickness",ThicknessEnum);
99  iomodel->FetchDataToInput(inputs2,elements,"md.geometry.surface",SurfaceEnum);
100  iomodel->FetchDataToInput(inputs2,elements,"md.geometry.base",BaseEnum);
101  iomodel->FetchDataToInput(inputs2,elements,"md.solidearth.sealevel",SealevelEnum,0);
102  iomodel->FetchDataToInput(inputs2,elements,"md.mask.ice_levelset",MaskIceLevelsetEnum);
103  iomodel->FetchDataToInput(inputs2,elements,"md.initialization.vx",VxEnum);
104  iomodel->FetchDataToInput(inputs2,elements,"md.initialization.vy",VyEnum);
105  iomodel->FetchDataToInput(inputs2,elements,"md.basalforcings.groundedice_melting_rate",BasalforcingsGroundediceMeltingRateEnum);
106  iomodel->FetchDataToInput(inputs2,elements,"md.smb.mass_balance",SmbMassBalanceEnum);
107  iomodel->FetchDataToInput(inputs2,elements,"md.balancethickness.thickening_rate",BalancethicknessThickeningRateEnum);
108 
109  if(iomodel->domaintype!=Domain2DhorizontalEnum){
110  iomodel->FetchDataToInput(inputs2,elements,"md.mesh.vertexonbase",MeshVertexonbaseEnum);
111  iomodel->FetchDataToInput(inputs2,elements,"md.mesh.vertexonsurface",MeshVertexonsurfaceEnum);
112  }
113 }/*}}}*/
114 void BalancethicknessAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
115  parameters->AddObject(iomodel->CopyConstantObject("md.balancethickness.stabilization",BalancethicknessStabilizationEnum));
116 }/*}}}*/
117 
118 /*Finite Element Analysis*/
119 void BalancethicknessAnalysis::Core(FemModel* femmodel){/*{{{*/
120  _error_("not implemented");
121 }/*}}}*/
122 ElementVector* BalancethicknessAnalysis::CreateDVector(Element* element){/*{{{*/
123  /*Default, return NULL*/
124  return NULL;
125 }/*}}}*/
126 ElementMatrix* BalancethicknessAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/
127 _error_("Not implemented");
128 }/*}}}*/
129 ElementMatrix* BalancethicknessAnalysis::CreateKMatrix(Element* element){/*{{{*/
130 
131  if(!element->IsOnBase()) return NULL;
132  Element* basalelement = element->SpawnBasalElement();
133 
134  ElementMatrix* Ke = NULL;
135  switch(element->FiniteElement()){
136  case P1Enum: case P2Enum:
137  Ke = CreateKMatrixCG(basalelement);
138  break;
139  case P1DGEnum:
140  Ke = CreateKMatrixDG(basalelement);
141  break;
142  default:
143  _error_("Element type " << EnumToStringx(element->FiniteElement()) << " not supported yet");
144  }
145 
146  int domaintype;
147  element->FindParam(&domaintype,DomainTypeEnum);
148  if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
149  return Ke;
150 }/*}}}*/
151 ElementMatrix* BalancethicknessAnalysis::CreateKMatrixCG(Element* element){/*{{{*/
152 
153  /*Intermediaries */
154  int stabilization;
155  int domaintype;
156  IssmDouble Jdet,D_scalar,h;
157  IssmDouble vel,vx,vy,dvxdx,dvydy;
158  IssmDouble dvx[2],dvy[2];
159  IssmDouble* xyz_list = NULL;
160 
161  /*Fetch number of nodes and dof for this finite element*/
162  int numnodes = element->GetNumberOfNodes();
163 
164  /*Initialize Element vector and other vectors*/
165  ElementMatrix* Ke = element->NewElementMatrix();
166  IssmDouble* basis = xNew<IssmDouble>(numnodes);
167  IssmDouble* dbasis = xNew<IssmDouble>(2*numnodes);
168  IssmDouble D[2][2];
169 
170  /*Retrieve all inputs and parameters*/
171  element->GetVerticesCoordinates(&xyz_list);
172  element->FindParam(&domaintype,DomainTypeEnum);
173  element->FindParam(&stabilization,BalancethicknessStabilizationEnum);
174  Input2* vxaverage_input=NULL;
175  Input2* vyaverage_input=NULL;
176  if(domaintype==Domain2DhorizontalEnum){
177  vxaverage_input=element->GetInput2(VxEnum); _assert_(vxaverage_input);
178  vyaverage_input=element->GetInput2(VyEnum); _assert_(vyaverage_input);
179  }
180  else{
181  vxaverage_input=element->GetInput2(VxAverageEnum); _assert_(vxaverage_input);
182  vyaverage_input=element->GetInput2(VyAverageEnum); _assert_(vyaverage_input);
183  }
184  h = element->CharacteristicLength();
185 
186  /* Start looping on the number of gaussian points: */
187  Gauss* gauss=element->NewGauss(2);
188  for(int ig=gauss->begin();ig<gauss->end();ig++){
189  gauss->GaussPoint(ig);
190 
191  element->JacobianDeterminant(&Jdet,xyz_list,gauss);
192  element->NodalFunctions(basis,gauss);
193  element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
194 
195  vxaverage_input->GetInputValue(&vx,gauss);
196  vyaverage_input->GetInputValue(&vy,gauss);
197  vxaverage_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
198  vyaverage_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
199  dvxdx=dvx[0];
200  dvydy=dvy[1];
201  D_scalar=gauss->weight*Jdet;
202 
203  for(int i=0;i<numnodes;i++){
204  for(int j=0;j<numnodes;j++){
205  /*\phi_i \phi_j \nabla\cdot v*/
206  Ke->values[i*numnodes+j] += D_scalar*basis[i]*basis[j]*(dvxdx+dvydy);
207  /*\phi_i v\cdot\nabla\phi_j*/
208  Ke->values[i*numnodes+j] += D_scalar*basis[i]*(vx*dbasis[0*numnodes+j] + vy*dbasis[1*numnodes+j]);
209  }
210  }
211 
212  if(stabilization==1){
213  /*Streamline upwinding*/
214  vel=sqrt(vx*vx+vy*vy);
215  D[0][0]=h/(2*vel)*vx*vx;
216  D[1][0]=h/(2*vel)*vy*vx;
217  D[0][1]=h/(2*vel)*vx*vy;
218  D[1][1]=h/(2*vel)*vy*vy;
219  }
220  else if(stabilization==2){
221  /*SSA*/
222  vxaverage_input->GetInputAverage(&vx);
223  vyaverage_input->GetInputAverage(&vy);
224  D[0][0]=h/2.0*fabs(vx);
225  D[0][1]=0.;
226  D[1][0]=0.;
227  D[1][1]=h/2.0*fabs(vy);
228  }
229  if(stabilization==1 || stabilization==2){
230  D[0][0]=D_scalar*D[0][0];
231  D[1][0]=D_scalar*D[1][0];
232  D[0][1]=D_scalar*D[0][1];
233  D[1][1]=D_scalar*D[1][1];
234 
235  for(int i=0;i<numnodes;i++){
236  for(int j=0;j<numnodes;j++){
237  Ke->values[i*numnodes+j] += (
238  dbasis[0*numnodes+i] *(D[0][0]*dbasis[0*numnodes+j] + D[0][1]*dbasis[1*numnodes+j]) +
239  dbasis[1*numnodes+i] *(D[1][0]*dbasis[0*numnodes+j] + D[1][1]*dbasis[1*numnodes+j])
240  );
241  }
242  }
243  }
244  }
245 
246  /*Clean up and return*/
247  xDelete<IssmDouble>(xyz_list);
248  xDelete<IssmDouble>(basis);
249  xDelete<IssmDouble>(dbasis);
250  delete gauss;
251  return Ke;
252 }/*}}}*/
253 ElementMatrix* BalancethicknessAnalysis::CreateKMatrixDG(Element* element){/*{{{*/
254 
255  /*Intermediaries */
256  int domaintype;
257  IssmDouble Jdet,D_scalar,vx,vy,dvxdx,dvydy,vel;
258  IssmDouble dvx[2],dvy[2];
259  IssmDouble* xyz_list = NULL;
260 
261  /*Fetch number of nodes and dof for this finite element*/
262  int numnodes = element->GetNumberOfNodes();
263 
264  /*Initialize Element vector and other vectors*/
265  ElementMatrix* Ke = element->NewElementMatrix();
266  IssmDouble* B = xNew<IssmDouble>(2*numnodes);
267  IssmDouble* Bprime = xNew<IssmDouble>(2*numnodes);
268  IssmDouble D[2][2];
269 
270  /*Retrieve all inputs and parameters*/
271  element->GetVerticesCoordinates(&xyz_list);
272  element->FindParam(&domaintype,DomainTypeEnum);
273  Input2* vxaverage_input=NULL;
274  Input2* vyaverage_input=NULL;
275  if(domaintype==Domain2DhorizontalEnum){
276  vxaverage_input=element->GetInput2(VxEnum); _assert_(vxaverage_input);
277  vyaverage_input=element->GetInput2(VyEnum); _assert_(vyaverage_input);
278  }
279  else{
280  vxaverage_input=element->GetInput2(VxAverageEnum); _assert_(vxaverage_input);
281  vyaverage_input=element->GetInput2(VyAverageEnum); _assert_(vyaverage_input);
282  }
283 
284  /* Start looping on the number of gaussian points: */
285  Gauss* gauss=element->NewGauss(2);
286  for(int ig=gauss->begin();ig<gauss->end();ig++){
287  gauss->GaussPoint(ig);
288 
289  element->JacobianDeterminant(&Jdet,xyz_list,gauss);
290  vxaverage_input->GetInputValue(&vx,gauss);
291  vyaverage_input->GetInputValue(&vy,gauss);
292  vxaverage_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
293  vyaverage_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
294  D_scalar=gauss->weight*Jdet;
295 
296  /*WARNING: B and Bprime are inverted compared to CG*/
297  GetB(Bprime,element,xyz_list,gauss);
298  GetBprime(B,element,xyz_list,gauss);
299 
300  D_scalar = - gauss->weight*Jdet;
301  D[0][0] = D_scalar*vx;
302  D[0][1] = 0.;
303  D[1][0] = 0.;
304  D[1][1] = D_scalar*vy;
305  TripleMultiply(B,2,numnodes,1,
306  &D[0][0],2,2,0,
307  Bprime,2,numnodes,0,
308  &Ke->values[0],1);
309 
310  }
311 
312  /*Clean up and return*/
313  xDelete<IssmDouble>(xyz_list);
314  xDelete<IssmDouble>(B);
315  xDelete<IssmDouble>(Bprime);
316  delete gauss;
317  return Ke;
318 }/*}}}*/
319 ElementVector* BalancethicknessAnalysis::CreatePVector(Element* element){/*{{{*/
320 
321  if(!element->IsOnBase()) return NULL;
322  Element* basalelement = element->SpawnBasalElement();
323 
324  ElementVector* pe = NULL;
325  switch(element->FiniteElement()){
326  case P1Enum: case P2Enum:
327  pe = CreatePVectorCG(basalelement);
328  break;
329  case P1DGEnum:
330  pe = CreatePVectorDG(basalelement);
331  break;
332  default:
333  _error_("Element type " << EnumToStringx(element->FiniteElement()) << " not supported yet");
334  }
335 
336  int domaintype;
337  element->FindParam(&domaintype,DomainTypeEnum);
338  if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
339  return pe;
340 }/*}}}*/
341 ElementVector* BalancethicknessAnalysis::CreatePVectorCG(Element* element){/*{{{*/
342 
343  /*Intermediaries */
344  IssmDouble dhdt,mb,ms,Jdet;
345  IssmDouble* xyz_list = NULL;
346 
347  /*Fetch number of nodes and dof for this finite element*/
348  int numnodes = element->GetNumberOfNodes();
349 
350  /*Initialize Element vector and other vectors*/
351  ElementVector* pe = element->NewElementVector();
352  IssmDouble* basis = xNew<IssmDouble>(numnodes);
353 
354  /*Retrieve all inputs and parameters*/
355  element->GetVerticesCoordinates(&xyz_list);
356  Input2* mb_input = element->GetInput2(BasalforcingsGroundediceMeltingRateEnum); _assert_(mb_input);
357  Input2* ms_input = element->GetInput2(SmbMassBalanceEnum); _assert_(ms_input);
358  Input2* dhdt_input = element->GetInput2(BalancethicknessThickeningRateEnum); _assert_(dhdt_input);
359 
360  /*Initialize mb_correction to 0, do not forget!:*/
361  /* Start looping on the number of gaussian points: */
362  Gauss* gauss=element->NewGauss(2);
363  for(int ig=gauss->begin();ig<gauss->end();ig++){
364  gauss->GaussPoint(ig);
365 
366  element->JacobianDeterminant(&Jdet,xyz_list,gauss);
367  element->NodalFunctions(basis,gauss);
368 
369  ms_input->GetInputValue(&ms,gauss);
370  mb_input->GetInputValue(&mb,gauss);
371  dhdt_input->GetInputValue(&dhdt,gauss);
372 
373  for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(ms-mb-dhdt)*basis[i];
374  }
375 
376  /*Clean up and return*/
377  xDelete<IssmDouble>(xyz_list);
378  xDelete<IssmDouble>(basis);
379  delete gauss;
380  return pe;
381 }/*}}}*/
382 ElementVector* BalancethicknessAnalysis::CreatePVectorDG(Element* element){/*{{{*/
383 
384  /*Intermediaries */
385  IssmDouble dhdt,mb,ms,Jdet;
386  IssmDouble* xyz_list = NULL;
387 
388  /*Fetch number of nodes and dof for this finite element*/
389  int numnodes = element->GetNumberOfNodes();
390 
391  /*Initialize Element vector and other vectors*/
392  ElementVector* pe = element->NewElementVector();
393  IssmDouble* basis = xNew<IssmDouble>(numnodes);
394 
395  /*Retrieve all inputs and parameters*/
396  element->GetVerticesCoordinates(&xyz_list);
397  Input2* mb_input = element->GetInput2(BasalforcingsGroundediceMeltingRateEnum); _assert_(mb_input);
398  Input2* ms_input = element->GetInput2(SmbMassBalanceEnum); _assert_(ms_input);
399  Input2* dhdt_input = element->GetInput2(BalancethicknessThickeningRateEnum); _assert_(dhdt_input);
400 
401  /*Initialize mb_correction to 0, do not forget!:*/
402  /* Start looping on the number of gaussian points: */
403  Gauss* gauss=element->NewGauss(2);
404  for(int ig=gauss->begin();ig<gauss->end();ig++){
405  gauss->GaussPoint(ig);
406 
407  element->JacobianDeterminant(&Jdet,xyz_list,gauss);
408  element->NodalFunctions(basis,gauss);
409 
410  ms_input->GetInputValue(&ms,gauss);
411  mb_input->GetInputValue(&mb,gauss);
412  dhdt_input->GetInputValue(&dhdt,gauss);
413 
414  for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(ms-mb-dhdt)*basis[i];
415  }
416 
417  /*Clean up and return*/
418  xDelete<IssmDouble>(xyz_list);
419  xDelete<IssmDouble>(basis);
420  delete gauss;
421  return pe;
422 }/*}}}*/
423 void BalancethicknessAnalysis::GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
424  /*Compute B matrix. B=[B1 B2 B3] where Bi is of size 3*2.
425  * For node i, Bi can be expressed in the actual coordinate system
426  * by:
427  * Bi=[ N ]
428  * [ N ]
429  * where N is the finiteelement function for node i.
430  *
431  * We assume B_prog has been allocated already, of size: 2x(1*numnodes)
432  */
433 
434  /*Fetch number of nodes for this finite element*/
435  int numnodes = element->GetNumberOfNodes();
436 
437  /*Get nodal functions*/
438  IssmDouble* basis=xNew<IssmDouble>(numnodes);
439  element->NodalFunctions(basis,gauss);
440 
441  /*Build B: */
442  for(int i=0;i<numnodes;i++){
443  B[numnodes*0+i] = basis[i];
444  B[numnodes*1+i] = basis[i];
445  }
446 
447  /*Clean-up*/
448  xDelete<IssmDouble>(basis);
449 }/*}}}*/
450 void BalancethicknessAnalysis::GetBprime(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
451  /*Compute B' matrix. B'=[B1' B2' B3'] where Bi' is of size 3*2.
452  * For node i, Bi' can be expressed in the actual coordinate system
453  * by:
454  * Bi_prime=[ dN/dx ]
455  * [ dN/dy ]
456  * where N is the finiteelement function for node i.
457  *
458  * We assume B' has been allocated already, of size: 3x(2*numnodes)
459  */
460 
461  /*Fetch number of nodes for this finite element*/
462  int numnodes = element->GetNumberOfNodes();
463 
464  /*Get nodal functions derivatives*/
465  IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
466  element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
467 
468  /*Build B': */
469  for(int i=0;i<numnodes;i++){
470  Bprime[numnodes*0+i] = dbasis[0*numnodes+i];
471  Bprime[numnodes*1+i] = dbasis[1*numnodes+i];
472  }
473 
474  /*Clean-up*/
475  xDelete<IssmDouble>(dbasis);
476 
477 }/*}}}*/
478 void BalancethicknessAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
479  _error_("not implemented yet");
480 }/*}}}*/
481 void BalancethicknessAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
482  /* WARNING: this gradient is valid for Soft balance thickness only */
483 
484  /*If on water, grad = 0: */
485  if(!element->IsIceInElement()) return;
486 
487  /*Intermediaries*/
488  IssmDouble Jdet,weight;
489  IssmDouble thickness,thicknessobs,dH[3],dp[3];
490  IssmDouble vx,vy,vel,dvx[2],dvy[2],dhdt,basal_melting,surface_mass_balance;
491  IssmDouble *xyz_list= NULL;
492 
493  /*Get list of cost functions*/
494  int *responses = NULL;
495  int num_responses,resp,solution;
496  element->FindParam(&num_responses,InversionNumCostFunctionsEnum);
497  element->FindParam(&responses,NULL,InversionCostFunctionsEnum);
498  element->FindParam(&solution,SolutionTypeEnum);
499  if(solution!=BalancethicknessSoftSolutionEnum) _error_("not implemented yet");
500  if(control_type!=ThicknessEnum) _error_("Control "<<EnumToStringx(control_type)<<" not supported");
501 
502  /*Fetch number of vertices for this finite element*/
503  int numvertices = element->GetNumberOfVertices();
504 
505  /*Initialize some vectors*/
506  IssmDouble* basis = xNew<IssmDouble>(numvertices);
507  IssmDouble* dbasis = xNew<IssmDouble>(2*numvertices);
508  IssmDouble* ge = xNewZeroInit<IssmDouble>(numvertices);
509  int* vertexpidlist = xNew<int>(numvertices);
510 
511  /*Retrieve all inputs we will be needing: */
512  element->GetVerticesCoordinates(&xyz_list);
513  element->GradientIndexing(&vertexpidlist[0],control_index);
514  DatasetInput2* weights_input = element->GetDatasetInput2(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
515  Input2* thickness_input = element->GetInput2(ThicknessEnum); _assert_(thickness_input);
516  Input2* thicknessobs_input = element->GetInput2(InversionThicknessObsEnum); _assert_(thicknessobs_input);
517  Input2* vx_input = element->GetInput2(VxEnum); _assert_(vx_input);
518  Input2* vy_input = element->GetInput2(VyEnum); _assert_(vy_input);
519  Input2* surface_mass_balance_input = element->GetInput2(SmbMassBalanceEnum); _assert_(surface_mass_balance_input);
520  Input2* basal_melting_input = element->GetInput2(BasalforcingsGroundediceMeltingRateEnum); _assert_(basal_melting_input);
521  Input2* dhdt_input = element->GetInput2(BalancethicknessThickeningRateEnum); _assert_(dhdt_input);
522 
523  /* Start looping on the number of gaussian points: */
524  Gauss* gauss=element->NewGauss(4);
525  for(int ig=gauss->begin();ig<gauss->end();ig++){
526  gauss->GaussPoint(ig);
527 
528  thickness_input->GetInputValue(&thickness, gauss);
529  thickness_input->GetInputDerivativeValue(&dH[0],xyz_list,gauss);
530  thicknessobs_input->GetInputValue(&thicknessobs, gauss);
531 
532  element->JacobianDeterminant(&Jdet,xyz_list,gauss);
533  element->NodalFunctionsP1(basis,gauss);
534  element->NodalFunctionsP1Derivatives(dbasis,xyz_list,gauss);
535 
536  /*Deal with first part (partial derivative a J with respect to k)*/
537  for(resp=0;resp<num_responses;resp++){
538 
539  weights_input->GetInputValue(&weight,gauss,responses[resp]);
540 
541  switch(responses[resp]){
542  case ThicknessAbsMisfitEnum:
543  for(int i=0;i<numvertices;i++) ge[i]+= (thicknessobs-thickness)*weight*Jdet*gauss->weight*basis[i];
544  break;
545  case ThicknessAbsGradientEnum:
546  for(int i=0;i<numvertices;i++) ge[i]+= - weight*dH[0]*dbasis[0*numvertices+i]*Jdet*gauss->weight;
547  for(int i=0;i<numvertices;i++) ge[i]+= - weight*dH[1]*dbasis[1*numvertices+i]*Jdet*gauss->weight;
548  break;
549  case ThicknessAlongGradientEnum:
550  vx_input->GetInputValue(&vx,gauss);
551  vy_input->GetInputValue(&vy,gauss);
552  vel = sqrt(vx*vx+vy*vy);
553  vx = vx/(vel+1.e-9);
554  vy = vy/(vel+1.e-9);
555  for(int i=0;i<numvertices;i++) ge[i]+= - weight*(dH[0]*vx+dH[1]*vy)*(dbasis[0*numvertices+i]*vx+dbasis[1*numvertices+i]*vy)*Jdet*gauss->weight;
556  break;
557  case ThicknessAcrossGradientEnum:
558  vx_input->GetInputValue(&vx,gauss);
559  vy_input->GetInputValue(&vy,gauss);
560  vel = sqrt(vx*vx+vy*vy);
561  vx = vx/(vel+1.e-9);
562  vy = vy/(vel+1.e-9);
563  for(int i=0;i<numvertices;i++) ge[i]+= - weight*(dH[0]*(-vy)+dH[1]*vx)*(dbasis[0*numvertices+i]*(-vy)+dbasis[1*numvertices+i]*vx)*Jdet*gauss->weight;
564  break;
565  case BalancethicknessMisfitEnum:
566  surface_mass_balance_input->GetInputValue(&surface_mass_balance,gauss);
567  basal_melting_input->GetInputValue(&basal_melting,gauss);
568  dhdt_input->GetInputValue(&dhdt,gauss);
569  vx_input->GetInputValue(&vx,gauss);
570  vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
571  vy_input->GetInputValue(&vy,gauss);
572  vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
573  for(int i=0;i<numvertices;i++){
574  ge[i]+= - weight*Jdet*gauss->weight*(
575  (vx*dH[0]+vy*dH[1] + thickness*(dvx[0]+dvy[1]))*(vx*dbasis[0*numvertices+i]+ vy*dbasis[1*numvertices+i] + basis[i]*(dvx[0]+dvy[1]))
576  -(surface_mass_balance-basal_melting-dhdt)*(vx*dbasis[0*numvertices+i]+ vy*dbasis[1*numvertices+i] + basis[i]*(dvx[0]+dvy[1]))
577  );
578  }
579  break;
580  default:
581  _error_("response " << EnumToStringx(responses[resp]) << " not supported yet");
582  }
583  }
584  }
585  gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL);
586 
587  /*Clean up and return*/
588  xDelete<IssmDouble>(xyz_list);
589  xDelete<IssmDouble>(basis);
590  xDelete<IssmDouble>(ge);
591  xDelete<int>(vertexpidlist);
592  xDelete<int>(responses);
593  delete gauss;
594 
595 }/*}}}*/
596 void BalancethicknessAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
597 
598  int domaintype;
599  element->FindParam(&domaintype,DomainTypeEnum);
600  switch(domaintype){
601  case Domain2DhorizontalEnum:
602  element->InputUpdateFromSolutionOneDof(solution,ThicknessEnum);
603  break;
604  case Domain3DEnum:
605  element->InputUpdateFromSolutionOneDofCollapsed(solution,ThicknessEnum);
606  break;
607  default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
608  }
609 }/*}}}*/
610 void BalancethicknessAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
611  /*Default, do nothing*/
612  return;
613 }/*}}}*/
BalancethicknessAnalysisEnum
@ BalancethicknessAnalysisEnum
Definition: EnumDefinitions.h:981
InversionThicknessObsEnum
@ InversionThicknessObsEnum
Definition: EnumDefinitions.h:631
BaseEnum
@ BaseEnum
Definition: EnumDefinitions.h:495
SmbMassBalanceEnum
@ SmbMassBalanceEnum
Definition: EnumDefinitions.h:748
_assert_
#define _assert_(ignore)
Definition: exceptions.h:37
IssmDouble
double IssmDouble
Definition: types.h:37
Element::IsOnBase
bool IsOnBase()
Definition: Element.cpp:1984
Nodes
Declaration of Nodes class.
Definition: Nodes.h:19
DatasetInput2
Definition: DatasetInput2.h:14
Element::GetNumberOfNodes
virtual int GetNumberOfNodes(void)=0
Element::FindParam
void FindParam(bool *pvalue, int paramenum)
Definition: Element.cpp:933
BalancethicknessThickeningRateEnum
@ BalancethicknessThickeningRateEnum
Definition: EnumDefinitions.h:474
DataSet::AddObject
int AddObject(Object *object)
Definition: DataSet.cpp:252
BalancethicknessMisfitEnum
@ BalancethicknessMisfitEnum
Definition: EnumDefinitions.h:471
InversionCostFunctionsCoefficientsEnum
@ InversionCostFunctionsCoefficientsEnum
Definition: EnumDefinitions.h:629
Parameters
Declaration of Parameters class.
Definition: Parameters.h:18
MaskIceLevelsetEnum
@ MaskIceLevelsetEnum
Definition: EnumDefinitions.h:641
BalancethicknessSoftSolutionEnum
@ BalancethicknessSoftSolutionEnum
Definition: EnumDefinitions.h:984
ADD_VAL
@ ADD_VAL
Definition: toolkitsenums.h:14
BalancethicknessAnalysis.h
: header file for generic external result object
Constraints
Declaration of Constraints class.
Definition: Constraints.h:13
Element::SpawnBasalElement
virtual Element * SpawnBasalElement(void)=0
BalancethicknessAnalysis::UpdateElements
void UpdateElements(Elements *elements, Inputs2 *inputs2, IoModel *iomodel, int analysis_counter, int analysis_type)
Definition: BalancethicknessAnalysis.cpp:75
MeshVertexonbaseEnum
@ MeshVertexonbaseEnum
Definition: EnumDefinitions.h:653
BalancethicknessAnalysis::CreateKMatrixDG
ElementMatrix * CreateKMatrixDG(Element *element)
Definition: BalancethicknessAnalysis.cpp:253
Elements
Declaration of Elements class.
Definition: Elements.h:17
TripleMultiply
int TripleMultiply(IssmDouble *a, int nrowa, int ncola, int itrna, IssmDouble *b, int nrowb, int ncolb, int itrnb, IssmDouble *c, int nrowc, int ncolc, int itrnc, IssmDouble *d, int iaddd)
Definition: MatrixUtils.cpp:20
IoModel::numberoffaces
int numberoffaces
Definition: IoModel.h:97
P1DGEnum
@ P1DGEnum
Definition: EnumDefinitions.h:1215
ElementVector::values
IssmDouble * values
Definition: ElementVector.h:24
IoModel::my_elements
bool * my_elements
Definition: IoModel.h:66
DatasetInput2::GetInputValue
void GetInputValue(IssmDouble *pvalue, Gauss *gauss, int index)
Definition: DatasetInput2.cpp:199
Element::NodalFunctionsP1Derivatives
virtual void NodalFunctionsP1Derivatives(IssmDouble *dbasis, IssmDouble *xyz_list, Gauss *gauss)=0
BalancethicknessAnalysis::GetSolutionFromInputs
void GetSolutionFromInputs(Vector< IssmDouble > *solution, Element *element)
Definition: BalancethicknessAnalysis.cpp:478
BalancethicknessAnalysis::GetBprime
void GetBprime(IssmDouble *B, Element *element, IssmDouble *xyz_list, Gauss *gauss)
Definition: BalancethicknessAnalysis.cpp:450
VyEnum
@ VyEnum
Definition: EnumDefinitions.h:850
Parameters::AddObject
void AddObject(Param *newparam)
Definition: Parameters.cpp:67
BalancethicknessAnalysis::GradientJ
void GradientJ(Vector< IssmDouble > *gradient, Element *element, int control_type, int control_index)
Definition: BalancethicknessAnalysis.cpp:481
Input2::GetInputDerivativeValue
virtual void GetInputDerivativeValue(IssmDouble *derivativevalues, IssmDouble *xyz_list, Gauss *gauss)
Definition: Input2.h:37
Element::GetInput2
virtual Input2 * GetInput2(int inputenum)=0
Element::DeleteMaterials
void DeleteMaterials(void)
Definition: Element.cpp:429
Element
Definition: Element.h:41
Element::NodalFunctions
virtual void NodalFunctions(IssmDouble *basis, Gauss *gauss)=0
P1Enum
@ P1Enum
Definition: EnumDefinitions.h:662
Element::GetDatasetInput2
virtual DatasetInput2 * GetDatasetInput2(int inputenum)
Definition: Element.h:250
Domain2DhorizontalEnum
@ Domain2DhorizontalEnum
Definition: EnumDefinitions.h:534
BalancethicknessAnalysis::CreateLoads
void CreateLoads(Loads *loads, IoModel *iomodel)
Definition: BalancethicknessAnalysis.cpp:21
Element::NodalFunctionsP1
virtual void NodalFunctionsP1(IssmDouble *basis, Gauss *gauss)=0
Element::NewElementVector
ElementVector * NewElementVector(int approximation_enum=NoneApproximationEnum)
Definition: Element.cpp:2505
IoModel::DeleteData
void DeleteData(int num,...)
Definition: IoModel.cpp:500
IoModel::numberofelements
int numberofelements
Definition: IoModel.h:96
IoModel::CopyConstantObject
Param * CopyConstantObject(const char *constant_name, int param_enum)
Definition: IoModel.cpp:418
BasalforcingsGroundediceMeltingRateEnum
@ BasalforcingsGroundediceMeltingRateEnum
Definition: EnumDefinitions.h:478
Vector::SetValues
void SetValues(int ssize, int *list, doubletype *values, InsMode mode)
Definition: Vector.h:153
DomainTypeEnum
@ DomainTypeEnum
Definition: EnumDefinitions.h:124
VxAverageEnum
@ VxAverageEnum
Definition: EnumDefinitions.h:845
BalancethicknessAnalysis::UpdateParameters
void UpdateParameters(Parameters *parameters, IoModel *iomodel, int solution_enum, int analysis_enum)
Definition: BalancethicknessAnalysis.cpp:114
Element::NewGauss
virtual Gauss * NewGauss(void)=0
InversionCostFunctionsEnum
@ InversionCostFunctionsEnum
Definition: EnumDefinitions.h:211
BalancethicknessAnalysis::CreateDVector
ElementVector * CreateDVector(Element *element)
Definition: BalancethicknessAnalysis.cpp:122
BalancethicknessAnalysis::CreateKMatrixCG
ElementMatrix * CreateKMatrixCG(Element *element)
Definition: BalancethicknessAnalysis.cpp:151
Element::InputUpdateFromSolutionOneDof
virtual void InputUpdateFromSolutionOneDof(IssmDouble *solution, int inputenum)=0
EnumToStringx
const char * EnumToStringx(int enum_in)
Definition: EnumToStringx.cpp:15
BalancethicknessAnalysis::CreatePVector
ElementVector * CreatePVector(Element *element)
Definition: BalancethicknessAnalysis.cpp:319
ThicknessAlongGradientEnum
@ ThicknessAlongGradientEnum
Definition: EnumDefinitions.h:839
IoModel::FindConstant
void FindConstant(bool *pvalue, const char *constant_name)
Definition: IoModel.cpp:2362
BalancethicknessAnalysis::CreateConstraints
void CreateConstraints(Constraints *constraints, IoModel *iomodel)
Definition: BalancethicknessAnalysis.cpp:9
Element::GetVerticesCoordinates
void GetVerticesCoordinates(IssmDouble **xyz_list)
Definition: Element.cpp:1446
SolutionTypeEnum
@ SolutionTypeEnum
Definition: EnumDefinitions.h:398
ThicknessAbsMisfitEnum
@ ThicknessAbsMisfitEnum
Definition: EnumDefinitions.h:837
BalancethicknessAnalysis::GetB
void GetB(IssmDouble *B, Element *element, IssmDouble *xyz_list, Gauss *gauss)
Definition: BalancethicknessAnalysis.cpp:423
BalancethicknessAnalysis::DofsPerNode
int DofsPerNode(int **doflist, int domaintype, int approximation)
Definition: BalancethicknessAnalysis.cpp:72
SurfaceEnum
@ SurfaceEnum
Definition: EnumDefinitions.h:823
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
Input2
Definition: Input2.h:18
FemModel
Definition: FemModel.h:31
BalancethicknessAnalysis::Core
void Core(FemModel *femmodel)
Definition: BalancethicknessAnalysis.cpp:119
SealevelEnum
@ SealevelEnum
Definition: EnumDefinitions.h:675
Element::IsIceInElement
bool IsIceInElement()
Definition: Element.cpp:2021
Loads
Declaration of Loads class.
Definition: Loads.h:16
BalancethicknessAnalysis::CreateNodes
void CreateNodes(Nodes *nodes, IoModel *iomodel, bool isamr=false)
Definition: BalancethicknessAnalysis.cpp:54
Numericalflux
Definition: Numericalflux.h:18
_error_
#define _error_(StreamArgs)
Definition: exceptions.h:49
BalancethicknessAnalysis::UpdateConstraints
void UpdateConstraints(FemModel *femmodel)
Definition: BalancethicknessAnalysis.cpp:610
IoModel::faces
int * faces
Definition: IoModel.h:88
Gauss::begin
virtual int begin(void)=0
DataSet::GetObjectByOffset
Object * GetObjectByOffset(int offset)
Definition: DataSet.cpp:334
Element::CharacteristicLength
virtual IssmDouble CharacteristicLength(void)=0
VxEnum
@ VxEnum
Definition: EnumDefinitions.h:846
MeshVertexonsurfaceEnum
@ MeshVertexonsurfaceEnum
Definition: EnumDefinitions.h:655
InversionNumCostFunctionsEnum
@ InversionNumCostFunctionsEnum
Definition: EnumDefinitions.h:224
Element::JacobianDeterminant
virtual void JacobianDeterminant(IssmDouble *Jdet, IssmDouble *xyz_list, Gauss *gauss)=0
Gauss::GaussPoint
virtual void GaussPoint(int ig)=0
Element::Update
virtual void Update(Inputs2 *inputs2, int index, IoModel *iomodel, int analysis_counter, int analysis_type, int finite_element)=0
VyAverageEnum
@ VyAverageEnum
Definition: EnumDefinitions.h:849
Element::FiniteElement
virtual int FiniteElement(void)=0
ThicknessEnum
@ ThicknessEnum
Definition: EnumDefinitions.h:840
BalancethicknessAnalysis::CreateJacobianMatrix
ElementMatrix * CreateJacobianMatrix(Element *element)
Definition: BalancethicknessAnalysis.cpp:126
BalancethicknessAnalysis::InputUpdateFromSolution
void InputUpdateFromSolution(IssmDouble *solution, Element *element)
Definition: BalancethicknessAnalysis.cpp:596
BalancethicknessAnalysis::CreatePVectorDG
ElementVector * CreatePVectorDG(Element *element)
Definition: BalancethicknessAnalysis.cpp:382
IoModel::FetchDataToInput
void FetchDataToInput(Inputs2 *inputs2, Elements *elements, const char *vector_name, int input_enum)
Definition: IoModel.cpp:1651
P2Enum
@ P2Enum
Definition: EnumDefinitions.h:1223
IoModelToConstraintsx
void IoModelToConstraintsx(Constraints *constraints, IoModel *iomodel, const char *spc_name, int analysis_type, int finite_element, int dof)
Definition: IoModelToConstraintsx.cpp:10
ElementVector
Definition: ElementVector.h:20
Input2::GetInputValue
virtual void GetInputValue(IssmDouble *pvalue, Gauss *gauss)
Definition: Input2.h:38
Gauss::weight
IssmDouble weight
Definition: Gauss.h:11
IoModel
Definition: IoModel.h:48
CreateFaces
void CreateFaces(IoModel *iomodel)
Definition: CreateFaces.cpp:9
BalancethicknessAnalysis::CreatePVectorCG
ElementVector * CreatePVectorCG(Element *element)
Definition: BalancethicknessAnalysis.cpp:341
BalancethicknessAnalysis::CreateKMatrix
ElementMatrix * CreateKMatrix(Element *element)
Definition: BalancethicknessAnalysis.cpp:129
ThicknessAbsGradientEnum
@ ThicknessAbsGradientEnum
Definition: EnumDefinitions.h:836
Element::GetNumberOfVertices
virtual int GetNumberOfVertices(void)=0
IoModel::domaintype
int domaintype
Definition: IoModel.h:78
ElementMatrix
Definition: ElementMatrix.h:19
BalancethicknessStabilizationEnum
@ BalancethicknessStabilizationEnum
Definition: EnumDefinitions.h:58
Vector< IssmDouble >
Element::GradientIndexing
void GradientIndexing(int *indexing, int control_index)
Definition: Element.cpp:1510
Input2::GetInputAverage
virtual void GetInputAverage(IssmDouble *pvalue)
Definition: Input2.h:33
Element::InputUpdateFromSolutionOneDofCollapsed
virtual void InputUpdateFromSolutionOneDofCollapsed(IssmDouble *solution, int inputenum)=0
Gauss
Definition: Gauss.h:8
Element::NodalFunctionsDerivatives
virtual void NodalFunctionsDerivatives(IssmDouble *dbasis, IssmDouble *xyz_list, Gauss *gauss)=0
ElementMatrix::values
IssmDouble * values
Definition: ElementMatrix.h:26
ThicknessAcrossGradientEnum
@ ThicknessAcrossGradientEnum
Definition: EnumDefinitions.h:838
Element::NewElementMatrix
ElementMatrix * NewElementMatrix(int approximation_enum=NoneApproximationEnum)
Definition: Element.cpp:2497
femmodel
FemModel * femmodel
Definition: esmfbinders.cpp:16
Generated on Thu Jul 2 2020 08:09:14 for Ice Sheet System Model by   doxygen 1.8.19