Index: /issm/trunk-jpl/src/c/analyses/SealevelriseAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/SealevelriseAnalysis.cpp (revision 20027) +++ /issm/trunk-jpl/src/c/analyses/SealevelriseAnalysis.cpp (revision 20028) @@ -43,6 +43,5 @@ IssmDouble* love_k=NULL; - bool legendre_precompute=false; - IssmDouble* legendre_coefficients=NULL; + IssmDouble* G_elastic = NULL; int M; @@ -54,35 +53,42 @@ 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(180/degacc+1); + G_elastic=xNew(M); + + /*compute combined legendre + love number (elastic green function:*/ + for(int i=0;i(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(180/degacc+1); - IssmDouble* x=xNew(M); - for(int i=0;i<=M;i++){ //watch out the <= - x[i]=-cos( (reCast(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(x); - } + /*free ressources:*/ + xDelete(G_elastic); /*free ressources: */ xDelete(love_h); xDelete(love_k); - xDelete(legendre_coefficients); + xDelete(G_elastic); }/*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 20027) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 20028) @@ -3515,13 +3515,6 @@ 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,15 +3543,13 @@ /*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(this->parameters->FindParamObject(SealevelriseLegendreCoefficientsEnum)); _assert_(parameter); - parameter->GetParameterValueByPointer(&legendre_coefficients,&M,NULL); + /*recover elastic green function:*/ + if(computeelastic){ + DoubleVecParam* parameter = static_cast(this->parameters->FindParamObject(SealevelriseGElasticEnum)); + _assert_(parameter); + parameter->GetParameterValueByPointer(&G_elastic_precomputed,&M); } @@ -3614,10 +3605,4 @@ if(computeeustatic) eustatic += rho_ice*area*I/(oceanarea*rho_water); - /*Speed up of love number comutation: */ - if(computeelastic){ - deltalove=xNew(nl); - for (int n=0;n(gsize); @@ -3629,6 +3614,4 @@ 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,5 +3621,4 @@ 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,18 +3627,6 @@ //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((M-1)*(1-alpha/PI)); - for(int n=0;n(alpha/PI*(M-1)); + G_elastic += G_elastic_precomputed[index]; } @@ -3673,11 +3643,4 @@ /*Assign output pointer:*/ *peustatic=eustatic; - - /*Free ressources:*/ - if(computeelastic){ - xDelete(love_h); - xDelete(love_k); - xDelete(deltalove); - } return; } @@ -3696,14 +3659,11 @@ 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,12 +3678,4 @@ 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,24 +3685,8 @@ 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,28 +3737,15 @@ 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(this->parameters->FindParamObject(SealevelriseLegendreCoefficientsEnum)); _assert_(parameter); - parameter->GetParameterValueByPointer(&legendre_coefficients,&M,NULL); - } - - /*Speed up of love number comutation for elastic mode: */ - deltalove=xNew(nl); - for (int n=0;n(this->parameters->FindParamObject(SealevelriseGElasticEnum)); _assert_(parameter); + parameter->GetParameterValueByPointer(&G_elastic_precomputed,&M); /*initialize G_elastic:*/ G_elastic=xNewZeroInit(gsize); } - if(compute_G_rigid){ - /*Initialize G_rigid and G_elastic:*/ - G_rigid=xNewZeroInit(gsize); - } + if(computerigid) G_rigid=xNewZeroInit(gsize); int* indices=xNew(gsize); @@ -3832,8 +3755,7 @@ 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,25 +3765,12 @@ 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((M-1)*(1-alpha/PI)); - for(int n=0;n(alpha/PI*(M-1)); + G_elastic[i] += G_elastic_precomputed[index]; } } @@ -3878,18 +3787,7 @@ xDelete(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(love_h); - xDelete(love_k); - xDelete(deltalove); - xDelete(G_elastic); - } - if(compute_G_rigid) xDelete(G_rigid); + if(computeelastic) xDelete(G_elastic); + if(computerigid) xDelete(G_rigid); return; Index: /issm/trunk-jpl/src/c/classes/Params/DoubleVecParam.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Params/DoubleVecParam.cpp (revision 20027) +++ /issm/trunk-jpl/src/c/classes/Params/DoubleVecParam.cpp (revision 20028) @@ -121,2 +121,11 @@ } /*}}}*/ + +/*DoubleVecParam specific routines:*/ +void DoubleVecParam::GetParameterValueByPointer(IssmDouble** pIssmDoublearray,int* pM){/*{{{*/ + + /*Assign output pointers:*/ + if(pM) *pM=M; + *pIssmDoublearray=values; +} +/*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Params/DoubleVecParam.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Params/DoubleVecParam.h (revision 20027) +++ /issm/trunk-jpl/src/c/classes/Params/DoubleVecParam.h (revision 20028) @@ -71,4 +71,8 @@ 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 */ Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 20027) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 20028) @@ -1047,9 +1047,5 @@ SealevelriseElasticEnum, SealevelriseEustaticEnum, - SealevelriseLegendrePrecomputeEnum, - SealevelriseLegendreCoefficientsEnum, - SealevelriseGRigidEnum, SealevelriseGElasticEnum, - SealevelriseStoreGreenFunctionsEnum, /*}}}*/ MaximumNumberOfDefinitionsEnum Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 20027) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 20028) @@ -1001,9 +1001,5 @@ 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"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 20027) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 20028) @@ -1025,9 +1025,5 @@ 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; Index: /issm/trunk-jpl/src/m/classes/slr.m =================================================================== --- /issm/trunk-jpl/src/m/classes/slr.m (revision 20027) +++ /issm/trunk-jpl/src/m/classes/slr.m (revision 20028) @@ -15,6 +15,4 @@ elastic = 0; eustatic = 0; - legendre_precompute = 0; - store_green_functions = 0; end methods @@ -41,10 +39,4 @@ self.eustatic=1; - %legendre coefficients: - self.legendre_precompute = 1; - - %optimization: - self.store_green_functions=1; - end % }}} function md = checkconsistency(self,md,solution,analyses) % {{{ @@ -57,5 +49,4 @@ 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,6 +68,4 @@ 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,6 +80,4 @@ 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,8 +89,6 @@ 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 % }}} Index: /issm/trunk-jpl/src/m/enum/EnumDefinitions.js =================================================================== --- /issm/trunk-jpl/src/m/enum/EnumDefinitions.js (revision 20027) +++ /issm/trunk-jpl/src/m/enum/EnumDefinitions.js (revision 20028) @@ -986,8 +986,4 @@ 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;} Index: /issm/trunk-jpl/src/m/enum/EnumDefinitions.py =================================================================== --- /issm/trunk-jpl/src/m/enum/EnumDefinitions.py (revision 20027) +++ /issm/trunk-jpl/src/m/enum/EnumDefinitions.py (revision 20028) @@ -993,8 +993,4 @@ 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] Index: sm/trunk-jpl/src/m/enum/SealevelriseGRigidEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/SealevelriseGRigidEnum.m (revision 20027) +++ (revision ) @@ -1,11 +1,0 @@ -function macro=SealevelriseGRigidEnum() -%SEALEVELRISEGRIGIDENUM - Enum of SealevelriseGRigid -% -% WARNING: DO NOT MODIFY THIS FILE -% this file has been automatically generated by src/c/shared/Enum/Synchronize.sh -% Please read src/c/shared/Enum/README for more information -% -% Usage: -% macro=SealevelriseGRigidEnum() - -macro=StringToEnum('SealevelriseGRigid'); Index: sm/trunk-jpl/src/m/enum/SealevelriseLegendreCoefficientsEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/SealevelriseLegendreCoefficientsEnum.m (revision 20027) +++ (revision ) @@ -1,11 +1,0 @@ -function macro=SealevelriseLegendreCoefficientsEnum() -%SEALEVELRISELEGENDRECOEFFICIENTSENUM - Enum of SealevelriseLegendreCoefficients -% -% WARNING: DO NOT MODIFY THIS FILE -% this file has been automatically generated by src/c/shared/Enum/Synchronize.sh -% Please read src/c/shared/Enum/README for more information -% -% Usage: -% macro=SealevelriseLegendreCoefficientsEnum() - -macro=StringToEnum('SealevelriseLegendreCoefficients'); Index: sm/trunk-jpl/src/m/enum/SealevelriseLegendrePrecomputeEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/SealevelriseLegendrePrecomputeEnum.m (revision 20027) +++ (revision ) @@ -1,11 +1,0 @@ -function macro=SealevelriseLegendrePrecomputeEnum() -%SEALEVELRISELEGENDREPRECOMPUTEENUM - Enum of SealevelriseLegendrePrecompute -% -% WARNING: DO NOT MODIFY THIS FILE -% this file has been automatically generated by src/c/shared/Enum/Synchronize.sh -% Please read src/c/shared/Enum/README for more information -% -% Usage: -% macro=SealevelriseLegendrePrecomputeEnum() - -macro=StringToEnum('SealevelriseLegendrePrecompute'); Index: sm/trunk-jpl/src/m/enum/SealevelriseStoreGreenFunctionsEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/SealevelriseStoreGreenFunctionsEnum.m (revision 20027) +++ (revision ) @@ -1,11 +1,0 @@ -function macro=SealevelriseStoreGreenFunctionsEnum() -%SEALEVELRISESTOREGREENFUNCTIONSENUM - Enum of SealevelriseStoreGreenFunctions -% -% WARNING: DO NOT MODIFY THIS FILE -% this file has been automatically generated by src/c/shared/Enum/Synchronize.sh -% Please read src/c/shared/Enum/README for more information -% -% Usage: -% macro=SealevelriseStoreGreenFunctionsEnum() - -macro=StringToEnum('SealevelriseStoreGreenFunctions');