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

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

CHG: took out UpdateStericRate and UpdateDynamicRate methods, not needed anymore.
Using new Inputs approach to speed up the retrieval of steric and dynamic sea-level change rates.
The global mean steric change rate is now a transient input (P0) instead of a transient parameter. This
will allow for future dsl class where we have a multi-model ensemble that dakota can play we at run time!
File size: 13.4 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.groundedice_levelset",MaskGroundediceLevelsetEnum);
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){
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
97 iomodel->FetchDataToInput(inputs2,elements,"md.dsl.sea_surface_height_change_above_geoid", DslSeaSurfaceHeightChangeAboveGeoidEnum);
98 iomodel->FetchDataToInput(inputs2,elements,"md.dsl.sea_water_pressure_change_at_sea_floor", DslSeaWaterPressureChangeAtSeaFloor);
99
100 //TransientInput2* transientinput=inputs2->SetDatasetTransientInput(DslGlobalAverageThermostericSeaLevelChangeEnum,counter, times,N);
101 }
102 else _error_("Dsl model " << dslmodel << " not implemented yet!");
103
104 /*Initialize cumdeltalthickness and sealevel rise rate input*/
105 InputUpdateFromConstantx(inputs2,elements,0.,SealevelriseCumDeltathicknessEnum);
106 InputUpdateFromConstantx(inputs2,elements,0.,SealevelNEsaRateEnum);
107 InputUpdateFromConstantx(inputs2,elements,0.,SealevelUEsaRateEnum);
108 InputUpdateFromConstantx(inputs2,elements,0.,SealevelRSLRateEnum);
109 InputUpdateFromConstantx(inputs2,elements,0.,SealevelEustaticMaskEnum);
110 InputUpdateFromConstantx(inputs2,elements,0.,SealevelEustaticOceanMaskEnum);
111
112}/*}}}*/
113void SealevelriseAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
114
115 int nl;
116 IssmDouble* love_h=NULL;
117 IssmDouble* love_k=NULL;
118 IssmDouble* love_l=NULL;
119 int dslmodel=0;
120
121 bool elastic=false;
122 IssmDouble* G_elastic = NULL;
123 IssmDouble* G_elastic_local = NULL;
124 IssmDouble* U_elastic = NULL;
125 IssmDouble* U_elastic_local = NULL;
126 IssmDouble* H_elastic = NULL;
127 IssmDouble* H_elastic_local = NULL;
128 int M,m,lower_row,upper_row;
129 IssmDouble degacc=.01;
130
131 int numoutputs;
132 char** requestedoutputs = NULL;
133
134 /*transition vectors: */
135 IssmDouble **transitions = NULL;
136 int *transitions_M = NULL;
137 int *transitions_N = NULL;
138 int ntransitions;
139
140 /*some constant parameters: */
141 parameters->AddObject(iomodel->CopyConstantObject("md.slr.geodetic_run_frequency",SealevelriseGeodeticRunFrequencyEnum));
142 parameters->AddObject(iomodel->CopyConstantObject("md.slr.reltol",SealevelriseReltolEnum));
143 parameters->AddObject(iomodel->CopyConstantObject("md.slr.abstol",SealevelriseAbstolEnum));
144 parameters->AddObject(iomodel->CopyConstantObject("md.slr.maxiter",SealevelriseMaxiterEnum));
145 parameters->AddObject(iomodel->CopyConstantObject("md.slr.loop_increment",SealevelriseLoopIncrementEnum));
146 parameters->AddObject(iomodel->CopyConstantObject("md.slr.rigid",SealevelriseRigidEnum));
147 parameters->AddObject(iomodel->CopyConstantObject("md.slr.horiz",SealevelriseHorizEnum));
148 parameters->AddObject(iomodel->CopyConstantObject("md.slr.elastic",SealevelriseElasticEnum));
149 parameters->AddObject(iomodel->CopyConstantObject("md.slr.rotation",SealevelriseRotationEnum));
150 parameters->AddObject(iomodel->CopyConstantObject("md.slr.tide_love_h",SealevelriseTidalLoveHEnum));
151 parameters->AddObject(iomodel->CopyConstantObject("md.slr.tide_love_k",SealevelriseTidalLoveKEnum));
152 parameters->AddObject(iomodel->CopyConstantObject("md.slr.fluid_love",SealevelriseFluidLoveEnum));
153 parameters->AddObject(iomodel->CopyConstantObject("md.slr.equatorial_moi",SealevelriseEquatorialMoiEnum));
154 parameters->AddObject(iomodel->CopyConstantObject("md.slr.polar_moi",SealevelrisePolarMoiEnum));
155 parameters->AddObject(iomodel->CopyConstantObject("md.slr.angular_velocity",SealevelriseAngularVelocityEnum));
156 parameters->AddObject(iomodel->CopyConstantObject("md.slr.ocean_area_scaling",SealevelriseOceanAreaScalingEnum));
157 parameters->AddObject(iomodel->CopyConstantObject("md.slr.geodetic",SealevelriseGeodeticEnum));
158
159
160 /*Deal with elasticity {{{*/
161 iomodel->FetchData(&elastic,"md.slr.elastic");
162 if(elastic){
163
164 /*love numbers: */
165 iomodel->FetchData(&love_h,&nl,NULL,"md.slr.love_h");
166 iomodel->FetchData(&love_k,&nl,NULL,"md.slr.love_k");
167 iomodel->FetchData(&love_l,&nl,NULL,"md.slr.love_l");
168
169 /*compute elastic green function for a range of angles*/
170 iomodel->FetchData(&degacc,"md.slr.degacc");
171 M=reCast<int,IssmDouble>(180./degacc+1.);
172
173 // AD performance is sensitive to calls to ensurecontiguous.
174 // // Providing "t" will cause ensurecontiguous to be called.
175 #ifdef _HAVE_AD_
176 G_elastic=xNew<IssmDouble>(M,"t");
177 U_elastic=xNew<IssmDouble>(M,"t");
178 H_elastic=xNew<IssmDouble>(M,"t");
179 #else
180 G_elastic=xNew<IssmDouble>(M);
181 U_elastic=xNew<IssmDouble>(M);
182 H_elastic=xNew<IssmDouble>(M);
183 #endif
184
185 /*compute combined legendre + love number (elastic green function:*/
186 m=DetermineLocalSize(M,IssmComm::GetComm());
187 GetOwnershipBoundariesFromRange(&lower_row,&upper_row,m,IssmComm::GetComm());
188 // AD performance is sensitive to calls to ensurecontiguous.
189 // // Providing "t" will cause ensurecontiguous to be called.
190 #ifdef _HAVE_AD_
191 G_elastic_local=xNew<IssmDouble>(m,"t");
192 U_elastic_local=xNew<IssmDouble>(m,"t");
193 H_elastic_local=xNew<IssmDouble>(m,"t");
194 #else
195 G_elastic_local=xNew<IssmDouble>(m);
196 U_elastic_local=xNew<IssmDouble>(m);
197 H_elastic_local=xNew<IssmDouble>(m);
198 #endif
199
200 for(int i=lower_row;i<upper_row;i++){
201 IssmDouble alpha,x;
202 alpha= reCast<IssmDouble>(i)*degacc * PI / 180.0;
203
204 G_elastic_local[i-lower_row]= (love_k[nl-1]-love_h[nl-1])/2.0/sin(alpha/2.0);
205 U_elastic_local[i-lower_row]= (love_h[nl-1])/2.0/sin(alpha/2.0);
206 H_elastic_local[i-lower_row]= 0;
207 IssmDouble Pn = 0.;
208 IssmDouble Pn1 = 0.;
209 IssmDouble Pn2 = 0.;
210 IssmDouble Pn_p = 0.;
211 IssmDouble Pn_p1 = 0.;
212 IssmDouble Pn_p2 = 0.;
213
214 for (int n=0;n<nl;n++) {
215 IssmDouble deltalove_G;
216 IssmDouble deltalove_U;
217
218 deltalove_G = (love_k[n]-love_k[nl-1]-love_h[n]+love_h[nl-1]);
219 deltalove_U = (love_h[n]-love_h[nl-1]);
220
221 /*compute legendre polynomials: P_n(cos\theta) & d P_n(cos\theta)/ d\theta: */
222 if(n==0){
223 Pn=1;
224 Pn_p=0;
225 }
226 else if(n==1){
227 Pn = cos(alpha);
228 Pn_p = 1;
229 }
230 else{
231 Pn = ( (2*n-1)*cos(alpha)*Pn1 - (n-1)*Pn2 ) /n;
232 Pn_p = ( (2*n-1)*(Pn1+cos(alpha)*Pn_p1) - (n-1)*Pn_p2 ) /n;
233 }
234 Pn2=Pn1; Pn1=Pn;
235 Pn_p2=Pn_p1; Pn_p1=Pn_p;
236
237 G_elastic_local[i-lower_row] += deltalove_G*Pn; // gravitational potential
238 U_elastic_local[i-lower_row] += deltalove_U*Pn; // vertical (up) displacement
239 H_elastic_local[i-lower_row] += sin(alpha)*love_l[n]*Pn_p; // horizontal displacements
240 }
241 }
242
243 /*merge G_elastic_local into G_elastic; U_elastic_local into U_elastic; H_elastic_local to H_elastic:{{{*/
244 int* recvcounts=xNew<int>(IssmComm::GetSize());
245 int* displs=xNew<int>(IssmComm::GetSize());
246
247 //recvcounts:
248 ISSM_MPI_Allgather(&m,1,ISSM_MPI_INT,recvcounts,1,ISSM_MPI_INT,IssmComm::GetComm());
249
250 /*displs: */
251 ISSM_MPI_Allgather(&lower_row,1,ISSM_MPI_INT,displs,1,ISSM_MPI_INT,IssmComm::GetComm());
252
253 /*All gather:*/
254 ISSM_MPI_Allgatherv(G_elastic_local, m, ISSM_MPI_DOUBLE, G_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
255 ISSM_MPI_Allgatherv(U_elastic_local, m, ISSM_MPI_DOUBLE, U_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
256 ISSM_MPI_Allgatherv(H_elastic_local, m, ISSM_MPI_DOUBLE, H_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
257 /*free ressources: */
258 xDelete<int>(recvcounts);
259 xDelete<int>(displs);
260
261 /*}}}*/
262
263 /*Avoid singularity at 0: */
264 G_elastic[0]=G_elastic[1];
265 parameters->AddObject(new DoubleVecParam(SealevelriseGElasticEnum,G_elastic,M));
266 U_elastic[0]=U_elastic[1];
267 parameters->AddObject(new DoubleVecParam(SealevelriseUElasticEnum,U_elastic,M));
268 H_elastic[0]=H_elastic[1];
269 parameters->AddObject(new DoubleVecParam(SealevelriseHElasticEnum,H_elastic,M));
270
271 /*free ressources: */
272 xDelete<IssmDouble>(love_h);
273 xDelete<IssmDouble>(love_k);
274 xDelete<IssmDouble>(love_l);
275 xDelete<IssmDouble>(G_elastic);
276 xDelete<IssmDouble>(G_elastic_local);
277 xDelete<IssmDouble>(U_elastic);
278 xDelete<IssmDouble>(U_elastic_local);
279 xDelete<IssmDouble>(H_elastic);
280 xDelete<IssmDouble>(H_elastic_local);
281 } /*}}}*/
282 /*Transitions:{{{ */
283 iomodel->FetchData(&transitions,&transitions_M,&transitions_N,&ntransitions,"md.slr.transitions");
284 if(transitions){
285 parameters->AddObject(new DoubleMatArrayParam(SealevelriseTransitionsEnum,transitions,ntransitions,transitions_M,transitions_N));
286
287 for(int i=0;i<ntransitions;i++){
288 IssmDouble* transition=transitions[i];
289 xDelete<IssmDouble>(transition);
290 }
291 xDelete<IssmDouble*>(transitions);
292 xDelete<int>(transitions_M);
293 xDelete<int>(transitions_N);
294 } /*}}}*/
295 /*Requested outputs {{{*/
296 iomodel->FindConstant(&requestedoutputs,&numoutputs,"md.slr.requested_outputs");
297 if(numoutputs)parameters->AddObject(new StringArrayParam(SealevelriseRequestedOutputsEnum,requestedoutputs,numoutputs));
298 iomodel->DeleteData(&requestedoutputs,numoutputs,"md.slr.requested_outputs");
299 /*}}}*/
300
301}/*}}}*/
302
303/*Finite Element Analysis*/
304void SealevelriseAnalysis::Core(FemModel* femmodel){/*{{{*/
305 _error_("not implemented");
306}/*}}}*/
307ElementVector* SealevelriseAnalysis::CreateDVector(Element* element){/*{{{*/
308 /*Default, return NULL*/
309 return NULL;
310}/*}}}*/
311ElementMatrix* SealevelriseAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/
312_error_("Not implemented");
313}/*}}}*/
314ElementMatrix* SealevelriseAnalysis::CreateKMatrix(Element* element){/*{{{*/
315 _error_("not implemented yet");
316}/*}}}*/
317ElementVector* SealevelriseAnalysis::CreatePVector(Element* element){/*{{{*/
318_error_("not implemented yet");
319}/*}}}*/
320void SealevelriseAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
321 _error_("not implemented yet");
322}/*}}}*/
323void SealevelriseAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
324 _error_("Not implemented yet");
325}/*}}}*/
326void SealevelriseAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
327 _error_("not implemeneted yet!");
328
329}/*}}}*/
330void SealevelriseAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
331 /*Default, do nothing*/
332 return;
333}/*}}}*/
Note: See TracBrowser for help on using the repository browser.