Changeset 20028

Timestamp:

01/30/16 13:45:48 (9 years ago)

Author:

Eric.Larour

Message:

CHG: faster implementation of the precomputation of Green functions for elastic case.

Location:

issm/trunk-jpl/src

Files:

• c/analyses/SealevelriseAnalysis.cpp (modified) (2 diffs)
• c/classes/Elements/Tria.cpp (modified) (14 diffs)
• c/classes/Params/DoubleVecParam.cpp (modified) (1 diff)
• c/classes/Params/DoubleVecParam.h (modified) (1 diff)
• c/shared/Enum/EnumDefinitions.h (modified) (1 diff)
• c/shared/Enum/EnumToStringx.cpp (modified) (1 diff)
• c/shared/Enum/StringToEnumx.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/SealevelriseGRigidEnum.m (deleted)
• m/enum/SealevelriseLegendreCoefficientsEnum.m (deleted)
• m/enum/SealevelriseLegendrePrecomputeEnum.m (deleted)
• m/enum/SealevelriseStoreGreenFunctionsEnum.m (deleted)

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

@@ -43 +43 @@
         IssmDouble* love_k=NULL;
         
-        bool        legendre_precompute=false;
-        IssmDouble* legendre_coefficients=NULL;
+        IssmDouble* G_elastic = NULL;
         int         M;
 
@@ -54 +53 @@
         parameters->AddObject(iomodel->CopyConstantObject(SealevelriseElasticEnum));
         parameters->AddObject(iomodel->CopyConstantObject(SealevelriseEustaticEnum));
-        parameters->AddObject(iomodel->CopyConstantObject(SealevelriseLegendrePrecomputeEnum));
-        parameters->AddObject(iomodel->CopyConstantObject(SealevelriseStoreGreenFunctionsEnum));
 
         /*love numbers: */
         iomodel->FetchData(&love_h,&nl,NULL,SealevelriseLoveHEnum);
         iomodel->FetchData(&love_k,&nl,NULL,SealevelriseLoveKEnum);
-        parameters->AddObject(new DoubleVecParam(SealevelriseLoveHEnum,love_h,nl));
-        parameters->AddObject(new DoubleVecParam(SealevelriseLoveKEnum,love_k,nl));
 
-        /*legendre coefficients: */
-        iomodel->Constant(&legendre_precompute,SealevelriseLegendrePrecomputeEnum);
-        if(legendre_precompute){
+        /*compute elastic green function for a range of angles*/
+        const IssmDouble degacc=.01; 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 <=
+                IssmDouble alpha,x;
+                alpha= reCast<IssmDouble>(i)*degacc * PI / 180.0;
+                
+                G_elastic[i]= (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]);
 
-                /*compute values at which to evaluate the legendre polynomical values:*/
-                const IssmDouble degacc=.01; M=reCast<int>(180/degacc+1);
-                IssmDouble* x=xNew<IssmDouble>(M);
-                for(int i=0;i<=M;i++){ //watch out the <=
-                        x[i]=-cos( (reCast<IssmDouble>(i)*degacc) * PI / 180.0);
+                        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[i] += deltalove*Pn;
                 }
-                /*compute legendre coefficient matrix:*/
-                legendre_coefficients=p_polynomial_value(M,nl-1,x);
-                parameters->AddObject(new DoubleMatParam(SealevelriseLegendreCoefficientsEnum,legendre_coefficients,M,nl));
+        }
+        parameters->AddObject(new DoubleVecParam(SealevelriseGElasticEnum,G_elastic,M));
 
-                /*free ressources:*/
-                xDelete<IssmDouble>(x);
-        }
+        /*free ressources:*/
+        xDelete<IssmDouble>(G_elastic);
 
         /*free ressources: */
         xDelete<IssmDouble>(love_h);
         xDelete<IssmDouble>(love_k);
-        xDelete<IssmDouble>(legendre_coefficients);
+        xDelete<IssmDouble>(G_elastic);
 
 }/*}}}*/

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

@@ -3515 +3515 @@
         IssmDouble lati,longi,ri;
 
-        /*love numbers:*/
-        IssmDouble* love_h=NULL;
-        IssmDouble* love_k=NULL; 
-        IssmDouble* deltalove=NULL;
-        int nl;
-
-        /*legendre coefficients:*/
-        bool        legendre_precompute=false;
-        IssmDouble* legendre_coefficients=NULL;
+        /*elastic green function:*/
+        IssmDouble* G_elastic_precomputed=NULL;
         int         M;
 
@@ -3550 +3543 @@
 
         /*recover love numbers and computational flags: */
-        this->parameters->FindParam(&love_h,&nl,SealevelriseLoveHEnum);
-        this->parameters->FindParam(&love_k,&nl,SealevelriseLoveKEnum);
         this->parameters->FindParam(&computerigid,SealevelriseRigidEnum);
         this->parameters->FindParam(&computeelastic,SealevelriseElasticEnum);
         this->parameters->FindParam(&computeeustatic,SealevelriseEustaticEnum);
 
-        /*recover legendre coefficients if they have been precomputed:*/
-        this->parameters->FindParam(&legendre_precompute,SealevelriseLegendrePrecomputeEnum);
-        if(legendre_precompute){
-                DoubleMatParam* parameter = static_cast<DoubleMatParam*>(this->parameters->FindParamObject(SealevelriseLegendreCoefficientsEnum)); _assert_(parameter);
-                parameter->GetParameterValueByPointer(&legendre_coefficients,&M,NULL);
+        /*recover elastic green function:*/
+        if(computeelastic){
+                DoubleVecParam* parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelriseGElasticEnum)); 
+                _assert_(parameter);
+                parameter->GetParameterValueByPointer(&G_elastic_precomputed,&M);
         }
 
@@ -3614 +3605 @@
         if(computeeustatic) eustatic += rho_ice*area*I/(oceanarea*rho_water); 
 
-        /*Speed up of love number comutation: */
-        if(computeelastic){
-                deltalove=xNew<IssmDouble>(nl);
-                for (int n=0;n<nl;n++) deltalove[n] = (love_k[n]-love_k[nl-1]-love_h[n]+love_h[nl-1]);
-        }
-
         if(computeelastic | computerigid){
                 int* indices=xNew<int>(gsize);
@@ -3629 +3614 @@
                         IssmDouble G_rigid=0;  //do not remove =0!
                         IssmDouble G_elastic=0;  //do not remove =0!
-                        IssmDouble Pn1,Pn2,Pn; //first two legendre polynomials
-                        IssmDouble cosalpha;
                         IssmDouble delPhi,delLambda;
 
@@ -3638 +3621 @@
                     delPhi=fabs(lati-late); delLambda=fabs(longi-longe);
                         alpha=2.*asin(sqrt(pow(sin(delPhi/2),2.0)+cos(lati)*cos(late)*pow(sin(delLambda/2),2)));
-                        cosalpha=cos(alpha); //compute this to speed up the elastic component.
 
                         //Rigid earth gravitational perturbation:
@@ -3645 +3627 @@
                         //Elastic component  (from Eq 17 in Adhikari et al, GMD 2015)
                         if(computeelastic){
-                                G_elastic = (love_k[nl-1]-love_h[nl-1])/2.0/sin(alpha/2.0);
-                                if(legendre_precompute){
-                                        int index=reCast<int,IssmDouble>((M-1)*(1-alpha/PI));
-                                        for(int n=0;n<nl;n++) G_elastic += deltalove[n]*legendre_coefficients[index*nl+n];
-                                }
-                                else{
-                                        for(int n=0;n<nl;n++){
-                                                if(n==0)Pn=1;
-                                                else if(n==1)Pn=cosalpha;
-                                                else Pn= ( (2*n-1)*cosalpha*Pn1 - (n-1)*Pn2 ) /n;
-                                                Pn2=Pn1; Pn1=Pn;
-                                                G_elastic += deltalove[n]*Pn;
-                                        }
-                                }
+                                int index=reCast<int,IssmDouble>(alpha/PI*(M-1));
+                                G_elastic += G_elastic_precomputed[index];
                         }
 
@@ -3673 +3643 @@
         /*Assign output pointer:*/
         *peustatic=eustatic;
-
-        /*Free ressources:*/
-        if(computeelastic){
-                xDelete<IssmDouble>(love_h);
-                xDelete<IssmDouble>(love_k);
-                xDelete<IssmDouble>(deltalove);
-        }
         return;
 }
@@ -3696 +3659 @@
         IssmDouble maxlong=-20;
 
-        /*love numbers:*/
-        IssmDouble* love_h=NULL;
-        IssmDouble* love_k=NULL; 
-        IssmDouble* deltalove=NULL;
-        int nl;
-
-        /*legendre coefficients:*/
-        bool        legendre_precompute=false;
-        IssmDouble* legendre_coefficients=NULL;
+        /*precomputed elastic green functions:*/
+        IssmDouble* G_elastic_precomputed = NULL;
         int         M;
+        
+        /*computation of Green functions:*/
+        IssmDouble* G_elastic= NULL;
+        IssmDouble* G_rigid= NULL;
 
         /*optimization:*/
@@ -3718 +3678 @@
         bool computeeustatic = true;
 
-        /*G_rigid and G_elasti variables, to speed up the computations: */
-        DoubleArrayInput*      G_rigid_input=NULL;
-        IssmDouble* G_rigid = NULL; int G_rigid_M;
-        bool        compute_G_rigid=false;
-        DoubleArrayInput*      G_elastic_input=NULL;
-        IssmDouble* G_elastic = NULL; int G_elastic_M;
-        bool        compute_G_elastic=false;
-
         /*early return if we are not on the ocean:*/
         if (!IsWaterInElement())return;
@@ -3733 +3685 @@
         this->parameters->FindParam(&computeelastic,SealevelriseElasticEnum);
         this->parameters->FindParam(&computeeustatic,SealevelriseEustaticEnum);
-        this->parameters->FindParam(&store_green_functions,SealevelriseStoreGreenFunctionsEnum);
         
         /*early return if rigid or elastic not requested:*/
         if(!computerigid && !computeelastic) return;
 
-        /*Try and recover, if it exists, G_rigid and G_elastic:*/
-        G_rigid_input=(DoubleArrayInput*)this->inputs->GetInput(SealevelriseGRigidEnum);
-        if(G_rigid_input){
-                compute_G_rigid=false;
-                G_rigid_input->GetValues(&G_rigid,&G_rigid_M);
-        }
-        else if(computerigid)compute_G_rigid=true;
-
-        G_elastic_input=(DoubleArrayInput*)this->inputs->GetInput(SealevelriseGElasticEnum);
-        if(G_elastic_input){
-                compute_G_elastic=false;
-                G_elastic_input->GetValues(&G_elastic,&G_elastic_M);
-        }
-        else if(computeelastic)compute_G_elastic=true;
-        
         /*recover material parameters: */
         rho_water=matpar->GetMaterialParameter(MaterialsRhoSeawaterEnum);
@@ -3801 +3737 @@
         longe=longe/180*PI;
         /*}}}*/
-        if(compute_G_elastic){
-                /*recover love numbers and legendre_coefficients:{{{*/
-                this->parameters->FindParam(&love_h,&nl,SealevelriseLoveHEnum);
-                this->parameters->FindParam(&love_k,&nl,SealevelriseLoveKEnum);
-
-                /*recover legendre coefficients if they have been precomputed:*/
-                this->parameters->FindParam(&legendre_precompute,SealevelriseLegendrePrecomputeEnum);
-                if(legendre_precompute){
-                        DoubleMatParam* parameter = static_cast<DoubleMatParam*>(this->parameters->FindParamObject(SealevelriseLegendreCoefficientsEnum)); _assert_(parameter);
-                        parameter->GetParameterValueByPointer(&legendre_coefficients,&M,NULL);
-                }
-
-                /*Speed up of love number comutation for elastic mode: */
-                deltalove=xNew<IssmDouble>(nl);
-                for (int n=0;n<nl;n++) deltalove[n] = (love_k[n]-love_k[nl-1]-love_h[n]+love_h[nl-1]);
-                //}}}
+        
+        if(computeelastic){
+        
+                /*recover elastic green function:*/
+                DoubleVecParam* parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelriseGElasticEnum)); _assert_(parameter);
+                parameter->GetParameterValueByPointer(&G_elastic_precomputed,&M);
 
                 /*initialize G_elastic:*/
                 G_elastic=xNewZeroInit<IssmDouble>(gsize);
         }
-        if(compute_G_rigid){
-                /*Initialize G_rigid and G_elastic:*/
-                G_rigid=xNewZeroInit<IssmDouble>(gsize);
-        }
+        if(computerigid) G_rigid=xNewZeroInit<IssmDouble>(gsize);
 
         int* indices=xNew<int>(gsize);
@@ -3832 +3755 @@
 
                 indices[i]=i; 
-                if(compute_G_elastic | compute_G_rigid){
+                if(computeelastic | computerigid){
+        
                         IssmDouble alpha;
-                        IssmDouble Pn1,Pn2,Pn; //first two legendre polynomials
-                        IssmDouble cosalpha;
                         IssmDouble delPhi,delLambda;
 
@@ -3843 +3765 @@
                         delPhi=fabs(lati-late); delLambda=fabs(longi-longe);
                         alpha=2.*asin(sqrt(pow(sin(delPhi/2),2.0)+cos(lati)*cos(late)*pow(sin(delLambda/2),2)));
-                        cosalpha=cos(alpha); //compute this to speed up the elastic component.
 
                         //Rigid earth gravitational perturbation:
-                        if(compute_G_rigid)G_rigid[i]=1.0/2.0/sin(alpha/2.0);
+                        if(computerigid)G_rigid[i]=1.0/2.0/sin(alpha/2.0);
 
                         //Elastic component  (from Eq 17 in Adhikari et al, GMD 2015)
-                        if(compute_G_elastic){
-                                G_elastic[i] = (love_k[nl-1]-love_h[nl-1])/2.0/sin(alpha/2.0);
-                                if(legendre_precompute){
-                                        int index=reCast<int,IssmDouble>((M-1)*(1-alpha/PI));
-                                        for(int n=0;n<nl;n++) G_elastic[i] += deltalove[n]*legendre_coefficients[index*nl+n];
-                                }
-                                else{
-                                        for(int n=0;n<nl;n++){
-                                                if(n==0)Pn=1;
-                                                else if(n==1)Pn=cosalpha;
-                                                else Pn= ( (2*n-1)*cosalpha*Pn1 - (n-1)*Pn2 ) /n;
-                                                Pn2=Pn1; Pn1=Pn;
-                                                G_elastic[i] += deltalove[n]*Pn;
-                                        }
-                                }
+                        if(computeelastic){
+                                int index=reCast<int,IssmDouble>(alpha/PI*(M-1));
+                                G_elastic[i] += G_elastic_precomputed[index];
                         }
                 }
@@ -3878 +3787 @@
         xDelete<int>(indices);
 
-        /*Save G_rigid and G_elastic if computed:*/
-        if(store_green_functions){
-                if(G_rigid)this->inputs->AddInput(new DoubleArrayInput(SealevelriseGRigidEnum,G_rigid,gsize));
-                if(G_elastic)this->inputs->AddInput(new DoubleArrayInput(SealevelriseGElasticEnum,G_elastic,gsize));
-        }
-
         /*Free ressources:*/
-        if(compute_G_elastic){
-                xDelete<IssmDouble>(love_h);
-                xDelete<IssmDouble>(love_k);
-                xDelete<IssmDouble>(deltalove);
-                xDelete<IssmDouble>(G_elastic);
-        }
-        if(compute_G_rigid) xDelete<IssmDouble>(G_rigid);
+        if(computeelastic) xDelete<IssmDouble>(G_elastic);
+        if(computerigid) xDelete<IssmDouble>(G_rigid);
 
         return;

issm/trunk-jpl/src/c/classes/Params/DoubleVecParam.cpp

@@ -121 +121 @@
 }
 /*}}}*/
+
+/*DoubleVecParam specific routines:*/
+void  DoubleVecParam::GetParameterValueByPointer(IssmDouble** pIssmDoublearray,int* pM){/*{{{*/
+        
+        /*Assign output pointers:*/
+        if(pM) *pM=M;
+        *pIssmDoublearray=values;
+}
+/*}}}*/

issm/trunk-jpl/src/c/classes/Params/DoubleVecParam.h

@@ -71 +71 @@
                 void  SetValue(IssmDouble** array, int M, int* mdim_array, int* ndim_array){_error_("Param "<< EnumToStringx(enum_type) << " cannot hold an array of matrices");}
                 /*}}}*/
+                /*DoubleVecParam specific routines:{{{*/
+                void  GetParameterValueByPointer(IssmDouble** pIssmDoublearray,int* pM);
+                /*}}}*/
+
 };
 #endif  /* _DOUBLEVECPARAM_H */

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

@@ -1047 +1047 @@
         SealevelriseElasticEnum,
         SealevelriseEustaticEnum,
-        SealevelriseLegendrePrecomputeEnum,
-        SealevelriseLegendreCoefficientsEnum,
-        SealevelriseGRigidEnum,
         SealevelriseGElasticEnum,
-        SealevelriseStoreGreenFunctionsEnum,
         /*}}}*/
         MaximumNumberOfDefinitionsEnum

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

@@ -1001 +1001 @@
                 case SealevelriseElasticEnum : return "SealevelriseElastic";
                 case SealevelriseEustaticEnum : return "SealevelriseEustatic";
-                case SealevelriseLegendrePrecomputeEnum : return "SealevelriseLegendrePrecompute";
-                case SealevelriseLegendreCoefficientsEnum : return "SealevelriseLegendreCoefficients";
-                case SealevelriseGRigidEnum : return "SealevelriseGRigid";
                 case SealevelriseGElasticEnum : return "SealevelriseGElastic";
-                case SealevelriseStoreGreenFunctionsEnum : return "SealevelriseStoreGreenFunctions";
                 case MaximumNumberOfDefinitionsEnum : return "MaximumNumberOfDefinitions";
                 default : return "unknown";

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

@@ -1025 +1025 @@
               else if (strcmp(name,"SealevelriseElastic")==0) return SealevelriseElasticEnum;
               else if (strcmp(name,"SealevelriseEustatic")==0) return SealevelriseEustaticEnum;
-              else if (strcmp(name,"SealevelriseLegendrePrecompute")==0) return SealevelriseLegendrePrecomputeEnum;
-              else if (strcmp(name,"SealevelriseLegendreCoefficients")==0) return SealevelriseLegendreCoefficientsEnum;
-              else if (strcmp(name,"SealevelriseGRigid")==0) return SealevelriseGRigidEnum;
               else if (strcmp(name,"SealevelriseGElastic")==0) return SealevelriseGElasticEnum;
-              else if (strcmp(name,"SealevelriseStoreGreenFunctions")==0) return SealevelriseStoreGreenFunctionsEnum;
               else if (strcmp(name,"MaximumNumberOfDefinitions")==0) return MaximumNumberOfDefinitionsEnum;
          else stage=10;

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

@@ -15 +15 @@
                 elastic         = 0;
                 eustatic         = 0;
-                legendre_precompute = 0;
-                store_green_functions = 0;
         end
         methods
@@ -41 +39 @@
                 self.eustatic=1;
                 
-                %legendre coefficients: 
-                self.legendre_precompute = 1;
-
-                %optimization: 
-                self.store_green_functions=1;
-
                 end % }}}
                 function md = checkconsistency(self,md,solution,analyses) % {{{
@@ -57 +49 @@
                         md = checkfield(md,'fieldname','slr.abstol','size',[1 1]);
                         md = checkfield(md,'fieldname','slr.maxiter','size',[1 1],'>=',1);
-                        md = checkfield(md,'fieldname','slr.store_green_functions','value',[0 1]);
 
                         %check that love numbers are provided at the same level of accuracy: 
@@ -77 +68 @@
                         fielddisplay(self,'elastic','elastic earth graviational potential perturbation');
                         fielddisplay(self,'eustatic','eustatic sea level rise');
-                        fielddisplay(self,'store_green_functions','store green functions (default 1) to speed up solutoin (though memory intense)');
-                        fielddisplay(self,'legendre_precompute','precompute legendre coefficients? (default is 0)');
 
                 end % }}}
@@ -91 +80 @@
                         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','store_green_functions','format','Boolean');
-                        WriteData(fid,'object',self,'class','sealevelrise','fieldname','legendre_precompute','format','Boolean');
                 end % }}}
                 function savemodeljs(self,fid,modelname) % {{{
@@ -102 +89 @@
                         writejs1Darray(fid,[modelname '.srl.love_h'],self.love_h);
                         writejs1Darray(fid,[modelname '.srl.love_k'],self.love_k);
-                        writejsdouble(fid,[modelname '.slr.store_green_functions'],self.store_green_functions);
                         writejsdouble(fid,[modelname '.slr.rigid'],self.rigid);
                         writejsdouble(fid,[modelname '.slr.eustatic'],self.eustatic);
-                        writejsdouble(fid,[modelname '.slr.legendre_precompute'],self.legendre_precompute);
 
                 end % }}}

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

@@ -986 +986 @@
 function SealevelriseElasticEnum(){ return 982;}
 function SealevelriseEustaticEnum(){ return 983;}
-function SealevelriseLegendrePrecomputeEnum(){ return 984;}
-function SealevelriseLegendreCoefficientsEnum(){ return 985;}
-function SealevelriseGRigidEnum(){ return 986;}
-function SealevelriseGElasticEnum(){ return 987;}
-function SealevelriseStoreGreenFunctionsEnum(){ return 988;}
-function MaximumNumberOfDefinitionsEnum(){ return 989;}
+function SealevelriseGElasticEnum(){ return 984;}
+function MaximumNumberOfDefinitionsEnum(){ return 985;}

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

@@ -993 +993 @@
 def SealevelriseElasticEnum(): return StringToEnum("SealevelriseElastic")[0]
 def SealevelriseEustaticEnum(): return StringToEnum("SealevelriseEustatic")[0]
-def SealevelriseLegendrePrecomputeEnum(): return StringToEnum("SealevelriseLegendrePrecompute")[0]
-def SealevelriseLegendreCoefficientsEnum(): return StringToEnum("SealevelriseLegendreCoefficients")[0]
-def SealevelriseGRigidEnum(): return StringToEnum("SealevelriseGRigid")[0]
 def SealevelriseGElasticEnum(): return StringToEnum("SealevelriseGElastic")[0]
-def SealevelriseStoreGreenFunctionsEnum(): return StringToEnum("SealevelriseStoreGreenFunctions")[0]
 def MaximumNumberOfDefinitionsEnum(): return StringToEnum("MaximumNumberOfDefinitions")[0]