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

Last change on this file since 22609 was 22609, checked in by Mathieu Morlighem, 7 years ago

CHG: fixing some enums
File size: 10.4 KB
Line 
1#include "./SealevelriseAnalysis.h"
2#include "../toolkits/toolkits.h"
3#include "../classes/classes.h"
4#include "../shared/shared.h"
5#include "../modules/modules.h"
6
7/*Model processing*/
8void SealevelriseAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
9 /*No constraints*/
10}/*}}}*/
11void SealevelriseAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
12 /*No loads*/
13}/*}}}*/
14void SealevelriseAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
15 ::CreateNodes(nodes,iomodel,SealevelriseAnalysisEnum,P1Enum);
16}/*}}}*/
17int SealevelriseAnalysis::DofsPerNode(int** doflist,int domaintype,int approximation){/*{{{*/
18 return 1;
19}/*}}}*/
20void SealevelriseAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
21
22 /*Update elements: */
23 int counter=0;
24 for(int i=0;i<iomodel->numberofelements;i++){
25 if(iomodel->my_elements[i]){
26 Element* element=(Element*)elements->GetObjectByOffset(counter);
27 element->Update(i,iomodel,analysis_counter,analysis_type,P1Enum);
28 counter++;
29 }
30 }
31
32 /*Create inputs: */
33 iomodel->FetchDataToInput(elements,"md.mask.ice_levelset",MaskIceLevelsetEnum);
34 iomodel->FetchDataToInput(elements,"md.mask.ocean_levelset",MaskOceanLevelsetEnum);
35 iomodel->FetchDataToInput(elements,"md.mask.land_levelset",MaskLandLevelsetEnum);
36 iomodel->FetchDataToInput(elements,"md.slr.deltathickness",SealevelriseDeltathicknessEnum);
37 iomodel->FetchDataToInput(elements,"md.slr.sealevel",SealevelEnum,0);
38 iomodel->FetchDataToInput(elements,"md.slr.steric_rate",SealevelriseStericRateEnum);
39
40}/*}}}*/
41void SealevelriseAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
42
43 int nl;
44 IssmDouble* love_h=NULL;
45 IssmDouble* love_k=NULL;
46 IssmDouble* love_l=NULL;
47
48 bool elastic=false;
49 IssmDouble* G_elastic = NULL;
50 IssmDouble* G_elastic_local = NULL;
51 IssmDouble* U_elastic = NULL;
52 IssmDouble* U_elastic_local = NULL;
53 IssmDouble* H_elastic = NULL;
54 IssmDouble* H_elastic_local = NULL;
55 int M,m,lower_row,upper_row;
56 IssmDouble degacc=.01;
57
58 int numoutputs;
59 char** requestedoutputs = NULL;
60
61 /*transition vectors: */
62 IssmDouble **transitions = NULL;
63 int *transitions_M = NULL;
64 int *transitions_N = NULL;
65 int ntransitions;
66
67 /*some constant parameters: */
68 parameters->AddObject(iomodel->CopyConstantObject("md.slr.reltol",SealevelriseReltolEnum));
69 parameters->AddObject(iomodel->CopyConstantObject("md.slr.abstol",SealevelriseAbstolEnum));
70 parameters->AddObject(iomodel->CopyConstantObject("md.slr.maxiter",SealevelriseMaxiterEnum));
71 parameters->AddObject(iomodel->CopyConstantObject("md.slr.rigid",SealevelriseRigidEnum));
72 parameters->AddObject(iomodel->CopyConstantObject("md.slr.elastic",SealevelriseElasticEnum));
73 parameters->AddObject(iomodel->CopyConstantObject("md.slr.rotation",SealevelriseRotationEnum));
74 parameters->AddObject(iomodel->CopyConstantObject("md.slr.tide_love_h",SealevelriseTidalLoveHEnum));
75 parameters->AddObject(iomodel->CopyConstantObject("md.slr.tide_love_k",SealevelriseTidalLoveKEnum));
76 parameters->AddObject(iomodel->CopyConstantObject("md.slr.fluid_love",SealevelriseFluidLoveEnum));
77 parameters->AddObject(iomodel->CopyConstantObject("md.slr.equatorial_moi",SealevelriseEquatorialMoiEnum));
78 parameters->AddObject(iomodel->CopyConstantObject("md.slr.polar_moi",SealevelrisePolarMoiEnum));
79 parameters->AddObject(iomodel->CopyConstantObject("md.slr.angular_velocity",SealevelriseAngularVelocityEnum));
80 parameters->AddObject(iomodel->CopyConstantObject("md.slr.ocean_area_scaling",SealevelriseOceanAreaScalingEnum));
81
82 iomodel->FetchData(&elastic,"md.slr.elastic");
83 if(elastic){
84
85 /*love numbers: */
86 iomodel->FetchData(&love_h,&nl,NULL,"md.slr.love_h");
87 iomodel->FetchData(&love_k,&nl,NULL,"md.slr.love_k");
88 iomodel->FetchData(&love_l,&nl,NULL,"md.slr.love_l");
89
90 /*compute elastic green function for a range of angles*/
91 iomodel->FetchData(&degacc,"md.slr.degacc");
92 M=reCast<int,IssmDouble>(180./degacc+1.);
93
94 // AD performance is sensitive to calls to ensurecontiguous.
95 // // Providing "t" will cause ensurecontiguous to be called.
96 #ifdef _HAVE_AD_
97 G_elastic=xNew<IssmDouble>(M,"t");
98 U_elastic=xNew<IssmDouble>(M,"t");
99 H_elastic=xNew<IssmDouble>(M,"t");
100 #else
101 G_elastic=xNew<IssmDouble>(M);
102 U_elastic=xNew<IssmDouble>(M);
103 H_elastic=xNew<IssmDouble>(M);
104 #endif
105
106
107 /*compute combined legendre + love number (elastic green function:*/
108 m=DetermineLocalSize(M,IssmComm::GetComm());
109 GetOwnershipBoundariesFromRange(&lower_row,&upper_row,m,IssmComm::GetComm());
110 // AD performance is sensitive to calls to ensurecontiguous.
111 // // Providing "t" will cause ensurecontiguous to be called.
112 #ifdef _HAVE_AD_
113 G_elastic_local=xNew<IssmDouble>(m,"t");
114 U_elastic_local=xNew<IssmDouble>(m,"t");
115 H_elastic_local=xNew<IssmDouble>(m,"t");
116 #else
117 G_elastic_local=xNew<IssmDouble>(m);
118 U_elastic_local=xNew<IssmDouble>(m);
119 H_elastic_local=xNew<IssmDouble>(m);
120 #endif
121
122 for(int i=lower_row;i<upper_row;i++){
123 IssmDouble alpha,x;
124 alpha= reCast<IssmDouble>(i)*degacc * PI / 180.0;
125
126 G_elastic_local[i-lower_row]= (love_k[nl-1]-love_h[nl-1])/2.0/sin(alpha/2.0);
127 U_elastic_local[i-lower_row]= (love_h[nl-1])/2.0/sin(alpha/2.0);
128 H_elastic_local[i-lower_row]= 0;
129 IssmDouble Pn = 0.;
130 IssmDouble Pn1 = 0.;
131 IssmDouble Pn2 = 0.;
132 IssmDouble Pn_p = 0.;
133 IssmDouble Pn_p1 = 0.;
134 IssmDouble Pn_p2 = 0.;
135
136 for (int n=0;n<nl;n++) {
137 IssmDouble deltalove_G;
138 IssmDouble deltalove_U;
139
140 deltalove_G = (love_k[n]-love_k[nl-1]-love_h[n]+love_h[nl-1]);
141 deltalove_U = (love_h[n]-love_h[nl-1]);
142
143 /*compute legendre polynomials: P_n(cos\theta) & d P_n(cos\theta)/ d\theta: */
144 if(n==0){
145 Pn=1;
146 Pn_p=0;
147 }
148 else if(n==1){
149 Pn = cos(alpha);
150 Pn_p = 1;
151 }
152 else{
153 Pn = ( (2*n-1)*cos(alpha)*Pn1 - (n-1)*Pn2 ) /n;
154 Pn_p = ( (2*n-1)*(Pn1+cos(alpha)*Pn_p1) - (n-1)*Pn_p2 ) /n;
155 }
156 Pn2=Pn1; Pn1=Pn;
157 Pn_p2=Pn_p1; Pn_p1=Pn_p;
158
159 G_elastic_local[i-lower_row] += deltalove_G*Pn; // gravitational potential
160 U_elastic_local[i-lower_row] += deltalove_U*Pn; // vertical (up) displacement
161 H_elastic_local[i-lower_row] += sin(alpha)*love_l[n]*Pn_p; // horizontal displacements
162 }
163 }
164
165 /*merge G_elastic_local into G_elastic; U_elastic_local into U_elastic; H_elastic_local to H_elastic:{{{*/
166 int* recvcounts=xNew<int>(IssmComm::GetSize());
167 int* displs=xNew<int>(IssmComm::GetSize());
168
169 //recvcounts:
170 ISSM_MPI_Allgather(&m,1,ISSM_MPI_INT,recvcounts,1,ISSM_MPI_INT,IssmComm::GetComm());
171
172 /*displs: */
173 ISSM_MPI_Allgather(&lower_row,1,ISSM_MPI_INT,displs,1,ISSM_MPI_INT,IssmComm::GetComm());
174
175 /*All gather:*/
176 ISSM_MPI_Allgatherv(G_elastic_local, m, ISSM_MPI_DOUBLE, G_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
177 ISSM_MPI_Allgatherv(U_elastic_local, m, ISSM_MPI_DOUBLE, U_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
178 ISSM_MPI_Allgatherv(H_elastic_local, m, ISSM_MPI_DOUBLE, H_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
179 /*free ressources: */
180 xDelete<int>(recvcounts);
181 xDelete<int>(displs);
182
183 /*}}}*/
184
185 /*Avoid singularity at 0: */
186 G_elastic[0]=G_elastic[1];
187 parameters->AddObject(new DoubleVecParam(SealevelriseGElasticEnum,G_elastic,M));
188 U_elastic[0]=U_elastic[1];
189 parameters->AddObject(new DoubleVecParam(SealevelriseUElasticEnum,U_elastic,M));
190 H_elastic[0]=H_elastic[1];
191 parameters->AddObject(new DoubleVecParam(SealevelriseHElasticEnum,H_elastic,M));
192
193 /*free ressources: */
194 xDelete<IssmDouble>(love_h);
195 xDelete<IssmDouble>(love_k);
196 xDelete<IssmDouble>(love_l);
197 xDelete<IssmDouble>(G_elastic);
198 xDelete<IssmDouble>(G_elastic_local);
199 xDelete<IssmDouble>(U_elastic);
200 xDelete<IssmDouble>(U_elastic_local);
201 xDelete<IssmDouble>(H_elastic);
202 xDelete<IssmDouble>(H_elastic_local);
203 }
204
205 /*Transitions: */
206 iomodel->FetchData(&transitions,&transitions_M,&transitions_N,&ntransitions,"md.slr.transitions");
207 if(transitions){
208 parameters->AddObject(new DoubleMatArrayParam(SealevelriseTransitionsEnum,transitions,ntransitions,transitions_M,transitions_N));
209
210 for(int i=0;i<ntransitions;i++){
211 IssmDouble* transition=transitions[i];
212 xDelete<IssmDouble>(transition);
213 }
214 xDelete<IssmDouble*>(transitions);
215 xDelete<int>(transitions_M);
216 xDelete<int>(transitions_N);
217 }
218
219 /*Requested outputs*/
220 iomodel->FindConstant(&requestedoutputs,&numoutputs,"md.slr.requested_outputs");
221 if(numoutputs)parameters->AddObject(new StringArrayParam(SealevelriseRequestedOutputsEnum,requestedoutputs,numoutputs));
222 iomodel->DeleteData(&requestedoutputs,numoutputs,"md.slr.requested_outputs");
223
224
225}/*}}}*/
226
227/*Finite Element Analysis*/
228void SealevelriseAnalysis::Core(FemModel* femmodel){/*{{{*/
229 _error_("not implemented");
230}/*}}}*/
231ElementVector* SealevelriseAnalysis::CreateDVector(Element* element){/*{{{*/
232 /*Default, return NULL*/
233 return NULL;
234}/*}}}*/
235ElementMatrix* SealevelriseAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/
236_error_("Not implemented");
237}/*}}}*/
238ElementMatrix* SealevelriseAnalysis::CreateKMatrix(Element* element){/*{{{*/
239 _error_("not implemented yet");
240}/*}}}*/
241ElementVector* SealevelriseAnalysis::CreatePVector(Element* element){/*{{{*/
242_error_("not implemented yet");
243}/*}}}*/
244void SealevelriseAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
245 _error_("not implemented yet");
246}/*}}}*/
247void SealevelriseAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
248 _error_("Not implemented yet");
249}/*}}}*/
250void SealevelriseAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
251
252 IssmDouble *deltaS = NULL;
253 IssmDouble *S = NULL;
254 int* sidlist = NULL;
255 int numvertices;
256
257 numvertices= element->GetNumberOfVertices();
258 sidlist=xNew<int>(numvertices);
259
260 element->GetVerticesSidList(sidlist);
261
262 deltaS = xNew<IssmDouble>(numvertices);
263 for(int i=0;i<numvertices;i++){
264 deltaS[i]=solution[sidlist[i]];
265 }
266
267 S = xNew<IssmDouble>(numvertices);
268 element->GetInputListOnVertices(S,SealevelEnum,0);
269
270 /*Add deltaS to S:*/
271 for (int i=0;i<numvertices;i++)S[i]+=deltaS[i];
272
273 /*Add S back into inputs: */
274 element->AddInput(SealevelEnum,S,P1Enum);
275
276 /*Free ressources:*/
277 xDelete<int>(sidlist);
278 xDelete<IssmDouble>(deltaS);
279 xDelete<IssmDouble>(S);
280
281}/*}}}*/
282void SealevelriseAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
283 /*Default, do nothing*/
284 return;
285}/*}}}*/
Note: See TracBrowser for help on using the repository browser.