Changeset 19613

Timestamp:

10/01/15 10:45:44 (10 years ago)

Author:

Eric.Larour

Message:

CHG: applying patches for the AD mode to compile.

Location:

issm/trunk-jpl/src/c

Files:

• classes/Elements/Element.cpp (modified) (3 diffs)
• modules/SurfaceMassBalancex/Gembx.cpp (modified) (13 diffs)

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

@@ -2375 +2375 @@
                 /*Calculate total system mass:*/
         sumMass=0; for(int i=0;i<m;i++) sumMass += dz[i]*d[i];
+
+        #ifndef _HAVE_ADOLC_ //we want to avoid the round operation at all cost. Not differentiable.
         dMass = sumMass + sumR + sumW - sumP - sumEC - initMass - sumMassAdd;
-        dMass = round(dMass * 100.0)/100.0;
+                dMass = round(dMass * 100.0)/100.0;
 
                 /*Check mass conservation:*/
         if (dMass != 0.0) _printf_("total system mass not conserved in MB function");
+                #endif
                 
                 /*Check bottom grid cell T is unchanged:*/
         if (T[m-1]!=T_bottom) _printf_("T(end)~=T_bottom" << "\n");
                 
+                /*Free ressources: */
+                xDelete<IssmDouble>(swf);
+
+                /*increase counter:*/
                 count++;
 
@@ -2399 +2406 @@
         this->AddInput(new DoubleArrayInput(SmbWEnum,W,m));
         this->AddInput(new DoubleArrayInput(SmbAEnum,a,m));
-        this->AddInput(new DoubleArrayInput(SmbSwfEnum,swf,m));
         this->AddInput(new DoubleInput(SmbMassBalanceEnum,(sumP + sumEC -sumR)/rho_water/dt));
         this->AddInput(new DoubleInput(SmbRunoffEnum,sumR/rho_water/dt));
@@ -2416 +2422 @@
         xDelete<IssmDouble>(a);
         xDelete<IssmDouble>(T);
-        xDelete<IssmDouble>(swf);
         delete gauss;
         /*}}}*/

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

@@ -43 +43 @@
         //into specified top structure depth (zTop). Also make sure top grid cell
         //structure length (dzTop) is greater than 5 cm
-        if (dgpTop != round(dgpTop)){
+        #ifndef _HAVE_ADOLC_  //avoid the round operation check!
+        if (dgpTop != round(dgpTop)){ 
                 _error_("top grid cell structure length does not go evenly into specified top structure depth, adjust dzTop or zTop");
         }
-        else if(dzTop < 0.05){
+        #endif
+        if(dzTop < 0.05){
                 _printf_("initial top grid cell length (dzTop) is < 0.05 m");
         }
-        gpTop=round(dgpTop);
+        gpTop=reCast<int,IssmDouble>(dgpTop);
 
         //initialize top grid depth vector
@@ -212 +214 @@
 
         // take absolute value of temperature gradient
-        for(int i=0;i<m;i++)dT[i]=abs(dT[i]);
+        for(int i=0;i<m;i++)dT[i]=fabs(dT[i]);
         
         /*Snow metamorphism. Depends on value of snow dendricity and wetness of the snowpack: */
@@ -418 +420 @@
                 for(int i=0;i<m;i++)if(W[i]>0)t0[i]=t0wet; // wet snow timescale
                 for(int i=0;i<m;i++)T[i]=TK[i] - 273.15; // change T from K to degC
-                for(int i=0;i<m;i++) t0warm[i]= abs(T[i]) * K + t0dry; //// 'warm' snow timescale
+                for(int i=0;i<m;i++) t0warm[i]= fabs(T[i]) * K + t0dry; //// 'warm' snow timescale
         for(int i=0;i<m;i++)if(W[i]==0.0 && T[i]>=-10)t0[i]= t0warm[i];
         for(int i=0;i<m;i++)if(T[i]<-10) t0[i] =  10 * K + t0dry; // 'cold' snow timescale
@@ -450 +452 @@
         if (a[0] > 1) _printf_("albedo > 1.0\n");
         else if (a[0] < 0) _printf_("albedo is negative\n");
-        else if (isnan(a[0])) _error_("albedo == NAN\n");
+        else if (xIsNan(a[0])) _error_("albedo == NAN\n");
 }  /*}}}*/
 void thermo(IssmDouble* pEC, IssmDouble* T, IssmDouble* dz, IssmDouble* d, IssmDouble* swf, IssmDouble dlwrf, IssmDouble Ta, IssmDouble V, IssmDouble eAir, IssmDouble pAir, IssmDouble Ws, IssmDouble dt0, int m, IssmDouble Vz, IssmDouble Tz,int sid) { /*{{{*/
@@ -692 +694 @@
 
         /* CALCULATE ENERGY SOURCES AND DIFFUSION FOR EVERY TIME STEP [dt]*/
-        for (int i=1;i<=dt0;i+=dt){
+        for (IssmDouble i=1;i<=dt0;i+=dt){
 
                 // PART OF ENERGY CONSERVATION CHECK
@@ -964 +966 @@
                         xDelete<IssmDouble>(Qs1);
                         xDelete<IssmDouble>(Qs2);
+                        
+                        
                 }
                 else{  //function of grid cell density
@@ -1158 +1162 @@
 
                 mass_diff = mass - massinit - P;
+                
+                #ifndef _HAVE_ADOLC_  //avoid round operation. only check in forward mode.
                 mass_diff = round(mass_diff * 100)/100;
-
                 if (mass_diff > 0) _error_("mass not conserved in accumulation function");
+                #endif
 
         }
@@ -1345 +1351 @@
                 X = 0;
                 for(int i=n-1;i>=0;i--){
-                        if(M[i]>0 || exsW[i]){
+                        if(M[i]>0 || reCast<int,IssmDouble>(exsW[i])){
                                 X=i;
                                 break;
@@ -1451 +1457 @@
                 celldelete(&R,n,D,D_size);
                 n=D_size;
+                xDelete<int>(D);
         
                 // calculate new grid lengths
@@ -1457 +1464 @@
                 //calculate Rsum: 
                 Rsum=cellsum(R,n);
-        
+
+                /*Free ressources:*/
+                xDelete<IssmDouble>(F);
+                xDelete<IssmDouble>(R);
         }
 
@@ -1492 +1502 @@
 
         // delete combined cells
-        xDelete<int>(D); D_size=0; for(int i=0;i<n;i++)if(m[i]!=99999)D_size++; D=xNew<int>(D_size); 
+        D_size=0; for(int i=0;i<n;i++)if(m[i]!=99999)D_size++; D=xNew<int>(D_size); 
         D_size=0; for(int i=0;i<n;i++)if(m[i]!=99999){ D[D_size] = i; D_size++;}
 
@@ -1507 +1517 @@
         celldelete(&EW,n,D,D_size);
         n=D_size;
+        xDelete<int>(D);
 
         // check if any of the top 10 cell depths are too large
@@ -1586 +1597 @@
         sumE1 = cellsum(EI,n) + cellsum(EW,n);
 
+        /*checks: */
+        for(int i=0;i<n;i++) if (W[i]<0) _error_("negative pore water generated in melt equations\n");
+        
+        /*only in forward mode! avoid round in AD mode as it is not differentiable: */
+        #ifndef _HAVE_ADOLC_
         dm = round(mSum0 - mSum1 + mAdd);
         dE = round(sumE0 - sumE1 - sumER +  addE);
-        
-        for(int i=0;i<n;i++) if (W[i]<0) _error_("negative pore water generated in melt equations\n");
-
         if (dm !=0  || dE !=0) _error_("mass or energy are not conserved in melt equations\n" 
                         << "dm: " << dm << " dE: " << dE << "\n");
-
+        #endif
         
         /*Free ressources:*/