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

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

CHG: chaning love class to lovenumbers.
File size: 20.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);
[25118]40 iomodel->FetchData(&geodetic,"md.solidearth.settings.computesealevelchange");
41 iomodel->FetchDataToInput(inputs2,elements,"md.solidearth.surfaceload.icethicknesschange",SurfaceloadIceThicknessChangeEnum);
42 iomodel->FetchDataToInput(inputs2,elements,"md.solidearth.sealevel",SealevelEnum,0);
[24335]43 iomodel->FetchDataToInput(inputs2,elements,"md.geometry.bed",BedEnum);
[25118]44 iomodel->FetchDataToInput(inputs2,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: */
63 inputs2->SetTransientInput(DslGlobalAverageThermostericSeaLevelChangeEnum,times,N);
64 TransientInput2* transientinput = inputs2->GetTransientInput(DslGlobalAverageThermostericSeaLevelChangeEnum);
65
66
67 for(int i=0;i<elements->Size();i++){
68 Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
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: */
[24469]84 iomodel->FetchDataToInput(inputs2,elements,"md.dsl.sea_surface_height_change_above_geoid", DslSeaSurfaceHeightChangeAboveGeoidEnum);
85 iomodel->FetchDataToInput(inputs2,elements,"md.dsl.sea_water_pressure_change_at_sea_floor", DslSeaWaterPressureChangeAtSeaFloor);
[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
99 /*deal with dsl.sea_surface_height_change_above_geoid {{{*/
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
113 TransientInput2* transientinput=inputs2->SetDatasetTransientInput(DslGlobalAverageThermostericSeaLevelChangeEnum,i, times,N);
114
115 for(int j=0;j<elements->Size();j++){
116 Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(j));
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 /*}}}*/
137 /*now do the same with the dsl.sea_surface_height_change_above_geoid field:{{{ */
138 iomodel->FetchData(&pstr,&pM,&pN,&nummodels,"md.dsl.sea_surface_height_change_above_geoid");
139
140 for (int i=0;i<nummodels;i++){
141 M=pM[i];
142 N=pN[i];
143 str=pstr[i];
144
145
146 //recover time vector:
147 times=xNew<IssmDouble>(N);
148 for(int t=0;t<N;t++) times[t] = str[(M-1)*N+t];
149
150 TransientInput2* transientinput=inputs2->SetDatasetTransientInput(DslSeaSurfaceHeightChangeAboveGeoidEnum,i, times,N);
151
152 for(int j=0;j<elements->Size();j++){
153
154 /*Get the right transient input*/
155 Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(j));
156
157 /*Get values and lid list*/
158 const int numvertices = element->GetNumberOfVertices();
159 int *vertexlids = xNew<int>(numvertices);
160 int *vertexsids = xNew<int>(numvertices);
161
162
163 /*Recover vertices ids needed to initialize inputs*/
164 _assert_(iomodel->elements);
165 for(int k=0;k<numvertices;k++){
166 vertexsids[k] =reCast<int>(iomodel->elements[numvertices*element->Sid()+k]); //ids for vertices are in the elements array from Matlab
167 vertexlids[k]=iomodel->my_vertices_lids[vertexsids[k]-1];
168 }
169
170 //element->GetVerticesLidList(vertexlids);
171 //element->GetVerticesSidList(vertexsids);
172 IssmDouble* values=xNew<IssmDouble>(numvertices);
173
174 for(int t=0;t<N;t++){
175 for (int k=0;k<numvertices;k++)values[k]=str[N*vertexsids[k]+t];
176
177 switch(element->ObjectEnum()){
178 case TriaEnum: transientinput->AddTriaTimeInput( t,numvertices,vertexlids,values,P1Enum); break;
179 case PentaEnum: transientinput->AddPentaTimeInput(t,numvertices,vertexlids,values,P1Enum); break;
180 default: _error_("Not implemented yet");
181 }
182 }
183 xDelete<IssmDouble>(values);
184 xDelete<int>(vertexlids);
185 xDelete<int>(vertexsids);
186 }
187
188 xDelete<IssmDouble>(times);
189 }
190
191 /*Delete data:*/
192 for(int i=0;i<nummodels;i++){
193 IssmDouble* str=pstr[i];
194 xDelete<IssmDouble>(str);
195 }
196 xDelete<IssmDouble*>(pstr);
197 xDelete<int>(pM);
198 xDelete<int>(pN);
199 /*}}}*/
200 /*now do the same with the dsl.sea_water_pressure_change_at_sea_floor field:{{{ */
201 iomodel->FetchData(&pstr,&pM,&pN,&nummodels,"md.dsl.sea_water_pressure_change_at_sea_floor");
202
203 for (int i=0;i<nummodels;i++){
204 M=pM[i];
205 N=pN[i];
206 str=pstr[i];
207
208 //recover time vector:
209 times=xNew<IssmDouble>(N);
210 for(int t=0;t<N;t++) times[t] = str[(M-1)*N+t];
211
212 TransientInput2* transientinput=inputs2->SetDatasetTransientInput(DslSeaWaterPressureChangeAtSeaFloor,i, times,N);
213
214 for(int j=0;j<elements->Size();j++){
215
216 /*Get the right transient input*/
217 Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(j));
218
219 /*Get values and lid list*/
220 const int numvertices = element->GetNumberOfVertices();
221 int *vertexlids = xNew<int>(numvertices);
222 int *vertexsids = xNew<int>(numvertices);
223
224
225 /*Recover vertices ids needed to initialize inputs*/
226 _assert_(iomodel->elements);
227 for(int k=0;k<numvertices;k++){
228 vertexsids[k] =reCast<int>(iomodel->elements[numvertices*element->Sid()+k]); //ids for vertices are in the elements array from Matlab
229 vertexlids[k]=iomodel->my_vertices_lids[vertexsids[k]-1];
230 }
231 //element->GetVerticesLidList(vertexlids);
232 //element->GetVerticesSidList(vertexsids);
233
234 IssmDouble* values=xNew<IssmDouble>(numvertices);
235
236 for(int t=0;t<N;t++){
237 for (int k=0;k<numvertices;k++)values[k]=str[N*vertexsids[k]+t];
238
239 switch(element->ObjectEnum()){
240 case TriaEnum: transientinput->AddTriaTimeInput( t,numvertices,vertexlids,values,P1Enum); break;
241 case PentaEnum: transientinput->AddPentaTimeInput(t,numvertices,vertexlids,values,P1Enum); break;
242 default: _error_("Not implemented yet");
243 }
244 }
245 xDelete<IssmDouble>(values);
246 xDelete<int>(vertexlids);
247 xDelete<int>(vertexsids);
248 }
249 xDelete<IssmDouble>(times);
250 }
251
252 /*Delete data:*/
253 for(int i=0;i<nummodels;i++){
254 IssmDouble* str=pstr[i];
255 xDelete<IssmDouble>(str);
256 }
257 xDelete<IssmDouble*>(pstr);
258 xDelete<int>(pM);
259 xDelete<int>(pN);
260 /*}}}*/
261
262 } /*}}}*/
[24469]263 else _error_("Dsl model " << dslmodel << " not implemented yet!");
[25060]264
[25066]265
[25060]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;
[25118]279 int lovemodel=1;
[19984]280 IssmDouble* love_h=NULL;
281 IssmDouble* love_k=NULL;
[20999]282 IssmDouble* love_l=NULL;
[25139]283 IssmDouble* love_th=NULL;
284 IssmDouble* love_tk=NULL;
285 IssmDouble* love_tl=NULL;
[24469]286 int dslmodel=0;
[23066]287
[24921]288 IssmDouble* G_rigid = NULL;
289 IssmDouble* G_rigid_local = NULL;
[20028]290 IssmDouble* G_elastic = NULL;
[20033]291 IssmDouble* G_elastic_local = NULL;
[20999]292 IssmDouble* U_elastic = NULL;
293 IssmDouble* U_elastic_local = NULL;
294 IssmDouble* H_elastic = NULL;
295 IssmDouble* H_elastic_local = NULL;
[20033]296 int M,m,lower_row,upper_row;
297 IssmDouble degacc=.01;
[24959]298 IssmDouble planetradius=0;
299 IssmDouble planetarea=0;
[19984]300
[20036]301 int numoutputs;
302 char** requestedoutputs = NULL;
303
[20136]304 /*transition vectors: */
305 IssmDouble **transitions = NULL;
306 int *transitions_M = NULL;
307 int *transitions_N = NULL;
308 int ntransitions;
309
[19986]310 /*some constant parameters: */
[24531]311 parameters->AddObject(iomodel->CopyConstantObject("md.dsl.model",DslModelEnum));
[25118]312 parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.runfrequency",SolidearthSettingsRunFrequencyEnum));
313 parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.reltol",SolidearthSettingsReltolEnum));
314 parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.abstol",SolidearthSettingsAbstolEnum));
315 parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.maxiter",SolidearthSettingsMaxiterEnum));
316 parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.rigid",SolidearthSettingsRigidEnum));
317 parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.horiz",SolidearthSettingsHorizEnum));
318 parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.elastic",SolidearthSettingsElasticEnum));
319 parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.rotation",SolidearthSettingsRotationEnum));
320 parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.rotational.equatorialmoi",RotationalEquatorialMoiEnum));
321 parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.rotational.polarmoi",RotationalPolarMoiEnum));
322 parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.rotational.angularvelocity",RotationalAngularVelocityEnum));
323 parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.ocean_area_scaling",SolidearthSettingsOceanAreaScalingEnum));
324 parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.computesealevelchange",SolidearthSettingsComputesealevelchangeEnum));
325 parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.planetradius",SolidearthPlanetRadiusEnum));
[19986]326
[24959]327 /*compute planet area and plug into parameters:*/
[25118]328 iomodel->FetchData(&planetradius,"md.solidearth.planetradius");
[24959]329 planetarea=4*PI*planetradius*planetradius;
[25118]330 parameters->AddObject(new DoubleParam(SolidearthPlanetAreaEnum,planetarea));
[24959]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: */
[25146]350 iomodel->FetchData(&love_h,&nl,NULL,"md.solidearth.lovenumbers.h");
351 iomodel->FetchData(&love_k,&nl,NULL,"md.solidearth.lovenumbers.k");
352 iomodel->FetchData(&love_l,&nl,NULL,"md.solidearth.lovenumbers.l");
353 iomodel->FetchData(&love_th,&nl,NULL,"md.solidearth.lovenumbers.th");
354 iomodel->FetchData(&love_tk,&nl,NULL,"md.solidearth.lovenumbers.tk");
355 iomodel->FetchData(&love_tl,&nl,NULL,"md.solidearth.lovenumbers.tl");
356 parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.lovenumbers.tk2secular",TidalLoveK2SecularEnum));
[19984]357
[25139]358 parameters->AddObject(new DoubleMatParam(LoadLoveHEnum,love_h,nl,1));
359 parameters->AddObject(new DoubleMatParam(LoadLoveKEnum,love_k,nl,1));
360 parameters->AddObject(new DoubleMatParam(LoadLoveLEnum,love_l,nl,1));
361 parameters->AddObject(new DoubleMatParam(TidalLoveHEnum,love_th,nl,1));
362 parameters->AddObject(new DoubleMatParam(TidalLoveKEnum,love_tk,nl,1));
363 parameters->AddObject(new DoubleMatParam(TidalLoveLEnum,love_tl,nl,1));
364
[24921]365 /*compute elastic green function for a range of angles*/
[25118]366 iomodel->FetchData(&degacc,"md.solidearth.settings.degacc");
[24921]367 M=reCast<int,IssmDouble>(180./degacc+1.);
[20022]368
[24921]369 // AD performance is sensitive to calls to ensurecontiguous.
370 // // Providing "t" will cause ensurecontiguous to be called.
371 #ifdef _HAVE_AD_
372 G_rigid=xNew<IssmDouble>(M,"t");
373 G_elastic=xNew<IssmDouble>(M,"t");
374 U_elastic=xNew<IssmDouble>(M,"t");
375 H_elastic=xNew<IssmDouble>(M,"t");
376 #else
377 G_rigid=xNew<IssmDouble>(M);
378 G_elastic=xNew<IssmDouble>(M);
379 U_elastic=xNew<IssmDouble>(M);
380 H_elastic=xNew<IssmDouble>(M);
381 #endif
[21742]382
[24921]383 /*compute combined legendre + love number (elastic green function:*/
384 m=DetermineLocalSize(M,IssmComm::GetComm());
385 GetOwnershipBoundariesFromRange(&lower_row,&upper_row,m,IssmComm::GetComm());
386 // AD performance is sensitive to calls to ensurecontiguous.
387 // // Providing "t" will cause ensurecontiguous to be called.
388 #ifdef _HAVE_AD_
389 G_elastic_local=xNew<IssmDouble>(m,"t");
390 G_rigid_local=xNew<IssmDouble>(m,"t");
391 U_elastic_local=xNew<IssmDouble>(m,"t");
392 H_elastic_local=xNew<IssmDouble>(m,"t");
393 #else
394 G_elastic_local=xNew<IssmDouble>(m);
395 G_rigid_local=xNew<IssmDouble>(m);
396 U_elastic_local=xNew<IssmDouble>(m);
397 H_elastic_local=xNew<IssmDouble>(m);
398 #endif
[21742]399
[24921]400 for(int i=lower_row;i<upper_row;i++){
401 IssmDouble alpha,x;
402 alpha= reCast<IssmDouble>(i)*degacc * PI / 180.0;
[20033]403
[24921]404 G_rigid_local[i-lower_row]= .5/sin(alpha/2.0);
405 G_elastic_local[i-lower_row]= (love_k[nl-1]-love_h[nl-1])*G_rigid_local[i-lower_row];
406 U_elastic_local[i-lower_row]= (love_h[nl-1])*G_rigid_local[i-lower_row];
407 H_elastic_local[i-lower_row]= 0;
408 IssmDouble Pn = 0.;
409 IssmDouble Pn1 = 0.;
410 IssmDouble Pn2 = 0.;
411 IssmDouble Pn_p = 0.;
412 IssmDouble Pn_p1 = 0.;
413 IssmDouble Pn_p2 = 0.;
[20029]414
[24921]415 for (int n=0;n<nl;n++) {
416 IssmDouble deltalove_G;
417 IssmDouble deltalove_U;
[21908]418
[24921]419 deltalove_G = (love_k[n]-love_k[nl-1]-love_h[n]+love_h[nl-1]);
420 deltalove_U = (love_h[n]-love_h[nl-1]);
[20033]421
[24921]422 /*compute legendre polynomials: P_n(cos\theta) & d P_n(cos\theta)/ d\theta: */
423 if(n==0){
424 Pn=1;
425 Pn_p=0;
426 }
427 else if(n==1){
428 Pn = cos(alpha);
429 Pn_p = 1;
430 }
431 else{
432 Pn = ( (2*n-1)*cos(alpha)*Pn1 - (n-1)*Pn2 ) /n;
433 Pn_p = ( (2*n-1)*(Pn1+cos(alpha)*Pn_p1) - (n-1)*Pn_p2 ) /n;
434 }
435 Pn2=Pn1; Pn1=Pn;
436 Pn_p2=Pn_p1; Pn_p1=Pn_p;
[23066]437
[24921]438 G_elastic_local[i-lower_row] += deltalove_G*Pn; // gravitational potential
439 U_elastic_local[i-lower_row] += deltalove_U*Pn; // vertical (up) displacement
440 H_elastic_local[i-lower_row] += sin(alpha)*love_l[n]*Pn_p; // horizontal displacements
[20033]441 }
[24921]442 }
[20033]443
[24921]444 /*merge G_elastic_local into G_elastic; U_elastic_local into U_elastic; H_elastic_local to H_elastic:{{{*/
445 int* recvcounts=xNew<int>(IssmComm::GetSize());
446 int* displs=xNew<int>(IssmComm::GetSize());
[20033]447
[24921]448 //recvcounts:
449 ISSM_MPI_Allgather(&m,1,ISSM_MPI_INT,recvcounts,1,ISSM_MPI_INT,IssmComm::GetComm());
[20033]450
[24921]451 /*displs: */
452 ISSM_MPI_Allgather(&lower_row,1,ISSM_MPI_INT,displs,1,ISSM_MPI_INT,IssmComm::GetComm());
[20033]453
[24921]454 /*All gather:*/
455 ISSM_MPI_Allgatherv(G_rigid_local, m, ISSM_MPI_DOUBLE, G_rigid, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
456 ISSM_MPI_Allgatherv(G_elastic_local, m, ISSM_MPI_DOUBLE, G_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
457 ISSM_MPI_Allgatherv(U_elastic_local, m, ISSM_MPI_DOUBLE, U_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
458 ISSM_MPI_Allgatherv(H_elastic_local, m, ISSM_MPI_DOUBLE, H_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
459 /*free ressources: */
460 xDelete<int>(recvcounts);
461 xDelete<int>(displs);
[20033]462
[24921]463 /*}}}*/
[20033]464
[24921]465 /*Avoid singularity at 0: */
466 G_rigid[0]=G_rigid[1];
467 parameters->AddObject(new DoubleVecParam(SealevelriseGRigidEnum,G_rigid,M));
468 G_elastic[0]=G_elastic[1];
469 parameters->AddObject(new DoubleVecParam(SealevelriseGElasticEnum,G_elastic,M));
470 U_elastic[0]=U_elastic[1];
471 parameters->AddObject(new DoubleVecParam(SealevelriseUElasticEnum,U_elastic,M));
472 H_elastic[0]=H_elastic[1];
473 parameters->AddObject(new DoubleVecParam(SealevelriseHElasticEnum,H_elastic,M));
[20028]474
[24921]475 /*free ressources: */
476 xDelete<IssmDouble>(love_h);
477 xDelete<IssmDouble>(love_k);
478 xDelete<IssmDouble>(love_l);
[25139]479 xDelete<IssmDouble>(love_th);
480 xDelete<IssmDouble>(love_tk);
481 xDelete<IssmDouble>(love_tl);
[24921]482 xDelete<IssmDouble>(G_rigid);
483 xDelete<IssmDouble>(G_rigid_local);
484 xDelete<IssmDouble>(G_elastic);
485 xDelete<IssmDouble>(G_elastic_local);
486 xDelete<IssmDouble>(U_elastic);
487 xDelete<IssmDouble>(U_elastic_local);
488 xDelete<IssmDouble>(H_elastic);
489 xDelete<IssmDouble>(H_elastic_local);
[24977]490
491 /*Indicate we have not yet ran the Geometry Core module: */
492 parameters->AddObject(new BoolParam(SealevelriseGeometryDoneEnum,false));
493
[24921]494/*}}}*/
[24481]495 /*Transitions:{{{ */
[25118]496 iomodel->FetchData(&transitions,&transitions_M,&transitions_N,&ntransitions,"md.solidearth.transitions");
[20136]497 if(transitions){
498 parameters->AddObject(new DoubleMatArrayParam(SealevelriseTransitionsEnum,transitions,ntransitions,transitions_M,transitions_N));
499
500 for(int i=0;i<ntransitions;i++){
501 IssmDouble* transition=transitions[i];
502 xDelete<IssmDouble>(transition);
503 }
504 xDelete<IssmDouble*>(transitions);
[20138]505 xDelete<int>(transitions_M);
506 xDelete<int>(transitions_N);
[24481]507 } /*}}}*/
508 /*Requested outputs {{{*/
[25118]509 iomodel->FindConstant(&requestedoutputs,&numoutputs,"md.solidearth.requested_outputs");
[20036]510 if(numoutputs)parameters->AddObject(new StringArrayParam(SealevelriseRequestedOutputsEnum,requestedoutputs,numoutputs));
[25118]511 iomodel->DeleteData(&requestedoutputs,numoutputs,"md.solidearth.requested_outputs");
[24481]512 /*}}}*/
[20029]513
[19984]514}/*}}}*/
515
516/*Finite Element Analysis*/
517void SealevelriseAnalysis::Core(FemModel* femmodel){/*{{{*/
518 _error_("not implemented");
519}/*}}}*/
520ElementVector* SealevelriseAnalysis::CreateDVector(Element* element){/*{{{*/
521 /*Default, return NULL*/
522 return NULL;
523}/*}}}*/
524ElementMatrix* SealevelriseAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/
525_error_("Not implemented");
526}/*}}}*/
527ElementMatrix* SealevelriseAnalysis::CreateKMatrix(Element* element){/*{{{*/
528 _error_("not implemented yet");
529}/*}}}*/
530ElementVector* SealevelriseAnalysis::CreatePVector(Element* element){/*{{{*/
531_error_("not implemented yet");
532}/*}}}*/
533void SealevelriseAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
534 _error_("not implemented yet");
535}/*}}}*/
536void SealevelriseAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
537 _error_("Not implemented yet");
538}/*}}}*/
539void SealevelriseAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
[22955]540 _error_("not implemeneted yet!");
[20153]541
[19984]542}/*}}}*/
543void SealevelriseAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
544 /*Default, do nothing*/
545 return;
546}/*}}}*/
Note: See TracBrowser for help on using the repository browser.