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

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

CHG: parallelized computation of G_elastic in the Sealevelrise UpdateParameters analysis phase.

File size: 5.9 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: */
33 iomodel->FetchDataToInput(elements,MaskIceLevelsetEnum);
34 iomodel->FetchDataToInput(elements,MaskOceanLevelsetEnum);
[19987]35 iomodel->FetchDataToInput(elements,MaskLandLevelsetEnum);
[19984]36 iomodel->FetchDataToInput(elements,SealevelriseDeltathicknessEnum);
37
38}/*}}}*/
39void SealevelriseAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
40
[19986]41 int nl;
[19984]42 IssmDouble* love_h=NULL;
43 IssmDouble* love_k=NULL;
[19986]44
[20029]45 bool elastic=false;
[20028]46 IssmDouble* G_elastic = NULL;
[20033]47 IssmDouble* G_elastic_local = NULL;
48 int M,m,lower_row,upper_row;
49 IssmDouble degacc=.01;
[19984]50
[19986]51 /*some constant parameters: */
[19984]52 parameters->AddObject(iomodel->CopyConstantObject(SealevelriseReltolEnum));
53 parameters->AddObject(iomodel->CopyConstantObject(SealevelriseAbstolEnum));
54 parameters->AddObject(iomodel->CopyConstantObject(SealevelriseMaxiterEnum));
55 parameters->AddObject(iomodel->CopyConstantObject(SealevelriseRigidEnum));
56 parameters->AddObject(iomodel->CopyConstantObject(SealevelriseElasticEnum));
57 parameters->AddObject(iomodel->CopyConstantObject(SealevelriseEustaticEnum));
[19986]58
[20029]59 iomodel->FetchData(&elastic,SealevelriseElasticEnum);
60 if(elastic){
[19984]61
[20029]62 /*love numbers: */
63 iomodel->FetchData(&love_h,&nl,NULL,SealevelriseLoveHEnum);
64 iomodel->FetchData(&love_k,&nl,NULL,SealevelriseLoveKEnum);
[20022]65
[20029]66 /*compute elastic green function for a range of angles*/
[20033]67 iomodel->FetchData(&degacc,SealevelriseDegaccEnum);
68 M=reCast<int>(180/degacc+1);
[20029]69 G_elastic=xNew<IssmDouble>(M);
[20028]70
[20029]71 /*compute combined legendre + love number (elastic green function:*/
[20033]72 m=DetermineLocalSize(M,IssmComm::GetComm());
73 GetOwnershipBoundariesFromRange(&lower_row,&upper_row,m,IssmComm::GetComm());
74 G_elastic_local=xNew<IssmDouble>(m);
75
76 for(int i=lower_row;i<upper_row;i++){
[20029]77 IssmDouble alpha,x;
78 alpha= reCast<IssmDouble>(i)*degacc * PI / 180.0;
79
[20033]80 G_elastic_local[i-lower_row]= (love_k[nl-1]-love_h[nl-1])/2.0/sin(alpha/2.0);
81 for (int n=0;n<nl;n++) {
82 IssmDouble Pn,Pn1,Pn2;
83 IssmDouble deltalove;
84
85 deltalove = (love_k[n]-love_k[nl-1]-love_h[n]+love_h[nl-1]);
86
87 if(n==0)Pn=1;
88 else if(n==1)Pn=cos(alpha);
89 else Pn= ( (2*n-1)*cos(alpha)*Pn1 - (n-1)*Pn2 ) /n;
90 Pn2=Pn1; Pn1=Pn;
91
92 G_elastic_local[i-lower_row] += deltalove*Pn;
93 }
94 }
95
96 /*merge G_elastic_local into G_elastic:{{{*/
97 int* recvcounts=xNew<int>(IssmComm::GetSize());
98 int* displs=xNew<int>(IssmComm::GetSize());
99
100 //recvcounts:
101 ISSM_MPI_Allgather(&m,1,ISSM_MPI_INT,recvcounts,1,ISSM_MPI_INT,IssmComm::GetComm());
102
103 /*displs: */
104 ISSM_MPI_Allgather(&lower_row,1,ISSM_MPI_INT,displs,1,ISSM_MPI_INT,IssmComm::GetComm());
105
106 /*All gather:*/
107 ISSM_MPI_Allgatherv(G_elastic_local, m, ISSM_MPI_DOUBLE, G_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
108 /*free ressources: */
109 xDelete<int>(recvcounts);
110 xDelete<int>(displs);
111
112 /*}}}*/
113
114 /*Old code: {{{*/
115 /*for(int i=0;i<M;i++){
116 IssmDouble alpha,x;
117 alpha= reCast<IssmDouble>(i)*degacc * PI / 180.0;
118
[20029]119 G_elastic[i]= (love_k[nl-1]-love_h[nl-1])/2.0/sin(alpha/2.0);
120 for (int n=0;n<nl;n++) {
121 IssmDouble Pn,Pn1,Pn2;
122 IssmDouble deltalove;
123
124 deltalove = (love_k[n]-love_k[nl-1]-love_h[n]+love_h[nl-1]);
125
126 if(n==0)Pn=1;
127 else if(n==1)Pn=cos(alpha);
128 else Pn= ( (2*n-1)*cos(alpha)*Pn1 - (n-1)*Pn2 ) /n;
129 Pn2=Pn1; Pn1=Pn;
130
131 G_elastic[i] += deltalove*Pn;
132 }
[20033]133 }*/
134 /*}}}*/
[19986]135
[20029]136 /*Avoid singularity at 0: */
137 G_elastic[0]=G_elastic[1];
138 parameters->AddObject(new DoubleVecParam(SealevelriseGElasticEnum,G_elastic,M));
[20028]139
[20029]140 /*free ressources: */
141 xDelete<IssmDouble>(love_h);
142 xDelete<IssmDouble>(love_k);
143 xDelete<IssmDouble>(G_elastic);
[20033]144 xDelete<IssmDouble>(G_elastic_local);
[20029]145 }
146
[19984]147}/*}}}*/
148
149/*Finite Element Analysis*/
150void SealevelriseAnalysis::Core(FemModel* femmodel){/*{{{*/
151 _error_("not implemented");
152}/*}}}*/
153ElementVector* SealevelriseAnalysis::CreateDVector(Element* element){/*{{{*/
154 /*Default, return NULL*/
155 return NULL;
156}/*}}}*/
157ElementMatrix* SealevelriseAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/
158_error_("Not implemented");
159}/*}}}*/
160ElementMatrix* SealevelriseAnalysis::CreateKMatrix(Element* element){/*{{{*/
161 _error_("not implemented yet");
162}/*}}}*/
163ElementVector* SealevelriseAnalysis::CreatePVector(Element* element){/*{{{*/
164_error_("not implemented yet");
165}/*}}}*/
166void SealevelriseAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
167 _error_("not implemented yet");
168}/*}}}*/
169void SealevelriseAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
170 _error_("Not implemented yet");
171}/*}}}*/
172void SealevelriseAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
173 _error_("not implemented yet");
174}/*}}}*/
175void SealevelriseAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
176 /*Default, do nothing*/
177 return;
178}/*}}}*/
Note: See TracBrowser for help on using the repository browser.