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

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

CHG: moving some stuff to analyses
File size: 11.7 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 int geodetic=0;
23
24 /*Update elements: */
25 int counter=0;
26 for(int i=0;i<iomodel->numberofelements;i++){
27 if(iomodel->my_elements[i]){
28 Element* element=(Element*)elements->GetObjectByOffset(counter);
29 element->Update(i,iomodel,analysis_counter,analysis_type,P1Enum);
30 counter++;
31 }
32 }
33
34 /*Create inputs: */
35 iomodel->FetchDataToInput(elements,"md.mask.groundedice_levelset",MaskGroundediceLevelsetEnum);
36 iomodel->FetchDataToInput(elements,"md.mask.ice_levelset",MaskIceLevelsetEnum);
37 //those only if we have requested geodetic computations:
38 iomodel->FetchData(&geodetic,"md.slr.geodetic");
39 if (geodetic){
40 char* masktype=NULL;
41 iomodel->FetchData(&masktype,"md.mask.type");
42 if (strcmp(masktype,"maskpsl")==0){
43 iomodel->FetchDataToInput(elements,"md.mask.ocean_levelset",MaskOceanLevelsetEnum);
44 iomodel->FetchDataToInput(elements,"md.mask.land_levelset",MaskLandLevelsetEnum);
45 }
46 xDelete<char>(masktype);
47 }
48 iomodel->FetchDataToInput(elements,"md.slr.deltathickness",SealevelriseDeltathicknessEnum);
49 iomodel->FetchDataToInput(elements,"md.slr.spcthickness",SealevelriseSpcthicknessEnum);
50 iomodel->FetchDataToInput(elements,"md.slr.sealevel",SealevelEnum,0);
51 iomodel->FetchDataToInput(elements,"md.geometry.bed",BedEnum);
52 iomodel->FetchDataToInput(elements,"md.slr.Ngia",SealevelNGiaRateEnum);
53 iomodel->FetchDataToInput(elements,"md.slr.Ugia",SealevelUGiaRateEnum);
54
55 /*Initialize cumdeltalthickness and sealevel rise rate input: unfortunately, we don't have femmodel, so we
56 * need to iterate on elements, otherwise we would have used InputUpdateFromConstantx: */
57 counter=0;
58 for(int i=0;i<iomodel->numberofelements;i++){
59 if(iomodel->my_elements[i]){
60 Element* element=(Element*)elements->GetObjectByOffset(counter);
61 element->InputUpdateFromConstant(0.0,SealevelriseCumDeltathicknessEnum);
62 element->InputUpdateFromConstant(0.0,SealevelNEsaRateEnum);
63 element->InputUpdateFromConstant(0.0,SealevelUEsaRateEnum);
64 element->InputUpdateFromConstant(0.0,SealevelRSLRateEnum);
65 element->InputUpdateFromConstant(0.0,SealevelEustaticMaskEnum);
66 element->InputUpdateFromConstant(0.0,SealevelEustaticOceanMaskEnum);
67 counter++;
68 }
69 }
70
71 iomodel->FetchDataToInput(elements,"md.slr.steric_rate",SealevelriseStericRateEnum);
72
73}/*}}}*/
74void SealevelriseAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
75
76 int nl;
77 IssmDouble* love_h=NULL;
78 IssmDouble* love_k=NULL;
79 IssmDouble* love_l=NULL;
80
81 bool elastic=false;
82 IssmDouble* G_elastic = NULL;
83 IssmDouble* G_elastic_local = NULL;
84 IssmDouble* U_elastic = NULL;
85 IssmDouble* U_elastic_local = NULL;
86 IssmDouble* H_elastic = NULL;
87 IssmDouble* H_elastic_local = NULL;
88 int M,m,lower_row,upper_row;
89 IssmDouble degacc=.01;
90
91 int numoutputs;
92 char** requestedoutputs = NULL;
93
94 /*transition vectors: */
95 IssmDouble **transitions = NULL;
96 int *transitions_M = NULL;
97 int *transitions_N = NULL;
98 int ntransitions;
99
100 /*some constant parameters: */
101 parameters->AddObject(iomodel->CopyConstantObject("md.slr.geodetic_run_frequency",SealevelriseGeodeticRunFrequencyEnum));
102 parameters->AddObject(iomodel->CopyConstantObject("md.slr.reltol",SealevelriseReltolEnum));
103 parameters->AddObject(iomodel->CopyConstantObject("md.slr.abstol",SealevelriseAbstolEnum));
104 parameters->AddObject(iomodel->CopyConstantObject("md.slr.maxiter",SealevelriseMaxiterEnum));
105 parameters->AddObject(iomodel->CopyConstantObject("md.slr.loop_increment",SealevelriseLoopIncrementEnum));
106 parameters->AddObject(iomodel->CopyConstantObject("md.slr.rigid",SealevelriseRigidEnum));
107 parameters->AddObject(iomodel->CopyConstantObject("md.slr.horiz",SealevelriseHorizEnum));
108 parameters->AddObject(iomodel->CopyConstantObject("md.slr.elastic",SealevelriseElasticEnum));
109 parameters->AddObject(iomodel->CopyConstantObject("md.slr.rotation",SealevelriseRotationEnum));
110 parameters->AddObject(iomodel->CopyConstantObject("md.slr.tide_love_h",SealevelriseTidalLoveHEnum));
111 parameters->AddObject(iomodel->CopyConstantObject("md.slr.tide_love_k",SealevelriseTidalLoveKEnum));
112 parameters->AddObject(iomodel->CopyConstantObject("md.slr.fluid_love",SealevelriseFluidLoveEnum));
113 parameters->AddObject(iomodel->CopyConstantObject("md.slr.equatorial_moi",SealevelriseEquatorialMoiEnum));
114 parameters->AddObject(iomodel->CopyConstantObject("md.slr.polar_moi",SealevelrisePolarMoiEnum));
115 parameters->AddObject(iomodel->CopyConstantObject("md.slr.angular_velocity",SealevelriseAngularVelocityEnum));
116 parameters->AddObject(iomodel->CopyConstantObject("md.slr.ocean_area_scaling",SealevelriseOceanAreaScalingEnum));
117 parameters->AddObject(iomodel->CopyConstantObject("md.slr.geodetic",SealevelriseGeodeticEnum));
118
119 iomodel->FetchData(&elastic,"md.slr.elastic");
120 if(elastic){
121
122 /*love numbers: */
123 iomodel->FetchData(&love_h,&nl,NULL,"md.slr.love_h");
124 iomodel->FetchData(&love_k,&nl,NULL,"md.slr.love_k");
125 iomodel->FetchData(&love_l,&nl,NULL,"md.slr.love_l");
126
127 /*compute elastic green function for a range of angles*/
128 iomodel->FetchData(&degacc,"md.slr.degacc");
129 M=reCast<int,IssmDouble>(180./degacc+1.);
130
131 // AD performance is sensitive to calls to ensurecontiguous.
132 // // Providing "t" will cause ensurecontiguous to be called.
133 #ifdef _HAVE_AD_
134 G_elastic=xNew<IssmDouble>(M,"t");
135 U_elastic=xNew<IssmDouble>(M,"t");
136 H_elastic=xNew<IssmDouble>(M,"t");
137 #else
138 G_elastic=xNew<IssmDouble>(M);
139 U_elastic=xNew<IssmDouble>(M);
140 H_elastic=xNew<IssmDouble>(M);
141 #endif
142
143
144 /*compute combined legendre + love number (elastic green function:*/
145 m=DetermineLocalSize(M,IssmComm::GetComm());
146 GetOwnershipBoundariesFromRange(&lower_row,&upper_row,m,IssmComm::GetComm());
147 // AD performance is sensitive to calls to ensurecontiguous.
148 // // Providing "t" will cause ensurecontiguous to be called.
149 #ifdef _HAVE_AD_
150 G_elastic_local=xNew<IssmDouble>(m,"t");
151 U_elastic_local=xNew<IssmDouble>(m,"t");
152 H_elastic_local=xNew<IssmDouble>(m,"t");
153 #else
154 G_elastic_local=xNew<IssmDouble>(m);
155 U_elastic_local=xNew<IssmDouble>(m);
156 H_elastic_local=xNew<IssmDouble>(m);
157 #endif
158
159 for(int i=lower_row;i<upper_row;i++){
160 IssmDouble alpha,x;
161 alpha= reCast<IssmDouble>(i)*degacc * PI / 180.0;
162
163 G_elastic_local[i-lower_row]= (love_k[nl-1]-love_h[nl-1])/2.0/sin(alpha/2.0);
164 U_elastic_local[i-lower_row]= (love_h[nl-1])/2.0/sin(alpha/2.0);
165 H_elastic_local[i-lower_row]= 0;
166 IssmDouble Pn = 0.;
167 IssmDouble Pn1 = 0.;
168 IssmDouble Pn2 = 0.;
169 IssmDouble Pn_p = 0.;
170 IssmDouble Pn_p1 = 0.;
171 IssmDouble Pn_p2 = 0.;
172
173 for (int n=0;n<nl;n++) {
174 IssmDouble deltalove_G;
175 IssmDouble deltalove_U;
176
177 deltalove_G = (love_k[n]-love_k[nl-1]-love_h[n]+love_h[nl-1]);
178 deltalove_U = (love_h[n]-love_h[nl-1]);
179
180 /*compute legendre polynomials: P_n(cos\theta) & d P_n(cos\theta)/ d\theta: */
181 if(n==0){
182 Pn=1;
183 Pn_p=0;
184 }
185 else if(n==1){
186 Pn = cos(alpha);
187 Pn_p = 1;
188 }
189 else{
190 Pn = ( (2*n-1)*cos(alpha)*Pn1 - (n-1)*Pn2 ) /n;
191 Pn_p = ( (2*n-1)*(Pn1+cos(alpha)*Pn_p1) - (n-1)*Pn_p2 ) /n;
192 }
193 Pn2=Pn1; Pn1=Pn;
194 Pn_p2=Pn_p1; Pn_p1=Pn_p;
195
196 G_elastic_local[i-lower_row] += deltalove_G*Pn; // gravitational potential
197 U_elastic_local[i-lower_row] += deltalove_U*Pn; // vertical (up) displacement
198 H_elastic_local[i-lower_row] += sin(alpha)*love_l[n]*Pn_p; // horizontal displacements
199 }
200 }
201
202 /*merge G_elastic_local into G_elastic; U_elastic_local into U_elastic; H_elastic_local to H_elastic:{{{*/
203 int* recvcounts=xNew<int>(IssmComm::GetSize());
204 int* displs=xNew<int>(IssmComm::GetSize());
205
206 //recvcounts:
207 ISSM_MPI_Allgather(&m,1,ISSM_MPI_INT,recvcounts,1,ISSM_MPI_INT,IssmComm::GetComm());
208
209 /*displs: */
210 ISSM_MPI_Allgather(&lower_row,1,ISSM_MPI_INT,displs,1,ISSM_MPI_INT,IssmComm::GetComm());
211
212 /*All gather:*/
213 ISSM_MPI_Allgatherv(G_elastic_local, m, ISSM_MPI_DOUBLE, G_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
214 ISSM_MPI_Allgatherv(U_elastic_local, m, ISSM_MPI_DOUBLE, U_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
215 ISSM_MPI_Allgatherv(H_elastic_local, m, ISSM_MPI_DOUBLE, H_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
216 /*free ressources: */
217 xDelete<int>(recvcounts);
218 xDelete<int>(displs);
219
220 /*}}}*/
221
222 /*Avoid singularity at 0: */
223 G_elastic[0]=G_elastic[1];
224 parameters->AddObject(new DoubleVecParam(SealevelriseGElasticEnum,G_elastic,M));
225 U_elastic[0]=U_elastic[1];
226 parameters->AddObject(new DoubleVecParam(SealevelriseUElasticEnum,U_elastic,M));
227 H_elastic[0]=H_elastic[1];
228 parameters->AddObject(new DoubleVecParam(SealevelriseHElasticEnum,H_elastic,M));
229
230 /*free ressources: */
231 xDelete<IssmDouble>(love_h);
232 xDelete<IssmDouble>(love_k);
233 xDelete<IssmDouble>(love_l);
234 xDelete<IssmDouble>(G_elastic);
235 xDelete<IssmDouble>(G_elastic_local);
236 xDelete<IssmDouble>(U_elastic);
237 xDelete<IssmDouble>(U_elastic_local);
238 xDelete<IssmDouble>(H_elastic);
239 xDelete<IssmDouble>(H_elastic_local);
240 }
241
242 /*Transitions: */
243 iomodel->FetchData(&transitions,&transitions_M,&transitions_N,&ntransitions,"md.slr.transitions");
244 if(transitions){
245 parameters->AddObject(new DoubleMatArrayParam(SealevelriseTransitionsEnum,transitions,ntransitions,transitions_M,transitions_N));
246
247 for(int i=0;i<ntransitions;i++){
248 IssmDouble* transition=transitions[i];
249 xDelete<IssmDouble>(transition);
250 }
251 xDelete<IssmDouble*>(transitions);
252 xDelete<int>(transitions_M);
253 xDelete<int>(transitions_N);
254 }
255
256 /*Requested outputs*/
257 iomodel->FindConstant(&requestedoutputs,&numoutputs,"md.slr.requested_outputs");
258 if(numoutputs)parameters->AddObject(new StringArrayParam(SealevelriseRequestedOutputsEnum,requestedoutputs,numoutputs));
259 iomodel->DeleteData(&requestedoutputs,numoutputs,"md.slr.requested_outputs");
260
261
262}/*}}}*/
263
264/*Finite Element Analysis*/
265void SealevelriseAnalysis::Core(FemModel* femmodel){/*{{{*/
266 _error_("not implemented");
267}/*}}}*/
268ElementVector* SealevelriseAnalysis::CreateDVector(Element* element){/*{{{*/
269 /*Default, return NULL*/
270 return NULL;
271}/*}}}*/
272ElementMatrix* SealevelriseAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/
273_error_("Not implemented");
274}/*}}}*/
275ElementMatrix* SealevelriseAnalysis::CreateKMatrix(Element* element){/*{{{*/
276 _error_("not implemented yet");
277}/*}}}*/
278ElementVector* SealevelriseAnalysis::CreatePVector(Element* element){/*{{{*/
279_error_("not implemented yet");
280}/*}}}*/
281void SealevelriseAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
282 _error_("not implemented yet");
283}/*}}}*/
284void SealevelriseAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
285 _error_("Not implemented yet");
286}/*}}}*/
287void SealevelriseAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
288 _error_("not implemeneted yet!");
289
290}/*}}}*/
291void SealevelriseAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
292 /*Default, do nothing*/
293 return;
294}/*}}}*/
Note: See TracBrowser for help on using the repository browser.