Index: ../trunk-jpl/src/c/modules/OceanExchangeDatax/OceanExchangeDatax.cpp =================================================================== --- ../trunk-jpl/src/c/modules/OceanExchangeDatax/OceanExchangeDatax.cpp (revision 26481) +++ ../trunk-jpl/src/c/modules/OceanExchangeDatax/OceanExchangeDatax.cpp (revision 26482) @@ -59,10 +59,10 @@ ISSM_MPI_Bcast(&oceantime,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); femmodel->parameters->SetParam(oceangridnxsize,OceanGridNxEnum); femmodel->parameters->SetParam(oceangridnysize,OceanGridNyEnum); + oceangridx=xNew(ngrids_ocean); + oceangridy=xNew(ngrids_ocean); if(my_rank==0){ - oceangridx = xNew(ngrids_ocean); ISSM_MPI_Recv(oceangridx,ngrids_ocean,ISSM_MPI_DOUBLE,0,10001005,tomitgcmcomm,&status); - oceangridy = xNew(ngrids_ocean); ISSM_MPI_Recv(oceangridy,ngrids_ocean,ISSM_MPI_DOUBLE,0,10001006,tomitgcmcomm,&status); /*Exchange varying parameters for the initialization*/ @@ -69,10 +69,7 @@ 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(ngrids_ocean); - oceangridy=xNew(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); Index: ../trunk-jpl/src/c/shared/Matrix/MatrixUtils.cpp =================================================================== --- ../trunk-jpl/src/c/shared/Matrix/MatrixUtils.cpp (revision 26481) +++ ../trunk-jpl/src/c/shared/Matrix/MatrixUtils.cpp (revision 26482) @@ -678,20 +678,15 @@ 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;iiGetNumberOfVertices(); int basinid; - IssmDouble tspin,beta0_basin,beta1_basin,noisespin_basin; //initialize scalars + IssmDouble tspin,beta0_basin,beta1_basin,noisespin_basin; IssmDouble* phi_basin = xNew(arorder); IssmDouble* smbspin = xNew(numvertices*arorder); - /*Get Basin ID*/ + /*Get Basin ID and Basin coefficients*/ this->GetInputValue(&basinid,SmbBasinsIdEnum); - - for(int ii=0;ii(numvertices*arorder); - for(int ii=0;ii(oldsmbspin); } @@ -3643,46 +3641,45 @@ xDelete(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(arorder); IssmDouble* smblist = xNew(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;iiinputs->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*arorder); - for(int ii=0;iiinputs->SetArrayInput(SmbValuesAutoregressionEnum,this->lid,temparray,numvertices*arorder); //save updated autoregression values - xDelete(temparray); //cleanup + /*Assign newest values and shift older values*/ + for(int i=0;iinputs->SetArrayInput(SmbValuesAutoregressionEnum,this->lid,temparray,numvertices*arorder); + xDelete(temparray); } /*Add input to element*/ this->AddInput(SmbMassBalanceEnum,smblist,P1Enum); Index: ../trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp =================================================================== --- ../trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp (revision 26481) +++ ../trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp (revision 26482) @@ -147,9 +147,12 @@ }/*}}}*/ 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(numbasins); IssmDouble* beta1 = xNew(numbasins); IssmDouble* phi = xNew(numbasins*arorder); @@ -157,26 +160,31 @@ 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(numbasins*nspin); //sample of basin-specific noise at each spinup step - for(int ii{0};ii(numbasins); - multivariateNormal(&temparray,numbasins,0.0,covmat); - for(int jj{0};jj(temparray); + int nspin{2*arorder+5}; + IssmDouble* noisespin = xNewZeroInit(numbasins*nspin); + my_rank=IssmComm::GetRank(); + if(my_rank==0){ + for(int i=0;i(numbasins); + multivariateNormal(&temparray,numbasins,0.0,covmat); + for(int j=0;j(temparray); + } } + ISSM_MPI_Bcast(noisespin,numbasins*nspin,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); + /*Initialize SmbValuesAutoregressionEnum*/ for(Object* &object:femmodel->elements->objects){ Element* element = xDynamicCast(object); //generate element object element->SmbautoregressionInit(numbasins,arorder,nspin,starttime,tstep_ar,tinit_ar,beta0,beta1,phi,noisespin); } + /*Cleanup*/ xDelete(beta0); xDelete(beta1); @@ -191,19 +199,20 @@ 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)parameters->FindParam(&numbasins,SmbNumBasinsEnum); femmodel->parameters->FindParam(&arorder,SmbAutoregressiveOrderEnum); IssmDouble tinit_ar; @@ -218,11 +227,15 @@ 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){ Index: ../trunk-jpl/src/c/shared/Random/random.cpp =================================================================== --- ../trunk-jpl/src/c/shared/Random/random.cpp (revision 26481) +++ ../trunk-jpl/src/c/shared/Random/random.cpp (revision 26482) @@ -20,26 +20,29 @@ /*}}}*/ 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 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 normdistri(mean,sdev); *prand = normdistri(generator); } /*}}}*/ void multivariateNormal(IssmDouble** prand, int dim, IssmDouble mean, IssmDouble* covariancematrix) { /*{{{*/ - IssmPDouble* sampleStandardNormal = xNew(dim); + + IssmPDouble* sampleStandardNormal = xNew(dim); IssmDouble* sampleMultivariateNormal = xNew(dim); IssmDouble* Lchol = xNewZeroInit(dim*dim); - for(int ii=0;ii(Lchol); } /*}}}*/ void multivariateNormal(IssmDouble** prand, int dim, IssmDouble* mean, IssmDouble* covariancematrix) { /*{{{*/ - IssmPDouble* sampleStandardNormal = xNew(dim); + + IssmPDouble* sampleStandardNormal = xNew(dim); IssmDouble* sampleMultivariateNormal = xNew(dim); IssmDouble* Lchol = xNewZeroInit(dim*dim); - for(int ii=0;ii(sampleStandardNormal); xDelete(Lchol);