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

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

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