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

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

CHG: new partition vector for solidearthsettings. New handling of rigid and
elastic table precomputations. New Europa planet radius. New earth default in
model constructor.
File size: 19.4 KB
RevLine 
[19984]1#include "./SealevelriseAnalysis.h"
[24959]2#include <math.h>
[19984]3#include "../toolkits/toolkits.h"
4#include "../classes/classes.h"
[25379]5#include "../classes/Inputs/TransientInput.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}/*}}}*/
[25379]22void SealevelriseAnalysis::UpdateElements(Elements* elements,Inputs* inputs,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);
[25379]32 element->Update(inputs,i,iomodel,analysis_counter,analysis_type,P1Enum);
[19984]33 counter++;
34 }
35 }
36
37 /*Create inputs: */
[25379]38 iomodel->FetchDataToInput(inputs,elements,"md.mask.ocean_levelset",MaskOceanLevelsetEnum);
39 iomodel->FetchDataToInput(inputs,elements,"md.mask.ice_levelset",MaskIceLevelsetEnum);
[25118]40 iomodel->FetchData(&geodetic,"md.solidearth.settings.computesealevelchange");
[25379]41 iomodel->FetchDataToInput(inputs,elements,"md.solidearth.surfaceload.icethicknesschange",SurfaceloadIceThicknessChangeEnum);
42 iomodel->FetchDataToInput(inputs,elements,"md.solidearth.sealevel",SealevelEnum,0);
43 iomodel->FetchDataToInput(inputs,elements,"md.geometry.bed",BedEnum);
44 iomodel->FetchDataToInput(inputs,elements,"md.solidearth.surfaceload.waterheightchange",SurfaceloadWaterHeightChangeEnum);
[24469]45
46 /*dynamic sea level: */
47 iomodel->FetchData(&dslmodel,"md.dsl.model");
[24505]48 if (dslmodel==1){ /*standard dsl model:{{{*/
[24481]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: */
[25379]63 inputs->SetTransientInput(DslGlobalAverageThermostericSeaLevelChangeEnum,times,N);
64 TransientInput* transientinput = inputs->GetTransientInput(DslGlobalAverageThermostericSeaLevelChangeEnum);
[24481]65
66
[25539]67 for(Object* & object : elements->objects){
68 Element* element=xDynamicCast<Element*>(object);
[24481]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: */
[25379]84 iomodel->FetchDataToInput(inputs,elements,"md.dsl.sea_surface_height_change_above_geoid", DslSeaSurfaceHeightChangeAboveGeoidEnum);
85 iomodel->FetchDataToInput(inputs,elements,"md.dsl.sea_water_pressure_change_at_sea_floor", DslSeaWaterPressureChangeAtSeaFloorEnum);
[24481]86
[24505]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
[25275]99 /*deal with dsl.global_average_thermosteric_sea_level_change {{{*/
[24505]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
[25379]113 TransientInput* transientinput=inputs->SetDatasetTransientInput(DslGlobalAverageThermostericSeaLevelChangeEnum,i, times,N);
[24505]114
[25539]115 for(Object* & object : elements->objects){
116 Element* element=xDynamicCast<Element*>(object);
[24505]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 /*}}}*/
[25379]137 iomodel->FetchDataToInput(inputs,elements,"md.dsl.sea_surface_height_change_above_geoid",DslSeaSurfaceHeightChangeAboveGeoidEnum);
138 iomodel->FetchDataToInput(inputs,elements,"md.dsl.sea_water_pressure_change_at_sea_floor",DslSeaWaterPressureChangeAtSeaFloorEnum);
[24505]139
140 } /*}}}*/
[24469]141 else _error_("Dsl model " << dslmodel << " not implemented yet!");
[25060]142
[25066]143
[25060]144
[22968]145 /*Initialize cumdeltalthickness and sealevel rise rate input*/
[25379]146 InputUpdateFromConstantx(inputs,elements,0.,SealevelriseCumDeltathicknessEnum);
147 InputUpdateFromConstantx(inputs,elements,0.,SealevelNEsaRateEnum);
148 InputUpdateFromConstantx(inputs,elements,0.,SealevelUEsaRateEnum);
149 InputUpdateFromConstantx(inputs,elements,0.,SealevelRSLRateEnum);
150 InputUpdateFromConstantx(inputs,elements,0.,SealevelEustaticMaskEnum);
151 InputUpdateFromConstantx(inputs,elements,0.,SealevelEustaticOceanMaskEnum);
[22968]152
[19984]153}/*}}}*/
154void SealevelriseAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
155
[19986]156 int nl;
[19984]157 IssmDouble* love_h=NULL;
158 IssmDouble* love_k=NULL;
[20999]159 IssmDouble* love_l=NULL;
[25139]160 IssmDouble* love_th=NULL;
161 IssmDouble* love_tk=NULL;
162 IssmDouble* love_tl=NULL;
[24469]163 int dslmodel=0;
[23066]164
[24921]165 IssmDouble* G_rigid = NULL;
166 IssmDouble* G_rigid_local = NULL;
[20028]167 IssmDouble* G_elastic = NULL;
[20033]168 IssmDouble* G_elastic_local = NULL;
[20999]169 IssmDouble* U_elastic = NULL;
170 IssmDouble* U_elastic_local = NULL;
171 IssmDouble* H_elastic = NULL;
172 IssmDouble* H_elastic_local = NULL;
[20033]173 int M,m,lower_row,upper_row;
174 IssmDouble degacc=.01;
[24959]175 IssmDouble planetradius=0;
176 IssmDouble planetarea=0;
[25655]177 bool rigid=false;
178 bool elastic=false;
[25751]179 bool rotation=false;
[19984]180
[20036]181 int numoutputs;
182 char** requestedoutputs = NULL;
183
[20136]184 /*transition vectors: */
185 IssmDouble **transitions = NULL;
186 int *transitions_M = NULL;
187 int *transitions_N = NULL;
188 int ntransitions;
189
[19986]190 /*some constant parameters: */
[24531]191 parameters->AddObject(iomodel->CopyConstantObject("md.dsl.model",DslModelEnum));
[25118]192 parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.runfrequency",SolidearthSettingsRunFrequencyEnum));
193 parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.reltol",SolidearthSettingsReltolEnum));
194 parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.abstol",SolidearthSettingsAbstolEnum));
195 parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.maxiter",SolidearthSettingsMaxiterEnum));
196 parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.rigid",SolidearthSettingsRigidEnum));
197 parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.horiz",SolidearthSettingsHorizEnum));
198 parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.elastic",SolidearthSettingsElasticEnum));
199 parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.rotation",SolidearthSettingsRotationEnum));
200 parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.rotational.equatorialmoi",RotationalEquatorialMoiEnum));
201 parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.rotational.polarmoi",RotationalPolarMoiEnum));
202 parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.rotational.angularvelocity",RotationalAngularVelocityEnum));
203 parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.ocean_area_scaling",SolidearthSettingsOceanAreaScalingEnum));
204 parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.computesealevelchange",SolidearthSettingsComputesealevelchangeEnum));
205 parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.planetradius",SolidearthPlanetRadiusEnum));
[25751]206 parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.glfraction",SolidearthSettingsGlfractionEnum));
[25627]207 parameters->AddObject(new DoubleParam(CumBslrEnum,0.0));
208 parameters->AddObject(new DoubleParam(CumBslrIceEnum,0.0));
209 parameters->AddObject(new DoubleParam(CumBslrHydroEnum,0.0));
210 parameters->AddObject(new DoubleParam(CumGmtslrEnum,0.0));
[19986]211
[24959]212 /*compute planet area and plug into parameters:*/
[25118]213 iomodel->FetchData(&planetradius,"md.solidearth.planetradius");
[24959]214 planetarea=4*PI*planetradius*planetradius;
[25118]215 parameters->AddObject(new DoubleParam(SolidearthPlanetAreaEnum,planetarea));
[24959]216
[25751]217 /*Deal with partition of the barystatic contribution:*/
218 iomodel->FetchData(&npartice,"md.solidearth.npartice");
219 parameters->AddObject(new IntParam(SolidearthNpartIceEnum,npartice));
220 if(npartice){
221 iomodel->FetchData(&partitionice,&nel,NULL,"md.solidearth.partitionice");
222 parameters->AddObject(new DoubleMatParam(SolidearthPartitionIceEnum,partitionice,nel,1));
223 parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.npartice",SolidearthNpartIceEnum));
224 bslrice_partition=xNewZeroInit<IssmDouble>(npartice);
225 parameters->AddObject(new DoubleMatParam(CumBslrIcePartitionEnum,bslrice_partition,npartice,1));
226 xDelete<IssmDouble>(partitionice);
227 }
228 iomodel->FetchData(&nparthydro,"md.solidearth.nparthydro");
229 parameters->AddObject(new IntParam(SolidearthNpartHydroEnum,nparthydro));
230 if(nparthydro){
231 iomodel->FetchData(&partitionhydro,&nel,NULL,"md.solidearth.partitionhydro");
232 parameters->AddObject(new DoubleMatParam(SolidearthPartitionHydroEnum,partitionhydro,nel,1));
233 parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.nparthydro",SolidearthNpartHydroEnum));
234 bslrhydro_partition=xNewZeroInit<IssmDouble>(nparthydro);
235 parameters->AddObject(new DoubleMatParam(CumBslrHydroPartitionEnum,bslrhydro_partition,nparthydro,1));
236 xDelete<IssmDouble>(partitionhydro);
237 }
238
239
[24505]240 /*Deal with dsl multi-model ensembles: {{{*/
241 iomodel->FetchData(&dslmodel,"md.dsl.model");
242 parameters->AddObject(iomodel->CopyConstantObject("md.dsl.compute_fingerprints",DslComputeFingerprintsEnum));
243 if(dslmodel==2){
[25273]244 IssmDouble modelid;
[24505]245 int nummodels;
[25273]246
247 /*create double param, not int param, because Dakota will be updating it as
248 * a double potentially: */
249 iomodel->FetchData(&modelid,"md.dsl.modelid");
[25274]250 parameters->AddObject(new DoubleParam(DslModelidEnum,modelid));
[24505]251 parameters->AddObject(iomodel->CopyConstantObject("md.dsl.nummodels",DslNummodelsEnum));
252 iomodel->FetchData(&nummodels,"md.dsl.nummodels");
253
254 /*quick checks: */
255 if(nummodels<=0)_error_("dslmme object in md.dsl field should contain at least 1 ensemble model!");
256 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!");
257 } /*}}}*/
[24481]258 /*Deal with elasticity {{{*/
[25655]259 iomodel->FetchData(&rigid,"md.solidearth.settings.rigid");
[25627]260 iomodel->FetchData(&elastic,"md.solidearth.settings.elastic");
[25751]261 iomodel->FetchData(&rotation,"md.solidearth.settings.rotation");
[19984]262
[25751]263 if(elastic | rigid){
264 /*compute green functions for a range of angles*/
265 iomodel->FetchData(&degacc,"md.solidearth.settings.degacc");
266 M=reCast<int,IssmDouble>(180./degacc+1.);
267 }
[25139]268
[25751]269 /*love numbers: */
270 if(elastic){
271 iomodel->FetchData(&love_h,&nl,NULL,"md.solidearth.lovenumbers.h");
272 iomodel->FetchData(&love_k,&nl,NULL,"md.solidearth.lovenumbers.k");
273 iomodel->FetchData(&love_l,&nl,NULL,"md.solidearth.lovenumbers.l");
274 iomodel->FetchData(&love_th,&nl,NULL,"md.solidearth.lovenumbers.th");
275 iomodel->FetchData(&love_tk,&nl,NULL,"md.solidearth.lovenumbers.tk");
276 iomodel->FetchData(&love_tl,&nl,NULL,"md.solidearth.lovenumbers.tl");
[20022]277
[25751]278 parameters->AddObject(new DoubleMatParam(LoadLoveHEnum,love_h,nl,1));
279 parameters->AddObject(new DoubleMatParam(LoadLoveKEnum,love_k,nl,1));
280 parameters->AddObject(new DoubleMatParam(LoadLoveLEnum,love_l,nl,1));
281 parameters->AddObject(new DoubleMatParam(TidalLoveHEnum,love_th,nl,1));
282 parameters->AddObject(new DoubleMatParam(TidalLoveKEnum,love_tk,nl,1));
283 parameters->AddObject(new DoubleMatParam(TidalLoveLEnum,love_tl,nl,1));
[21742]284
[25751]285 // AD performance is sensitive to calls to ensurecontiguous.
286 // // Providing "t" will cause ensurecontiguous to be called.
287 #ifdef _HAVE_AD_
288 G_elastic=xNew<IssmDouble>(M,"t");
289 U_elastic=xNew<IssmDouble>(M,"t");
290 H_elastic=xNew<IssmDouble>(M,"t");
291 #else
292 G_elastic=xNew<IssmDouble>(M);
293 U_elastic=xNew<IssmDouble>(M);
294 H_elastic=xNew<IssmDouble>(M);
295 #endif
296 }
297 if(rigid){
298 #ifdef _HAVE_AD_
299 G_rigid=xNew<IssmDouble>(M,"t");
300 #else
301 G_rigid=xNew<IssmDouble>(M);
302 #endif
303 }
304
305 if(rotation)parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.lovenumbers.tk2secular",TidalLoveK2SecularEnum));
[21742]306
[25751]307 if(rigid | elastic){
[20033]308
[25751]309 /*compute combined legendre + love number (elastic green function:*/
310 m=DetermineLocalSize(M,IssmComm::GetComm());
311 GetOwnershipBoundariesFromRange(&lower_row,&upper_row,m,IssmComm::GetComm());
312 }
313 if(elastic){
314 #ifdef _HAVE_AD_
315 G_elastic_local=xNew<IssmDouble>(m,"t");
316 U_elastic_local=xNew<IssmDouble>(m,"t");
317 H_elastic_local=xNew<IssmDouble>(m,"t");
318 #else
319 G_elastic_local=xNew<IssmDouble>(m);
320 U_elastic_local=xNew<IssmDouble>(m);
321 H_elastic_local=xNew<IssmDouble>(m);
322 #endif
323 }
324 if(rigid){
325 #ifdef _HAVE_AD_
326 G_rigid_local=xNew<IssmDouble>(m,"t");
327 #else
328 G_rigid_local=xNew<IssmDouble>(m);
329 #endif
330 }
[20029]331
[25751]332 if(rigid){
333 for(int i=lower_row;i<upper_row;i++){
334 IssmDouble alpha,x;
335 alpha= reCast<IssmDouble>(i)*degacc * PI / 180.0;
336 G_rigid_local[i-lower_row]= .5/sin(alpha/2.0);
337 }
338 }
339 if(elastic){
340 for(int i=lower_row;i<upper_row;i++){
341 IssmDouble alpha,x;
342 alpha= reCast<IssmDouble>(i)*degacc * PI / 180.0;
[21908]343
[25751]344 G_elastic_local[i-lower_row]= (love_k[nl-1]-love_h[nl-1])*G_rigid_local[i-lower_row];
345 U_elastic_local[i-lower_row]= (love_h[nl-1])*G_rigid_local[i-lower_row];
346 H_elastic_local[i-lower_row]= 0;
347 IssmDouble Pn = 0.;
348 IssmDouble Pn1 = 0.;
349 IssmDouble Pn2 = 0.;
350 IssmDouble Pn_p = 0.;
351 IssmDouble Pn_p1 = 0.;
352 IssmDouble Pn_p2 = 0.;
[20033]353
[25751]354 for (int n=0;n<nl;n++) {
355 IssmDouble deltalove_G;
356 IssmDouble deltalove_U;
[25655]357
[25751]358 deltalove_G = (love_k[n]-love_k[nl-1]-love_h[n]+love_h[nl-1]);
359 deltalove_U = (love_h[n]-love_h[nl-1]);
360
361 /*compute legendre polynomials: P_n(cos\theta) & d P_n(cos\theta)/ d\theta: */
362 if(n==0){
363 Pn=1;
364 Pn_p=0;
[25627]365 }
[25751]366 else if(n==1){
367 Pn = cos(alpha);
368 Pn_p = 1;
369 }
370 else{
371 Pn = ( (2*n-1)*cos(alpha)*Pn1 - (n-1)*Pn2 ) /n;
372 Pn_p = ( (2*n-1)*(Pn1+cos(alpha)*Pn_p1) - (n-1)*Pn_p2 ) /n;
373 }
374 Pn2=Pn1; Pn1=Pn;
375 Pn_p2=Pn_p1; Pn_p1=Pn_p;
376
377 G_elastic_local[i-lower_row] += deltalove_G*Pn; // gravitational potential
378 U_elastic_local[i-lower_row] += deltalove_U*Pn; // vertical (up) displacement
379 H_elastic_local[i-lower_row] += sin(alpha)*love_l[n]*Pn_p; // horizontal displacements
[25655]380 }
[25751]381 }
382 }
383 if(rigid){
[25627]384
[25751]385 /*merge G_elastic_local into G_elastic; U_elastic_local into U_elastic; H_elastic_local to H_elastic:{{{*/
386 int* recvcounts=xNew<int>(IssmComm::GetSize());
387 int* displs=xNew<int>(IssmComm::GetSize());
[25655]388
[25751]389 //recvcounts:
390 ISSM_MPI_Allgather(&m,1,ISSM_MPI_INT,recvcounts,1,ISSM_MPI_INT,IssmComm::GetComm());
[25655]391
[25751]392 /*displs: */
393 ISSM_MPI_Allgather(&lower_row,1,ISSM_MPI_INT,displs,1,ISSM_MPI_INT,IssmComm::GetComm());
[25655]394
[25751]395 /*All gather:*/
396 ISSM_MPI_Allgatherv(G_rigid_local, m, ISSM_MPI_DOUBLE, G_rigid, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
397 if(elastic){
[25655]398 ISSM_MPI_Allgatherv(G_elastic_local, m, ISSM_MPI_DOUBLE, G_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
399 ISSM_MPI_Allgatherv(U_elastic_local, m, ISSM_MPI_DOUBLE, U_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
400 ISSM_MPI_Allgatherv(H_elastic_local, m, ISSM_MPI_DOUBLE, H_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
[25751]401 }
402
403 /*free resources: */
404 xDelete<int>(recvcounts);
405 xDelete<int>(displs);
[25655]406
[25751]407 /*Avoid singularity at 0: */
408 G_rigid[0]=G_rigid[1];
409 G_elastic[0]=G_elastic[1];
410 U_elastic[0]=U_elastic[1];
411 H_elastic[0]=H_elastic[1];
[25655]412
[25751]413 /*Save our precomputed tables into parameters*/
414 parameters->AddObject(new DoubleVecParam(SealevelriseGRigidEnum,G_rigid,M));
415 if(elastic){
[25655]416 parameters->AddObject(new DoubleVecParam(SealevelriseGElasticEnum,G_elastic,M));
417 parameters->AddObject(new DoubleVecParam(SealevelriseUElasticEnum,U_elastic,M));
418 parameters->AddObject(new DoubleVecParam(SealevelriseHElasticEnum,H_elastic,M));
[25751]419 }
[25655]420
[25751]421 /*free resources: */
422 xDelete<IssmDouble>(G_rigid);
423 xDelete<IssmDouble>(G_rigid_local);
424 if(elastic){
[25655]425 xDelete<IssmDouble>(love_h);
426 xDelete<IssmDouble>(love_k);
427 xDelete<IssmDouble>(love_l);
428 xDelete<IssmDouble>(love_th);
429 xDelete<IssmDouble>(love_tk);
430 xDelete<IssmDouble>(love_tl);
431 xDelete<IssmDouble>(G_elastic);
432 xDelete<IssmDouble>(G_elastic_local);
433 xDelete<IssmDouble>(U_elastic);
434 xDelete<IssmDouble>(U_elastic_local);
435 xDelete<IssmDouble>(H_elastic);
436 xDelete<IssmDouble>(H_elastic_local);
437 }
438 }
[25751]439 /*}}}*/
440
[25655]441 /*Indicate we have not yet run the Geometry Core module: */
[24977]442 parameters->AddObject(new BoolParam(SealevelriseGeometryDoneEnum,false));
443
[24481]444 /*Transitions:{{{ */
[25118]445 iomodel->FetchData(&transitions,&transitions_M,&transitions_N,&ntransitions,"md.solidearth.transitions");
[20136]446 if(transitions){
447 parameters->AddObject(new DoubleMatArrayParam(SealevelriseTransitionsEnum,transitions,ntransitions,transitions_M,transitions_N));
448
449 for(int i=0;i<ntransitions;i++){
450 IssmDouble* transition=transitions[i];
451 xDelete<IssmDouble>(transition);
452 }
453 xDelete<IssmDouble*>(transitions);
[20138]454 xDelete<int>(transitions_M);
455 xDelete<int>(transitions_N);
[24481]456 } /*}}}*/
457 /*Requested outputs {{{*/
[25118]458 iomodel->FindConstant(&requestedoutputs,&numoutputs,"md.solidearth.requested_outputs");
[20036]459 if(numoutputs)parameters->AddObject(new StringArrayParam(SealevelriseRequestedOutputsEnum,requestedoutputs,numoutputs));
[25118]460 iomodel->DeleteData(&requestedoutputs,numoutputs,"md.solidearth.requested_outputs");
[24481]461 /*}}}*/
[20029]462
[19984]463}/*}}}*/
464
465/*Finite Element Analysis*/
466void SealevelriseAnalysis::Core(FemModel* femmodel){/*{{{*/
467 _error_("not implemented");
468}/*}}}*/
469ElementVector* SealevelriseAnalysis::CreateDVector(Element* element){/*{{{*/
470 /*Default, return NULL*/
471 return NULL;
472}/*}}}*/
473ElementMatrix* SealevelriseAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/
474_error_("Not implemented");
475}/*}}}*/
476ElementMatrix* SealevelriseAnalysis::CreateKMatrix(Element* element){/*{{{*/
477 _error_("not implemented yet");
478}/*}}}*/
479ElementVector* SealevelriseAnalysis::CreatePVector(Element* element){/*{{{*/
480_error_("not implemented yet");
481}/*}}}*/
482void SealevelriseAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
483 _error_("not implemented yet");
484}/*}}}*/
[25317]485void SealevelriseAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_interp,int control_index){/*{{{*/
[19984]486 _error_("Not implemented yet");
487}/*}}}*/
488void SealevelriseAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
[22955]489 _error_("not implemeneted yet!");
[20153]490
[19984]491}/*}}}*/
492void SealevelriseAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
493 /*Default, do nothing*/
494 return;
495}/*}}}*/
Note: See TracBrowser for help on using the repository browser.