source: issm/oecreview/Archive/19101-20495/ISSM-20135-20136.diff@ 20498

Last change on this file since 20498 was 20498, checked in by Mathieu Morlighem, 9 years ago

CHG: done with Archive/19101-20495
File size: 15.4 KB
RevLine 
[20498]1Index: ../trunk-jpl/src/c/analyses/SealevelriseAnalysis.cpp
2===================================================================
3--- ../trunk-jpl/src/c/analyses/SealevelriseAnalysis.cpp (revision 20135)
4+++ ../trunk-jpl/src/c/analyses/SealevelriseAnalysis.cpp (revision 20136)
5@@ -51,6 +51,12 @@
6 int numoutputs;
7 char** requestedoutputs = NULL;
8
9+ /*transition vectors: */
10+ IssmDouble **transitions = NULL;
11+ int *transitions_M = NULL;
12+ int *transitions_N = NULL;
13+ int ntransitions;
14+
15 /*some constant parameters: */
16 parameters->AddObject(iomodel->CopyConstantObject(SealevelriseReltolEnum));
17 parameters->AddObject(iomodel->CopyConstantObject(SealevelriseAbstolEnum));
18@@ -125,6 +131,20 @@
19 xDelete<IssmDouble>(G_elastic_local);
20 }
21
22+ /*Transitions: */
23+ iomodel->FetchData(&transitions,&transitions_M,&transitions_N,&ntransitions,SealevelriseTransitionsEnum);
24+ if(transitions){
25+ parameters->AddObject(new DoubleMatArrayParam(SealevelriseTransitionsEnum,transitions,ntransitions,transitions_M,transitions_N));
26+
27+ for(int i=0;i<ntransitions;i++){
28+ IssmDouble* transition=transitions[i];
29+ xDelete<IssmDouble>(transition);
30+ }
31+ xDelete<IssmDouble*>(transitions);
32+ xDelete<IssmDouble*>(transitions_M);
33+ xDelete<IssmDouble*>(transitions_N);
34+ }
35+
36 /*Requested outputs*/
37 iomodel->FetchData(&requestedoutputs,&numoutputs,SealevelriseRequestedOutputsEnum);
38 parameters->AddObject(new IntParam(SealevelriseNumRequestedOutputsEnum,numoutputs));
39Index: ../trunk-jpl/src/c/cores/transient_core.cpp
40===================================================================
41--- ../trunk-jpl/src/c/cores/transient_core.cpp (revision 20135)
42+++ ../trunk-jpl/src/c/cores/transient_core.cpp (revision 20136)
43@@ -20,7 +20,7 @@
44
45 /*parameters: */
46 IssmDouble starttime,finaltime,dt,yts;
47- bool isstressbalance,ismasstransport,issmb,isFS,isthermal,isgroundingline,isgia,islevelset,isdamageevolution,ishydrology,iscalving;
48+ bool isstressbalance,ismasstransport,issmb,isFS,isthermal,isgroundingline,isgia,isslr,iscoupler,islevelset,isdamageevolution,ishydrology,iscalving;
49 bool save_results,dakota_analysis;
50 bool time_adapt;
51 int output_frequency;
52@@ -50,6 +50,8 @@
53 femmodel->parameters->FindParam(&issmb,TransientIssmbEnum);
54 femmodel->parameters->FindParam(&isthermal,TransientIsthermalEnum);
55 femmodel->parameters->FindParam(&isgia,TransientIsgiaEnum);
56+ femmodel->parameters->FindParam(&isslr,TransientIsslrEnum);
57+ femmodel->parameters->FindParam(&iscoupler,TransientIscouplerEnum);
58 femmodel->parameters->FindParam(&isgroundingline,TransientIsgroundinglineEnum);
59 femmodel->parameters->FindParam(&islevelset,TransientIslevelsetEnum);
60 femmodel->parameters->FindParam(&isdamageevolution,TransientIsdamageevolutionEnum);
61@@ -140,6 +142,9 @@
62 #endif
63 }
64
65+ /*Sea level rise: */
66+ if(isslr | iscoupler) sealevelrise_core(femmodel);
67+
68 /*unload results*/
69 if(VerboseSolution()) _printf0_(" computing requested outputs\n");
70 femmodel->RequestedOutputsx(&femmodel->results,requested_outputs,numoutputs,save_results);
71Index: ../trunk-jpl/src/c/cores/cores.h
72===================================================================
73--- ../trunk-jpl/src/c/cores/cores.h (revision 20135)
74+++ ../trunk-jpl/src/c/cores/cores.h (revision 20136)
75@@ -49,8 +49,8 @@
76 void smb_core(FemModel* femmodel);
77 void damage_core(FemModel* femmodel);
78 void sealevelrise_core(FemModel* femmodel);
79-Vector<IssmDouble> * sealevelrise_core_eustatic(FemModel* femmodel);
80-Vector<IssmDouble> * sealevelrise_core_noneustatic(FemModel* femmodel,Vector<IssmDouble>* Sg_eustatic);
81+Vector<IssmDouble>* sealevelrise_core_eustatic(FemModel* femmodel);
82+Vector<IssmDouble>* sealevelrise_core_noneustatic(FemModel* femmodel,Vector<IssmDouble>* Sg_eustatic);
83 IssmDouble objectivefunction(IssmDouble search_scalar,FemModel* femmodel);
84
85 //optimization
86@@ -61,6 +61,8 @@
87 void WriteLockFile(char* filename);
88 void ResetBoundaryConditions(FemModel* femmodel, int analysis_type);
89 void PrintBanner(void);
90+void TransferForcing(FemModel* femmodel,int forcingenum);
91+void TransferSealevel(FemModel* femmodel,int forcingenum);
92
93 //solution configuration
94 void CorePointerFromSolutionEnum(void (**psolutioncore)(FemModel*),Parameters* parameters,int solutiontype);
95Index: ../trunk-jpl/src/c/cores/sealevelrise_core.cpp
96===================================================================
97--- ../trunk-jpl/src/c/cores/sealevelrise_core.cpp (revision 20135)
98+++ ../trunk-jpl/src/c/cores/sealevelrise_core.cpp (revision 20136)
99@@ -9,43 +9,281 @@
100 #include "../modules/modules.h"
101 #include "../solutionsequences/solutionsequences.h"
102
103-void sealevelrise_core(FemModel* femmodel){
104+void sealevelrise_core(FemModel* femmodel){ /*{{{*/
105
106 Vector<IssmDouble> *Sg = NULL;
107 Vector<IssmDouble> *Sg_eustatic = NULL;
108- bool save_results;
109+ bool save_results,isslr,iscoupler;
110 int configuration_type;
111+ int solution_type;
112 int numoutputs = 0;
113 char **requested_outputs = NULL;
114-
115- if(VerboseSolution()) _printf0_(" computing sea level rise\n");
116-
117+
118 /*Recover some parameters: */
119 femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
120+ femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
121 femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
122- femmodel->parameters->FindParam(&numoutputs,SealevelriseNumRequestedOutputsEnum);
123- if(numoutputs) femmodel->parameters->FindParam(&requested_outputs,&numoutputs,SealevelriseRequestedOutputsEnum);
124+ femmodel->parameters->FindParam(&isslr,TransientIsslrEnum);
125+ femmodel->parameters->FindParam(&iscoupler,TransientIscouplerEnum);
126
127+ /*several cases here, depending on value of iscoupler and isslr:
128+ ( !iscoupler & !isslr) we are not interested in being here :)
129+ ( !iscoupler & isslr) we are running in uncoupled mode
130+ ( iscoupler & isslr) we are running in coupled mode, and better be earth
131+ ( iscoupler & !isslr) we are running in coupled mode, and better be an ice cap
132+ */
133+
134+ /*early return: */
135+ if( !iscoupler & !isslr) return; //we are not interested in being here :)
136+
137+ /*In what follows we assume we are all running slr, either in coupled, or uncoupled mode:*/
138+ if(VerboseSolution()) _printf0_(" computing sea level rise\n");
139+
140 /*set configuration: */
141- femmodel->SetCurrentConfiguration(SealevelriseAnalysisEnum);
142+ if(isslr)femmodel->SetCurrentConfiguration(SealevelriseAnalysisEnum);
143
144- /*call two sub cores:*/
145- Sg_eustatic=sealevelrise_core_eustatic(femmodel); //generalized eustatic (Farrel and Clark, Eq 4, 1st, 3rd and 4rd terms on the RHS.
146+ /*transfer deltathickness forcing from ice caps to earth model: */
147+ if(iscoupler) TransferForcing(femmodel,SealevelriseDeltathicknessEnum);
148
149- Sg=sealevelrise_core_noneustatic(femmodel,Sg_eustatic); //sea-level rise dependent terms (2nd and 5th terms on the RHS)
150+ /*call sea-level rise sub cores:*/
151+ if(isslr){
152+ Sg_eustatic=sealevelrise_core_eustatic(femmodel); //generalized eustatic (Farrel and Clark, Eq 4, 1st, 3rd and 4rd terms on the RHS.
153
154- /*get results out:*/
155- InputUpdateFromVectorx(femmodel,Sg,SealevelriseSEnum,VertexSIdEnum);
156+ Sg=sealevelrise_core_noneustatic(femmodel,Sg_eustatic); //sea-level rise dependent terms (2nd and 5th terms on the RHS)
157
158- if(save_results){
159- if(VerboseSolution()) _printf0_(" saving results\n");
160- femmodel->RequestedOutputsx(&femmodel->results,requested_outputs,numoutputs);
161+ /*get results into elements:*/
162+ InputUpdateFromVectorx(femmodel,Sg,SealevelriseSEnum,VertexSIdEnum);
163+
164+ if(save_results){
165+ if(VerboseSolution()) _printf0_(" saving results\n");
166+ femmodel->parameters->FindParam(&requested_outputs,&numoutputs,SealevelriseRequestedOutputsEnum);
167+ femmodel->RequestedOutputsx(&femmodel->results,requested_outputs,numoutputs);
168+ }
169+
170+ if(solution_type==SealevelriseSolutionEnum)femmodel->RequestedDependentsx();
171+
172+ /*Free ressources:*/
173+ delete Sg;
174+ delete Sg_eustatic;
175+ if(numoutputs){for(int i=0;i<numoutputs;i++){xDelete<char>(requested_outputs[i]);} xDelete<char*>(requested_outputs);}
176 }
177+
178+ /*transfer sea-level back to ice caps: */
179+ if(iscoupler)TransferSealevel(femmodel,SealevelriseSEnum);
180+}
181+/*}}}*/
182+void TransferForcing(FemModel* femmodel,int forcingenum){ /*{{{*/
183+
184+ /*forcing being transferred from models to earth: */
185+ IssmDouble** forcings=NULL;
186+ IssmDouble* forcing=NULL;
187+ Vector<IssmDouble>* forcingglobal=NULL;
188+ int* nvs=NULL;
189
190- /*Free ressources:*/
191- if(numoutputs){for(int i=0;i<numoutputs;i++){xDelete<char>(requested_outputs[i]);} xDelete<char*>(requested_outputs);}
192- delete Sg;
193- delete Sg_eustatic;
194+ /*transition vectors:*/
195+ IssmDouble** transitions=NULL;
196+ int ntransitions;
197+ int* transitions_m=NULL;
198+ int* transitions_n=NULL;
199+ int nv;
200+
201+ /*communicators:*/
202+ ISSM_MPI_Comm tocomm;
203+ ISSM_MPI_Comm* fromcomms=NULL;
204+ ISSM_MPI_Status status;
205+ int my_rank;
206+ int modelid,earthid;
207+ int nummodels;
208
209+ /*Recover some parameters: */
210+ femmodel->parameters->FindParam(&modelid,ModelIdEnum);
211+ femmodel->parameters->FindParam(&earthid,EarthIdEnum);
212+ femmodel->parameters->FindParam(&nummodels,NumModelsEnum);
213+ my_rank=IssmComm::GetRank();
214+
215+ /*retrieve the inter communicators that will be used to send data from each ice cap to the earth: */
216+ if(modelid==earthid)femmodel->parameters->FindParam(&fromcomms,&nv,IcecapToEarthCommEnum);
217+ else femmodel->parameters->FindParam(&tocomm, IcecapToEarthCommEnum);
218
219-}
220+ /*For each icecap, retrieve the forcing vector that will be sent to the earth model: */
221+ if(modelid!=earthid){
222+ nv=femmodel->vertices->NumberOfVertices();
223+ GetVectorFromInputsx(&forcing,femmodel,forcingenum,VertexSIdEnum);
224+ }
225+
226+ /*Send the forcing to the earth model:{{{*/
227+ if(my_rank==0){
228+ if(modelid==earthid){
229+ forcings=xNew<IssmDouble*>(nummodels-1);
230+ nvs=xNew<int>(nummodels-1);
231+ for(int i=0;i<earthid;i++){
232+ MPI_Recv(nvs+i, 1, ISSM_MPI_INT, 0,i, fromcomms[i], &status);
233+ forcings[i]=xNew<IssmDouble>(nvs[i]);
234+ MPI_Recv(forcings[i], nvs[i], ISSM_MPI_DOUBLE, 0,i, fromcomms[i], &status);
235+ }
236+
237+ }
238+ else{
239+ ISSM_MPI_Send(&nv, 1, ISSM_MPI_INT, 0, modelid, tocomm);
240+ ISSM_MPI_Send(forcing, nv, ISSM_MPI_DOUBLE, 0, modelid, tocomm);
241+ }
242+ }
243+ /*}}}*/
244+
245+ /*On the earth model, consolidate all the forcings into one, and update the elements dataset accordingly: {{{*/
246+ if(modelid==earthid){
247+
248+ /*Out of all the delta thicknesses, build one delta thickness vector made of all the ice cap contributions.
249+ *First, build the global delta thickness vector in the earth model: */
250+ nv=femmodel->vertices->NumberOfVertices();
251+ forcingglobal= new Vector<IssmDouble>(nv);
252+
253+ /*Retrieve transition vectors, used to plug from each ice cap into the global forcing:*/
254+ femmodel->parameters->FindParam(&transitions,&ntransitions,&transitions_m,&transitions_n,SealevelriseTransitionsEnum);
255+
256+ if(ntransitions!=earthid)_error_("TransferForcing error message: number of transition vectors is not equal to the number of icecaps!");
257+
258+ /*Go through all the delta thicknesses coming from each ice cap: */
259+ if(my_rank==0){
260+ for(int i=0;i<earthid;i++){
261+
262+ IssmDouble* forcingfromcap= forcings[i]; //careful, this only exists on rank 0 of the earth model!
263+ IssmDouble* transition=transitions[i];
264+ int M=transitions_m[i];
265+
266+ /*build index to plug values: */
267+ int* index=xNew<int>(M); for(int i=0;i<M;i++)index[i]=reCast<int>(transition[i])-1; //matlab indexing!
268+
269+
270+ /*We are going to plug this vector into the earth model, at the right vertices corresponding to this particular
271+ * ice cap: */
272+ forcingglobal->SetValues(M,index,forcingfromcap,INS_VAL);
273+ xDelete<int>(index);
274+ }
275+ }
276+
277+ /*Assemble vector:*/
278+ forcingglobal->Assemble();
279+
280+ /*Plug into elements:*/
281+ InputUpdateFromVectorx(femmodel,forcingglobal,forcingenum,VertexSIdEnum);
282+ }
283+ /*}}}*/
284+
285+ /*Free ressources:{{{*/
286+ if(forcings){
287+ for(int i=0;i<nummodels-1;i++){
288+ IssmDouble* temp=forcings[i]; xDelete<IssmDouble>(temp);
289+ }
290+ xDelete<IssmDouble*>(forcings);
291+ }
292+ if(forcing)xDelete<IssmDouble>(forcing);
293+ if(forcingglobal)delete forcingglobal;
294+ if(transitions){
295+ for(int i=0;i<earthid;i++){
296+ IssmDouble* temp=transitions[i];
297+ xDelete<IssmDouble>(temp);
298+ }
299+ xDelete<IssmDouble*>(transitions);
300+ xDelete<int>(transitions_m);
301+ xDelete<int>(transitions_n);
302+ }
303+ if(nvs)xDelete<int>(nvs);
304+ /*}}}*/
305+
306+} /*}}}*/
307+void TransferSealevel(FemModel* femmodel,int forcingenum){ /*{{{*/
308+
309+ /*forcing being transferred from earth to ice caps: */
310+ IssmDouble* forcing=NULL;
311+ IssmDouble* forcingglobal=NULL;
312+
313+ /*transition vectors:*/
314+ IssmDouble** transitions=NULL;
315+ int ntransitions;
316+ int* transitions_m=NULL;
317+ int* transitions_n=NULL;
318+ int nv;
319+
320+ /*communicators:*/
321+ ISSM_MPI_Comm fromcomm;
322+ ISSM_MPI_Comm* tocomms=NULL;
323+ ISSM_MPI_Status status;
324+ int my_rank;
325+ int modelid,earthid;
326+ int nummodels;
327+ int numcoms;
328+
329+ /*Recover some parameters: */
330+ femmodel->parameters->FindParam(&modelid,ModelIdEnum);
331+ femmodel->parameters->FindParam(&earthid,EarthIdEnum);
332+ femmodel->parameters->FindParam(&nummodels,NumModelsEnum);
333+ my_rank=IssmComm::GetRank();
334+
335+ /*retrieve the inter communicators that will be used to send data from earth to ice caps:*/
336+ if(modelid==earthid)femmodel->parameters->FindParam(&tocomms,&numcoms,IcecapToEarthCommEnum);
337+ else femmodel->parameters->FindParam(&fromcomm, IcecapToEarthCommEnum);
338+
339+
340+ /*Retrieve sea-level on earth model: */
341+ if(modelid==earthid){
342+ nv=femmodel->vertices->NumberOfVertices();
343+ GetVectorFromInputsx(&forcingglobal,femmodel,forcingenum,VertexSIdEnum);
344+ }
345+
346+ /*Send the forcing to the ice caps:{{{*/
347+ if(my_rank==0){
348+
349+ if(modelid==earthid){
350+
351+ /*Retrieve transition vectors, used to figure out global forcing contribution to each ice cap's own elements: */
352+ femmodel->parameters->FindParam(&transitions,&ntransitions,&transitions_m,&transitions_n,SealevelriseTransitionsEnum);
353+
354+ if(ntransitions!=earthid)_error_("TransferSeaLevel error message: number of transition vectors is not equal to the number of icecaps!");
355+
356+ for(int i=0;i<earthid;i++){
357+ nv=transitions_m[i];
358+ forcing=xNew<IssmDouble>(nv);
359+ IssmDouble* transition=transitions[i];
360+ for(int j=0;j<nv;j++){
361+ forcing[j]=forcingglobal[reCast<int>(transition[j])-1];
362+ }
363+ ISSM_MPI_Send(&nv, 1, ISSM_MPI_INT, 0, i, tocomms[i]);
364+ ISSM_MPI_Send(forcing, nv, ISSM_MPI_DOUBLE, 0, i, tocomms[i]);
365+ }
366+ }
367+ else{
368+ ISSM_MPI_Recv(&nv, 1, ISSM_MPI_INT, 0, modelid, fromcomm, &status);
369+ forcing=xNew<IssmDouble>(nv);
370+ ISSM_MPI_Recv(forcing, nv, ISSM_MPI_DOUBLE, 0, modelid, fromcomm, &status);
371+ }
372+ }
373+ /*}}}*/
374+
375+ /*On each ice cap, spread the forcing across cpus, and update the elements dataset accordingly: {{{*/
376+ if(modelid!=earthid){
377+
378+ ISSM_MPI_Bcast(&nv,1,ISSM_MPI_INT,0,IssmComm::GetComm());
379+ if(my_rank!=0)forcing=xNew<IssmDouble>(nv);
380+ ISSM_MPI_Bcast(forcing,nv,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
381+
382+ /*Plug into elements:*/
383+ InputUpdateFromVectorx(femmodel,forcing,forcingenum,VertexSIdEnum);
384+ }
385+ /*}}}*/
386+
387+ /*Free ressources:{{{*/
388+ if(forcingglobal)xDelete<IssmDouble>(forcingglobal);
389+ if(forcing)xDelete<IssmDouble>(forcing);
390+ if(transitions){
391+ for(int i=0;i<ntransitions;i++){
392+ IssmDouble* temp=transitions[i];
393+ xDelete<IssmDouble>(temp);
394+ }
395+ xDelete<IssmDouble*>(transitions);
396+ xDelete<int>(transitions_m);
397+ xDelete<int>(transitions_n);
398+ }
399+ /*}}}*/
400+
401+} /*}}}*/
Note: See TracBrowser for help on using the repository browser.