Changeset 20007


Ignore:
Timestamp:
01/27/16 15:25:57 (9 years ago)
Author:
Eric.Larour
Message:

CHG: validation of the sealevelrise_core against Surendra's matlab code.
Also, improvements in the treatment of the Green function, which can be reused
in the non eustatic part of the solution.

Location:
issm/trunk-jpl/src/c
Files:
2 added
11 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/Makefile.am

    r19984 r20007  
    468468if SEALEVELRISE
    469469issm_sources +=  ./cores/sealevelrise_core.cpp\
    470                                         ./analyses/SealevelriseAnalysis.cpp
     470                                 ./cores/sealevelrise_core_eustatic.cpp\
     471                                 ./cores/sealevelrise_core_noneustatic.cpp\
     472                                 ./analyses/SealevelriseAnalysis.cpp
    471473endif
    472474#}}}

  • issm/trunk-jpl/src/c/classes/Elements/Element.h

    r19987 r20007  
    302302                #endif
    303303                #ifdef _HAVE_SEALEVELRISE_
    304                 virtual void    Sealevelrise(Vector<IssmDouble>* pSgi,Vector<IssmDouble>* pSgo,IssmDouble* peustatic,IssmDouble* Sg_old,IssmDouble* Sgi_old,IssmDouble* Sgo_old,IssmDouble* x,IssmDouble* y,IssmDouble* z,IssmDouble oceanarea)=0;
     304                virtual void          SealevelriseEustatic(Vector<IssmDouble>* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea)=0;
     305                virtual void          SealevelriseNonEustatic(Vector<IssmDouble>* pSgo,IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea)=0;
     306                virtual IssmDouble    OceanAverage(IssmDouble* Sg)=0;
    305307                virtual IssmDouble    OceanArea(void)=0;
     308                virtual IssmDouble    GetArea3D(void)=0;
    306309                #endif
    307310

  • issm/trunk-jpl/src/c/classes/Elements/Penta.h

    r19984 r20007  
    171171                void           ViscousHeating(IssmDouble* pphi,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input);
    172172                void           ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum);
     173                IssmDouble     GetArea3D(void){_error_("not implemented yet!");};
    173174
    174175                #ifdef _HAVE_DAKOTA_
    … …  
    181182                #endif
    182183                #ifdef _HAVE_SEALEVELRISE_
    183                 void    Sealevelrise(Vector<IssmDouble>* pSgi,Vector<IssmDouble>* pSgo,IssmDouble* peustatic,IssmDouble* Sg_old,IssmDouble* Sgi_old,IssmDouble* Sgo_old,IssmDouble* x,IssmDouble* y,IssmDouble* z, IssmDouble oceanarea){_error_("not implemented yet!");};
     184                void    SealevelriseEustatic(Vector<IssmDouble>* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){_error_("not implemented yet!");};
     185                void    SealevelriseNonEustatic(Vector<IssmDouble>* pSgo,IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){_error_("not implemented yet!");};
    184186                IssmDouble    OceanArea(void){_error_("not implemented yet!");};
     187                IssmDouble    OceanAverage(IssmDouble* Sg){_error_("not implemented yet!");};
    185188                #endif
    186189

  • issm/trunk-jpl/src/c/classes/Elements/Seg.h

    r19984 r20007  
    161161                void        ViscousHeating(IssmDouble* pphi,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){_error_("not implemented yet");};
    162162                void        ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum){_error_("not implemented");};
     163                IssmDouble     GetArea3D(void){_error_("not implemented yet!");};
    163164
    164165#ifdef _HAVE_GIA_
    … …  
    167168
    168169#ifdef _HAVE_SEALEVELRISE_
    169                 void    Sealevelrise(Vector<IssmDouble>* pSgi,Vector<IssmDouble>* pSgo,IssmDouble* peustatic,IssmDouble* Sg_old,IssmDouble* Sgi_old,IssmDouble* Sgo_old,IssmDouble* x,IssmDouble* y,IssmDouble* z,IssmDouble oceanarea){_error_("not implemented yet!");};
     170                void    SealevelriseEustatic(Vector<IssmDouble>* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){_error_("not implemented yet!");};
     171                void    SealevelriseNonEustatic(Vector<IssmDouble>* pSgo,IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){_error_("not implemented yet!");};
    170172                IssmDouble    OceanArea(void){_error_("not implemented yet!");};
     173                IssmDouble    OceanAverage(IssmDouble* Sg){_error_("not implemented yet!");};
    171174#endif
    172175

  • issm/trunk-jpl/src/c/classes/Elements/Tetra.h

    r19984 r20007  
    168168                void        ViscousHeating(IssmDouble* pphi,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input);
    169169                void        ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum);
     170                IssmDouble     GetArea3D(void){_error_("not implemented yet!");};
    170171
    171172#ifdef _HAVE_GIA_
    … …  
    173174#endif
    174175#ifdef _HAVE_SEALEVELRISE_
    175                 void    Sealevelrise(Vector<IssmDouble>* pSgi,Vector<IssmDouble>* pSgo,IssmDouble* peustatic,IssmDouble* Sg_old,IssmDouble* Sgi_old,IssmDouble* Sgo_old,IssmDouble* x,IssmDouble* y,IssmDouble* z,IssmDouble oceanarea){_error_("not implemented yet!");};
     176                void    SealevelriseEustatic(Vector<IssmDouble>* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){_error_("not implemented yet!");};
     177                void    SealevelriseNonEustatic(Vector<IssmDouble>* pSgo,IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){_error_("not implemented yet!");};
    176178                IssmDouble    OceanArea(void){_error_("not implemented yet!");};
     179                IssmDouble    OceanAverage(IssmDouble* Sg){_error_("not implemented yet!");};
    177180#endif
    178181

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

    r19991 r20007  
    35033503
    35043504#ifdef _HAVE_SEALEVELRISE_
    3505 void    Tria::Sealevelrise(Vector<IssmDouble>* pSgi,Vector<IssmDouble>* pSgo,IssmDouble* peustatic,IssmDouble* Sg_old,IssmDouble* Sgi_old,IssmDouble* Sgo_old,IssmDouble* x,IssmDouble* y,IssmDouble* z,IssmDouble oceanarea){ /*{{{*/
     3505void    Tria::SealevelriseEustatic(Vector<IssmDouble>* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){ /*{{{*/
    35063506
    35073507        /*diverse:*/
    35083508        int gsize;
    35093509        bool spherical=true;
    3510         IssmDouble x0,y0,z0,a;
    3511         IssmDouble xyz_list[NUMVERTICES][3];
     3510        IssmDouble llr_list[NUMVERTICES][3];
    35123511        IssmDouble area;
    35133512        IssmDouble I;  //change in ice thickness or water level(Farrel and Clarke, Equ. 4)
    3514         IssmDouble Me;
    35153513        IssmDouble rho;
    3516        
     3514        IssmDouble late,longe,re;
     3515        IssmDouble lati,longi,ri;
     3516
    35173517        /*love numbers:*/
    35183518        IssmDouble* love_h=NULL;
    … …  
    35293529        IssmDouble rho_ice,rho_water,rho_earth;
    35303530
    3531         /*Solution outputs: */
    3532         Vector<IssmDouble>* pSolution=NULL;
    3533        
    35343531        /*Initialize eustatic component: do not skip this step :):*/
    35353532        IssmDouble eustatic = 0;
    … …  
    35393536        bool computeelastic= true;
    35403537        bool computeeustatic = true;
     3538
    35413539       
     3540        /*early return if we are not on an ice cap:*/
     3541        if(IsWaterInElement() | !IsIceInElement()){
     3542                *peustatic=0; //do not forget to assign this pointer, otherwise, global eustatic will be garbage!
     3543                return;
     3544        }
     3545
    35423546        /*recover material parameters: */
    35433547        rho_ice=matpar->GetMaterialParameter(MaterialsRhoIceEnum);
    … …  
    35543558        /*recover legendre coefficients if they have been precomputed:*/
    35553559        this->parameters->FindParam(&legendre_precompute,SealevelriseLegendrePrecomputeEnum);
    3556         this->parameters->FindParam(&legendre_coefficients,&M,NULL,SealevelriseLegendreCoefficientsEnum);
     3560        if(legendre_precompute)this->parameters->FindParam(&legendre_coefficients,&M,NULL,SealevelriseLegendreCoefficientsEnum);
    35573561
    35583562        /*how many dofs are we working with here? */
    … …  
    35603564
    35613565        /*retrieve coordinates: */
    3562         ::GetVerticesCoordinates(&xyz_list[0][0],this->vertices,NUMVERTICES);
    3563 
    3564         /* Where is the centroid of this element?. To avoid issues with lat,long
    3565          * being between [0,180] and [0 360], and issues with jumps at the central
    3566          * meridian and poles, we do everything in cartesian coordinate system: */
    3567         x0=(xyz_list[0][0]+xyz_list[1][0]+xyz_list[2][0])/3;
    3568         y0=(xyz_list[0][1]+xyz_list[1][1]+xyz_list[2][1])/3;
    3569         z0=(xyz_list[0][2]+xyz_list[1][2]+xyz_list[2][2])/3;
    3570 
    3571         /*radius at this location?:*/
    3572         a=sqrt(pow(xyz_list[0][0],2)+pow(xyz_list[0][1],2)+pow(xyz_list[0][2],2)); //a from Farrel and Clark, Eq 4.
    3573 
    3574         /*Compute mass of the earth:*/
    3575         Me= rho_earth*4/3*PI*pow(a,3);
     3566        ::GetVerticesCoordinates(&llr_list[0][0],this->vertices,NUMVERTICES,spherical);
     3567       
     3568        /* Where is the centroid of this element?:*/
     3569        IssmDouble minlong=400;
     3570        IssmDouble maxlong=-20;
     3571        for (int i=0;i<NUMVERTICES;i++){
     3572                llr_list[i][0]=(90-llr_list[i][0]);
     3573                if(llr_list[i][1]<0)llr_list[i][1]=180+(180+llr_list[i][1]);
     3574                if(llr_list[i][1]>maxlong)maxlong=llr_list[i][1];
     3575                if(llr_list[i][1]<minlong)minlong=llr_list[i][1];
     3576        }
     3577        if(minlong==0 && maxlong>180){
     3578                if (llr_list[0][1]==0)llr_list[0][1]=360;
     3579                if (llr_list[1][1]==0)llr_list[1][1]=360;
     3580                if (llr_list[2][1]==0)llr_list[2][1]=360;
     3581        }
     3582       
     3583        // correction at the north pole
     3584        if(llr_list[0][0]==0)llr_list[0][1]=(llr_list[1][1]+llr_list[2][1])/2.0;
     3585        if(llr_list[1][0]==0)llr_list[1][1]=(llr_list[0][1]+llr_list[2][1])/2.0;
     3586        if(llr_list[2][0]==0)llr_list[2][1]=(llr_list[0][1]+llr_list[1][1])/2.0;
     3587                       
     3588        //correction at the south pole
     3589        if(llr_list[0][0]==180)llr_list[0][1]=(llr_list[1][1]+llr_list[2][1])/2.0;
     3590        if(llr_list[1][0]==180)llr_list[1][1]=(llr_list[0][1]+llr_list[2][1])/2.0;
     3591        if(llr_list[2][0]==180)llr_list[2][1]=(llr_list[0][1]+llr_list[1][1])/2.0;
     3592
     3593        late=(llr_list[0][0]+llr_list[1][0]+llr_list[2][0])/3.0;
     3594        longe=(llr_list[0][1]+llr_list[1][1]+llr_list[2][1])/3.0;
     3595
     3596        late=90-late;
     3597        if(longe>180)longe=(longe-180)-180;
     3598
     3599        late=late/180*PI;
     3600        longe=longe/180*PI;
    35763601
    35773602        /*Compute area of element:*/
    35783603        area=GetArea3D();
    35793604
    3580         /*Compute ice thickness or sea level change: */
    3581         if(IsWaterInElement()){
    3582 
    3583                 IssmDouble phi_I_I_O,phi_O_O_O;
    3584                
    3585                 /*From Sg_old, recover water sea level rise:*/
    3586                 I=0; for(int i=0;i<NUMVERTICES;i++) I+=Sg_old[this->vertices[i]->Sid()]/NUMVERTICES;
    3587                 rho=rho_water;
    3588                 pSolution=pSgo;
    3589 
    3590                 /*Compute eustatic component: */
    3591                 if(computeeustatic){
    3592                         phi_I_I_O=0; for(int i=0;i<NUMVERTICES;i++) phi_I_I_O+=area/oceanarea * Sgi_old[this->vertices[i]->Sid()]/NUMVERTICES;
    3593                         phi_O_O_O=0; for(int i=0;i<NUMVERTICES;i++) phi_O_O_O+=area/oceanarea * Sgo_old[this->vertices[i]->Sid()]/NUMVERTICES;
    3594                         eustatic -= (phi_I_I_O + phi_O_O_O); //Watch out the sign "-"!
    3595                 }
    3596 
    3597         }
    3598         else if(IsIceInElement()){
    3599                 Input*  deltathickness_input=inputs->GetInput(SealevelriseDeltathicknessEnum);
    3600                 if (!deltathickness_input)_error_("delta thickness input needed to compute sea level rise!");
    3601                 deltathickness_input->GetInputAverage(&I);
    3602                 rho=rho_ice;
    3603                 pSolution=pSgi;
    3604 
    3605                 /*Compute eustatic compoent:*/
    3606                 if(computeeustatic)eustatic -= rho_ice*area*I/(oceanarea*rho_water); //Watch out the sign "-"!
    3607         }
    3608         else if(IsLandInElement()){
    3609                 //do nothing
    3610         }
    3611         else _error_("Tria::Sealevelrise error: partition of land, ice and ocean is not complete!");
     3605        /*Compute ice thickness change: */
     3606        Input*  deltathickness_input=inputs->GetInput(SealevelriseDeltathicknessEnum);
     3607        if (!deltathickness_input)_error_("delta thickness input needed to compute sea level rise!");
     3608        deltathickness_input->GetInputAverage(&I);
     3609
     3610        /*Compute eustatic compoent:*/
     3611        if(computeeustatic) eustatic += rho_ice*area*I/(oceanarea*rho_water);
    36123612
    36133613        /*Speed up of love number comutation: */
    … …  
    36253625                        IssmDouble Pn1,Pn2,Pn; //first two legendre polynomials
    36263626                        IssmDouble cosalpha;
    3627 
    3628                         /*Compute alpha angle between centroid and current vertex and cosine of this angle: */
    3629                         alpha = acos(
    3630                                         (x[i]*x0+y[i]*y0+z[i]*z0)  //scalar product  of two position vectors
    3631                                         / sqrt(pow(x0,2.0)+pow(y0,2.0)+pow(z0,2.0))  //norm of first position vector
    3632                                         / sqrt(pow(x[i],2.0)+pow(y[i],2.0)+pow(z[i],2.0))  //norm of second position vector
    3633                                         );
     3627                        IssmDouble delPhi,delLambda;
     3628
     3629                        /*Compute alpha angle between centroid and current vertex : */
     3630                        lati=latitude[i]/180*PI; longi=longitude[i]/180*PI;
     3631
     3632                    delPhi=fabs(lati-late); delLambda=fabs(longi-longe);
     3633                        alpha=2.*asin(sqrt(pow(sin(delPhi/2),2.0)+cos(lati)*cos(late)*pow(sin(delLambda/2),2)));
    36343634                        cosalpha=cos(alpha); //compute this to speed up the elastic component.
    36353635
    … …  
    36523652
    36533653                        /*Add all components to the pSgi or pSgo solution vectors:*/
    3654                         pSolution->SetValue(i,rho*a/Me*I*area*(G_rigid+G_elastic),ADD_VAL);
     3654                        pSgi->SetValue(i,3*rho_ice/rho_earth*area/eartharea*I*(G_rigid+G_elastic),ADD_VAL);
     3655                       
    36553656                }
    36563657        }
    … …  
    36603661
    36613662        /*Free ressources:*/
    3662         xDelete<IssmDouble>(love_h);
    3663         xDelete<IssmDouble>(love_k);
    3664         xDelete<IssmDouble>(deltalove);
     3663        if(computeelastic){
     3664                xDelete<IssmDouble>(love_h);
     3665                xDelete<IssmDouble>(love_k);
     3666                xDelete<IssmDouble>(deltalove);
     3667        }
    36653668        return;
     3669}
     3670/*}}}*/
     3671void    Tria::SealevelriseNonEustatic(Vector<IssmDouble>* pSgo,IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){ /*{{{*/
     3672
     3673        /*diverse:*/
     3674        int gsize;
     3675        bool spherical=true;
     3676        IssmDouble llr_list[NUMVERTICES][3];
     3677        IssmDouble area;
     3678        IssmDouble I;  //change in water water level(Farrel and Clarke, Equ. 4)
     3679        IssmDouble late,longe,re;
     3680        IssmDouble lati,longi,ri;
     3681        IssmDouble minlong=400;
     3682        IssmDouble maxlong=-20;
     3683
     3684        /*love numbers:*/
     3685        IssmDouble* love_h=NULL;
     3686        IssmDouble* love_k=NULL;
     3687        IssmDouble* deltalove=NULL;
     3688        int nl;
     3689
     3690        /*legendre coefficients:*/
     3691        bool        legendre_precompute=false;
     3692        IssmDouble* legendre_coefficients=NULL;
     3693        int         M;
     3694
     3695        /*ice properties: */
     3696        IssmDouble rho_ice,rho_water,rho_earth;
     3697
     3698        /*Computational flags:*/
     3699        bool computerigid = true;
     3700        bool computeelastic= true;
     3701        bool computeeustatic = true;
     3702
     3703        /*G_rigid and G_elasti variables, to speed up the computations: */
     3704        DoubleArrayInput*      G_rigid_input=NULL;
     3705        IssmDouble* G_rigid = NULL; int G_rigid_M;
     3706        bool        compute_G_rigid=false;
     3707        DoubleArrayInput*      G_elastic_input=NULL;
     3708        IssmDouble* G_elastic = NULL; int G_elastic_M;
     3709        bool        compute_G_elastic=false;
     3710
     3711        /*early return if we are not on the ocean:*/
     3712        if (!IsWaterInElement())return;
     3713
     3714        /*recover computational flags: */
     3715        this->parameters->FindParam(&computerigid,SealevelriseRigidEnum);
     3716        this->parameters->FindParam(&computeelastic,SealevelriseElasticEnum);
     3717        this->parameters->FindParam(&computeeustatic,SealevelriseEustaticEnum);
     3718       
     3719        /*early return if rigid or elastic not requested:*/
     3720        if(!computerigid && !computeelastic) return;
     3721
     3722        /*Try and recover, if it exists, G_rigid and G_elastic:*/
     3723        G_rigid_input=(DoubleArrayInput*)this->inputs->GetInput(SealevelriseGRigidEnum);
     3724        if(G_rigid_input){
     3725                compute_G_rigid=false;
     3726                G_rigid_input->GetValues(&G_rigid,&G_rigid_M);
     3727        }
     3728        else if(computerigid)compute_G_rigid=true;
     3729
     3730        G_elastic_input=(DoubleArrayInput*)this->inputs->GetInput(SealevelriseGElasticEnum);
     3731        if(G_elastic_input){
     3732                compute_G_elastic=false;
     3733                G_elastic_input->GetValues(&G_elastic,&G_elastic_M);
     3734        }
     3735        else if(computeelastic)compute_G_elastic=true;
     3736       
     3737        /*recover material parameters: */
     3738        rho_water=matpar->GetMaterialParameter(MaterialsRhoSeawaterEnum);
     3739        rho_earth=matpar->GetMaterialParameter(MaterialsEarthDensityEnum);
     3740
     3741        /*how many dofs are we working with here? */
     3742        this->parameters->FindParam(&gsize,MeshNumberofverticesEnum);
     3743
     3744        /*From Sg_old, recover water sea level rise:*/
     3745        I=0; for(int i=0;i<NUMVERTICES;i++) I+=Sg_old[this->vertices[i]->Sid()]/NUMVERTICES;
     3746
     3747        /*Compute area of element:*/
     3748        area=GetArea3D();
     3749
     3750        /* Where is the centroid of this element?:{{{*/
     3751        ::GetVerticesCoordinates(&llr_list[0][0],this->vertices,NUMVERTICES,spherical);
     3752
     3753        minlong=400; maxlong=-20;
     3754        for (int i=0;i<NUMVERTICES;i++){
     3755                llr_list[i][0]=(90-llr_list[i][0]);
     3756                if(llr_list[i][1]<0)llr_list[i][1]=180+(180+llr_list[i][1]);
     3757                if(llr_list[i][1]>maxlong)maxlong=llr_list[i][1];
     3758                if(llr_list[i][1]<minlong)minlong=llr_list[i][1];
     3759        }
     3760        if(minlong==0 && maxlong>180){
     3761                if (llr_list[0][1]==0)llr_list[0][1]=360;
     3762                if (llr_list[1][1]==0)llr_list[1][1]=360;
     3763                if (llr_list[2][1]==0)llr_list[2][1]=360;
     3764        }
     3765
     3766        // correction at the north pole
     3767        if(llr_list[0][0]==0)llr_list[0][1]=(llr_list[1][1]+llr_list[2][1])/2.0;
     3768        if(llr_list[1][0]==0)llr_list[1][1]=(llr_list[0][1]+llr_list[2][1])/2.0;
     3769        if(llr_list[2][0]==0)llr_list[2][1]=(llr_list[0][1]+llr_list[1][1])/2.0;
     3770
     3771        //correction at the south pole
     3772        if(llr_list[0][0]==180)llr_list[0][1]=(llr_list[1][1]+llr_list[2][1])/2.0;
     3773        if(llr_list[1][0]==180)llr_list[1][1]=(llr_list[0][1]+llr_list[2][1])/2.0;
     3774        if(llr_list[2][0]==180)llr_list[2][1]=(llr_list[0][1]+llr_list[1][1])/2.0;
     3775
     3776        late=(llr_list[0][0]+llr_list[1][0]+llr_list[2][0])/3.0;
     3777        longe=(llr_list[0][1]+llr_list[1][1]+llr_list[2][1])/3.0;
     3778
     3779        late=90-late;
     3780        if(longe>180)longe=(longe-180)-180;
     3781
     3782        late=late/180*PI;
     3783        longe=longe/180*PI;
     3784        /*}}}*/
     3785        if(compute_G_elastic){
     3786                /*recover love numbers and legendre_coefficients:{{{*/
     3787                this->parameters->FindParam(&love_h,&nl,SealevelriseLoveHEnum);
     3788                this->parameters->FindParam(&love_k,&nl,SealevelriseLoveKEnum);
     3789
     3790                /*recover legendre coefficients if they have been precomputed:*/
     3791                this->parameters->FindParam(&legendre_precompute,SealevelriseLegendrePrecomputeEnum);
     3792                if(legendre_precompute)this->parameters->FindParam(&legendre_coefficients,&M,NULL,SealevelriseLegendreCoefficientsEnum);
     3793
     3794                /*Speed up of love number comutation for elastic mode: */
     3795                deltalove=xNew<IssmDouble>(nl);
     3796                for (int n=0;n<nl;n++) deltalove[n] = (love_k[n]-love_k[nl-1]-love_h[n]+love_h[nl-1]);
     3797                //}}}
     3798
     3799                /*initialize G_elastic:*/
     3800                G_elastic=xNew<IssmDouble>(gsize);
     3801        }
     3802        if(compute_G_rigid){
     3803                /*Initialize G_rigid and G_elastic:*/
     3804                G_rigid=xNew<IssmDouble>(gsize);
     3805        }
     3806
     3807        for(int i=0;i<gsize;i++){
     3808
     3809                if(compute_G_elastic | compute_G_rigid){
     3810                        IssmDouble alpha;
     3811                        IssmDouble Pn1,Pn2,Pn; //first two legendre polynomials
     3812                        IssmDouble cosalpha;
     3813                        IssmDouble delPhi,delLambda;
     3814
     3815                        /*Compute alpha angle between centroid and current vertex : */
     3816                        lati=latitude[i]/180*PI; longi=longitude[i]/180*PI;
     3817
     3818                        delPhi=fabs(lati-late); delLambda=fabs(longi-longe);
     3819                        alpha=2.*asin(sqrt(pow(sin(delPhi/2),2.0)+cos(lati)*cos(late)*pow(sin(delLambda/2),2)));
     3820                        cosalpha=cos(alpha); //compute this to speed up the elastic component.
     3821
     3822                        //Rigid earth gravitational perturbation:
     3823                        if(compute_G_rigid)G_rigid[i]=1.0/2.0/sin(alpha/2.0);
     3824
     3825                        //Elastic component  (from Eq 17 in Adhikari et al, GMD 2015)
     3826                        if(compute_G_elastic){
     3827                                G_elastic[i] = (love_k[nl-1]-love_h[nl-1])/2.0/sin(alpha/2.0);
     3828                                if(legendre_precompute){
     3829                                        for(int n=0;n<nl;n++) G_elastic[i] += deltalove[n]*legendre_coefficients[reCast<int,IssmDouble>((M-1)*(cosalpha+1.0)/2.0)*nl+n];
     3830                                }
     3831                                else{
     3832                                        for(int n=0;n<nl;n++){
     3833                                                Pn=legendre(Pn1,Pn2,reCast<int,IssmDouble>(cosalpha),n); Pn1=Pn2; Pn2=Pn;
     3834                                                G_elastic[i] += deltalove[n]*Pn;
     3835                                        }
     3836                                }
     3837                        }
     3838                }
     3839
     3840                /*Add all components to the pSgi or pSgo solution vectors:*/
     3841                if(computerigid)pSgo->SetValue(i,3*rho_water/rho_earth*area/eartharea*I*G_rigid[i],ADD_VAL);
     3842                if(computeelastic)pSgo->SetValue(i,3*rho_water/rho_earth*area/eartharea*I*G_elastic[i],ADD_VAL);
     3843        }
     3844
     3845        /*Save G_rigid and G_elastic if computed:*/
     3846        if(G_rigid)this->inputs->AddInput(new DoubleArrayInput(SealevelriseGRigidEnum,G_rigid,gsize));
     3847        if(G_elastic)this->inputs->AddInput(new DoubleArrayInput(SealevelriseGElasticEnum,G_elastic,gsize));
     3848
     3849        /*Free ressources:*/
     3850        if(compute_G_elastic){
     3851                xDelete<IssmDouble>(love_h);
     3852                xDelete<IssmDouble>(love_k);
     3853                xDelete<IssmDouble>(deltalove);
     3854                xDelete<IssmDouble>(G_elastic);
     3855        }
     3856        if(compute_G_rigid) xDelete<IssmDouble>(G_rigid);
     3857
     3858        return;
     3859}
     3860/*}}}*/
     3861IssmDouble Tria::OceanAverage(IssmDouble* Sg){ /*{{{*/
     3862
     3863        if(IsWaterInElement()){
     3864               
     3865                IssmDouble area;
     3866
     3867                /*Compute area of element:*/
     3868                area=GetArea3D();
     3869
     3870                /*Average Sg over vertices:*/
     3871                IssmDouble Sg_avg=0; for(int i=0;i<NUMVERTICES;i++) Sg_avg+=Sg[this->vertices[i]->Sid()]/NUMVERTICES;
     3872
     3873                /*return: */
     3874                return area*Sg_avg;
     3875        }
     3876        else return 0;
     3877
    36663878}
    36673879/*}}}*/

  • issm/trunk-jpl/src/c/classes/Elements/Tria.h

    r19984 r20007  
    144144                #endif
    145145                #ifdef _HAVE_SEALEVELRISE_
    146                 void    Sealevelrise(Vector<IssmDouble>* pSgi,Vector<IssmDouble>* pSgo,IssmDouble* peustatic,IssmDouble* Sg_old,IssmDouble* Sgi_old,IssmDouble* Sgo_old,IssmDouble* x,IssmDouble* y,IssmDouble* z,IssmDouble oceanarea);
     146                void    SealevelriseEustatic(Vector<IssmDouble>* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea);
     147                void    SealevelriseNonEustatic(Vector<IssmDouble>* pSgo,IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea);
     148                IssmDouble OceanAverage(IssmDouble* Sg);
    147149                IssmDouble OceanArea(void);
    148150                #endif

  • issm/trunk-jpl/src/c/classes/FemModel.cpp

    r19984 r20007  
    22062206#endif
    22072207#ifdef _HAVE_SEALEVELRISE_
    2208 void FemModel::Sealevelrise(Vector<IssmDouble>* pSgi,Vector<IssmDouble>* pSgo, IssmDouble* peustatic, Vector<IssmDouble>* pSg_old, Vector<IssmDouble>* pSgi_old,Vector<IssmDouble>* pSgo_old,IssmDouble* x, IssmDouble* y, IssmDouble* z) { /*{{{*/
     2208void FemModel::SealevelriseEustatic(Vector<IssmDouble>* pSgi, IssmDouble* peustatic, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius) { /*{{{*/
    22092209
    22102210        /*serialized vectors:*/
    2211         IssmDouble* Sg_old=NULL;
    2212         IssmDouble* Sgi_old=NULL;
    2213         IssmDouble* Sgo_old=NULL;
    22142211        IssmDouble  eustatic=0;
    22152212        IssmDouble  eustatic_cpu=0;
    … …  
    22172214        IssmDouble  oceanarea=0;
    22182215        IssmDouble  oceanarea_cpu=0;
     2216        IssmDouble  eartharea=0;
     2217        IssmDouble  eartharea_cpu=0;
    22192218        int         ns;
    22202219       
    2221         /*Serialize vectors from previous iteration:*/
    2222 
    2223         Sg_old=pSg_old->ToMPISerial();
    2224         Sgi_old=pSgi_old->ToMPISerial();
    2225         Sgo_old=pSgo_old->ToMPISerial();
    2226 
    22272220        /*Go through elements, and add contribution from each element to the deflection vector wg:*/
    22282221        ns = elements->Size();
    … …  
    22322225                Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
    22332226                oceanarea_cpu += element->OceanArea();
     2227                eartharea_cpu+= element->GetArea3D();
    22342228        }
    22352229        ISSM_MPI_Reduce (&oceanarea_cpu,&oceanarea,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
    22362230        ISSM_MPI_Bcast(&oceanarea,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
     2231       
     2232        ISSM_MPI_Reduce (&eartharea_cpu,&eartharea,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
     2233        ISSM_MPI_Bcast(&eartharea,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
    22372234
    22382235        /*Call the sea level rise core: */
    … …  
    22422239       
    22432240                Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
    2244                 element->Sealevelrise(pSgi,pSgo,&eustatic_cpu_e,Sg_old,Sgi_old,Sgo_old,x,y,z,oceanarea);
     2241                element->SealevelriseEustatic(pSgi,&eustatic_cpu_e,latitude,longitude,radius,oceanarea,eartharea);
    22452242                eustatic_cpu+=eustatic_cpu_e;
    22462243        }
    … …  
    22542251        *peustatic=eustatic;
    22552252
     2253}
     2254/*}}}*/
     2255void FemModel::SealevelriseNonEustatic(Vector<IssmDouble>* pSgo, Vector<IssmDouble>* pSg_old, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, bool verboseconvolution){/*{{{*/
     2256
     2257        /*serialized vectors:*/
     2258        IssmDouble* Sg_old=NULL;
     2259       
     2260        IssmDouble  oceanarea=0;
     2261        IssmDouble  oceanarea_cpu=0;
     2262        IssmDouble  eartharea=0;
     2263        IssmDouble  eartharea_cpu=0;
     2264
     2265        int         ns;
     2266       
     2267        /*Serialize vectors from previous iteration:*/
     2268        Sg_old=pSg_old->ToMPISerial();
     2269
     2270        /*Go through elements, and add contribution from each element to the deflection vector wg:*/
     2271        ns = elements->Size();
     2272
     2273        /*First, figure out the area of the ocean, which is needed to compute the eustatic component: */
     2274        for(int i=0;i<ns;i++){
     2275                Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
     2276                oceanarea_cpu += element->OceanArea();
     2277                eartharea_cpu+= element->GetArea3D();
     2278        }
     2279        ISSM_MPI_Reduce (&oceanarea_cpu,&oceanarea,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
     2280        ISSM_MPI_Bcast(&oceanarea,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
     2281       
     2282        ISSM_MPI_Reduce (&eartharea_cpu,&eartharea,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
     2283        ISSM_MPI_Bcast(&eartharea,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
     2284
     2285        /*Call the sea level rise core: */
     2286        for(int i=0;i<ns;i++){
     2287
     2288                if(verboseconvolution)if(IssmComm::GetRank()==0)if(VerboseConvergence())if(i%1000)_printf_("\r" << "      convolution progress: " << (float)i/(float)ns*100 << "\%");
     2289               
     2290                Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
     2291                element->SealevelriseNonEustatic(pSgo, Sg_old,latitude,longitude,radius,oceanarea,eartharea);
     2292        }
     2293        if(verboseconvolution)if(IssmComm::GetRank()==0)if(VerboseConvergence())_printf_("\n");
     2294
     2295
    22562296        /*Free ressources:*/
    22572297        xDelete<IssmDouble>(Sg_old);
    2258         xDelete<IssmDouble>(Sgi_old);
    2259         xDelete<IssmDouble>(Sgo_old);
     2298}
     2299/*}}}*/
     2300IssmDouble FemModel::SealevelriseOceanAverage(Vector<IssmDouble>* Sg) { /*{{{*/
     2301
     2302        IssmDouble* Sg_serial=NULL;
     2303        IssmDouble  oceanvalue,oceanvalue_cpu;
     2304        IssmDouble  oceanarea,oceanarea_cpu;
     2305
     2306        /*Serialize vectors from previous iteration:*/
     2307        Sg_serial=Sg->ToMPISerial();
     2308
     2309        /*Initialize:*/
     2310        oceanvalue_cpu=0;
     2311        oceanarea_cpu=0;
     2312
     2313        /*Go through elements, and add contribution from each element and divide by overall ocean area:*/
     2314        for(int i=0;i<elements->Size();i++){
     2315                Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
     2316                oceanarea_cpu += element->OceanArea();
     2317                oceanvalue_cpu += element->OceanAverage(Sg_serial);
     2318        }
     2319        ISSM_MPI_Reduce (&oceanarea_cpu,&oceanarea,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
     2320        ISSM_MPI_Bcast(&oceanarea,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
     2321       
     2322        ISSM_MPI_Reduce (&oceanvalue_cpu,&oceanvalue,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
     2323        ISSM_MPI_Bcast(&oceanvalue,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
     2324
     2325
     2326        /*Free ressources:*/
     2327        xDelete<IssmDouble>(Sg_serial);
     2328       
     2329        return oceanvalue/oceanarea;
    22602330}
    22612331/*}}}*/

  • issm/trunk-jpl/src/c/classes/FemModel.h

    r19984 r20007  
    113113                #endif
    114114                #ifdef _HAVE_SEALEVELRISE_
    115                 void Sealevelrise(Vector<IssmDouble>* pSgi,Vector<IssmDouble>* pSgo, IssmDouble* peustatic, Vector<IssmDouble>* pSg_old, Vector<IssmDouble>* pSgi_old,Vector<IssmDouble>* pSgo_old,IssmDouble* x, IssmDouble* y, IssmDouble* z);
     115                void SealevelriseEustatic(Vector<IssmDouble>* pSgi, IssmDouble* peustatic, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius);
     116                void SealevelriseNonEustatic(Vector<IssmDouble>* pSgo, Vector<IssmDouble>* pSg_old, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius,bool verboseconvolution);
     117                IssmDouble SealevelriseOceanAverage(Vector<IssmDouble>* Sg);
    116118                #endif
    117119                void TimeAdaptx(IssmDouble* pdt);

  • issm/trunk-jpl/src/c/cores/cores.h

    r19984 r20007  
    4646void dummy_core(FemModel* femmodel);
    4747void gia_core(FemModel* femmodel);
    48 void sealevelrise_core(FemModel* femmodel);
    4948void smb_core(FemModel* femmodel);
    5049void damage_core(FemModel* femmodel);
     50void sealevelrise_core(FemModel* femmodel);
     51Vector<IssmDouble> * sealevelrise_core_eustatic(FemModel* femmodel);
     52Vector<IssmDouble> * sealevelrise_core_noneustatic(FemModel* femmodel,Vector<IssmDouble>* Sg_eustatic);
    5153IssmDouble objectivefunction(IssmDouble search_scalar,FemModel* femmodel);
    5254

  • issm/trunk-jpl/src/c/cores/sealevelrise_core.cpp

    r19984 r20007  
    1010#include "../solutionsequences/solutionsequences.h"
    1111
    12 void slrconvergence(bool* pconverged, Vector<IssmDouble>* Sg,Vector<IssmDouble>* Sg_old,IssmDouble eps_rel,IssmDouble eps_abs);
    13 
    1412void sealevelrise_core(FemModel* femmodel){
    1513
    1614        Vector<IssmDouble> *Sg    = NULL;
    17         Vector<IssmDouble> *Sg_old    = NULL;
    18         Vector<IssmDouble> *Sgi    = NULL; //ice convolution
    19         Vector<IssmDouble> *Sgi_old    = NULL;
    20         Vector<IssmDouble> *Sgo    = NULL; //ocean convolution
    21         Vector<IssmDouble> *Sgo_old    = NULL;
    22 
    23         /*parameters: */
    24         int count;
     15        Vector<IssmDouble> *Sg_eustatic    = NULL;
    2516        bool save_results;
    26         int  gsize;
    27         int  configuration_type;
    28         bool spherical=true;
    29         bool converged=true;
    30         int max_nonlinear_iterations;
    31         IssmDouble           eps_rel;
    32         IssmDouble           eps_abs;
    33         IssmDouble          *x    = NULL;
    34         IssmDouble          *y    = NULL;
    35         IssmDouble          *z    = NULL;
    36         IssmDouble           eustatic;
    37 
    38 
    39         /*Recover some parameters: */
    40         femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
    41         femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
    42         femmodel->parameters->FindParam(&max_nonlinear_iterations,SealevelriseMaxiterEnum);
    43         femmodel->parameters->FindParam(&eps_rel,SealevelriseReltolEnum);
    44         femmodel->parameters->FindParam(&eps_abs,SealevelriseAbstolEnum);
     17        int configuration_type;
    4518
    4619        if(VerboseSolution()) _printf0_("   computing sea level rise\n");
    4720
    48         /*Call on core computations: */
     21        /*Recover some parameters: */
     22        femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
     23        femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
     24
     25        /*set configuration: */
    4926        femmodel->SetCurrentConfiguration(SealevelriseAnalysisEnum);
    5027
    51         /*first, recover lat,long and radius vectors from vertices: */
    52         VertexCoordinatesx(&x,&y,&z,femmodel->vertices); //no need for z coordinate
     28        /*call two sub cores:*/
     29        Sg_eustatic=sealevelrise_core_eustatic(femmodel); //generalized eustatic (Farrel and Clark, Eq 4, 1st, 3rd and 4rd terms on the RHS.
    5330
    54         /*Figure out size of g-set deflection vector and allocate solution vector: */
    55         gsize      = femmodel->nodes->NumberOfDofs(configuration_type,GsetEnum);
    56        
    57         /*Initialize:*/
    58         Sg = new Vector<IssmDouble>(gsize);
    59         Sg->Assemble();
    60         Sg_old = new Vector<IssmDouble>(gsize);
    61         Sg_old->Assemble();
     31        Sg=sealevelrise_core_noneustatic(femmodel,Sg_eustatic); //sea-level rise dependent terms (2nd and 5th terms on the RHS)
    6232
    63         Sgi = new Vector<IssmDouble>(gsize);
    64         Sgi->Assemble();
    65         Sgi_old = new Vector<IssmDouble>(gsize);
    66         Sgi_old->Assemble();
    67 
    68         Sgo = new Vector<IssmDouble>(gsize);
    69         Sgo->Assemble();
    70         Sgo_old = new Vector<IssmDouble>(gsize);
    71         Sgo_old->Assemble();
    72 
    73         count=1;
    74         converged=false;
    75        
    76         /*Start loop: */
    77         for(;;){
    78 
    79                 //save pointer to old sea level rise
    80                 delete Sgi_old; Sgi_old=Sgi;
    81                 delete Sgo_old; Sgo_old=Sgo;
    82                 delete Sg_old; Sg_old=Sg;
    83 
    84                 /*Initialize solution vector: */
    85                 Sg = new Vector<IssmDouble>(gsize); Sg->Assemble();
    86                 Sgi = new Vector<IssmDouble>(gsize); Sgi->Assemble();
    87                 Sgo = new Vector<IssmDouble>(gsize); Sgo->Assemble();
    88 
    89                 /*call the main module: */
    90         femmodel->Sealevelrise(Sgi,Sgo, &eustatic, Sg_old, Sgi_old, Sgo_old,x, y, z);
    91                
    92                 /*assemble solution vectors: */
    93                 Sgi->Assemble();
    94                 Sgo->Assemble();
    95                
    96                 /*Sg is the sum of the ice and ocean convolutions + eustatic component: (eq 4 in Farrel and Clark)*/
    97                 Sgi->Copy(Sg); Sg->AXPY(Sgo,1);
    98                 Sg->Shift(eustatic);
    99 
    100                 /*convergence criterion:*/
    101                 slrconvergence(&converged,Sg,Sg_old,eps_rel,eps_abs);
    102 
    103                 /*Increase count: */
    104                 count++;
    105                 if(converged==true){
    106                         break;
    107                 }
    108                 if(count>=max_nonlinear_iterations){
    109                         _printf0_("   maximum number of nonlinear iterations (" << max_nonlinear_iterations << ") exceeded\n");
    110                         converged=true;
    111                         break;
    112                 }       
    113                
    114         }
    115         if(VerboseConvergence()) _printf0_("\n   total number of iterations: " << count-1 << "\n");
    116 
     33        /*get results out:*/
    11734        InputUpdateFromVectorx(femmodel,Sg,SealevelriseSEnum,VertexSIdEnum);
    11835
    … …  
    12239                femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],1);
    12340        }
    124         xDelete<IssmDouble>(x);
    125         xDelete<IssmDouble>(y);
    126         xDelete<IssmDouble>(z);
    127         delete Sg_old;
    12841        delete Sg;
     42        delete Sg_eustatic;
    12943}
    130 
    131 void slrconvergence(bool* pconverged, Vector<IssmDouble>* Sg,Vector<IssmDouble>* Sg_old,IssmDouble eps_rel,IssmDouble eps_abs){ /*{{{*/
    132        
    133         bool converged=true;
    134         IssmDouble ndS,nS;
    135         Vector<IssmDouble> *dSg    = NULL;
    136 
    137         //compute norm(du) and norm(u) if requested
    138         dSg=Sg_old->Duplicate(); Sg_old->Copy(dSg); dSg->AYPX(Sg,-1.0);
    139         ndS=dSg->Norm(NORM_TWO);
    140        
    141         if(!xIsNan<IssmDouble>(eps_rel)){
    142                 nS=Sg_old->Norm(NORM_TWO);
    143         }
    144 
    145         if (xIsNan<IssmDouble>(ndS) || xIsNan<IssmDouble>(nS)) _error_("convergence criterion is NaN!");
    146 
    147         //clean up
    148         delete dSg;
    149 
    150         //print
    151         if(!xIsNan<IssmDouble>(eps_rel)){
    152                 if((ndS/nS)<eps_rel){
    153                         if(VerboseConvergence()) _printf0_(setw(50) << left << "      convergence criterion: norm(dS)/norm(S)" << ndS/nS*100 << " < " << eps_rel*100 << " %\n");
    154                 }
    155                 else{
    156                         if(VerboseConvergence()) _printf0_(setw(50) << left << "      convergence criterion: norm(dS)/norm(S)" << ndS/nS*100 << " > " << eps_rel*100 << " %\n");
    157                         converged=false;
    158                 }
    159         }
    160         if(!xIsNan<IssmDouble>(eps_abs)){
    161                 if(ndS<eps_abs){
    162                         if(VerboseConvergence()) _printf0_(setw(50) << left << "      convergence criterion: norm(dS)" << ndS << " < " << eps_abs << " \n");
    163                 }
    164                 else{
    165                         if(VerboseConvergence()) _printf0_(setw(50) << left << "      convergence criterion: norm(dS)" << ndS << " > " << eps_abs << " \n");
    166                         converged=false;
    167                 }
    168         }
    169 
    170         /*assign output*/
    171         *pconverged=converged;
    172 
    173 } /*}}}*/
Note: See TracChangeset for help on using the changeset viewer.