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

Last change on this file since 24861 was 24861, checked in by Mathieu Morlighem, 5 years ago

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