Changeset 26479

Timestamp:

10/13/21 09:57:34 (3 years ago)

Author:

Mathieu Morlighem

Message:

CHG: fixing style and AD after Vincent's commit

Location:

issm/trunk-jpl/src/c

Files:

• classes/Elements/Element.cpp (modified) (4 diffs)
• modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp (modified) (2 diffs)
• shared/Matrix/MatrixUtils.cpp (modified) (1 diff)
• shared/Random/random.cpp (modified) (1 diff)

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

@@ -3604 +3604 @@
 /*}}}*/
 void       Element::SmbautoregressionInit(int numbasins,int arorder,int nspin,IssmDouble starttime,IssmDouble tstep_ar,IssmDouble tinit_ar,IssmDouble* beta0,IssmDouble* beta1,IssmDouble* phi,IssmDouble* noisespin){/*{{{*/
+
         const int numvertices = this->GetNumberOfVertices();
         int         basinid;
         IssmDouble  tspin,beta0_basin,beta1_basin,noisespin_basin; //initialize scalars 
-        IssmDouble* phi_basin               = xNew<IssmDouble>(arorder);
-   IssmDouble* smbspin                 = xNew<IssmDouble>(numvertices*arorder);
+        IssmDouble* phi_basin   = xNew<IssmDouble>(arorder);
+   IssmDouble* smbspin     = xNew<IssmDouble>(numvertices*arorder);
+
+        /*Get Basin ID*/
         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];
-        for(int jj{0};jj<nspin;jj++){   
+        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];
-      IssmDouble* oldsmbspin = xNewZeroInit<IssmDouble>(numvertices*arorder);
-      for(int ii{0};ii<numvertices*arorder;ii++){oldsmbspin[ii]=smbspin[ii];} //copy smbspin in oldsmbspin
+      IssmDouble* oldsmbspin = xNew<IssmDouble>(numvertices*arorder);
+      for(int ii=0;ii<numvertices*arorder;ii++) oldsmbspin[ii]=smbspin[ii]; //copy smbspin in oldsmbspin
+
       for(int v=0;v<numvertices;v++){
-         double autoregressionterm{0.0};
-         for(int ii{0};ii<arorder;ii++){autoregressionterm += phi_basin[ii]*smbspin[v+ii*numvertices];}
+         IssmDouble autoregressionterm = 0.;
+         for(int ii=0;ii<arorder;ii++) autoregressionterm += phi_basin[ii]*smbspin[v+ii*numvertices];
          smbspin[v] = beta0_basin+beta1_basin*(tspin-tinit_ar)+autoregressionterm+noisespin_basin; //compute newest values in smbspin
       }
-      for(int ii{0};ii<(arorder-1)*numvertices;ii++){smbspin[ii+numvertices]=oldsmbspin[ii];} //correct older values in smbspin
-      xDelete<IssmDouble>(oldsmbspin); //cleanup
-        }
-        this->inputs->SetArrayInput(SmbValuesAutoregressionEnum,this->lid,smbspin,numvertices*arorder);//save spin-up autoregression values
-   xDelete<IssmDouble>(smbspin); //cleanup
-   xDelete<IssmDouble>(phi_basin); //cleanup
+
+                /*correct older values in smbspin*/
+      for(int ii=0;ii<(arorder-1)*numvertices;ii++){
+                        smbspin[ii+numvertices]=oldsmbspin[ii];
+                } 
+
+      xDelete<IssmDouble>(oldsmbspin); 
+        }
+        this->inputs->SetArrayInput(SmbValuesAutoregressionEnum,this->lid,smbspin,numvertices*arorder);
+
+        /*Cleanup and return*/
+   xDelete<IssmDouble>(smbspin); 
+   xDelete<IssmDouble>(phi_basin);
 }/*}}}*/
 void       Element::Smbautoregression(bool isstepforar,int arorder,IssmDouble telapsed_ar,IssmDouble* beta0,IssmDouble* beta1,IssmDouble* phi,IssmDouble* noiseterms){/*{{{*/
@@ -3634 +3647 @@
         int         basinid,M,N;
         IssmDouble  beta0_basin,beta1_basin,noise_basin; //initialize scalars 
-        IssmDouble* phi_basin               = xNew<IssmDouble>(arorder);
-   IssmDouble* smblist                 = xNew<IssmDouble>(numvertices);
+        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
+
         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];
@@ -3646 +3661 @@
         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];}
+      for(int ii=0;ii<numvertices;ii++){smblist[ii]=smbvaluesautoregression[ii];}
         }
         /*If AR model timestep: compute SMB values on vertices from AR*/
         else{
                 for(int v=0;v<numvertices;v++){
-                        double autoregressionterm{0.0};
-                        for(int ii{0};ii<arorder;ii++){autoregressionterm += phi_basin[ii]*smbvaluesautoregression[v+ii*numvertices];} //compute autoregressive term
-                        smblist[v] = beta0_basin+beta1_basin*telapsed_ar+autoregressionterm+noise_basin; //stochastic SMB value
+
+                        /*compute autoregressive term*/
+                        IssmDouble autoregressionterm=0.;
+                        for(int ii=0;ii<arorder;ii++){
+                                autoregressionterm += phi_basin[ii]*smbvaluesautoregression[v+ii*numvertices];
+                        }
+
+                        /*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
+                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
@@ -3664 +3687 @@
         /*Add input to element*/
         this->AddInput(SmbMassBalanceEnum,smblist,P1Enum);
+
         /*Cleanup*/
         xDelete<IssmDouble>(phi_basin);

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

@@ -3 +3 @@
  */
 
+#include <config.h>
 #include "./SurfaceMassBalancex.h"
 #include "../../shared/shared.h"
@@ -195 +196 @@
         /*Determine if this is a time step for the AR model*/
         bool isstepforar{false};
-   if((std::fmod(time,tstep_ar)<std::fmod((time-dt),tstep_ar)) || (time<=starttime+dt) || tstep_ar==dt){isstepforar = true;}
+
+        #ifndef _HAVE_AD_
+   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;

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

@@ -677 +677 @@
    Follows the Cholesky–Banachiewicz algorithm
    Lchol should point to an IssmDouble* of same dimensions as A*/
-   for(int ii{0};ii<(ndim*ndim);ii++){Lchol[ii]=0;}; //ensure zero-initialization
-   for(int ii{0};ii<ndim;ii++){
-                for(int jj{0};jj<=ii;jj++){
-                        double sum{0.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);}
-                }
-        }
+
+        /*ensure zero-initialization*/
+        for(int ii=0;ii<(ndim*ndim);ii++) Lchol[ii]=0;;
+
+        for(int ii=0;ii<ndim;ii++){
+                for(int jj{0};jj<=ii;jj++){
+                        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);
+                        }
+                }
+        }
 } /*}}}*/

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

@@ -20 +20 @@
 /*}}}*/
 
-void univariateNormal(IssmDouble* prand, IssmDouble mean, IssmDouble sdev) { /*{{{*/
+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<double> normdistri(mean,sdev); //Normal probability distribution
-        double tfunc;
+   std::normal_distribution<IssmPDouble> normdistri(mean,sdev); //Normal probability distribution
         *prand = normdistri(generator);
 } /*}}}*/
 void multivariateNormal(IssmDouble** prand, int dim, IssmDouble mean, IssmDouble* covariancematrix) { /*{{{*/
-   IssmDouble* sampleStandardNormal     = xNew<IssmDouble>(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 ii=0;ii<dim;ii++){univariateNormal(&(sampleStandardNormal[ii]),0.0,1.0);}
    CholeskyRealPositiveDefinite(Lchol,covariancematrix,dim);
-   for(int ii{0};ii<dim;ii++){ //matrix by vector multiplication
-      double sum{0.0};
-      for(int jj{0};jj<dim;jj++){sum += sampleStandardNormal[jj]*Lchol[ii*dim+jj];} //entry-by-entry multiplication along matrix row
-      sampleMultivariateNormal[ii] = mean+sum; //assign value
+   for(int ii=0;ii<dim;ii++){ //matrix by vector multiplication
+      /*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;
    }
+
+   /*Assign output pointer and cleanup*/
    *prand = sampleMultivariateNormal;
-   xDelete<IssmDouble>(sampleStandardNormal);
+   xDelete<IssmPDouble>(sampleStandardNormal);
    xDelete<IssmDouble>(Lchol);
 } /*}}}*/
 void multivariateNormal(IssmDouble** prand, int dim, IssmDouble* mean, IssmDouble* covariancematrix) { /*{{{*/
-        IssmDouble* sampleStandardNormal     = xNew<IssmDouble>(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 ii=0;ii<dim;ii++) univariateNormal(&(sampleStandardNormal[ii]),0.0,1.0);
+
         CholeskyRealPositiveDefinite(Lchol,covariancematrix,dim);
-        for(int ii{0};ii<dim;ii++){ //matrix by vector multiplication
-                double sum{0.0};
-      for(int jj{0};jj<dim;jj++){sum += sampleStandardNormal[jj]*Lchol[ii*dim+jj];} //entry-by-entry multiplication along matrix row
+
+        for(int ii=0;ii<dim;ii++){ //matrix by vector multiplication
+                IssmDouble sum = 0.;
+      for(int jj=0.;jj<dim;jj++){
+                        sum += sampleStandardNormal[jj]*Lchol[ii*dim+jj];
+                }
       sampleMultivariateNormal[ii] = mean[ii]+sum; //assign value
         }
         *prand = sampleMultivariateNormal;
-        xDelete<IssmDouble>(sampleStandardNormal);
+        xDelete<IssmPDouble>(sampleStandardNormal);
         xDelete<IssmDouble>(Lchol);
 } /*}}}*/