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

../trunk-jpl/src/c/modules/OceanExchangeDatax/OceanExchangeDatax.cpp

@@ -59 +59 @@
                 ISSM_MPI_Bcast(&oceantime,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
                 femmodel->parameters->SetParam(oceangridnxsize,OceanGridNxEnum);
                 femmodel->parameters->SetParam(oceangridnysize,OceanGridNyEnum);
+                oceangridx=xNew<IssmDouble>(ngrids_ocean);
+                oceangridy=xNew<IssmDouble>(ngrids_ocean);
                 if(my_rank==0){
-                        oceangridx = xNew<IssmDouble>(ngrids_ocean);
                         ISSM_MPI_Recv(oceangridx,ngrids_ocean,ISSM_MPI_DOUBLE,0,10001005,tomitgcmcomm,&status);
-                        oceangridy = xNew<IssmDouble>(ngrids_ocean);
                         ISSM_MPI_Recv(oceangridy,ngrids_ocean,ISSM_MPI_DOUBLE,0,10001006,tomitgcmcomm,&status);
 
                         /*Exchange varying parameters for the initialization*/
@@ -69 +69 @@
                         ISSM_MPI_Send(&time,1,ISSM_MPI_DOUBLE,0,10001001,tomitgcmcomm);
                         ISSM_MPI_Recv(&oceantime,1,ISSM_MPI_DOUBLE,0,10001002,tomitgcmcomm,&status);
                 }
-                if(my_rank!=0){
-                        oceangridx=xNew<IssmDouble>(ngrids_ocean);
-                        oceangridy=xNew<IssmDouble>(ngrids_ocean);
-                }
+                
                 ISSM_MPI_Bcast(oceangridx,ngrids_ocean,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
                 ISSM_MPI_Bcast(oceangridy,ngrids_ocean,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
                 femmodel->parameters->SetParam(oceangridx,ngrids_ocean,OceanGridXEnum);

../trunk-jpl/src/c/shared/Matrix/MatrixUtils.cpp

@@ -678 +678 @@
    Lchol should point to an IssmDouble* of same dimensions as A*/
 
         /*ensure zero-initialization*/
-        for(int ii=0;ii<(ndim*ndim);ii++) Lchol[ii]=0;;
+        for(int i=0;i<(ndim*ndim);i++) Lchol[i]=0;;
 
-        for(int ii=0;ii<ndim;ii++){
-                for(int jj{0};jj<=ii;jj++){
+        for(int i=0;i<ndim;i++){
+                for(int j=0;j<=i;j++){
                         IssmDouble sum=0.;
-                        for(int kk{0};kk<jj;kk++){
-                                sum += Lchol[ii*ndim+kk]*Lchol[jj*ndim+kk];
-                        }
-                        if(ii==jj){
-                                Lchol[ii*ndim+jj] = sqrt(A[ii*ndim+jj]-sum);
-                        }
-                        else{
-                                Lchol[ii*ndim+jj] = 1./Lchol[jj*ndim+jj] * (A[ii*ndim+jj]-sum);
-                        }
+                        for(int k=0;k<j;k++) sum += Lchol[i*ndim+k]*Lchol[j*ndim+k];
+                        
+                        if(i==j) Lchol[i*ndim+j] = sqrt(A[i*ndim+j]-sum);
+                        else Lchol[i*ndim+j]     = 1./Lchol[j*ndim+j] * (A[i*ndim+j]-sum);
                 }
         }
 } /*}}}*/

../trunk-jpl/src/c/classes/Elements/Element.cpp

@@ -3606 +3606 @@
 
         const int numvertices = this->GetNumberOfVertices();
         int         basinid;
-        IssmDouble  tspin,beta0_basin,beta1_basin,noisespin_basin; //initialize scalars 
+        IssmDouble  tspin,beta0_basin,beta1_basin,noisespin_basin;  
         IssmDouble* phi_basin   = xNew<IssmDouble>(arorder);
    IssmDouble* smbspin     = xNew<IssmDouble>(numvertices*arorder);
 
-        /*Get Basin ID*/
+        /*Get Basin ID and Basin coefficients*/
         this->GetInputValue(&basinid,SmbBasinsIdEnum);
-
-        for(int ii=0;ii<arorder;ii++) phi_basin[ii] = phi[basinid*arorder+ii];
-
+        for(int i=0;i<arorder;i++) phi_basin[i] = phi[basinid*arorder+i];
         beta0_basin   = beta0[basinid];
         beta1_basin   = beta1[basinid];
-        for(int jj=0;jj<nspin;jj++){    
-                tspin = starttime-((nspin-jj)*tstep_ar); //current time in spin-up loop
-      noisespin_basin = noisespin[jj*numbasins+basinid];
+
+        /*Loop over number of spin-up steps for all vertices*/  
+        for(int j=0;j<nspin;j++){       
+                tspin = starttime-((nspin-j)*tstep_ar); 
+      noisespin_basin = noisespin[j*numbasins+basinid];
       IssmDouble* oldsmbspin = xNew<IssmDouble>(numvertices*arorder);
-      for(int ii=0;ii<numvertices*arorder;ii++) oldsmbspin[ii]=smbspin[ii]; //copy smbspin in oldsmbspin
+      for(int i=0;i<numvertices*arorder;i++) oldsmbspin[i]=smbspin[i]; 
 
       for(int v=0;v<numvertices;v++){
          IssmDouble autoregressionterm = 0.;
-         for(int ii=0;ii<arorder;ii++) autoregressionterm += phi_basin[ii]*smbspin[v+ii*numvertices];
+         for(int i=0;i<arorder;i++) autoregressionterm += phi_basin[i]*smbspin[v+i*numvertices];
          smbspin[v] = beta0_basin+beta1_basin*(tspin-tinit_ar)+autoregressionterm+noisespin_basin; //compute newest values in smbspin
       }
 
-                /*correct older values in smbspin*/
-      for(int ii=0;ii<(arorder-1)*numvertices;ii++){
-                        smbspin[ii+numvertices]=oldsmbspin[ii];
-                } 
+                /*Correct older values in smbspin*/
+      for(int i=0;i<(arorder-1)*numvertices;i++) smbspin[i+numvertices]=oldsmbspin[i]; 
 
       xDelete<IssmDouble>(oldsmbspin); 
         }
@@ -3643 +3641 @@
    xDelete<IssmDouble>(phi_basin);
 }/*}}}*/
 void       Element::Smbautoregression(bool isstepforar,int arorder,IssmDouble telapsed_ar,IssmDouble* beta0,IssmDouble* beta1,IssmDouble* phi,IssmDouble* noiseterms){/*{{{*/
+        
         const int numvertices = this->GetNumberOfVertices();
         int         basinid,M,N;
-        IssmDouble  beta0_basin,beta1_basin,noise_basin; //initialize scalars 
+        IssmDouble  beta0_basin,beta1_basin,noise_basin;  
         IssmDouble* phi_basin  = xNew<IssmDouble>(arorder);
    IssmDouble* smblist    = xNew<IssmDouble>(numvertices);
-        IssmDouble* smbvaluesautoregression = NULL; //array for past SMB values that we are about to retrieve
+        IssmDouble* smbvaluesautoregression = NULL;
 
+        /*Get Basin ID and Basin coefficients*/
         this->GetInputValue(&basinid,SmbBasinsIdEnum);
-
-        for(int ii=0;ii<arorder;ii++){phi_basin[ii] = phi[basinid*arorder+ii];}
+        for(int ii=0;ii<arorder;ii++) phi_basin[ii] = phi[basinid*arorder+ii]; 
         beta0_basin   = beta0[basinid];
         beta1_basin   = beta1[basinid];
-        noise_basin   = noiseterms[basinid]; //note that noiseterms is computed at every timestep
-        this->inputs->GetArray(SmbValuesAutoregressionEnum,this->lid,&smbvaluesautoregression,&M); //get array of past SMB values to compute AR model
+        noise_basin   = noiseterms[basinid]; 
+        this->inputs->GetArray(SmbValuesAutoregressionEnum,this->lid,&smbvaluesautoregression,&M);
+
         /*If not AR model timestep: take the old SMB values*/
         if(isstepforar==false){
-                //VV do something with the lapse rate here if needed (12Oct2021)
-      for(int ii=0;ii<numvertices;ii++){smblist[ii]=smbvaluesautoregression[ii];}
+                /*VV do something with the lapse rate here if needed (12Oct2021)*/
+      for(int i=0;i<numvertices;i++) smblist[i]=smbvaluesautoregression[i]; 
         }
         /*If AR model timestep: compute SMB values on vertices from AR*/
         else{
                 for(int v=0;v<numvertices;v++){
 
-                        /*compute autoregressive term*/
+                        /*Compute autoregressive term*/
                         IssmDouble autoregressionterm=0.;
-                        for(int ii=0;ii<arorder;ii++){
-                                autoregressionterm += phi_basin[ii]*smbvaluesautoregression[v+ii*numvertices];
-                        }
+                        for(int ii=0;ii<arorder;ii++) autoregressionterm += phi_basin[ii]*smbvaluesautoregression[v+ii*numvertices];
 
-                        /*stochastic SMB value*/
+                        /*Stochastic SMB value*/
                         smblist[v] = beta0_basin+beta1_basin*telapsed_ar+autoregressionterm+noise_basin;
                 }
                 /*Update autoregression smb values*/
                 IssmDouble* temparray = xNew<IssmDouble>(numvertices*arorder);
-                for(int ii=0;ii<numvertices;ii++) temparray[ii] = smblist[ii]; //first store newly computed smb values
-                for(int ii=0;ii<(arorder-1)*numvertices;ii++){
-                        temparray[ii+numvertices] = smbvaluesautoregression[ii];
-                } //second shift older smb values
-                this->inputs->SetArrayInput(SmbValuesAutoregressionEnum,this->lid,temparray,numvertices*arorder); //save updated autoregression values
-                xDelete<IssmDouble>(temparray); //cleanup
+                /*Assign newest values and shift older values*/
+                for(int i=0;i<numvertices;i++) temparray[i] = smblist[i]; 
+                for(int i=0;i<(arorder-1)*numvertices;i++) temparray[i+numvertices] = smbvaluesautoregression[i];
+                this->inputs->SetArrayInput(SmbValuesAutoregressionEnum,this->lid,temparray,numvertices*arorder); 
+                xDelete<IssmDouble>(temparray);
         }
         /*Add input to element*/
         this->AddInput(SmbMassBalanceEnum,smblist,P1Enum);

../trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp

@@ -147 +147 @@
 
 }/*}}}*/
 void SmbautoregressionInitx(FemModel* femmodel){/*{{{*/
+        
         /*Initialization step of Smbautoregressionx*/
-        int M,N,Nphi,arorder,numbasins;
+        int M,N,Nphi,arorder,numbasins,my_rank;
         IssmDouble starttime,tstep_ar,tinit_ar;
+        femmodel->parameters->FindParam(&numbasins,SmbNumBasinsEnum);
+        femmodel->parameters->FindParam(&arorder,SmbAutoregressiveOrderEnum);
         IssmDouble* beta0    = xNew<IssmDouble>(numbasins);
         IssmDouble* beta1    = xNew<IssmDouble>(numbasins);
         IssmDouble* phi      = xNew<IssmDouble>(numbasins*arorder);
@@ -157 +160 @@
         femmodel->parameters->FindParam(&starttime,TimesteppingStartTimeEnum);
    femmodel->parameters->FindParam(&tstep_ar,SmbAutoregressionTimestepEnum);
         femmodel->parameters->FindParam(&tinit_ar,SmbAutoregressionInitialTimeEnum);
-        femmodel->parameters->FindParam(&numbasins,SmbNumBasinsEnum);
-        femmodel->parameters->FindParam(&arorder,SmbAutoregressiveOrderEnum);
         femmodel->parameters->FindParam(&beta0,&M,SmbBeta0Enum);    _assert_(M==numbasins);
         femmodel->parameters->FindParam(&beta1,&M,SmbBeta1Enum);    _assert_(M==numbasins);
         femmodel->parameters->FindParam(&phi,&M,&Nphi,SmbPhiEnum);  _assert_(M==numbasins); _assert_(Nphi==arorder);
         femmodel->parameters->FindParam(&covmat,&M,&N,SmbCovmatEnum); _assert_(M==numbasins); _assert_(N==numbasins);
+        
         /*AR model spin-up*/
-        int nspin{2*arorder+5}; //number of spin-up steps to be executed
-        IssmDouble* noisespin = xNewZeroInit<IssmDouble>(numbasins*nspin); //sample of basin-specific noise at each spinup step
-        for(int ii{0};ii<nspin;ii++){
-                IssmDouble* temparray = xNew<IssmDouble>(numbasins);
-                multivariateNormal(&temparray,numbasins,0.0,covmat);
-                for(int jj{0};jj<numbasins;jj++){noisespin[ii*numbasins+jj]=temparray[jj];}
-                xDelete<IssmDouble>(temparray);
+        int nspin{2*arorder+5}; 
+        IssmDouble* noisespin = xNewZeroInit<IssmDouble>(numbasins*nspin); 
+        my_rank=IssmComm::GetRank();
+        if(my_rank==0){
+                for(int i=0;i<nspin;i++){
+                        IssmDouble* temparray = xNew<IssmDouble>(numbasins);
+                        multivariateNormal(&temparray,numbasins,0.0,covmat);
+                        for(int j=0;j<numbasins;j++){noisespin[i*numbasins+j]=temparray[j];}
+                        xDelete<IssmDouble>(temparray);
+                }
         }
+        ISSM_MPI_Bcast(noisespin,numbasins*nspin,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
+        
         /*Initialize SmbValuesAutoregressionEnum*/
         for(Object* &object:femmodel->elements->objects){
       Element* element      = xDynamicCast<Element*>(object); //generate element object
                 element->SmbautoregressionInit(numbasins,arorder,nspin,starttime,tstep_ar,tinit_ar,beta0,beta1,phi,noisespin);
         }
+        
         /*Cleanup*/
         xDelete<IssmDouble>(beta0);
         xDelete<IssmDouble>(beta1);
@@ -191 +199 @@
         femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
         femmodel->parameters->FindParam(&starttime,TimesteppingStartTimeEnum);
    femmodel->parameters->FindParam(&tstep_ar,SmbAutoregressionTimestepEnum);
+
         /*Initialize module at first time step*/
         if(time<=starttime+dt){SmbautoregressionInitx(femmodel);}
         /*Determine if this is a time step for the AR model*/
-        bool isstepforar{false};
+        bool isstepforar = false;
 
         #ifndef _HAVE_AD_
-   if((fmod(time,tstep_ar)<fmod((time-dt),tstep_ar)) || (time<=starttime+dt) || tstep_ar==dt){isstepforar = true;}
+   if((fmod(time,tstep_ar)<fmod((time-dt),tstep_ar)) || (time<=starttime+dt) || tstep_ar==dt) isstepforar = true;
         #else
         _error_("not implemented yet");
         #endif
 
         /*Load parameters*/
-        int M,N,Nphi,arorder,numbasins;
+        int M,N,Nphi,arorder,numbasins,my_rank;
         femmodel->parameters->FindParam(&numbasins,SmbNumBasinsEnum);
         femmodel->parameters->FindParam(&arorder,SmbAutoregressiveOrderEnum);
         IssmDouble tinit_ar;
@@ -218 +227 @@
         femmodel->parameters->FindParam(&phi,&M,&Nphi,SmbPhiEnum);  _assert_(M==numbasins); _assert_(Nphi==arorder);
         femmodel->parameters->FindParam(&covmat,&M,&N,SmbCovmatEnum); _assert_(M==numbasins); _assert_(N==numbasins);
 
-        /*time elapsed with respect to AR model initial time*/
+        /*Time elapsed with respect to AR model initial time*/
         IssmDouble telapsed_ar = time-tinit_ar; 
 
         /*Before looping through elements: compute noise term specific to each basin from covmat*/
-        multivariateNormal(&noiseterms,numbasins,0.0,covmat);
+        my_rank=IssmComm::GetRank();
+        if(my_rank==0){
+                multivariateNormal(&noiseterms,numbasins,0.0,covmat);
+        }
+        ISSM_MPI_Bcast(noiseterms,numbasins,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
 
         /*Loop over each element to compute SMB at vertices*/
         for(Object* &object:femmodel->elements->objects){

../trunk-jpl/src/c/shared/Random/random.cpp

@@ -20 +20 @@
 /*}}}*/
 
 void univariateNormal(IssmPDouble* prand, IssmPDouble mean, IssmPDouble sdev) { /*{{{*/
-        /*univariateNormal      generates a random value follwoing Normal distribution*/
-        unsigned seed = std::chrono::steady_clock::now().time_since_epoch().count(); //random seed using time_since_epoch
-   std::default_random_engine generator(seed); //generator of random numbers
-   std::normal_distribution<IssmPDouble> normdistri(mean,sdev); //Normal probability distribution
+
+        /*Random seed using time_since_epoch*/
+        unsigned seed = std::chrono::steady_clock::now().time_since_epoch().count(); 
+   std::default_random_engine generator(seed);
+        /*Normal Probability Distribution*/
+   std::normal_distribution<IssmPDouble> normdistri(mean,sdev); 
         *prand = normdistri(generator);
 } /*}}}*/
 void multivariateNormal(IssmDouble** prand, int dim, IssmDouble mean, IssmDouble* covariancematrix) { /*{{{*/
-   IssmPDouble* sampleStandardNormal     = xNew<IssmPDouble>(dim);
+   
+        IssmPDouble* sampleStandardNormal    = xNew<IssmPDouble>(dim);
    IssmDouble* sampleMultivariateNormal = xNew<IssmDouble>(dim);
    IssmDouble* Lchol                    = xNewZeroInit<IssmDouble>(dim*dim);
 
-   for(int ii=0;ii<dim;ii++){univariateNormal(&(sampleStandardNormal[ii]),0.0,1.0);}
+   for(int i=0;i<dim;i++) univariateNormal(&(sampleStandardNormal[i]),0.0,1.0); 
    CholeskyRealPositiveDefinite(Lchol,covariancematrix,dim);
-   for(int ii=0;ii<dim;ii++){ //matrix by vector multiplication
-      /*entry-by-entry multiplication along matrix row*/
+   
+        /*Matrix by vector multiplication*/
+        for(int i=0;i<dim;i++){ 
+      /*Entry-by-entry multiplication along matrix row*/
       IssmDouble sum=0.;
-      for(int jj{0};jj<dim;jj++){
-         sum += sampleStandardNormal[jj]*Lchol[ii*dim+jj];
-      }
-      sampleMultivariateNormal[ii] = mean+sum;
+      for(int j=0;j<dim;j++) sum += sampleStandardNormal[j]*Lchol[i*dim+j]; 
+      sampleMultivariateNormal[i] = mean+sum;
    }
 
    /*Assign output pointer and cleanup*/
@@ -48 +51 @@
    xDelete<IssmDouble>(Lchol);
 } /*}}}*/
 void multivariateNormal(IssmDouble** prand, int dim, IssmDouble* mean, IssmDouble* covariancematrix) { /*{{{*/
-        IssmPDouble* sampleStandardNormal     = xNew<IssmPDouble>(dim);
+        
+        IssmPDouble* sampleStandardNormal    = xNew<IssmPDouble>(dim);
         IssmDouble* sampleMultivariateNormal = xNew<IssmDouble>(dim);
         IssmDouble* Lchol                    = xNewZeroInit<IssmDouble>(dim*dim);
-        for(int ii=0;ii<dim;ii++) univariateNormal(&(sampleStandardNormal[ii]),0.0,1.0);
+        for(int i=0;i<dim;i++) univariateNormal(&(sampleStandardNormal[i]),0.0,1.0);
 
         CholeskyRealPositiveDefinite(Lchol,covariancematrix,dim);
 
-        for(int ii=0;ii<dim;ii++){ //matrix by vector multiplication
+        /*Matrix by vector multiplication*/
+        for(int i=0;i<dim;i++){
                 IssmDouble sum = 0.;
-      for(int jj=0.;jj<dim;jj++){
-                        sum += sampleStandardNormal[jj]*Lchol[ii*dim+jj];
-                }
-      sampleMultivariateNormal[ii] = mean[ii]+sum; //assign value
+      for(int j=0;j<dim;j++) sum += sampleStandardNormal[j]*Lchol[i*dim+j]; 
+      sampleMultivariateNormal[i] = mean[i]+sum;
         }
+   
+        /*Assign output pointer and cleanup*/
         *prand = sampleMultivariateNormal;
         xDelete<IssmPDouble>(sampleStandardNormal);
         xDelete<IssmDouble>(Lchol);