Changeset 19986

Timestamp:

01/22/16 10:47:46 (9 years ago)

Author:

Eric.Larour

Message:

CHG: added precomputation of legendre polynomials.

Location:

issm/trunk-jpl/src

Files:

• c/analyses/SealevelriseAnalysis.cpp (modified) (2 diffs)
• c/classes/Elements/Tria.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)
• m/classes/slr.m (modified) (6 diffs)
• m/enum/EnumDefinitions.js (modified) (1 diff)
• m/enum/EnumDefinitions.py (modified) (1 diff)

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

@@ -38 +38 @@
 void SealevelriseAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
 
-        int nl;
+        int         nl;
         IssmDouble* love_h=NULL;
         IssmDouble* love_k=NULL;
+        
+        bool        legendre_precompute=false;
+        IssmDouble* legendre_coefficients=NULL;
+        int         M;
 
+        /*some constant parameters: */
         parameters->AddObject(iomodel->CopyConstantObject(SealevelriseReltolEnum));
         parameters->AddObject(iomodel->CopyConstantObject(SealevelriseAbstolEnum));
@@ -48 +53 @@
         parameters->AddObject(iomodel->CopyConstantObject(SealevelriseElasticEnum));
         parameters->AddObject(iomodel->CopyConstantObject(SealevelriseEustaticEnum));
+        parameters->AddObject(iomodel->CopyConstantObject(SealevelriseLegendrePrecomputeEnum));
+
         
+        /*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){
+                iomodel->FetchData(&legendre_coefficients,&M,&nl,SealevelriseLoveKEnum);
+                parameters->AddObject(new DoubleMatParam(SealevelriseLegendreCoefficientsEnum,legendre_coefficients,M,nl));
+        }
 
         /*free ressources: */
         xDelete<IssmDouble>(love_h);
         xDelete<IssmDouble>(love_k);
+        xDelete<IssmDouble>(legendre_coefficients);
 
 }/*}}}*/

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

@@ -3521 +3521 @@
         int nl;
 
+        /*legendre coefficients:*/
+        bool        legendre_precompute=false;
+        IssmDouble* legendre_coefficients=NULL;
+        int         M;
+
         /*ice properties: */
         IssmDouble rho_ice,rho_water,rho_earth;
-
-        /*other constants: */
-        //IssmDouble g; 
 
         /*Solution outputs: */
@@ -3550 +3552 @@
         this->parameters->FindParam(&computeeustatic,SealevelriseEustaticEnum);
 
+        /*recover legendre coefficients if they have been precomputed:*/
+        this->parameters->FindParam(&legendre_precompute,SealevelriseLegendrePrecomputeEnum);
+        this->parameters->FindParam(&legendre_coefficients,&M,NULL,SealevelriseLegendrePrecomputeEnum);
+
         /*how many dofs are we working with here? */
         this->parameters->FindParam(&gsize,MeshNumberofverticesEnum);
@@ -3614 +3620 @@
                         IssmDouble G_rigid=0;  //do not remove =0!
                         IssmDouble G_elastic=0;  //do not remove =0!
-                        IssmDouble Pn1,Pn2; //first two legendre polynomials
+                        IssmDouble Pn1,Pn2,Pn; //first two legendre polynomials
                         IssmDouble cosalpha;
 
@@ -3631 +3637 @@
                         if(computeelastic){
                                 G_elastic = (love_k[nl-1]-love_h[nl-1])/2.0/sin(alpha/2.0);
-                                for(int n=0;n<nl;n++) G_elastic += deltalove[n]*legendre(Pn1,Pn2,n,cosalpha);
+                                if(legendre_precompute){
+                                        for(int n=0;n<nl;n++) G_elastic += deltalove[n]*legendre_coefficients[(int)((M-1)*(cosalpha+1.0)/2.0)*nl+n];
+                                }
+                                else{
+                                        for(int n=0;n<nl;n++){
+                                                Pn=legendre(Pn1,Pn2,n,cosalpha); Pn1=Pn2; Pn2=Pn;
+                                                G_elastic += deltalove[n]*Pn;
+                                        }
+                                }
                         }
 

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

@@ -1032 +1032 @@
         SealevelriseElasticEnum,
         SealevelriseEustaticEnum,
+        SealevelriseLegendrePrecomputeEnum,
+        SealevelriseLegendreCoefficientsEnum,
         /*}}}*/
         MaximumNumberOfDefinitionsEnum

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

@@ -986 +986 @@
                 case SealevelriseElasticEnum : return "SealevelriseElastic";
                 case SealevelriseEustaticEnum : return "SealevelriseEustatic";
+                case SealevelriseLegendrePrecomputeEnum : return "SealevelriseLegendrePrecompute";
+                case SealevelriseLegendreCoefficientsEnum : return "SealevelriseLegendreCoefficients";
                 case MaximumNumberOfDefinitionsEnum : return "MaximumNumberOfDefinitions";
                 default : return "unknown";

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

@@ -1010 +1010 @@
               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,"MaximumNumberOfDefinitions")==0) return MaximumNumberOfDefinitionsEnum;
          else stage=10;

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

@@ -15 +15 @@
                 elastic         = 0;
                 eustatic         = 0;
+                legendre_precompute = 0;
+                legendre_coefficients = NaN;
         end
         methods
@@ -38 +40 @@
                 self.elastic=1;
                 self.eustatic=1;
+                
+                %legendre coefficients: 
+                legendre_precompute = 0;
+                legendre_coefficients = NaN;
 
                 end % }}}
@@ -55 +61 @@
                         end
 
+                        %check that the legendre coefficients have indeed been computed if requested:
+                        if self.legendre_precompute,
+                                md = checkfield(md,'fieldname','slr.legendre_coefficients','NaN',1,'Inf',1,'size',[NaN length(self.love_h)]);
+                        end
+
                 end % }}}
                 function disp(self) % {{{
@@ -68 +79 @@
                         fielddisplay(self,'elastic','elastic earth graviational potential perturbation');
                         fielddisplay(self,'eustatic','eustatic sea level rise');
+                        fielddisplay(self,'legendre_precompute','precompute legendre coefficients? (default is 0)');
+                        if(self.legendre_precompute)
+                                fielddisplay(self,'legendre_coefficients','precomputed legendre coefficients');
+                        end
 
                 end % }}}
@@ -80 +95 @@
                         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','legendre_precompute','format','Boolean');
+                        if(self.legendre_precompute),
+                                WriteData(fid,'object',self,'class','sealevelrise','fieldname','legendre_coefficients','format','DoubleMat','mattype',1);
+                        end
                 end % }}}
                 function savemodeljs(self,fid,modelname) % {{{
@@ -91 +110 @@
                         writejsdouble(fid,[modelname '.slr.rigid'],self.rigid);
                         writejsdouble(fid,[modelname '.slr.eustatic'],self.eustatic);
+                        writejsdouble(fid,[modelname '.slr.legendre_precompute'],self.legendre_precompute);
+                        if self.legendre_precompute,
+                                writejs2Darray(fid,[modelname '.srl.legendre_coefficients'],self.legendre_coefficients);
+                        end
 
                 end % }}}

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

@@ -971 +971 @@
 function SealevelriseElasticEnum(){ return 967;}
 function SealevelriseEustaticEnum(){ return 968;}
-function MaximumNumberOfDefinitionsEnum(){ return 969;}
+function SealevelriseLegendrePrecomputeEnum(){ return 969;}
+function SealevelriseLegendreCoefficientsEnum(){ return 970;}
+function MaximumNumberOfDefinitionsEnum(){ return 971;}

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

@@ -978 +978 @@
 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 MaximumNumberOfDefinitionsEnum(): return StringToEnum("MaximumNumberOfDefinitions")[0]