source: issm/trunk-jpl/src/c/analyses/SealevelriseAnalysis.cpp@ 24940

Last change on this file since 24940 was 24940, checked in by Eric.Larour, 5 years ago

CHG: replacing earth area computations with a constant coming from the
slr class.
File size: 19.7 KB
RevLine 
[19984]1#include "./SealevelriseAnalysis.h"
2#include "../toolkits/toolkits.h"
3#include "../classes/classes.h"
[24481]4#include "../classes/Inputs2/TransientInput2.h"
[19984]5#include "../shared/shared.h"
6#include "../modules/modules.h"
7
8/*Model processing*/
9void SealevelriseAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
10 /*No constraints*/
11}/*}}}*/
12void SealevelriseAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
13 /*No loads*/
14}/*}}}*/
[23585]15void SealevelriseAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel,bool isamr){/*{{{*/
[19984]16 ::CreateNodes(nodes,iomodel,SealevelriseAnalysisEnum,P1Enum);
17}/*}}}*/
18int SealevelriseAnalysis::DofsPerNode(int** doflist,int domaintype,int approximation){/*{{{*/
19 return 1;
20}/*}}}*/
[24335]21void SealevelriseAnalysis::UpdateElements(Elements* elements,Inputs2* inputs2,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
[19984]22
[22955]23 int geodetic=0;
[24469]24 int dslmodel=0;
[22955]25
[19984]26 /*Update elements: */
27 int counter=0;
28 for(int i=0;i<iomodel->numberofelements;i++){
29 if(iomodel->my_elements[i]){
30 Element* element=(Element*)elements->GetObjectByOffset(counter);
[24335]31 element->Update(inputs2,i,iomodel,analysis_counter,analysis_type,P1Enum);
[19984]32 counter++;
33 }
34 }
35
36 /*Create inputs: */
[24861]37 iomodel->FetchDataToInput(inputs2,elements,"md.mask.ocean_levelset",MaskOceanLevelsetEnum);
[24335]38 iomodel->FetchDataToInput(inputs2,elements,"md.mask.ice_levelset",MaskIceLevelsetEnum);
[22955]39 iomodel->FetchData(&geodetic,"md.slr.geodetic");
[24335]40 iomodel->FetchDataToInput(inputs2,elements,"md.slr.deltathickness",SealevelriseDeltathicknessEnum);
41 iomodel->FetchDataToInput(inputs2,elements,"md.slr.spcthickness",SealevelriseSpcthicknessEnum);
42 iomodel->FetchDataToInput(inputs2,elements,"md.slr.sealevel",SealevelEnum,0);
43 iomodel->FetchDataToInput(inputs2,elements,"md.geometry.bed",BedEnum);
44 iomodel->FetchDataToInput(inputs2,elements,"md.slr.Ngia",SealevelNGiaRateEnum);
45 iomodel->FetchDataToInput(inputs2,elements,"md.slr.Ugia",SealevelUGiaRateEnum);
46 iomodel->FetchDataToInput(inputs2,elements,"md.slr.hydro_rate",SealevelriseHydroRateEnum);
[24469]47
48 /*dynamic sea level: */
49 iomodel->FetchData(&dslmodel,"md.dsl.model");
[24505]50 if (dslmodel==1){ /*standard dsl model:{{{*/
[24481]51
52 /*deal with global mean steric rate: */
53 IssmDouble* str=NULL;
54 IssmDouble* times = NULL;
55 int M,N;
56
57 /*fetch str vector:*/
58 iomodel->FetchData(&str,&M,&N,"md.dsl.global_average_thermosteric_sea_level_change"); _assert_(M==2);
59
60 //recover time vector:
61 times=xNew<IssmDouble>(N);
62 for(int t=0;t<N;t++) times[t] = str[N+t];
63
64 /*create transient input: */
65 inputs2->SetTransientInput(DslGlobalAverageThermostericSeaLevelChangeEnum,times,N);
66 TransientInput2* transientinput = inputs2->GetTransientInput(DslGlobalAverageThermostericSeaLevelChangeEnum);
67
68
69 for(int i=0;i<elements->Size();i++){
70 Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
71
72 for(int t=0;t<N;t++){
73 switch(element->ObjectEnum()){
74 case TriaEnum: transientinput->AddTriaTimeInput( t,1,&element->lid,&str[t],P0Enum); break;
75 case PentaEnum: transientinput->AddPentaTimeInput(t,1,&element->lid,&str[t],P0Enum); break;
76 default: _error_("Not implemented yet");
77 }
78 }
79 }
80
81 /*cleanup:*/
82 xDelete<IssmDouble>(times);
83 iomodel->DeleteData(str,"md.dsl.global_average_thermosteric_sea_level_change");
84
85 /*deal with dynamic sea level fields: */
[24469]86 iomodel->FetchDataToInput(inputs2,elements,"md.dsl.sea_surface_height_change_above_geoid", DslSeaSurfaceHeightChangeAboveGeoidEnum);
87 iomodel->FetchDataToInput(inputs2,elements,"md.dsl.sea_water_pressure_change_at_sea_floor", DslSeaWaterPressureChangeAtSeaFloor);
[24481]88
[24505]89 } /*}}}*/
90 else if (dslmodel==2){ /*multi-model ensemble dsl model:{{{*/
91
92 /*variables:*/
93 int nummodels;
94 IssmDouble** pstr=NULL;
95 IssmDouble* str=NULL;
96 IssmDouble* times = NULL;
97 int* pM = NULL;
98 int* pN = NULL;
99 int M,N;
100
101 /*deal with dsl.sea_surface_height_change_above_geoid {{{*/
102 iomodel->FetchData(&pstr,&pM,&pN,&nummodels,"md.dsl.global_average_thermosteric_sea_level_change");
103
104 /*go through the mat array and create a dataset of transient inputs:*/
105 for (int i=0;i<nummodels;i++){
106
107 M=pM[i];
108 N=pN[i];
109 str=pstr[i];
110
111 //recover time vector:
112 times=xNew<IssmDouble>(N);
113 for(int t=0;t<N;t++) times[t] = str[(M-1)*N+t];
114
115 TransientInput2* transientinput=inputs2->SetDatasetTransientInput(DslGlobalAverageThermostericSeaLevelChangeEnum,i, times,N);
116
117 for(int j=0;j<elements->Size();j++){
118 Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(j));
119
120 for(int t=0;t<N;t++){
121 switch(element->ObjectEnum()){
122 case TriaEnum: transientinput->AddTriaTimeInput( t,1,&element->lid,&str[t],P0Enum); break;
123 case PentaEnum: transientinput->AddPentaTimeInput(t,1,&element->lid,&str[t],P0Enum); break;
124 default: _error_("Not implemented yet");
125 }
126 }
127 }
128 xDelete<IssmDouble>(times);
129 }
130 /*Delete data:*/
131 for(int i=0;i<nummodels;i++){
132 IssmDouble* str=pstr[i];
133 xDelete<IssmDouble>(str);
134 }
135 xDelete<IssmDouble*>(pstr);
136 xDelete<int>(pM);
137 xDelete<int>(pN);
138 /*}}}*/
139 /*now do the same with the dsl.sea_surface_height_change_above_geoid field:{{{ */
140 iomodel->FetchData(&pstr,&pM,&pN,&nummodels,"md.dsl.sea_surface_height_change_above_geoid");
141
142 for (int i=0;i<nummodels;i++){
143 M=pM[i];
144 N=pN[i];
145 str=pstr[i];
146
147
148 //recover time vector:
149 times=xNew<IssmDouble>(N);
150 for(int t=0;t<N;t++) times[t] = str[(M-1)*N+t];
151
152 TransientInput2* transientinput=inputs2->SetDatasetTransientInput(DslSeaSurfaceHeightChangeAboveGeoidEnum,i, times,N);
153
154 for(int j=0;j<elements->Size();j++){
155
156 /*Get the right transient input*/
157 Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(j));
158
159 /*Get values and lid list*/
160 const int numvertices = element->GetNumberOfVertices();
161 int *vertexlids = xNew<int>(numvertices);
162 int *vertexsids = xNew<int>(numvertices);
163
164
165 /*Recover vertices ids needed to initialize inputs*/
166 _assert_(iomodel->elements);
167 for(int k=0;k<numvertices;k++){
168 vertexsids[k] =reCast<int>(iomodel->elements[numvertices*element->Sid()+k]); //ids for vertices are in the elements array from Matlab
169 vertexlids[k]=iomodel->my_vertices_lids[vertexsids[k]-1];
170 }
171
172 //element->GetVerticesLidList(vertexlids);
173 //element->GetVerticesSidList(vertexsids);
174 IssmDouble* values=xNew<IssmDouble>(numvertices);
175
176 for(int t=0;t<N;t++){
177 for (int k=0;k<numvertices;k++)values[k]=str[N*vertexsids[k]+t];
178
179 switch(element->ObjectEnum()){
180 case TriaEnum: transientinput->AddTriaTimeInput( t,numvertices,vertexlids,values,P1Enum); break;
181 case PentaEnum: transientinput->AddPentaTimeInput(t,numvertices,vertexlids,values,P1Enum); break;
182 default: _error_("Not implemented yet");
183 }
184 }
185 xDelete<IssmDouble>(values);
186 xDelete<int>(vertexlids);
187 xDelete<int>(vertexsids);
188 }
189
190 xDelete<IssmDouble>(times);
191 }
192
193 /*Delete data:*/
194 for(int i=0;i<nummodels;i++){
195 IssmDouble* str=pstr[i];
196 xDelete<IssmDouble>(str);
197 }
198 xDelete<IssmDouble*>(pstr);
199 xDelete<int>(pM);
200 xDelete<int>(pN);
201 /*}}}*/
202 /*now do the same with the dsl.sea_water_pressure_change_at_sea_floor field:{{{ */
203 iomodel->FetchData(&pstr,&pM,&pN,&nummodels,"md.dsl.sea_water_pressure_change_at_sea_floor");
204
205 for (int i=0;i<nummodels;i++){
206 M=pM[i];
207 N=pN[i];
208 str=pstr[i];
209
210 //recover time vector:
211 times=xNew<IssmDouble>(N);
212 for(int t=0;t<N;t++) times[t] = str[(M-1)*N+t];
213
214 TransientInput2* transientinput=inputs2->SetDatasetTransientInput(DslSeaWaterPressureChangeAtSeaFloor,i, times,N);
215
216 for(int j=0;j<elements->Size();j++){
217
218 /*Get the right transient input*/
219 Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(j));
220
221 /*Get values and lid list*/
222 const int numvertices = element->GetNumberOfVertices();
223 int *vertexlids = xNew<int>(numvertices);
224 int *vertexsids = xNew<int>(numvertices);
225
226
227 /*Recover vertices ids needed to initialize inputs*/
228 _assert_(iomodel->elements);
229 for(int k=0;k<numvertices;k++){
230 vertexsids[k] =reCast<int>(iomodel->elements[numvertices*element->Sid()+k]); //ids for vertices are in the elements array from Matlab
231 vertexlids[k]=iomodel->my_vertices_lids[vertexsids[k]-1];
232 }
233 //element->GetVerticesLidList(vertexlids);
234 //element->GetVerticesSidList(vertexsids);
235
236 IssmDouble* values=xNew<IssmDouble>(numvertices);
237
238 for(int t=0;t<N;t++){
239 for (int k=0;k<numvertices;k++)values[k]=str[N*vertexsids[k]+t];
240
241 switch(element->ObjectEnum()){
242 case TriaEnum: transientinput->AddTriaTimeInput( t,numvertices,vertexlids,values,P1Enum); break;
243 case PentaEnum: transientinput->AddPentaTimeInput(t,numvertices,vertexlids,values,P1Enum); break;
244 default: _error_("Not implemented yet");
245 }
246 }
247 xDelete<IssmDouble>(values);
248 xDelete<int>(vertexlids);
249 xDelete<int>(vertexsids);
250 }
251 xDelete<IssmDouble>(times);
252 }
253
254 /*Delete data:*/
255 for(int i=0;i<nummodels;i++){
256 IssmDouble* str=pstr[i];
257 xDelete<IssmDouble>(str);
258 }
259 xDelete<IssmDouble*>(pstr);
260 xDelete<int>(pM);
261 xDelete<int>(pN);
262 /*}}}*/
263
264 } /*}}}*/
[24469]265 else _error_("Dsl model " << dslmodel << " not implemented yet!");
266
[22968]267 /*Initialize cumdeltalthickness and sealevel rise rate input*/
[24335]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);
[22968]274
[19984]275}/*}}}*/
276void SealevelriseAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
277
[19986]278 int nl;
[19984]279 IssmDouble* love_h=NULL;
280 IssmDouble* love_k=NULL;
[20999]281 IssmDouble* love_l=NULL;
[24469]282 int dslmodel=0;
[23066]283
[24921]284 IssmDouble* G_rigid = NULL;
285 IssmDouble* G_rigid_local = NULL;
[20028]286 IssmDouble* G_elastic = NULL;
[20033]287 IssmDouble* G_elastic_local = NULL;
[20999]288 IssmDouble* U_elastic = NULL;
289 IssmDouble* U_elastic_local = NULL;
290 IssmDouble* H_elastic = NULL;
291 IssmDouble* H_elastic_local = NULL;
[20033]292 int M,m,lower_row,upper_row;
293 IssmDouble degacc=.01;
[19984]294
[20036]295 int numoutputs;
296 char** requestedoutputs = NULL;
297
[20136]298 /*transition vectors: */
299 IssmDouble **transitions = NULL;
300 int *transitions_M = NULL;
301 int *transitions_N = NULL;
302 int ntransitions;
303
[19986]304 /*some constant parameters: */
[24531]305 parameters->AddObject(iomodel->CopyConstantObject("md.dsl.model",DslModelEnum));
[22961]306 parameters->AddObject(iomodel->CopyConstantObject("md.slr.geodetic_run_frequency",SealevelriseGeodeticRunFrequencyEnum));
[20690]307 parameters->AddObject(iomodel->CopyConstantObject("md.slr.reltol",SealevelriseReltolEnum));
308 parameters->AddObject(iomodel->CopyConstantObject("md.slr.abstol",SealevelriseAbstolEnum));
309 parameters->AddObject(iomodel->CopyConstantObject("md.slr.maxiter",SealevelriseMaxiterEnum));
[22955]310 parameters->AddObject(iomodel->CopyConstantObject("md.slr.loop_increment",SealevelriseLoopIncrementEnum));
[20690]311 parameters->AddObject(iomodel->CopyConstantObject("md.slr.rigid",SealevelriseRigidEnum));
[22955]312 parameters->AddObject(iomodel->CopyConstantObject("md.slr.horiz",SealevelriseHorizEnum));
[20690]313 parameters->AddObject(iomodel->CopyConstantObject("md.slr.elastic",SealevelriseElasticEnum));
[21331]314 parameters->AddObject(iomodel->CopyConstantObject("md.slr.rotation",SealevelriseRotationEnum));
315 parameters->AddObject(iomodel->CopyConstantObject("md.slr.tide_love_h",SealevelriseTidalLoveHEnum));
316 parameters->AddObject(iomodel->CopyConstantObject("md.slr.tide_love_k",SealevelriseTidalLoveKEnum));
[21344]317 parameters->AddObject(iomodel->CopyConstantObject("md.slr.fluid_love",SealevelriseFluidLoveEnum));
318 parameters->AddObject(iomodel->CopyConstantObject("md.slr.equatorial_moi",SealevelriseEquatorialMoiEnum));
319 parameters->AddObject(iomodel->CopyConstantObject("md.slr.polar_moi",SealevelrisePolarMoiEnum));
[21331]320 parameters->AddObject(iomodel->CopyConstantObject("md.slr.angular_velocity",SealevelriseAngularVelocityEnum));
[21295]321 parameters->AddObject(iomodel->CopyConstantObject("md.slr.ocean_area_scaling",SealevelriseOceanAreaScalingEnum));
[22955]322 parameters->AddObject(iomodel->CopyConstantObject("md.slr.geodetic",SealevelriseGeodeticEnum));
[24940]323 parameters->AddObject(iomodel->CopyConstantObject("md.slr.eartharea",SealevelEarthAreaEnum));
[19986]324
[24505]325 /*Deal with dsl multi-model ensembles: {{{*/
326 iomodel->FetchData(&dslmodel,"md.dsl.model");
327 parameters->AddObject(iomodel->CopyConstantObject("md.dsl.compute_fingerprints",DslComputeFingerprintsEnum));
328 if(dslmodel==2){
329 int modelid;
330 int nummodels;
331
332 parameters->AddObject(iomodel->CopyConstantObject("md.dsl.modelid",DslModelidEnum));
333 parameters->AddObject(iomodel->CopyConstantObject("md.dsl.nummodels",DslNummodelsEnum));
334 iomodel->FetchData(&modelid,"md.dsl.modelid");
335 iomodel->FetchData(&nummodels,"md.dsl.nummodels");
336
337 /*quick checks: */
338 if(nummodels<=0)_error_("dslmme object in md.dsl field should contain at least 1 ensemble model!");
339 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!");
340 } /*}}}*/
[24481]341 /*Deal with elasticity {{{*/
[24921]342 /*love numbers: */
343 iomodel->FetchData(&love_h,&nl,NULL,"md.slr.love_h");
344 iomodel->FetchData(&love_k,&nl,NULL,"md.slr.love_k");
345 iomodel->FetchData(&love_l,&nl,NULL,"md.slr.love_l");
[19984]346
[24921]347 /*compute elastic green function for a range of angles*/
348 iomodel->FetchData(&degacc,"md.slr.degacc");
349 M=reCast<int,IssmDouble>(180./degacc+1.);
[20022]350
[24921]351 // AD performance is sensitive to calls to ensurecontiguous.
352 // // Providing "t" will cause ensurecontiguous to be called.
353 #ifdef _HAVE_AD_
354 G_rigid=xNew<IssmDouble>(M,"t");
355 G_elastic=xNew<IssmDouble>(M,"t");
356 U_elastic=xNew<IssmDouble>(M,"t");
357 H_elastic=xNew<IssmDouble>(M,"t");
358 #else
359 G_rigid=xNew<IssmDouble>(M);
360 G_elastic=xNew<IssmDouble>(M);
361 U_elastic=xNew<IssmDouble>(M);
362 H_elastic=xNew<IssmDouble>(M);
363 #endif
[21742]364
[24921]365 /*compute combined legendre + love number (elastic green function:*/
366 m=DetermineLocalSize(M,IssmComm::GetComm());
367 GetOwnershipBoundariesFromRange(&lower_row,&upper_row,m,IssmComm::GetComm());
368 // AD performance is sensitive to calls to ensurecontiguous.
369 // // Providing "t" will cause ensurecontiguous to be called.
370 #ifdef _HAVE_AD_
371 G_elastic_local=xNew<IssmDouble>(m,"t");
372 G_rigid_local=xNew<IssmDouble>(m,"t");
373 U_elastic_local=xNew<IssmDouble>(m,"t");
374 H_elastic_local=xNew<IssmDouble>(m,"t");
375 #else
376 G_elastic_local=xNew<IssmDouble>(m);
377 G_rigid_local=xNew<IssmDouble>(m);
378 U_elastic_local=xNew<IssmDouble>(m);
379 H_elastic_local=xNew<IssmDouble>(m);
380 #endif
[21742]381
[24921]382 for(int i=lower_row;i<upper_row;i++){
383 IssmDouble alpha,x;
384 alpha= reCast<IssmDouble>(i)*degacc * PI / 180.0;
[20033]385
[24921]386 G_rigid_local[i-lower_row]= .5/sin(alpha/2.0);
387 G_elastic_local[i-lower_row]= (love_k[nl-1]-love_h[nl-1])*G_rigid_local[i-lower_row];
388 U_elastic_local[i-lower_row]= (love_h[nl-1])*G_rigid_local[i-lower_row];
389 H_elastic_local[i-lower_row]= 0;
390 IssmDouble Pn = 0.;
391 IssmDouble Pn1 = 0.;
392 IssmDouble Pn2 = 0.;
393 IssmDouble Pn_p = 0.;
394 IssmDouble Pn_p1 = 0.;
395 IssmDouble Pn_p2 = 0.;
[20029]396
[24921]397 for (int n=0;n<nl;n++) {
398 IssmDouble deltalove_G;
399 IssmDouble deltalove_U;
[21908]400
[24921]401 deltalove_G = (love_k[n]-love_k[nl-1]-love_h[n]+love_h[nl-1]);
402 deltalove_U = (love_h[n]-love_h[nl-1]);
[20033]403
[24921]404 /*compute legendre polynomials: P_n(cos\theta) & d P_n(cos\theta)/ d\theta: */
405 if(n==0){
406 Pn=1;
407 Pn_p=0;
408 }
409 else if(n==1){
410 Pn = cos(alpha);
411 Pn_p = 1;
412 }
413 else{
414 Pn = ( (2*n-1)*cos(alpha)*Pn1 - (n-1)*Pn2 ) /n;
415 Pn_p = ( (2*n-1)*(Pn1+cos(alpha)*Pn_p1) - (n-1)*Pn_p2 ) /n;
416 }
417 Pn2=Pn1; Pn1=Pn;
418 Pn_p2=Pn_p1; Pn_p1=Pn_p;
[23066]419
[24921]420 G_elastic_local[i-lower_row] += deltalove_G*Pn; // gravitational potential
421 U_elastic_local[i-lower_row] += deltalove_U*Pn; // vertical (up) displacement
422 H_elastic_local[i-lower_row] += sin(alpha)*love_l[n]*Pn_p; // horizontal displacements
[20033]423 }
[24921]424 }
[20033]425
[24921]426 /*merge G_elastic_local into G_elastic; U_elastic_local into U_elastic; H_elastic_local to H_elastic:{{{*/
427 int* recvcounts=xNew<int>(IssmComm::GetSize());
428 int* displs=xNew<int>(IssmComm::GetSize());
[20033]429
[24921]430 //recvcounts:
431 ISSM_MPI_Allgather(&m,1,ISSM_MPI_INT,recvcounts,1,ISSM_MPI_INT,IssmComm::GetComm());
[20033]432
[24921]433 /*displs: */
434 ISSM_MPI_Allgather(&lower_row,1,ISSM_MPI_INT,displs,1,ISSM_MPI_INT,IssmComm::GetComm());
[20033]435
[24921]436 /*All gather:*/
437 ISSM_MPI_Allgatherv(G_rigid_local, m, ISSM_MPI_DOUBLE, G_rigid, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
438 ISSM_MPI_Allgatherv(G_elastic_local, m, ISSM_MPI_DOUBLE, G_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
439 ISSM_MPI_Allgatherv(U_elastic_local, m, ISSM_MPI_DOUBLE, U_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
440 ISSM_MPI_Allgatherv(H_elastic_local, m, ISSM_MPI_DOUBLE, H_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
441 /*free ressources: */
442 xDelete<int>(recvcounts);
443 xDelete<int>(displs);
[20033]444
[24921]445 /*}}}*/
[20033]446
[24921]447 /*Avoid singularity at 0: */
448 G_rigid[0]=G_rigid[1];
449 parameters->AddObject(new DoubleVecParam(SealevelriseGRigidEnum,G_rigid,M));
450 G_elastic[0]=G_elastic[1];
451 parameters->AddObject(new DoubleVecParam(SealevelriseGElasticEnum,G_elastic,M));
452 U_elastic[0]=U_elastic[1];
453 parameters->AddObject(new DoubleVecParam(SealevelriseUElasticEnum,U_elastic,M));
454 H_elastic[0]=H_elastic[1];
455 parameters->AddObject(new DoubleVecParam(SealevelriseHElasticEnum,H_elastic,M));
[20028]456
[24921]457 /*free ressources: */
458 xDelete<IssmDouble>(love_h);
459 xDelete<IssmDouble>(love_k);
460 xDelete<IssmDouble>(love_l);
461 xDelete<IssmDouble>(G_rigid);
462 xDelete<IssmDouble>(G_rigid_local);
463 xDelete<IssmDouble>(G_elastic);
464 xDelete<IssmDouble>(G_elastic_local);
465 xDelete<IssmDouble>(U_elastic);
466 xDelete<IssmDouble>(U_elastic_local);
467 xDelete<IssmDouble>(H_elastic);
468 xDelete<IssmDouble>(H_elastic_local);
469/*}}}*/
[24481]470 /*Transitions:{{{ */
[20690]471 iomodel->FetchData(&transitions,&transitions_M,&transitions_N,&ntransitions,"md.slr.transitions");
[20136]472 if(transitions){
473 parameters->AddObject(new DoubleMatArrayParam(SealevelriseTransitionsEnum,transitions,ntransitions,transitions_M,transitions_N));
474
475 for(int i=0;i<ntransitions;i++){
476 IssmDouble* transition=transitions[i];
477 xDelete<IssmDouble>(transition);
478 }
479 xDelete<IssmDouble*>(transitions);
[20138]480 xDelete<int>(transitions_M);
481 xDelete<int>(transitions_N);
[24481]482 } /*}}}*/
483 /*Requested outputs {{{*/
[20690]484 iomodel->FindConstant(&requestedoutputs,&numoutputs,"md.slr.requested_outputs");
[20036]485 if(numoutputs)parameters->AddObject(new StringArrayParam(SealevelriseRequestedOutputsEnum,requestedoutputs,numoutputs));
[20690]486 iomodel->DeleteData(&requestedoutputs,numoutputs,"md.slr.requested_outputs");
[24481]487 /*}}}*/
[20029]488
[19984]489}/*}}}*/
490
491/*Finite Element Analysis*/
492void SealevelriseAnalysis::Core(FemModel* femmodel){/*{{{*/
493 _error_("not implemented");
494}/*}}}*/
495ElementVector* SealevelriseAnalysis::CreateDVector(Element* element){/*{{{*/
496 /*Default, return NULL*/
497 return NULL;
498}/*}}}*/
499ElementMatrix* SealevelriseAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/
500_error_("Not implemented");
501}/*}}}*/
502ElementMatrix* SealevelriseAnalysis::CreateKMatrix(Element* element){/*{{{*/
503 _error_("not implemented yet");
504}/*}}}*/
505ElementVector* SealevelriseAnalysis::CreatePVector(Element* element){/*{{{*/
506_error_("not implemented yet");
507}/*}}}*/
508void SealevelriseAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
509 _error_("not implemented yet");
510}/*}}}*/
511void SealevelriseAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
512 _error_("Not implemented yet");
513}/*}}}*/
514void SealevelriseAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
[22955]515 _error_("not implemeneted yet!");
[20153]516
[19984]517}/*}}}*/
518void SealevelriseAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
519 /*Default, do nothing*/
520 return;
521}/*}}}*/
Note: See TracBrowser for help on using the repository browser.