source: issm/oecreview/Archive/25834-26739/ISSM-26481-26482.diff@ 26740

Last change on this file since 26740 was 26740, checked in by Mathieu Morlighem, 3 years ago

CHG: added 25834-26739
File size: 16.4 KB
RevLine 
[26740]1Index: ../trunk-jpl/src/c/modules/OceanExchangeDatax/OceanExchangeDatax.cpp
2===================================================================
3--- ../trunk-jpl/src/c/modules/OceanExchangeDatax/OceanExchangeDatax.cpp (revision 26481)
4+++ ../trunk-jpl/src/c/modules/OceanExchangeDatax/OceanExchangeDatax.cpp (revision 26482)
5@@ -59,10 +59,10 @@
6 ISSM_MPI_Bcast(&oceantime,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
7 femmodel->parameters->SetParam(oceangridnxsize,OceanGridNxEnum);
8 femmodel->parameters->SetParam(oceangridnysize,OceanGridNyEnum);
9+ oceangridx=xNew<IssmDouble>(ngrids_ocean);
10+ oceangridy=xNew<IssmDouble>(ngrids_ocean);
11 if(my_rank==0){
12- oceangridx = xNew<IssmDouble>(ngrids_ocean);
13 ISSM_MPI_Recv(oceangridx,ngrids_ocean,ISSM_MPI_DOUBLE,0,10001005,tomitgcmcomm,&status);
14- oceangridy = xNew<IssmDouble>(ngrids_ocean);
15 ISSM_MPI_Recv(oceangridy,ngrids_ocean,ISSM_MPI_DOUBLE,0,10001006,tomitgcmcomm,&status);
16
17 /*Exchange varying parameters for the initialization*/
18@@ -69,10 +69,7 @@
19 ISSM_MPI_Send(&time,1,ISSM_MPI_DOUBLE,0,10001001,tomitgcmcomm);
20 ISSM_MPI_Recv(&oceantime,1,ISSM_MPI_DOUBLE,0,10001002,tomitgcmcomm,&status);
21 }
22- if(my_rank!=0){
23- oceangridx=xNew<IssmDouble>(ngrids_ocean);
24- oceangridy=xNew<IssmDouble>(ngrids_ocean);
25- }
26+
27 ISSM_MPI_Bcast(oceangridx,ngrids_ocean,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
28 ISSM_MPI_Bcast(oceangridy,ngrids_ocean,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
29 femmodel->parameters->SetParam(oceangridx,ngrids_ocean,OceanGridXEnum);
30Index: ../trunk-jpl/src/c/shared/Matrix/MatrixUtils.cpp
31===================================================================
32--- ../trunk-jpl/src/c/shared/Matrix/MatrixUtils.cpp (revision 26481)
33+++ ../trunk-jpl/src/c/shared/Matrix/MatrixUtils.cpp (revision 26482)
34@@ -678,20 +678,15 @@
35 Lchol should point to an IssmDouble* of same dimensions as A*/
36
37 /*ensure zero-initialization*/
38- for(int ii=0;ii<(ndim*ndim);ii++) Lchol[ii]=0;;
39+ for(int i=0;i<(ndim*ndim);i++) Lchol[i]=0;;
40
41- for(int ii=0;ii<ndim;ii++){
42- for(int jj{0};jj<=ii;jj++){
43+ for(int i=0;i<ndim;i++){
44+ for(int j=0;j<=i;j++){
45 IssmDouble sum=0.;
46- for(int kk{0};kk<jj;kk++){
47- sum += Lchol[ii*ndim+kk]*Lchol[jj*ndim+kk];
48- }
49- if(ii==jj){
50- Lchol[ii*ndim+jj] = sqrt(A[ii*ndim+jj]-sum);
51- }
52- else{
53- Lchol[ii*ndim+jj] = 1./Lchol[jj*ndim+jj] * (A[ii*ndim+jj]-sum);
54- }
55+ for(int k=0;k<j;k++) sum += Lchol[i*ndim+k]*Lchol[j*ndim+k];
56+
57+ if(i==j) Lchol[i*ndim+j] = sqrt(A[i*ndim+j]-sum);
58+ else Lchol[i*ndim+j] = 1./Lchol[j*ndim+j] * (A[i*ndim+j]-sum);
59 }
60 }
61 } /*}}}*/
62Index: ../trunk-jpl/src/c/classes/Elements/Element.cpp
63===================================================================
64--- ../trunk-jpl/src/c/classes/Elements/Element.cpp (revision 26481)
65+++ ../trunk-jpl/src/c/classes/Elements/Element.cpp (revision 26482)
66@@ -3606,33 +3606,31 @@
67
68 const int numvertices = this->GetNumberOfVertices();
69 int basinid;
70- IssmDouble tspin,beta0_basin,beta1_basin,noisespin_basin; //initialize scalars
71+ IssmDouble tspin,beta0_basin,beta1_basin,noisespin_basin;
72 IssmDouble* phi_basin = xNew<IssmDouble>(arorder);
73 IssmDouble* smbspin = xNew<IssmDouble>(numvertices*arorder);
74
75- /*Get Basin ID*/
76+ /*Get Basin ID and Basin coefficients*/
77 this->GetInputValue(&basinid,SmbBasinsIdEnum);
78-
79- for(int ii=0;ii<arorder;ii++) phi_basin[ii] = phi[basinid*arorder+ii];
80-
81+ for(int i=0;i<arorder;i++) phi_basin[i] = phi[basinid*arorder+i];
82 beta0_basin = beta0[basinid];
83 beta1_basin = beta1[basinid];
84- for(int jj=0;jj<nspin;jj++){
85- tspin = starttime-((nspin-jj)*tstep_ar); //current time in spin-up loop
86- noisespin_basin = noisespin[jj*numbasins+basinid];
87+
88+ /*Loop over number of spin-up steps for all vertices*/
89+ for(int j=0;j<nspin;j++){
90+ tspin = starttime-((nspin-j)*tstep_ar);
91+ noisespin_basin = noisespin[j*numbasins+basinid];
92 IssmDouble* oldsmbspin = xNew<IssmDouble>(numvertices*arorder);
93- for(int ii=0;ii<numvertices*arorder;ii++) oldsmbspin[ii]=smbspin[ii]; //copy smbspin in oldsmbspin
94+ for(int i=0;i<numvertices*arorder;i++) oldsmbspin[i]=smbspin[i];
95
96 for(int v=0;v<numvertices;v++){
97 IssmDouble autoregressionterm = 0.;
98- for(int ii=0;ii<arorder;ii++) autoregressionterm += phi_basin[ii]*smbspin[v+ii*numvertices];
99+ for(int i=0;i<arorder;i++) autoregressionterm += phi_basin[i]*smbspin[v+i*numvertices];
100 smbspin[v] = beta0_basin+beta1_basin*(tspin-tinit_ar)+autoregressionterm+noisespin_basin; //compute newest values in smbspin
101 }
102
103- /*correct older values in smbspin*/
104- for(int ii=0;ii<(arorder-1)*numvertices;ii++){
105- smbspin[ii+numvertices]=oldsmbspin[ii];
106- }
107+ /*Correct older values in smbspin*/
108+ for(int i=0;i<(arorder-1)*numvertices;i++) smbspin[i+numvertices]=oldsmbspin[i];
109
110 xDelete<IssmDouble>(oldsmbspin);
111 }
112@@ -3643,46 +3641,45 @@
113 xDelete<IssmDouble>(phi_basin);
114 }/*}}}*/
115 void Element::Smbautoregression(bool isstepforar,int arorder,IssmDouble telapsed_ar,IssmDouble* beta0,IssmDouble* beta1,IssmDouble* phi,IssmDouble* noiseterms){/*{{{*/
116+
117 const int numvertices = this->GetNumberOfVertices();
118 int basinid,M,N;
119- IssmDouble beta0_basin,beta1_basin,noise_basin; //initialize scalars
120+ IssmDouble beta0_basin,beta1_basin,noise_basin;
121 IssmDouble* phi_basin = xNew<IssmDouble>(arorder);
122 IssmDouble* smblist = xNew<IssmDouble>(numvertices);
123- IssmDouble* smbvaluesautoregression = NULL; //array for past SMB values that we are about to retrieve
124+ IssmDouble* smbvaluesautoregression = NULL;
125
126+ /*Get Basin ID and Basin coefficients*/
127 this->GetInputValue(&basinid,SmbBasinsIdEnum);
128-
129- for(int ii=0;ii<arorder;ii++){phi_basin[ii] = phi[basinid*arorder+ii];}
130+ for(int ii=0;ii<arorder;ii++) phi_basin[ii] = phi[basinid*arorder+ii];
131 beta0_basin = beta0[basinid];
132 beta1_basin = beta1[basinid];
133- noise_basin = noiseterms[basinid]; //note that noiseterms is computed at every timestep
134- this->inputs->GetArray(SmbValuesAutoregressionEnum,this->lid,&smbvaluesautoregression,&M); //get array of past SMB values to compute AR model
135+ noise_basin = noiseterms[basinid];
136+ this->inputs->GetArray(SmbValuesAutoregressionEnum,this->lid,&smbvaluesautoregression,&M);
137+
138 /*If not AR model timestep: take the old SMB values*/
139 if(isstepforar==false){
140- //VV do something with the lapse rate here if needed (12Oct2021)
141- for(int ii=0;ii<numvertices;ii++){smblist[ii]=smbvaluesautoregression[ii];}
142+ /*VV do something with the lapse rate here if needed (12Oct2021)*/
143+ for(int i=0;i<numvertices;i++) smblist[i]=smbvaluesautoregression[i];
144 }
145 /*If AR model timestep: compute SMB values on vertices from AR*/
146 else{
147 for(int v=0;v<numvertices;v++){
148
149- /*compute autoregressive term*/
150+ /*Compute autoregressive term*/
151 IssmDouble autoregressionterm=0.;
152- for(int ii=0;ii<arorder;ii++){
153- autoregressionterm += phi_basin[ii]*smbvaluesautoregression[v+ii*numvertices];
154- }
155+ for(int ii=0;ii<arorder;ii++) autoregressionterm += phi_basin[ii]*smbvaluesautoregression[v+ii*numvertices];
156
157- /*stochastic SMB value*/
158+ /*Stochastic SMB value*/
159 smblist[v] = beta0_basin+beta1_basin*telapsed_ar+autoregressionterm+noise_basin;
160 }
161 /*Update autoregression smb values*/
162 IssmDouble* temparray = xNew<IssmDouble>(numvertices*arorder);
163- for(int ii=0;ii<numvertices;ii++) temparray[ii] = smblist[ii]; //first store newly computed smb values
164- for(int ii=0;ii<(arorder-1)*numvertices;ii++){
165- temparray[ii+numvertices] = smbvaluesautoregression[ii];
166- } //second shift older smb values
167- this->inputs->SetArrayInput(SmbValuesAutoregressionEnum,this->lid,temparray,numvertices*arorder); //save updated autoregression values
168- xDelete<IssmDouble>(temparray); //cleanup
169+ /*Assign newest values and shift older values*/
170+ for(int i=0;i<numvertices;i++) temparray[i] = smblist[i];
171+ for(int i=0;i<(arorder-1)*numvertices;i++) temparray[i+numvertices] = smbvaluesautoregression[i];
172+ this->inputs->SetArrayInput(SmbValuesAutoregressionEnum,this->lid,temparray,numvertices*arorder);
173+ xDelete<IssmDouble>(temparray);
174 }
175 /*Add input to element*/
176 this->AddInput(SmbMassBalanceEnum,smblist,P1Enum);
177Index: ../trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp
178===================================================================
179--- ../trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp (revision 26481)
180+++ ../trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp (revision 26482)
181@@ -147,9 +147,12 @@
182
183 }/*}}}*/
184 void SmbautoregressionInitx(FemModel* femmodel){/*{{{*/
185+
186 /*Initialization step of Smbautoregressionx*/
187- int M,N,Nphi,arorder,numbasins;
188+ int M,N,Nphi,arorder,numbasins,my_rank;
189 IssmDouble starttime,tstep_ar,tinit_ar;
190+ femmodel->parameters->FindParam(&numbasins,SmbNumBasinsEnum);
191+ femmodel->parameters->FindParam(&arorder,SmbAutoregressiveOrderEnum);
192 IssmDouble* beta0 = xNew<IssmDouble>(numbasins);
193 IssmDouble* beta1 = xNew<IssmDouble>(numbasins);
194 IssmDouble* phi = xNew<IssmDouble>(numbasins*arorder);
195@@ -157,26 +160,31 @@
196 femmodel->parameters->FindParam(&starttime,TimesteppingStartTimeEnum);
197 femmodel->parameters->FindParam(&tstep_ar,SmbAutoregressionTimestepEnum);
198 femmodel->parameters->FindParam(&tinit_ar,SmbAutoregressionInitialTimeEnum);
199- femmodel->parameters->FindParam(&numbasins,SmbNumBasinsEnum);
200- femmodel->parameters->FindParam(&arorder,SmbAutoregressiveOrderEnum);
201 femmodel->parameters->FindParam(&beta0,&M,SmbBeta0Enum); _assert_(M==numbasins);
202 femmodel->parameters->FindParam(&beta1,&M,SmbBeta1Enum); _assert_(M==numbasins);
203 femmodel->parameters->FindParam(&phi,&M,&Nphi,SmbPhiEnum); _assert_(M==numbasins); _assert_(Nphi==arorder);
204 femmodel->parameters->FindParam(&covmat,&M,&N,SmbCovmatEnum); _assert_(M==numbasins); _assert_(N==numbasins);
205+
206 /*AR model spin-up*/
207- int nspin{2*arorder+5}; //number of spin-up steps to be executed
208- IssmDouble* noisespin = xNewZeroInit<IssmDouble>(numbasins*nspin); //sample of basin-specific noise at each spinup step
209- for(int ii{0};ii<nspin;ii++){
210- IssmDouble* temparray = xNew<IssmDouble>(numbasins);
211- multivariateNormal(&temparray,numbasins,0.0,covmat);
212- for(int jj{0};jj<numbasins;jj++){noisespin[ii*numbasins+jj]=temparray[jj];}
213- xDelete<IssmDouble>(temparray);
214+ int nspin{2*arorder+5};
215+ IssmDouble* noisespin = xNewZeroInit<IssmDouble>(numbasins*nspin);
216+ my_rank=IssmComm::GetRank();
217+ if(my_rank==0){
218+ for(int i=0;i<nspin;i++){
219+ IssmDouble* temparray = xNew<IssmDouble>(numbasins);
220+ multivariateNormal(&temparray,numbasins,0.0,covmat);
221+ for(int j=0;j<numbasins;j++){noisespin[i*numbasins+j]=temparray[j];}
222+ xDelete<IssmDouble>(temparray);
223+ }
224 }
225+ ISSM_MPI_Bcast(noisespin,numbasins*nspin,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
226+
227 /*Initialize SmbValuesAutoregressionEnum*/
228 for(Object* &object:femmodel->elements->objects){
229 Element* element = xDynamicCast<Element*>(object); //generate element object
230 element->SmbautoregressionInit(numbasins,arorder,nspin,starttime,tstep_ar,tinit_ar,beta0,beta1,phi,noisespin);
231 }
232+
233 /*Cleanup*/
234 xDelete<IssmDouble>(beta0);
235 xDelete<IssmDouble>(beta1);
236@@ -191,19 +199,20 @@
237 femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
238 femmodel->parameters->FindParam(&starttime,TimesteppingStartTimeEnum);
239 femmodel->parameters->FindParam(&tstep_ar,SmbAutoregressionTimestepEnum);
240+
241 /*Initialize module at first time step*/
242 if(time<=starttime+dt){SmbautoregressionInitx(femmodel);}
243 /*Determine if this is a time step for the AR model*/
244- bool isstepforar{false};
245+ bool isstepforar = false;
246
247 #ifndef _HAVE_AD_
248- if((fmod(time,tstep_ar)<fmod((time-dt),tstep_ar)) || (time<=starttime+dt) || tstep_ar==dt){isstepforar = true;}
249+ if((fmod(time,tstep_ar)<fmod((time-dt),tstep_ar)) || (time<=starttime+dt) || tstep_ar==dt) isstepforar = true;
250 #else
251 _error_("not implemented yet");
252 #endif
253
254 /*Load parameters*/
255- int M,N,Nphi,arorder,numbasins;
256+ int M,N,Nphi,arorder,numbasins,my_rank;
257 femmodel->parameters->FindParam(&numbasins,SmbNumBasinsEnum);
258 femmodel->parameters->FindParam(&arorder,SmbAutoregressiveOrderEnum);
259 IssmDouble tinit_ar;
260@@ -218,11 +227,15 @@
261 femmodel->parameters->FindParam(&phi,&M,&Nphi,SmbPhiEnum); _assert_(M==numbasins); _assert_(Nphi==arorder);
262 femmodel->parameters->FindParam(&covmat,&M,&N,SmbCovmatEnum); _assert_(M==numbasins); _assert_(N==numbasins);
263
264- /*time elapsed with respect to AR model initial time*/
265+ /*Time elapsed with respect to AR model initial time*/
266 IssmDouble telapsed_ar = time-tinit_ar;
267
268 /*Before looping through elements: compute noise term specific to each basin from covmat*/
269- multivariateNormal(&noiseterms,numbasins,0.0,covmat);
270+ my_rank=IssmComm::GetRank();
271+ if(my_rank==0){
272+ multivariateNormal(&noiseterms,numbasins,0.0,covmat);
273+ }
274+ ISSM_MPI_Bcast(noiseterms,numbasins,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
275
276 /*Loop over each element to compute SMB at vertices*/
277 for(Object* &object:femmodel->elements->objects){
278Index: ../trunk-jpl/src/c/shared/Random/random.cpp
279===================================================================
280--- ../trunk-jpl/src/c/shared/Random/random.cpp (revision 26481)
281+++ ../trunk-jpl/src/c/shared/Random/random.cpp (revision 26482)
282@@ -20,26 +20,29 @@
283 /*}}}*/
284
285 void univariateNormal(IssmPDouble* prand, IssmPDouble mean, IssmPDouble sdev) { /*{{{*/
286- /*univariateNormal generates a random value follwoing Normal distribution*/
287- unsigned seed = std::chrono::steady_clock::now().time_since_epoch().count(); //random seed using time_since_epoch
288- std::default_random_engine generator(seed); //generator of random numbers
289- std::normal_distribution<IssmPDouble> normdistri(mean,sdev); //Normal probability distribution
290+
291+ /*Random seed using time_since_epoch*/
292+ unsigned seed = std::chrono::steady_clock::now().time_since_epoch().count();
293+ std::default_random_engine generator(seed);
294+ /*Normal Probability Distribution*/
295+ std::normal_distribution<IssmPDouble> normdistri(mean,sdev);
296 *prand = normdistri(generator);
297 } /*}}}*/
298 void multivariateNormal(IssmDouble** prand, int dim, IssmDouble mean, IssmDouble* covariancematrix) { /*{{{*/
299- IssmPDouble* sampleStandardNormal = xNew<IssmPDouble>(dim);
300+
301+ IssmPDouble* sampleStandardNormal = xNew<IssmPDouble>(dim);
302 IssmDouble* sampleMultivariateNormal = xNew<IssmDouble>(dim);
303 IssmDouble* Lchol = xNewZeroInit<IssmDouble>(dim*dim);
304
305- for(int ii=0;ii<dim;ii++){univariateNormal(&(sampleStandardNormal[ii]),0.0,1.0);}
306+ for(int i=0;i<dim;i++) univariateNormal(&(sampleStandardNormal[i]),0.0,1.0);
307 CholeskyRealPositiveDefinite(Lchol,covariancematrix,dim);
308- for(int ii=0;ii<dim;ii++){ //matrix by vector multiplication
309- /*entry-by-entry multiplication along matrix row*/
310+
311+ /*Matrix by vector multiplication*/
312+ for(int i=0;i<dim;i++){
313+ /*Entry-by-entry multiplication along matrix row*/
314 IssmDouble sum=0.;
315- for(int jj{0};jj<dim;jj++){
316- sum += sampleStandardNormal[jj]*Lchol[ii*dim+jj];
317- }
318- sampleMultivariateNormal[ii] = mean+sum;
319+ for(int j=0;j<dim;j++) sum += sampleStandardNormal[j]*Lchol[i*dim+j];
320+ sampleMultivariateNormal[i] = mean+sum;
321 }
322
323 /*Assign output pointer and cleanup*/
324@@ -48,20 +51,22 @@
325 xDelete<IssmDouble>(Lchol);
326 } /*}}}*/
327 void multivariateNormal(IssmDouble** prand, int dim, IssmDouble* mean, IssmDouble* covariancematrix) { /*{{{*/
328- IssmPDouble* sampleStandardNormal = xNew<IssmPDouble>(dim);
329+
330+ IssmPDouble* sampleStandardNormal = xNew<IssmPDouble>(dim);
331 IssmDouble* sampleMultivariateNormal = xNew<IssmDouble>(dim);
332 IssmDouble* Lchol = xNewZeroInit<IssmDouble>(dim*dim);
333- for(int ii=0;ii<dim;ii++) univariateNormal(&(sampleStandardNormal[ii]),0.0,1.0);
334+ for(int i=0;i<dim;i++) univariateNormal(&(sampleStandardNormal[i]),0.0,1.0);
335
336 CholeskyRealPositiveDefinite(Lchol,covariancematrix,dim);
337
338- for(int ii=0;ii<dim;ii++){ //matrix by vector multiplication
339+ /*Matrix by vector multiplication*/
340+ for(int i=0;i<dim;i++){
341 IssmDouble sum = 0.;
342- for(int jj=0.;jj<dim;jj++){
343- sum += sampleStandardNormal[jj]*Lchol[ii*dim+jj];
344- }
345- sampleMultivariateNormal[ii] = mean[ii]+sum; //assign value
346+ for(int j=0;j<dim;j++) sum += sampleStandardNormal[j]*Lchol[i*dim+j];
347+ sampleMultivariateNormal[i] = mean[i]+sum;
348 }
349+
350+ /*Assign output pointer and cleanup*/
351 *prand = sampleMultivariateNormal;
352 xDelete<IssmPDouble>(sampleStandardNormal);
353 xDelete<IssmDouble>(Lchol);
Note: See TracBrowser for help on using the repository browser.