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

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

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