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

Last change on this file since 21295 was 21295, checked in by adhikari, 8 years ago

CHG: ocean_area_scaling option added to capture true eustatic rate
File size: 9.1 KB
RevLine 
[19984]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: */
[20690]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);
[20704]36 iomodel->FetchDataToInput(elements,"md.slr.deltathickness",SealevelriseDeltathicknessEnum);
[20690]37 iomodel->FetchDataToInput(elements,"md.slr.sealevel",SealevelEnum,0);
[19984]38
39}/*}}}*/
40void SealevelriseAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
41
[19986]42 int nl;
[19984]43 IssmDouble* love_h=NULL;
44 IssmDouble* love_k=NULL;
[20999]45 IssmDouble* love_l=NULL;
[19986]46
[20029]47 bool elastic=false;
[20028]48 IssmDouble* G_elastic = NULL;
[20033]49 IssmDouble* G_elastic_local = NULL;
[20999]50 IssmDouble* U_elastic = NULL;
51 IssmDouble* U_elastic_local = NULL;
52 IssmDouble* H_elastic = NULL;
53 IssmDouble* H_elastic_local = NULL;
[20033]54 int M,m,lower_row,upper_row;
55 IssmDouble degacc=.01;
[19984]56
[20036]57 int numoutputs;
58 char** requestedoutputs = NULL;
59
[20136]60 /*transition vectors: */
61 IssmDouble **transitions = NULL;
62 int *transitions_M = NULL;
63 int *transitions_N = NULL;
64 int ntransitions;
65
[19986]66 /*some constant parameters: */
[20690]67 parameters->AddObject(iomodel->CopyConstantObject("md.slr.reltol",SealevelriseReltolEnum));
68 parameters->AddObject(iomodel->CopyConstantObject("md.slr.abstol",SealevelriseAbstolEnum));
69 parameters->AddObject(iomodel->CopyConstantObject("md.slr.maxiter",SealevelriseMaxiterEnum));
70 parameters->AddObject(iomodel->CopyConstantObject("md.slr.rigid",SealevelriseRigidEnum));
71 parameters->AddObject(iomodel->CopyConstantObject("md.slr.elastic",SealevelriseElasticEnum));
[21295]72 parameters->AddObject(iomodel->CopyConstantObject("md.slr.ocean_area_scaling",SealevelriseOceanAreaScalingEnum));
[19986]73
[20690]74 iomodel->FetchData(&elastic,"md.slr.elastic");
[20029]75 if(elastic){
[19984]76
[20029]77 /*love numbers: */
[20690]78 iomodel->FetchData(&love_h,&nl,NULL,"md.slr.love_h");
79 iomodel->FetchData(&love_k,&nl,NULL,"md.slr.love_k");
[20999]80 iomodel->FetchData(&love_l,&nl,NULL,"md.slr.love_l");
[20022]81
[20029]82 /*compute elastic green function for a range of angles*/
[20690]83 iomodel->FetchData(&degacc,"md.slr.degacc");
[20034]84 M=reCast<int,IssmDouble>(180./degacc+1.);
[20029]85 G_elastic=xNew<IssmDouble>(M);
[20999]86 U_elastic=xNew<IssmDouble>(M);
87 H_elastic=xNew<IssmDouble>(M);
[20044]88
89 /*compute combined legendre + love number (elastic green function:*/
[20033]90 m=DetermineLocalSize(M,IssmComm::GetComm());
91 GetOwnershipBoundariesFromRange(&lower_row,&upper_row,m,IssmComm::GetComm());
92 G_elastic_local=xNew<IssmDouble>(m);
[20999]93 U_elastic_local=xNew<IssmDouble>(m);
94 H_elastic_local=xNew<IssmDouble>(m);
[20033]95
96 for(int i=lower_row;i<upper_row;i++){
[20029]97 IssmDouble alpha,x;
98 alpha= reCast<IssmDouble>(i)*degacc * PI / 180.0;
99
[20033]100 G_elastic_local[i-lower_row]= (love_k[nl-1]-love_h[nl-1])/2.0/sin(alpha/2.0);
[20999]101 U_elastic_local[i-lower_row]= (love_h[nl-1])/2.0/sin(alpha/2.0);
102 H_elastic_local[i-lower_row]= 0;
[20054]103 IssmDouble Pn,Pn1,Pn2;
[20999]104 IssmDouble Pn_p,Pn_p1,Pn_p2;
[20033]105 for (int n=0;n<nl;n++) {
[20999]106 IssmDouble deltalove_G;
107 IssmDouble deltalove_U;
[20033]108
[20999]109 deltalove_G = (love_k[n]-love_k[nl-1]-love_h[n]+love_h[nl-1]);
110 deltalove_U = (love_h[n]-love_h[nl-1]);
111
112 /*compute legendre polynomials: P_n(cos\theta) & d P_n(cos\theta)/ d\theta: */
113 if(n==0){
114 Pn=1;
115 Pn_p=0;
116 }
117 else if(n==1){
118 Pn = cos(alpha);
119 Pn_p = 1;
120 }
121 else{
122 Pn = ( (2*n-1)*cos(alpha)*Pn1 - (n-1)*Pn2 ) /n;
123 Pn_p = ( (2*n-1)*(Pn1+cos(alpha)*Pn_p1) - (n-1)*Pn_p2 ) /n;
124 }
[20033]125 Pn2=Pn1; Pn1=Pn;
[20999]126 Pn_p2=Pn_p1; Pn_p1=Pn_p;
[20033]127
[20999]128 G_elastic_local[i-lower_row] += deltalove_G*Pn; // gravitational potential
129 U_elastic_local[i-lower_row] += deltalove_U*Pn; // vertical (up) displacement
130 H_elastic_local[i-lower_row] += sin(alpha)*love_l[n]*Pn_p; // horizontal displacements
[20033]131 }
132 }
133
[20999]134 /*merge G_elastic_local into G_elastic; U_elastic_local into U_elastic; H_elastic_local to H_elastic:{{{*/
[20033]135 int* recvcounts=xNew<int>(IssmComm::GetSize());
136 int* displs=xNew<int>(IssmComm::GetSize());
137
138 //recvcounts:
139 ISSM_MPI_Allgather(&m,1,ISSM_MPI_INT,recvcounts,1,ISSM_MPI_INT,IssmComm::GetComm());
140
141 /*displs: */
142 ISSM_MPI_Allgather(&lower_row,1,ISSM_MPI_INT,displs,1,ISSM_MPI_INT,IssmComm::GetComm());
143
144 /*All gather:*/
145 ISSM_MPI_Allgatherv(G_elastic_local, m, ISSM_MPI_DOUBLE, G_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
[20999]146 ISSM_MPI_Allgatherv(U_elastic_local, m, ISSM_MPI_DOUBLE, U_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
147 ISSM_MPI_Allgatherv(H_elastic_local, m, ISSM_MPI_DOUBLE, H_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
[20033]148 /*free ressources: */
149 xDelete<int>(recvcounts);
150 xDelete<int>(displs);
151
152 /*}}}*/
153
[20029]154 /*Avoid singularity at 0: */
155 G_elastic[0]=G_elastic[1];
156 parameters->AddObject(new DoubleVecParam(SealevelriseGElasticEnum,G_elastic,M));
[20999]157 U_elastic[0]=U_elastic[1];
158 parameters->AddObject(new DoubleVecParam(SealevelriseUElasticEnum,U_elastic,M));
159 H_elastic[0]=H_elastic[1];
160 parameters->AddObject(new DoubleVecParam(SealevelriseHElasticEnum,H_elastic,M));
[20028]161
[20029]162 /*free ressources: */
163 xDelete<IssmDouble>(love_h);
164 xDelete<IssmDouble>(love_k);
[20999]165 xDelete<IssmDouble>(love_l);
[20029]166 xDelete<IssmDouble>(G_elastic);
[20033]167 xDelete<IssmDouble>(G_elastic_local);
[20999]168 xDelete<IssmDouble>(U_elastic);
169 xDelete<IssmDouble>(U_elastic_local);
170 xDelete<IssmDouble>(H_elastic);
171 xDelete<IssmDouble>(H_elastic_local);
[20029]172 }
[20036]173
[20136]174 /*Transitions: */
[20690]175 iomodel->FetchData(&transitions,&transitions_M,&transitions_N,&ntransitions,"md.slr.transitions");
[20136]176 if(transitions){
177 parameters->AddObject(new DoubleMatArrayParam(SealevelriseTransitionsEnum,transitions,ntransitions,transitions_M,transitions_N));
178
179 for(int i=0;i<ntransitions;i++){
180 IssmDouble* transition=transitions[i];
181 xDelete<IssmDouble>(transition);
182 }
183 xDelete<IssmDouble*>(transitions);
[20138]184 xDelete<int>(transitions_M);
185 xDelete<int>(transitions_N);
[20136]186 }
187
[20036]188 /*Requested outputs*/
[20690]189 iomodel->FindConstant(&requestedoutputs,&numoutputs,"md.slr.requested_outputs");
[20036]190 parameters->AddObject(new IntParam(SealevelriseNumRequestedOutputsEnum,numoutputs));
191 if(numoutputs)parameters->AddObject(new StringArrayParam(SealevelriseRequestedOutputsEnum,requestedoutputs,numoutputs));
[20690]192 iomodel->DeleteData(&requestedoutputs,numoutputs,"md.slr.requested_outputs");
[20029]193
[20036]194
[19984]195}/*}}}*/
196
197/*Finite Element Analysis*/
198void SealevelriseAnalysis::Core(FemModel* femmodel){/*{{{*/
199 _error_("not implemented");
200}/*}}}*/
201ElementVector* SealevelriseAnalysis::CreateDVector(Element* element){/*{{{*/
202 /*Default, return NULL*/
203 return NULL;
204}/*}}}*/
205ElementMatrix* SealevelriseAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/
206_error_("Not implemented");
207}/*}}}*/
208ElementMatrix* SealevelriseAnalysis::CreateKMatrix(Element* element){/*{{{*/
209 _error_("not implemented yet");
210}/*}}}*/
211ElementVector* SealevelriseAnalysis::CreatePVector(Element* element){/*{{{*/
212_error_("not implemented yet");
213}/*}}}*/
214void SealevelriseAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
215 _error_("not implemented yet");
216}/*}}}*/
217void SealevelriseAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
218 _error_("Not implemented yet");
219}/*}}}*/
220void SealevelriseAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
[20153]221
222 IssmDouble *deltaS = NULL;
223 IssmDouble *S = NULL;
224 int* sidlist = NULL;
225 int numvertices;
226
227 numvertices= element->GetNumberOfVertices();
228 sidlist=xNew<int>(numvertices);
229
230 element->GetVerticesSidList(sidlist);
231
232 deltaS = xNew<IssmDouble>(numvertices);
233 for(int i=0;i<numvertices;i++){
234 deltaS[i]=solution[sidlist[i]];
235 }
236
237 S = xNew<IssmDouble>(numvertices);
[20154]238 element->GetInputListOnVertices(S,SealevelEnum,0);
[20153]239
240 /*Add deltaS to S:*/
241 for (int i=0;i<numvertices;i++)S[i]+=deltaS[i];
242
243 /*Add S back into inputs: */
[20154]244 element->AddInput(SealevelEnum,S,P1Enum);
[20153]245
246 /*Free ressources:*/
247 xDelete<int>(sidlist);
248 xDelete<IssmDouble>(deltaS);
249 xDelete<IssmDouble>(S);
250
[19984]251}/*}}}*/
252void SealevelriseAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
253 /*Default, do nothing*/
254 return;
255}/*}}}*/
Note: See TracBrowser for help on using the repository browser.