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

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

CHG: 10% increase in sea level solver speed, using a G_rigid kernel pre computed
in the SealevelriseAnalysis section, before G_elastic kernels.

File size: 19.6 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));
[19986]323
[24505]324 /*Deal with dsl multi-model ensembles: {{{*/
325 iomodel->FetchData(&dslmodel,"md.dsl.model");
326 parameters->AddObject(iomodel->CopyConstantObject("md.dsl.compute_fingerprints",DslComputeFingerprintsEnum));
327 if(dslmodel==2){
328 int modelid;
329 int nummodels;
330
331 parameters->AddObject(iomodel->CopyConstantObject("md.dsl.modelid",DslModelidEnum));
332 parameters->AddObject(iomodel->CopyConstantObject("md.dsl.nummodels",DslNummodelsEnum));
333 iomodel->FetchData(&modelid,"md.dsl.modelid");
334 iomodel->FetchData(&nummodels,"md.dsl.nummodels");
335
336 /*quick checks: */
337 if(nummodels<=0)_error_("dslmme object in md.dsl field should contain at least 1 ensemble model!");
338 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!");
339 } /*}}}*/
[24481]340 /*Deal with elasticity {{{*/
[24921]341 /*love numbers: */
342 iomodel->FetchData(&love_h,&nl,NULL,"md.slr.love_h");
343 iomodel->FetchData(&love_k,&nl,NULL,"md.slr.love_k");
344 iomodel->FetchData(&love_l,&nl,NULL,"md.slr.love_l");
[19984]345
[24921]346 /*compute elastic green function for a range of angles*/
347 iomodel->FetchData(&degacc,"md.slr.degacc");
348 M=reCast<int,IssmDouble>(180./degacc+1.);
[20022]349
[24921]350 // AD performance is sensitive to calls to ensurecontiguous.
351 // // Providing "t" will cause ensurecontiguous to be called.
352 #ifdef _HAVE_AD_
353 G_rigid=xNew<IssmDouble>(M,"t");
354 G_elastic=xNew<IssmDouble>(M,"t");
355 U_elastic=xNew<IssmDouble>(M,"t");
356 H_elastic=xNew<IssmDouble>(M,"t");
357 #else
358 G_rigid=xNew<IssmDouble>(M);
359 G_elastic=xNew<IssmDouble>(M);
360 U_elastic=xNew<IssmDouble>(M);
361 H_elastic=xNew<IssmDouble>(M);
362 #endif
[21742]363
[24921]364 /*compute combined legendre + love number (elastic green function:*/
365 m=DetermineLocalSize(M,IssmComm::GetComm());
366 GetOwnershipBoundariesFromRange(&lower_row,&upper_row,m,IssmComm::GetComm());
367 // AD performance is sensitive to calls to ensurecontiguous.
368 // // Providing "t" will cause ensurecontiguous to be called.
369 #ifdef _HAVE_AD_
370 G_elastic_local=xNew<IssmDouble>(m,"t");
371 G_rigid_local=xNew<IssmDouble>(m,"t");
372 U_elastic_local=xNew<IssmDouble>(m,"t");
373 H_elastic_local=xNew<IssmDouble>(m,"t");
374 #else
375 G_elastic_local=xNew<IssmDouble>(m);
376 G_rigid_local=xNew<IssmDouble>(m);
377 U_elastic_local=xNew<IssmDouble>(m);
378 H_elastic_local=xNew<IssmDouble>(m);
379 #endif
[21742]380
[24921]381 for(int i=lower_row;i<upper_row;i++){
382 IssmDouble alpha,x;
383 alpha= reCast<IssmDouble>(i)*degacc * PI / 180.0;
[20033]384
[24921]385 G_rigid_local[i-lower_row]= .5/sin(alpha/2.0);
386 G_elastic_local[i-lower_row]= (love_k[nl-1]-love_h[nl-1])*G_rigid_local[i-lower_row];
387 U_elastic_local[i-lower_row]= (love_h[nl-1])*G_rigid_local[i-lower_row];
388 H_elastic_local[i-lower_row]= 0;
389 IssmDouble Pn = 0.;
390 IssmDouble Pn1 = 0.;
391 IssmDouble Pn2 = 0.;
392 IssmDouble Pn_p = 0.;
393 IssmDouble Pn_p1 = 0.;
394 IssmDouble Pn_p2 = 0.;
[20029]395
[24921]396 for (int n=0;n<nl;n++) {
397 IssmDouble deltalove_G;
398 IssmDouble deltalove_U;
[21908]399
[24921]400 deltalove_G = (love_k[n]-love_k[nl-1]-love_h[n]+love_h[nl-1]);
401 deltalove_U = (love_h[n]-love_h[nl-1]);
[20033]402
[24921]403 /*compute legendre polynomials: P_n(cos\theta) & d P_n(cos\theta)/ d\theta: */
404 if(n==0){
405 Pn=1;
406 Pn_p=0;
407 }
408 else if(n==1){
409 Pn = cos(alpha);
410 Pn_p = 1;
411 }
412 else{
413 Pn = ( (2*n-1)*cos(alpha)*Pn1 - (n-1)*Pn2 ) /n;
414 Pn_p = ( (2*n-1)*(Pn1+cos(alpha)*Pn_p1) - (n-1)*Pn_p2 ) /n;
415 }
416 Pn2=Pn1; Pn1=Pn;
417 Pn_p2=Pn_p1; Pn_p1=Pn_p;
[23066]418
[24921]419 G_elastic_local[i-lower_row] += deltalove_G*Pn; // gravitational potential
420 U_elastic_local[i-lower_row] += deltalove_U*Pn; // vertical (up) displacement
421 H_elastic_local[i-lower_row] += sin(alpha)*love_l[n]*Pn_p; // horizontal displacements
[20033]422 }
[24921]423 }
[20033]424
[24921]425 /*merge G_elastic_local into G_elastic; U_elastic_local into U_elastic; H_elastic_local to H_elastic:{{{*/
426 int* recvcounts=xNew<int>(IssmComm::GetSize());
427 int* displs=xNew<int>(IssmComm::GetSize());
[20033]428
[24921]429 //recvcounts:
430 ISSM_MPI_Allgather(&m,1,ISSM_MPI_INT,recvcounts,1,ISSM_MPI_INT,IssmComm::GetComm());
[20033]431
[24921]432 /*displs: */
433 ISSM_MPI_Allgather(&lower_row,1,ISSM_MPI_INT,displs,1,ISSM_MPI_INT,IssmComm::GetComm());
[20033]434
[24921]435 /*All gather:*/
436 ISSM_MPI_Allgatherv(G_rigid_local, m, ISSM_MPI_DOUBLE, G_rigid, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
437 ISSM_MPI_Allgatherv(G_elastic_local, m, ISSM_MPI_DOUBLE, G_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
438 ISSM_MPI_Allgatherv(U_elastic_local, m, ISSM_MPI_DOUBLE, U_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
439 ISSM_MPI_Allgatherv(H_elastic_local, m, ISSM_MPI_DOUBLE, H_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
440 /*free ressources: */
441 xDelete<int>(recvcounts);
442 xDelete<int>(displs);
[20033]443
[24921]444 /*}}}*/
[20033]445
[24921]446 /*Avoid singularity at 0: */
447 G_rigid[0]=G_rigid[1];
448 parameters->AddObject(new DoubleVecParam(SealevelriseGRigidEnum,G_rigid,M));
449 G_elastic[0]=G_elastic[1];
450 parameters->AddObject(new DoubleVecParam(SealevelriseGElasticEnum,G_elastic,M));
451 U_elastic[0]=U_elastic[1];
452 parameters->AddObject(new DoubleVecParam(SealevelriseUElasticEnum,U_elastic,M));
453 H_elastic[0]=H_elastic[1];
454 parameters->AddObject(new DoubleVecParam(SealevelriseHElasticEnum,H_elastic,M));
[20028]455
[24921]456 /*free ressources: */
457 xDelete<IssmDouble>(love_h);
458 xDelete<IssmDouble>(love_k);
459 xDelete<IssmDouble>(love_l);
460 xDelete<IssmDouble>(G_rigid);
461 xDelete<IssmDouble>(G_rigid_local);
462 xDelete<IssmDouble>(G_elastic);
463 xDelete<IssmDouble>(G_elastic_local);
464 xDelete<IssmDouble>(U_elastic);
465 xDelete<IssmDouble>(U_elastic_local);
466 xDelete<IssmDouble>(H_elastic);
467 xDelete<IssmDouble>(H_elastic_local);
468/*}}}*/
[24481]469 /*Transitions:{{{ */
[20690]470 iomodel->FetchData(&transitions,&transitions_M,&transitions_N,&ntransitions,"md.slr.transitions");
[20136]471 if(transitions){
472 parameters->AddObject(new DoubleMatArrayParam(SealevelriseTransitionsEnum,transitions,ntransitions,transitions_M,transitions_N));
473
474 for(int i=0;i<ntransitions;i++){
475 IssmDouble* transition=transitions[i];
476 xDelete<IssmDouble>(transition);
477 }
478 xDelete<IssmDouble*>(transitions);
[20138]479 xDelete<int>(transitions_M);
480 xDelete<int>(transitions_N);
[24481]481 } /*}}}*/
482 /*Requested outputs {{{*/
[20690]483 iomodel->FindConstant(&requestedoutputs,&numoutputs,"md.slr.requested_outputs");
[20036]484 if(numoutputs)parameters->AddObject(new StringArrayParam(SealevelriseRequestedOutputsEnum,requestedoutputs,numoutputs));
[20690]485 iomodel->DeleteData(&requestedoutputs,numoutputs,"md.slr.requested_outputs");
[24481]486 /*}}}*/
[20029]487
[19984]488}/*}}}*/
489
490/*Finite Element Analysis*/
491void SealevelriseAnalysis::Core(FemModel* femmodel){/*{{{*/
492 _error_("not implemented");
493}/*}}}*/
494ElementVector* SealevelriseAnalysis::CreateDVector(Element* element){/*{{{*/
495 /*Default, return NULL*/
496 return NULL;
497}/*}}}*/
498ElementMatrix* SealevelriseAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/
499_error_("Not implemented");
500}/*}}}*/
501ElementMatrix* SealevelriseAnalysis::CreateKMatrix(Element* element){/*{{{*/
502 _error_("not implemented yet");
503}/*}}}*/
504ElementVector* SealevelriseAnalysis::CreatePVector(Element* element){/*{{{*/
505_error_("not implemented yet");
506}/*}}}*/
507void SealevelriseAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
508 _error_("not implemented yet");
509}/*}}}*/
510void SealevelriseAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
511 _error_("Not implemented yet");
512}/*}}}*/
513void SealevelriseAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
[22955]514 _error_("not implemeneted yet!");
[20153]515
[19984]516}/*}}}*/
517void SealevelriseAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
518 /*Default, do nothing*/
519 return;
520}/*}}}*/
Note: See TracBrowser for help on using the repository browser.