Changeset 20033

Timestamp:

01/30/16 19:12:44 (9 years ago)

Author:

Eric.Larour

Message:

CHG: parallelized computation of G_elastic in the Sealevelrise UpdateParameters analysis phase.

Location:

issm/trunk-jpl/src

Files:

• c/analyses/SealevelriseAnalysis.cpp (modified) (4 diffs)
• c/shared/Enum/EnumDefinitions.h (modified) (1 diff)
• c/shared/Enum/EnumToStringx.cpp (modified) (1 diff)
• c/shared/Enum/StringToEnumx.cpp (modified) (1 diff)
• c/shared/Numerics/legendre.cpp (modified) (1 diff)
• m/classes/slr.m (modified) (6 diffs)
• m/enum/EnumDefinitions.js (modified) (1 diff)
• m/enum/EnumDefinitions.py (modified) (1 diff)
• m/enum/SealevelriseDegaccEnum.m (added)

issm/trunk-jpl/src/c/analyses/SealevelriseAnalysis.cpp

@@ -45 +45 @@
         bool elastic=false;
         IssmDouble* G_elastic = NULL;
-        int         M;
+        IssmDouble* G_elastic_local = NULL;
+        int         M,m,lower_row,upper_row;
+        IssmDouble  degacc=.01;
 
         /*some constant parameters: */
@@ -63 +65 @@
 
                 /*compute elastic green function for a range of angles*/
-                IssmDouble degacc=.01; M=reCast<int,IssmDouble>(180./degacc+1.);
+                iomodel->FetchData(&degacc,SealevelriseDegaccEnum);
+                M=reCast<int>(180/degacc+1);
                 G_elastic=xNew<IssmDouble>(M);
 
                 /*compute combined legendre + love number (elastic green function:*/
-                for(int i=0;i<M;i++){ //watch out the <=
+                m=DetermineLocalSize(M,IssmComm::GetComm());
+                GetOwnershipBoundariesFromRange(&lower_row,&upper_row,m,IssmComm::GetComm());
+                G_elastic_local=xNew<IssmDouble>(m);
+
+                for(int i=lower_row;i<upper_row;i++){
+                        IssmDouble alpha,x;
+                        alpha= reCast<IssmDouble>(i)*degacc * PI / 180.0;
+
+                        G_elastic_local[i-lower_row]= (love_k[nl-1]-love_h[nl-1])/2.0/sin(alpha/2.0);
+                        for (int n=0;n<nl;n++) {
+                                IssmDouble Pn,Pn1,Pn2;
+                                IssmDouble deltalove;
+
+                                deltalove = (love_k[n]-love_k[nl-1]-love_h[n]+love_h[nl-1]);
+
+                                if(n==0)Pn=1;
+                                else if(n==1)Pn=cos(alpha);
+                                else Pn= ( (2*n-1)*cos(alpha)*Pn1 - (n-1)*Pn2 ) /n;
+                                Pn2=Pn1; Pn1=Pn;
+
+                                G_elastic_local[i-lower_row] += deltalove*Pn;
+                        }
+                }
+
+                /*merge G_elastic_local into G_elastic:{{{*/
+                int* recvcounts=xNew<int>(IssmComm::GetSize());
+                int* displs=xNew<int>(IssmComm::GetSize());
+
+                //recvcounts:
+                ISSM_MPI_Allgather(&m,1,ISSM_MPI_INT,recvcounts,1,ISSM_MPI_INT,IssmComm::GetComm());
+
+                /*displs: */
+                ISSM_MPI_Allgather(&lower_row,1,ISSM_MPI_INT,displs,1,ISSM_MPI_INT,IssmComm::GetComm());
+
+                /*All gather:*/
+                ISSM_MPI_Allgatherv(G_elastic_local, m, ISSM_MPI_DOUBLE, G_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
+                /*free ressources: */
+                xDelete<int>(recvcounts);
+                xDelete<int>(displs);
+
+                /*}}}*/
+
+                /*Old code: {{{*/
+                /*for(int i=0;i<M;i++){ 
                         IssmDouble alpha,x;
                         alpha= reCast<IssmDouble>(i)*degacc * PI / 180.0;
@@ -85 +131 @@
                                 G_elastic[i] += deltalove*Pn;
                         }
-                }
+                }*/ 
+                /*}}}*/
 
                 /*Avoid singularity at 0: */
@@ -91 +138 @@
                 parameters->AddObject(new DoubleVecParam(SealevelriseGElasticEnum,G_elastic,M));
 
-                /*free ressources:*/
-                xDelete<IssmDouble>(G_elastic);
-
                 /*free ressources: */
                 xDelete<IssmDouble>(love_h);
                 xDelete<IssmDouble>(love_k);
                 xDelete<IssmDouble>(G_elastic);
+                xDelete<IssmDouble>(G_elastic_local);
         }
 

issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h

@@ -1048 +1048 @@
         SealevelriseEustaticEnum,
         SealevelriseGElasticEnum,
+        SealevelriseDegaccEnum,
         /*}}}*/
         MaximumNumberOfDefinitionsEnum

issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp

@@ -1002 +1002 @@
                 case SealevelriseEustaticEnum : return "SealevelriseEustatic";
                 case SealevelriseGElasticEnum : return "SealevelriseGElastic";
+                case SealevelriseDegaccEnum : return "SealevelriseDegacc";
                 case MaximumNumberOfDefinitionsEnum : return "MaximumNumberOfDefinitions";
                 default : return "unknown";

issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp

@@ -1026 +1026 @@
               else if (strcmp(name,"SealevelriseEustatic")==0) return SealevelriseEustaticEnum;
               else if (strcmp(name,"SealevelriseGElastic")==0) return SealevelriseGElasticEnum;
+              else if (strcmp(name,"SealevelriseDegacc")==0) return SealevelriseDegaccEnum;
               else if (strcmp(name,"MaximumNumberOfDefinitions")==0) return MaximumNumberOfDefinitionsEnum;
          else stage=10;

issm/trunk-jpl/src/c/shared/Numerics/legendre.cpp

@@ -140 +140 @@
 #include "./recast.h"
 
-IssmDouble* p_polynomial_value ( int m, int n, IssmDouble* x){
+IssmDouble *p_polynomial_value ( int m, int n, IssmDouble* x){
 
         /******************************************************************************{{{/

issm/trunk-jpl/src/m/classes/slr.m

@@ -15 +15 @@
                 elastic         = 0;
                 eustatic         = 0;
+                degacc         = 0;
         end
         methods
@@ -38 +39 @@
                 self.elastic=1;
                 self.eustatic=1;
+
+                %numerical discretization accuracy
+                self.degacc=.01;
                 
                 end % }}}
@@ -49 +53 @@
                         md = checkfield(md,'fieldname','slr.abstol','size',[1 1]);
                         md = checkfield(md,'fieldname','slr.maxiter','size',[1 1],'>=',1);
+                        md = checkfield(md,'fieldname','slr.degacc','size',[1 1],'>=',1e-10);
 
                         %check that love numbers are provided at the same level of accuracy: 
@@ -68 +73 @@
                         fielddisplay(self,'elastic','elastic earth graviational potential perturbation');
                         fielddisplay(self,'eustatic','eustatic sea level rise');
+                        fielddisplay(self,'degacc','accuracy (default .01 deg) for numerical discretization of the Green''s functions');
 
                 end % }}}
@@ -80 +86 @@
                         WriteData(fid,'object',self,'class','sealevelrise','fieldname','elastic','format','Boolean');
                         WriteData(fid,'object',self,'class','sealevelrise','fieldname','eustatic','format','Boolean');
+                        WriteData(fid,'object',self,'class','sealevelrise','fieldname','degacc','format','Double');
                 end % }}}
                 function savemodeljs(self,fid,modelname) % {{{
@@ -91 +98 @@
                         writejsdouble(fid,[modelname '.slr.rigid'],self.rigid);
                         writejsdouble(fid,[modelname '.slr.eustatic'],self.eustatic);
+                        writejsdouble(fid,[modelname '.slr.degacc'],self.degacc);
 
                 end % }}}

issm/trunk-jpl/src/m/enum/EnumDefinitions.js

@@ -987 +987 @@
 function SealevelriseEustaticEnum(){ return 983;}
 function SealevelriseGElasticEnum(){ return 984;}
-function MaximumNumberOfDefinitionsEnum(){ return 985;}
+function SealevelriseDegaccEnum(){ return 985;}
+function MaximumNumberOfDefinitionsEnum(){ return 986;}

issm/trunk-jpl/src/m/enum/EnumDefinitions.py

@@ -994 +994 @@
 def SealevelriseEustaticEnum(): return StringToEnum("SealevelriseEustatic")[0]
 def SealevelriseGElasticEnum(): return StringToEnum("SealevelriseGElastic")[0]
+def SealevelriseDegaccEnum(): return StringToEnum("SealevelriseDegacc")[0]
 def MaximumNumberOfDefinitionsEnum(): return StringToEnum("MaximumNumberOfDefinitions")[0]