Ice Sheet System Model  4.18
Code documentation
SealevelriseAnalysis.cpp
Go to the documentation of this file.
1 #include "./SealevelriseAnalysis.h"
2 #include <math.h>
3 #include "../toolkits/toolkits.h"
4 #include "../classes/classes.h"
5 #include "../classes/Inputs2/TransientInput2.h"
6 #include "../shared/shared.h"
7 #include "../modules/modules.h"
8 
9 /*Model processing*/
10 void SealevelriseAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
11  /*No constraints*/
12 }/*}}}*/
13 void SealevelriseAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
14  /*No loads*/
15 }/*}}}*/
16 void SealevelriseAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel,bool isamr){/*{{{*/
17  ::CreateNodes(nodes,iomodel,SealevelriseAnalysisEnum,P1Enum);
18 }/*}}}*/
19 int SealevelriseAnalysis::DofsPerNode(int** doflist,int domaintype,int approximation){/*{{{*/
20  return 1;
21 }/*}}}*/
22 void SealevelriseAnalysis::UpdateElements(Elements* elements,Inputs2* inputs2,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
23 
24  int geodetic=0;
25  int dslmodel=0;
26 
27  /*Update elements: */
28  int counter=0;
29  for(int i=0;i<iomodel->numberofelements;i++){
30  if(iomodel->my_elements[i]){
31  Element* element=(Element*)elements->GetObjectByOffset(counter);
32  element->Update(inputs2,i,iomodel,analysis_counter,analysis_type,P1Enum);
33  counter++;
34  }
35  }
36 
37  /*Create inputs: */
38  iomodel->FetchDataToInput(inputs2,elements,"md.mask.ocean_levelset",MaskOceanLevelsetEnum);
39  iomodel->FetchDataToInput(inputs2,elements,"md.mask.ice_levelset",MaskIceLevelsetEnum);
40  iomodel->FetchData(&geodetic,"md.solidearth.settings.computesealevelchange");
41  iomodel->FetchDataToInput(inputs2,elements,"md.solidearth.surfaceload.icethicknesschange",SurfaceloadIceThicknessChangeEnum);
42  iomodel->FetchDataToInput(inputs2,elements,"md.solidearth.sealevel",SealevelEnum,0);
43  iomodel->FetchDataToInput(inputs2,elements,"md.geometry.bed",BedEnum);
44  iomodel->FetchDataToInput(inputs2,elements,"md.solidearth.surfaceload.waterheightchange",SurfaceloadWaterHeightChangeEnum);
45 
46  /*dynamic sea level: */
47  iomodel->FetchData(&dslmodel,"md.dsl.model");
48  if (dslmodel==1){ /*standard dsl model:{{{*/
49 
50  /*deal with global mean steric rate: */
51  IssmDouble* str=NULL;
52  IssmDouble* times = NULL;
53  int M,N;
54 
55  /*fetch str vector:*/
56  iomodel->FetchData(&str,&M,&N,"md.dsl.global_average_thermosteric_sea_level_change"); _assert_(M==2);
57 
58  //recover time vector:
59  times=xNew<IssmDouble>(N);
60  for(int t=0;t<N;t++) times[t] = str[N+t];
61 
62  /*create transient input: */
63  inputs2->SetTransientInput(DslGlobalAverageThermostericSeaLevelChangeEnum,times,N);
64  TransientInput2* transientinput = inputs2->GetTransientInput(DslGlobalAverageThermostericSeaLevelChangeEnum);
65 
66 
67  for(int i=0;i<elements->Size();i++){
68  Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
69 
70  for(int t=0;t<N;t++){
71  switch(element->ObjectEnum()){
72  case TriaEnum: transientinput->AddTriaTimeInput( t,1,&element->lid,&str[t],P0Enum); break;
73  case PentaEnum: transientinput->AddPentaTimeInput(t,1,&element->lid,&str[t],P0Enum); break;
74  default: _error_("Not implemented yet");
75  }
76  }
77  }
78 
79  /*cleanup:*/
80  xDelete<IssmDouble>(times);
81  iomodel->DeleteData(str,"md.dsl.global_average_thermosteric_sea_level_change");
82 
83  /*deal with dynamic sea level fields: */
84  iomodel->FetchDataToInput(inputs2,elements,"md.dsl.sea_surface_height_change_above_geoid", DslSeaSurfaceHeightChangeAboveGeoidEnum);
85  iomodel->FetchDataToInput(inputs2,elements,"md.dsl.sea_water_pressure_change_at_sea_floor", DslSeaWaterPressureChangeAtSeaFloor);
86 
87  } /*}}}*/
88  else if (dslmodel==2){ /*multi-model ensemble dsl model:{{{*/
89 
90  /*variables:*/
91  int nummodels;
92  IssmDouble** pstr=NULL;
93  IssmDouble* str=NULL;
94  IssmDouble* times = NULL;
95  int* pM = NULL;
96  int* pN = NULL;
97  int M,N;
98 
99  /*deal with dsl.sea_surface_height_change_above_geoid {{{*/
100  iomodel->FetchData(&pstr,&pM,&pN,&nummodels,"md.dsl.global_average_thermosteric_sea_level_change");
101 
102  /*go through the mat array and create a dataset of transient inputs:*/
103  for (int i=0;i<nummodels;i++){
104 
105  M=pM[i];
106  N=pN[i];
107  str=pstr[i];
108 
109  //recover time vector:
110  times=xNew<IssmDouble>(N);
111  for(int t=0;t<N;t++) times[t] = str[(M-1)*N+t];
112 
113  TransientInput2* transientinput=inputs2->SetDatasetTransientInput(DslGlobalAverageThermostericSeaLevelChangeEnum,i, times,N);
114 
115  for(int j=0;j<elements->Size();j++){
116  Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(j));
117 
118  for(int t=0;t<N;t++){
119  switch(element->ObjectEnum()){
120  case TriaEnum: transientinput->AddTriaTimeInput( t,1,&element->lid,&str[t],P0Enum); break;
121  case PentaEnum: transientinput->AddPentaTimeInput(t,1,&element->lid,&str[t],P0Enum); break;
122  default: _error_("Not implemented yet");
123  }
124  }
125  }
126  xDelete<IssmDouble>(times);
127  }
128  /*Delete data:*/
129  for(int i=0;i<nummodels;i++){
130  IssmDouble* str=pstr[i];
131  xDelete<IssmDouble>(str);
132  }
133  xDelete<IssmDouble*>(pstr);
134  xDelete<int>(pM);
135  xDelete<int>(pN);
136  /*}}}*/
137  /*now do the same with the dsl.sea_surface_height_change_above_geoid field:{{{ */
138  iomodel->FetchData(&pstr,&pM,&pN,&nummodels,"md.dsl.sea_surface_height_change_above_geoid");
139 
140  for (int i=0;i<nummodels;i++){
141  M=pM[i];
142  N=pN[i];
143  str=pstr[i];
144 
145 
146  //recover time vector:
147  times=xNew<IssmDouble>(N);
148  for(int t=0;t<N;t++) times[t] = str[(M-1)*N+t];
149 
150  TransientInput2* transientinput=inputs2->SetDatasetTransientInput(DslSeaSurfaceHeightChangeAboveGeoidEnum,i, times,N);
151 
152  for(int j=0;j<elements->Size();j++){
153 
154  /*Get the right transient input*/
155  Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(j));
156 
157  /*Get values and lid list*/
158  const int numvertices = element->GetNumberOfVertices();
159  int *vertexlids = xNew<int>(numvertices);
160  int *vertexsids = xNew<int>(numvertices);
161 
162 
163  /*Recover vertices ids needed to initialize inputs*/
164  _assert_(iomodel->elements);
165  for(int k=0;k<numvertices;k++){
166  vertexsids[k] =reCast<int>(iomodel->elements[numvertices*element->Sid()+k]); //ids for vertices are in the elements array from Matlab
167  vertexlids[k]=iomodel->my_vertices_lids[vertexsids[k]-1];
168  }
169 
170  //element->GetVerticesLidList(vertexlids);
171  //element->GetVerticesSidList(vertexsids);
172  IssmDouble* values=xNew<IssmDouble>(numvertices);
173 
174  for(int t=0;t<N;t++){
175  for (int k=0;k<numvertices;k++)values[k]=str[N*vertexsids[k]+t];
176 
177  switch(element->ObjectEnum()){
178  case TriaEnum: transientinput->AddTriaTimeInput( t,numvertices,vertexlids,values,P1Enum); break;
179  case PentaEnum: transientinput->AddPentaTimeInput(t,numvertices,vertexlids,values,P1Enum); break;
180  default: _error_("Not implemented yet");
181  }
182  }
183  xDelete<IssmDouble>(values);
184  xDelete<int>(vertexlids);
185  xDelete<int>(vertexsids);
186  }
187 
188  xDelete<IssmDouble>(times);
189  }
190 
191  /*Delete data:*/
192  for(int i=0;i<nummodels;i++){
193  IssmDouble* str=pstr[i];
194  xDelete<IssmDouble>(str);
195  }
196  xDelete<IssmDouble*>(pstr);
197  xDelete<int>(pM);
198  xDelete<int>(pN);
199  /*}}}*/
200  /*now do the same with the dsl.sea_water_pressure_change_at_sea_floor field:{{{ */
201  iomodel->FetchData(&pstr,&pM,&pN,&nummodels,"md.dsl.sea_water_pressure_change_at_sea_floor");
202 
203  for (int i=0;i<nummodels;i++){
204  M=pM[i];
205  N=pN[i];
206  str=pstr[i];
207 
208  //recover time vector:
209  times=xNew<IssmDouble>(N);
210  for(int t=0;t<N;t++) times[t] = str[(M-1)*N+t];
211 
212  TransientInput2* transientinput=inputs2->SetDatasetTransientInput(DslSeaWaterPressureChangeAtSeaFloor,i, times,N);
213 
214  for(int j=0;j<elements->Size();j++){
215 
216  /*Get the right transient input*/
217  Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(j));
218 
219  /*Get values and lid list*/
220  const int numvertices = element->GetNumberOfVertices();
221  int *vertexlids = xNew<int>(numvertices);
222  int *vertexsids = xNew<int>(numvertices);
223 
224 
225  /*Recover vertices ids needed to initialize inputs*/
226  _assert_(iomodel->elements);
227  for(int k=0;k<numvertices;k++){
228  vertexsids[k] =reCast<int>(iomodel->elements[numvertices*element->Sid()+k]); //ids for vertices are in the elements array from Matlab
229  vertexlids[k]=iomodel->my_vertices_lids[vertexsids[k]-1];
230  }
231  //element->GetVerticesLidList(vertexlids);
232  //element->GetVerticesSidList(vertexsids);
233 
234  IssmDouble* values=xNew<IssmDouble>(numvertices);
235 
236  for(int t=0;t<N;t++){
237  for (int k=0;k<numvertices;k++)values[k]=str[N*vertexsids[k]+t];
238 
239  switch(element->ObjectEnum()){
240  case TriaEnum: transientinput->AddTriaTimeInput( t,numvertices,vertexlids,values,P1Enum); break;
241  case PentaEnum: transientinput->AddPentaTimeInput(t,numvertices,vertexlids,values,P1Enum); break;
242  default: _error_("Not implemented yet");
243  }
244  }
245  xDelete<IssmDouble>(values);
246  xDelete<int>(vertexlids);
247  xDelete<int>(vertexsids);
248  }
249  xDelete<IssmDouble>(times);
250  }
251 
252  /*Delete data:*/
253  for(int i=0;i<nummodels;i++){
254  IssmDouble* str=pstr[i];
255  xDelete<IssmDouble>(str);
256  }
257  xDelete<IssmDouble*>(pstr);
258  xDelete<int>(pM);
259  xDelete<int>(pN);
260  /*}}}*/
261 
262  } /*}}}*/
263  else _error_("Dsl model " << dslmodel << " not implemented yet!");
264 
265 
266 
267  /*Initialize cumdeltalthickness and sealevel rise rate input*/
268  InputUpdateFromConstantx(inputs2,elements,0.,SealevelriseCumDeltathicknessEnum);
269  InputUpdateFromConstantx(inputs2,elements,0.,SealevelNEsaRateEnum);
270  InputUpdateFromConstantx(inputs2,elements,0.,SealevelUEsaRateEnum);
271  InputUpdateFromConstantx(inputs2,elements,0.,SealevelRSLRateEnum);
272  InputUpdateFromConstantx(inputs2,elements,0.,SealevelEustaticMaskEnum);
273  InputUpdateFromConstantx(inputs2,elements,0.,SealevelEustaticOceanMaskEnum);
274 
275 }/*}}}*/
276 void SealevelriseAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
277 
278  int nl;
279  IssmDouble* love_h=NULL;
280  IssmDouble* love_k=NULL;
281  IssmDouble* love_l=NULL;
282  IssmDouble* love_th=NULL;
283  IssmDouble* love_tk=NULL;
284  IssmDouble* love_tl=NULL;
285  int dslmodel=0;
286 
287  IssmDouble* G_rigid = NULL;
288  IssmDouble* G_rigid_local = NULL;
289  IssmDouble* G_elastic = NULL;
290  IssmDouble* G_elastic_local = NULL;
291  IssmDouble* U_elastic = NULL;
292  IssmDouble* U_elastic_local = NULL;
293  IssmDouble* H_elastic = NULL;
294  IssmDouble* H_elastic_local = NULL;
295  int M,m,lower_row,upper_row;
296  IssmDouble degacc=.01;
297  IssmDouble planetradius=0;
298  IssmDouble planetarea=0;
299 
300  int numoutputs;
301  char** requestedoutputs = NULL;
302 
303  /*transition vectors: */
304  IssmDouble **transitions = NULL;
305  int *transitions_M = NULL;
306  int *transitions_N = NULL;
307  int ntransitions;
308 
309  /*some constant parameters: */
310  parameters->AddObject(iomodel->CopyConstantObject("md.dsl.model",DslModelEnum));
311  parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.runfrequency",SolidearthSettingsRunFrequencyEnum));
312  parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.reltol",SolidearthSettingsReltolEnum));
313  parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.abstol",SolidearthSettingsAbstolEnum));
314  parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.maxiter",SolidearthSettingsMaxiterEnum));
315  parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.rigid",SolidearthSettingsRigidEnum));
316  parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.horiz",SolidearthSettingsHorizEnum));
317  parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.elastic",SolidearthSettingsElasticEnum));
318  parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.rotation",SolidearthSettingsRotationEnum));
319  parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.rotational.equatorialmoi",RotationalEquatorialMoiEnum));
320  parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.rotational.polarmoi",RotationalPolarMoiEnum));
321  parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.rotational.angularvelocity",RotationalAngularVelocityEnum));
322  parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.ocean_area_scaling",SolidearthSettingsOceanAreaScalingEnum));
323  parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.computesealevelchange",SolidearthSettingsComputesealevelchangeEnum));
324  parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.planetradius",SolidearthPlanetRadiusEnum));
325 
326  /*compute planet area and plug into parameters:*/
327  iomodel->FetchData(&planetradius,"md.solidearth.planetradius");
328  planetarea=4*PI*planetradius*planetradius;
329  parameters->AddObject(new DoubleParam(SolidearthPlanetAreaEnum,planetarea));
330 
331  /*Deal with dsl multi-model ensembles: {{{*/
332  iomodel->FetchData(&dslmodel,"md.dsl.model");
333  parameters->AddObject(iomodel->CopyConstantObject("md.dsl.compute_fingerprints",DslComputeFingerprintsEnum));
334  if(dslmodel==2){
335  int modelid;
336  int nummodels;
337 
338  parameters->AddObject(iomodel->CopyConstantObject("md.dsl.modelid",DslModelidEnum));
339  parameters->AddObject(iomodel->CopyConstantObject("md.dsl.nummodels",DslNummodelsEnum));
340  iomodel->FetchData(&modelid,"md.dsl.modelid");
341  iomodel->FetchData(&nummodels,"md.dsl.nummodels");
342 
343  /*quick checks: */
344  if(nummodels<=0)_error_("dslmme object in md.dsl field should contain at least 1 ensemble model!");
345  if(modelid<=0 || modelid>nummodels)_error_("modelid field in dslmme object of md.dsl field should be between 1 and the number of ensemble runs!");
346  } /*}}}*/
347  /*Deal with elasticity {{{*/
348  /*love numbers: */
349  iomodel->FetchData(&love_h,&nl,NULL,"md.solidearth.lovenumbers.h");
350  iomodel->FetchData(&love_k,&nl,NULL,"md.solidearth.lovenumbers.k");
351  iomodel->FetchData(&love_l,&nl,NULL,"md.solidearth.lovenumbers.l");
352  iomodel->FetchData(&love_th,&nl,NULL,"md.solidearth.lovenumbers.th");
353  iomodel->FetchData(&love_tk,&nl,NULL,"md.solidearth.lovenumbers.tk");
354  iomodel->FetchData(&love_tl,&nl,NULL,"md.solidearth.lovenumbers.tl");
355  parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.lovenumbers.tk2secular",TidalLoveK2SecularEnum));
356 
357  parameters->AddObject(new DoubleMatParam(LoadLoveHEnum,love_h,nl,1));
358  parameters->AddObject(new DoubleMatParam(LoadLoveKEnum,love_k,nl,1));
359  parameters->AddObject(new DoubleMatParam(LoadLoveLEnum,love_l,nl,1));
360  parameters->AddObject(new DoubleMatParam(TidalLoveHEnum,love_th,nl,1));
361  parameters->AddObject(new DoubleMatParam(TidalLoveKEnum,love_tk,nl,1));
362  parameters->AddObject(new DoubleMatParam(TidalLoveLEnum,love_tl,nl,1));
363 
364  /*compute elastic green function for a range of angles*/
365  iomodel->FetchData(&degacc,"md.solidearth.settings.degacc");
366  M=reCast<int,IssmDouble>(180./degacc+1.);
367 
368  // AD performance is sensitive to calls to ensurecontiguous.
369  // // Providing "t" will cause ensurecontiguous to be called.
370  #ifdef _HAVE_AD_
371  G_rigid=xNew<IssmDouble>(M,"t");
372  G_elastic=xNew<IssmDouble>(M,"t");
373  U_elastic=xNew<IssmDouble>(M,"t");
374  H_elastic=xNew<IssmDouble>(M,"t");
375  #else
376  G_rigid=xNew<IssmDouble>(M);
377  G_elastic=xNew<IssmDouble>(M);
378  U_elastic=xNew<IssmDouble>(M);
379  H_elastic=xNew<IssmDouble>(M);
380  #endif
381 
382  /*compute combined legendre + love number (elastic green function:*/
383  m=DetermineLocalSize(M,IssmComm::GetComm());
384  GetOwnershipBoundariesFromRange(&lower_row,&upper_row,m,IssmComm::GetComm());
385  // AD performance is sensitive to calls to ensurecontiguous.
386  // // Providing "t" will cause ensurecontiguous to be called.
387  #ifdef _HAVE_AD_
388  G_elastic_local=xNew<IssmDouble>(m,"t");
389  G_rigid_local=xNew<IssmDouble>(m,"t");
390  U_elastic_local=xNew<IssmDouble>(m,"t");
391  H_elastic_local=xNew<IssmDouble>(m,"t");
392  #else
393  G_elastic_local=xNew<IssmDouble>(m);
394  G_rigid_local=xNew<IssmDouble>(m);
395  U_elastic_local=xNew<IssmDouble>(m);
396  H_elastic_local=xNew<IssmDouble>(m);
397  #endif
398 
399  for(int i=lower_row;i<upper_row;i++){
400  IssmDouble alpha,x;
401  alpha= reCast<IssmDouble>(i)*degacc * PI / 180.0;
402 
403  G_rigid_local[i-lower_row]= .5/sin(alpha/2.0);
404  G_elastic_local[i-lower_row]= (love_k[nl-1]-love_h[nl-1])*G_rigid_local[i-lower_row];
405  U_elastic_local[i-lower_row]= (love_h[nl-1])*G_rigid_local[i-lower_row];
406  H_elastic_local[i-lower_row]= 0;
407  IssmDouble Pn = 0.;
408  IssmDouble Pn1 = 0.;
409  IssmDouble Pn2 = 0.;
410  IssmDouble Pn_p = 0.;
411  IssmDouble Pn_p1 = 0.;
412  IssmDouble Pn_p2 = 0.;
413 
414  for (int n=0;n<nl;n++) {
415  IssmDouble deltalove_G;
416  IssmDouble deltalove_U;
417 
418  deltalove_G = (love_k[n]-love_k[nl-1]-love_h[n]+love_h[nl-1]);
419  deltalove_U = (love_h[n]-love_h[nl-1]);
420 
421  /*compute legendre polynomials: P_n(cos\theta) & d P_n(cos\theta)/ d\theta: */
422  if(n==0){
423  Pn=1;
424  Pn_p=0;
425  }
426  else if(n==1){
427  Pn = cos(alpha);
428  Pn_p = 1;
429  }
430  else{
431  Pn = ( (2*n-1)*cos(alpha)*Pn1 - (n-1)*Pn2 ) /n;
432  Pn_p = ( (2*n-1)*(Pn1+cos(alpha)*Pn_p1) - (n-1)*Pn_p2 ) /n;
433  }
434  Pn2=Pn1; Pn1=Pn;
435  Pn_p2=Pn_p1; Pn_p1=Pn_p;
436 
437  G_elastic_local[i-lower_row] += deltalove_G*Pn; // gravitational potential
438  U_elastic_local[i-lower_row] += deltalove_U*Pn; // vertical (up) displacement
439  H_elastic_local[i-lower_row] += sin(alpha)*love_l[n]*Pn_p; // horizontal displacements
440  }
441  }
442 
443  /*merge G_elastic_local into G_elastic; U_elastic_local into U_elastic; H_elastic_local to H_elastic:{{{*/
444  int* recvcounts=xNew<int>(IssmComm::GetSize());
445  int* displs=xNew<int>(IssmComm::GetSize());
446 
447  //recvcounts:
448  ISSM_MPI_Allgather(&m,1,ISSM_MPI_INT,recvcounts,1,ISSM_MPI_INT,IssmComm::GetComm());
449 
450  /*displs: */
451  ISSM_MPI_Allgather(&lower_row,1,ISSM_MPI_INT,displs,1,ISSM_MPI_INT,IssmComm::GetComm());
452 
453  /*All gather:*/
454  ISSM_MPI_Allgatherv(G_rigid_local, m, ISSM_MPI_DOUBLE, G_rigid, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
455  ISSM_MPI_Allgatherv(G_elastic_local, m, ISSM_MPI_DOUBLE, G_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
456  ISSM_MPI_Allgatherv(U_elastic_local, m, ISSM_MPI_DOUBLE, U_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
457  ISSM_MPI_Allgatherv(H_elastic_local, m, ISSM_MPI_DOUBLE, H_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
458  /*free ressources: */
459  xDelete<int>(recvcounts);
460  xDelete<int>(displs);
461 
462  /*}}}*/
463 
464  /*Avoid singularity at 0: */
465  G_rigid[0]=G_rigid[1];
466  parameters->AddObject(new DoubleVecParam(SealevelriseGRigidEnum,G_rigid,M));
467  G_elastic[0]=G_elastic[1];
468  parameters->AddObject(new DoubleVecParam(SealevelriseGElasticEnum,G_elastic,M));
469  U_elastic[0]=U_elastic[1];
470  parameters->AddObject(new DoubleVecParam(SealevelriseUElasticEnum,U_elastic,M));
471  H_elastic[0]=H_elastic[1];
472  parameters->AddObject(new DoubleVecParam(SealevelriseHElasticEnum,H_elastic,M));
473 
474  /*free ressources: */
475  xDelete<IssmDouble>(love_h);
476  xDelete<IssmDouble>(love_k);
477  xDelete<IssmDouble>(love_l);
478  xDelete<IssmDouble>(love_th);
479  xDelete<IssmDouble>(love_tk);
480  xDelete<IssmDouble>(love_tl);
481  xDelete<IssmDouble>(G_rigid);
482  xDelete<IssmDouble>(G_rigid_local);
483  xDelete<IssmDouble>(G_elastic);
484  xDelete<IssmDouble>(G_elastic_local);
485  xDelete<IssmDouble>(U_elastic);
486  xDelete<IssmDouble>(U_elastic_local);
487  xDelete<IssmDouble>(H_elastic);
488  xDelete<IssmDouble>(H_elastic_local);
489 
490  /*Indicate we have not yet ran the Geometry Core module: */
491  parameters->AddObject(new BoolParam(SealevelriseGeometryDoneEnum,false));
492 
493 /*}}}*/
494  /*Transitions:{{{ */
495  iomodel->FetchData(&transitions,&transitions_M,&transitions_N,&ntransitions,"md.solidearth.transitions");
496  if(transitions){
497  parameters->AddObject(new DoubleMatArrayParam(SealevelriseTransitionsEnum,transitions,ntransitions,transitions_M,transitions_N));
498 
499  for(int i=0;i<ntransitions;i++){
500  IssmDouble* transition=transitions[i];
501  xDelete<IssmDouble>(transition);
502  }
503  xDelete<IssmDouble*>(transitions);
504  xDelete<int>(transitions_M);
505  xDelete<int>(transitions_N);
506  } /*}}}*/
507  /*Requested outputs {{{*/
508  iomodel->FindConstant(&requestedoutputs,&numoutputs,"md.solidearth.requested_outputs");
509  if(numoutputs)parameters->AddObject(new StringArrayParam(SealevelriseRequestedOutputsEnum,requestedoutputs,numoutputs));
510  iomodel->DeleteData(&requestedoutputs,numoutputs,"md.solidearth.requested_outputs");
511  /*}}}*/
512 
513 }/*}}}*/
514 
515 /*Finite Element Analysis*/
516 void SealevelriseAnalysis::Core(FemModel* femmodel){/*{{{*/
517  _error_("not implemented");
518 }/*}}}*/
519 ElementVector* SealevelriseAnalysis::CreateDVector(Element* element){/*{{{*/
520  /*Default, return NULL*/
521  return NULL;
522 }/*}}}*/
523 ElementMatrix* SealevelriseAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/
524 _error_("Not implemented");
525 }/*}}}*/
526 ElementMatrix* SealevelriseAnalysis::CreateKMatrix(Element* element){/*{{{*/
527  _error_("not implemented yet");
528 }/*}}}*/
529 ElementVector* SealevelriseAnalysis::CreatePVector(Element* element){/*{{{*/
530 _error_("not implemented yet");
531 }/*}}}*/
532 void SealevelriseAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
533  _error_("not implemented yet");
534 }/*}}}*/
535 void SealevelriseAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
536  _error_("Not implemented yet");
537 }/*}}}*/
538 void SealevelriseAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
539  _error_("not implemeneted yet!");
540 
541 }/*}}}*/
542 void SealevelriseAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
543  /*Default, do nothing*/
544  return;
545 }/*}}}*/
DataSet::Size
int Size()
Definition: DataSet.cpp:399
DslSeaWaterPressureChangeAtSeaFloor
@ DslSeaWaterPressureChangeAtSeaFloor
Definition: EnumDefinitions.h:542
DslSeaSurfaceHeightChangeAboveGeoidEnum
@ DslSeaSurfaceHeightChangeAboveGeoidEnum
Definition: EnumDefinitions.h:541
DslGlobalAverageThermostericSeaLevelChangeEnum
@ DslGlobalAverageThermostericSeaLevelChangeEnum
Definition: EnumDefinitions.h:540
Element::lid
int lid
Definition: Element.h:46
SealevelUEsaRateEnum
@ SealevelUEsaRateEnum
Definition: EnumDefinitions.h:686
SealevelriseAnalysis::CreateDVector
ElementVector * CreateDVector(Element *element)
Definition: SealevelriseAnalysis.cpp:519
SolidearthSettingsAbstolEnum
@ SolidearthSettingsAbstolEnum
Definition: EnumDefinitions.h:306
SealevelEustaticOceanMaskEnum
@ SealevelEustaticOceanMaskEnum
Definition: EnumDefinitions.h:677
RotationalPolarMoiEnum
@ RotationalPolarMoiEnum
Definition: EnumDefinitions.h:326
_assert_
#define _assert_(ignore)
Definition: exceptions.h:37
TidalLoveK2SecularEnum
@ TidalLoveK2SecularEnum
Definition: EnumDefinitions.h:314
IssmDouble
double IssmDouble
Definition: types.h:37
Nodes
Declaration of Nodes class.
Definition: Nodes.h:19
SealevelriseAnalysis::UpdateElements
void UpdateElements(Elements *elements, Inputs2 *inputs2, IoModel *iomodel, int analysis_counter, int analysis_type)
Definition: SealevelriseAnalysis.cpp:22
SealevelriseAnalysis::CreateLoads
void CreateLoads(Loads *loads, IoModel *iomodel)
Definition: SealevelriseAnalysis.cpp:13
DslModelidEnum
@ DslModelidEnum
Definition: EnumDefinitions.h:126
SealevelriseAnalysis::DofsPerNode
int DofsPerNode(int **doflist, int domaintype, int approximation)
Definition: SealevelriseAnalysis.cpp:19
SealevelriseHElasticEnum
@ SealevelriseHElasticEnum
Definition: EnumDefinitions.h:322
RotationalAngularVelocityEnum
@ RotationalAngularVelocityEnum
Definition: EnumDefinitions.h:307
MaskOceanLevelsetEnum
@ MaskOceanLevelsetEnum
Definition: EnumDefinitions.h:640
Parameters
Declaration of Parameters class.
Definition: Parameters.h:18
MaskIceLevelsetEnum
@ MaskIceLevelsetEnum
Definition: EnumDefinitions.h:641
IssmComm::GetComm
static ISSM_MPI_Comm GetComm(void)
Definition: IssmComm.cpp:30
Constraints
Declaration of Constraints class.
Definition: Constraints.h:13
IoModel::my_vertices_lids
int * my_vertices_lids
Definition: IoModel.h:73
BedEnum
@ BedEnum
Definition: EnumDefinitions.h:499
P0Enum
@ P0Enum
Definition: EnumDefinitions.h:661
TransientInput2
Definition: TransientInput2.h:13
Elements
Declaration of Elements class.
Definition: Elements.h:17
GetOwnershipBoundariesFromRange
void GetOwnershipBoundariesFromRange(int *plower_row, int *pupper_row, int range, ISSM_MPI_Comm comm)
Definition: GetOwnershipBoundariesFromRange.cpp:16
SealevelEustaticMaskEnum
@ SealevelEustaticMaskEnum
Definition: EnumDefinitions.h:676
SealevelRSLRateEnum
@ SealevelRSLRateEnum
Definition: EnumDefinitions.h:683
IoModel::my_elements
bool * my_elements
Definition: IoModel.h:66
LoadLoveHEnum
@ LoadLoveHEnum
Definition: EnumDefinitions.h:315
SolidearthSettingsRigidEnum
@ SolidearthSettingsRigidEnum
Definition: EnumDefinitions.h:329
Element::Sid
int Sid()
Definition: Element.cpp:3578
SolidearthSettingsHorizEnum
@ SolidearthSettingsHorizEnum
Definition: EnumDefinitions.h:323
TidalLoveKEnum
@ TidalLoveKEnum
Definition: EnumDefinitions.h:312
Parameters::AddObject
void AddObject(Param *newparam)
Definition: Parameters.cpp:67
SolidearthSettingsOceanAreaScalingEnum
@ SolidearthSettingsOceanAreaScalingEnum
Definition: EnumDefinitions.h:325
SealevelriseTransitionsEnum
@ SealevelriseTransitionsEnum
Definition: EnumDefinitions.h:332
DetermineLocalSize
int DetermineLocalSize(int global_size, ISSM_MPI_Comm comm)
Definition: DetermineLocalSize.cpp:9
Element
Definition: Element.h:41
P1Enum
@ P1Enum
Definition: EnumDefinitions.h:662
DslModelEnum
@ DslModelEnum
Definition: EnumDefinitions.h:125
SolidearthSettingsReltolEnum
@ SolidearthSettingsReltolEnum
Definition: EnumDefinitions.h:327
SealevelriseCumDeltathicknessEnum
@ SealevelriseCumDeltathicknessEnum
Definition: EnumDefinitions.h:688
ISSM_MPI_DOUBLE
#define ISSM_MPI_DOUBLE
Definition: issmmpi.h:125
DslComputeFingerprintsEnum
@ DslComputeFingerprintsEnum
Definition: EnumDefinitions.h:128
ISSM_MPI_INT
#define ISSM_MPI_INT
Definition: issmmpi.h:127
IoModel::DeleteData
void DeleteData(int num,...)
Definition: IoModel.cpp:500
SolidearthSettingsMaxiterEnum
@ SolidearthSettingsMaxiterEnum
Definition: EnumDefinitions.h:324
SealevelriseAnalysis::CreateConstraints
void CreateConstraints(Constraints *constraints, IoModel *iomodel)
Definition: SealevelriseAnalysis.cpp:10
IoModel::numberofelements
int numberofelements
Definition: IoModel.h:96
SolidearthPlanetRadiusEnum
@ SolidearthPlanetRadiusEnum
Definition: EnumDefinitions.h:303
DoubleParam
Definition: DoubleParam.h:20
IoModel::CopyConstantObject
Param * CopyConstantObject(const char *constant_name, int param_enum)
Definition: IoModel.cpp:418
IssmComm::GetSize
static int GetSize(void)
Definition: IssmComm.cpp:46
SealevelriseUElasticEnum
@ SealevelriseUElasticEnum
Definition: EnumDefinitions.h:333
alpha
IssmDouble alpha(IssmDouble x, IssmDouble y, IssmDouble z, int testid)
Definition: fsanalyticals.cpp:221
IoModel::FindConstant
void FindConstant(bool *pvalue, const char *constant_name)
Definition: IoModel.cpp:2362
SealevelriseAnalysis::UpdateConstraints
void UpdateConstraints(FemModel *femmodel)
Definition: SealevelriseAnalysis.cpp:542
TidalLoveHEnum
@ TidalLoveHEnum
Definition: EnumDefinitions.h:311
StringArrayParam
Definition: StringArrayParam.h:20
SealevelriseAnalysis::Core
void Core(FemModel *femmodel)
Definition: SealevelriseAnalysis.cpp:516
SurfaceloadWaterHeightChangeEnum
@ SurfaceloadWaterHeightChangeEnum
Definition: EnumDefinitions.h:692
SealevelriseAnalysis.h
: header file for generic external result object
SolidearthPlanetAreaEnum
@ SolidearthPlanetAreaEnum
Definition: EnumDefinitions.h:304
IoModel::FetchData
void FetchData(bool *pboolean, const char *data_name)
Definition: IoModel.cpp:933
Inputs2
Declaration of Inputs class.
Definition: Inputs2.h:23
SealevelriseAnalysisEnum
@ SealevelriseAnalysisEnum
Definition: EnumDefinitions.h:1266
InputUpdateFromConstantx
void InputUpdateFromConstantx(FemModel *femmodel, bool constant, int name)
Definition: InputUpdateFromConstantx.cpp:10
SealevelriseGeometryDoneEnum
@ SealevelriseGeometryDoneEnum
Definition: EnumDefinitions.h:309
PI
const double PI
Definition: constants.h:11
SealevelriseAnalysis::UpdateParameters
void UpdateParameters(Parameters *parameters, IoModel *iomodel, int solution_enum, int analysis_enum)
Definition: SealevelriseAnalysis.cpp:276
TransientInput2::AddPentaTimeInput
void AddPentaTimeInput(IssmDouble time, int numindices, int *indices, IssmDouble *values_in, int interp_in)
Definition: TransientInput2.cpp:180
Object::ObjectEnum
virtual int ObjectEnum()=0
FemModel
Definition: FemModel.h:31
SealevelEnum
@ SealevelEnum
Definition: EnumDefinitions.h:675
ISSM_MPI_Allgather
int ISSM_MPI_Allgather(void *sendbuf, int sendcount, ISSM_MPI_Datatype sendtype, void *recvbuf, int recvcount, ISSM_MPI_Datatype recvtype, ISSM_MPI_Comm comm)
Definition: issmmpi.cpp:45
Inputs2::GetTransientInput
TransientInput2 * GetTransientInput(int enum_type)
Definition: Inputs2.cpp:406
Loads
Declaration of Loads class.
Definition: Loads.h:16
SealevelriseRequestedOutputsEnum
@ SealevelriseRequestedOutputsEnum
Definition: EnumDefinitions.h:328
SealevelriseAnalysis::InputUpdateFromSolution
void InputUpdateFromSolution(IssmDouble *solution, Element *element)
Definition: SealevelriseAnalysis.cpp:538
_error_
#define _error_(StreamArgs)
Definition: exceptions.h:49
SealevelriseGElasticEnum
@ SealevelriseGElasticEnum
Definition: EnumDefinitions.h:319
SolidearthSettingsRotationEnum
@ SolidearthSettingsRotationEnum
Definition: EnumDefinitions.h:330
DataSet::GetObjectByOffset
Object * GetObjectByOffset(int offset)
Definition: DataSet.cpp:334
SurfaceloadIceThicknessChangeEnum
@ SurfaceloadIceThicknessChangeEnum
Definition: EnumDefinitions.h:691
SealevelriseAnalysis::CreateNodes
void CreateNodes(Nodes *nodes, IoModel *iomodel, bool isamr=false)
Definition: SealevelriseAnalysis.cpp:16
DoubleVecParam
Definition: DoubleVecParam.h:20
PentaEnum
@ PentaEnum
Definition: EnumDefinitions.h:1231
SolidearthSettingsRunFrequencyEnum
@ SolidearthSettingsRunFrequencyEnum
Definition: EnumDefinitions.h:321
Element::Update
virtual void Update(Inputs2 *inputs2, int index, IoModel *iomodel, int analysis_counter, int analysis_type, int finite_element)=0
SealevelriseAnalysis::CreateKMatrix
ElementMatrix * CreateKMatrix(Element *element)
Definition: SealevelriseAnalysis.cpp:526
ISSM_MPI_Allgatherv
int ISSM_MPI_Allgatherv(void *sendbuf, int sendcount, ISSM_MPI_Datatype sendtype, void *recvbuf, int *recvcounts, int *displs, ISSM_MPI_Datatype recvtype, ISSM_MPI_Comm comm)
Definition: issmmpi.cpp:79
SealevelriseAnalysis::GetSolutionFromInputs
void GetSolutionFromInputs(Vector< IssmDouble > *solution, Element *element)
Definition: SealevelriseAnalysis.cpp:532
Inputs2::SetTransientInput
void SetTransientInput(int enum_in, IssmDouble *times, int numtimes)
Definition: Inputs2.cpp:676
IoModel::FetchDataToInput
void FetchDataToInput(Inputs2 *inputs2, Elements *elements, const char *vector_name, int input_enum)
Definition: IoModel.cpp:1651
TidalLoveLEnum
@ TidalLoveLEnum
Definition: EnumDefinitions.h:313
LoadLoveKEnum
@ LoadLoveKEnum
Definition: EnumDefinitions.h:316
LoadLoveLEnum
@ LoadLoveLEnum
Definition: EnumDefinitions.h:317
DoubleMatParam
Definition: DoubleMatParam.h:20
ElementVector
Definition: ElementVector.h:20
IoModel::elements
int * elements
Definition: IoModel.h:79
IoModel
Definition: IoModel.h:48
BoolParam
Definition: BoolParam.h:20
Element::GetNumberOfVertices
virtual int GetNumberOfVertices(void)=0
SealevelriseGRigidEnum
@ SealevelriseGRigidEnum
Definition: EnumDefinitions.h:318
SolidearthSettingsElasticEnum
@ SolidearthSettingsElasticEnum
Definition: EnumDefinitions.h:308
ElementMatrix
Definition: ElementMatrix.h:19
Vector< IssmDouble >
SealevelNEsaRateEnum
@ SealevelNEsaRateEnum
Definition: EnumDefinitions.h:679
SealevelriseAnalysis::CreateJacobianMatrix
ElementMatrix * CreateJacobianMatrix(Element *element)
Definition: SealevelriseAnalysis.cpp:523
DoubleMatArrayParam
Definition: DoubleMatArrayParam.h:20
TriaEnum
@ TriaEnum
Definition: EnumDefinitions.h:1318
RotationalEquatorialMoiEnum
@ RotationalEquatorialMoiEnum
Definition: EnumDefinitions.h:310
SolidearthSettingsComputesealevelchangeEnum
@ SolidearthSettingsComputesealevelchangeEnum
Definition: EnumDefinitions.h:320
TransientInput2::AddTriaTimeInput
void AddTriaTimeInput(IssmDouble time, int numindices, int *indices, IssmDouble *values_in, int interp_in)
Definition: TransientInput2.cpp:138
SealevelriseAnalysis::GradientJ
void GradientJ(Vector< IssmDouble > *gradient, Element *element, int control_type, int control_index)
Definition: SealevelriseAnalysis.cpp:535
SealevelriseAnalysis::CreatePVector
ElementVector * CreatePVector(Element *element)
Definition: SealevelriseAnalysis.cpp:529
DslNummodelsEnum
@ DslNummodelsEnum
Definition: EnumDefinitions.h:127
Inputs2::SetDatasetTransientInput
TransientInput2 * SetDatasetTransientInput(int enum_in, int id, IssmDouble *times, int numtimes)
Definition: Inputs2.cpp:649
femmodel
FemModel * femmodel
Definition: esmfbinders.cpp:16
Generated on Thu Jul 2 2020 08:09:15 for Ice Sheet System Model by   doxygen 1.8.19