Changeset 20999


Ignore:
Timestamp:
07/26/16 14:42:18 (9 years ago)
Author:
adhikari
Message:

CHG: added capabilities to compute absolute sea level and 3D crustal motion

Location:
issm/trunk-jpl/src/c
Files:
16 edited

Legend:

Unmodified
Added
Removed

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

    r20704 r20999  
    4343        IssmDouble* love_h=NULL;
    4444        IssmDouble* love_k=NULL;
     45        IssmDouble* love_l=NULL;
    4546       
    4647        bool elastic=false;
    4748        IssmDouble* G_elastic = NULL;
    4849        IssmDouble* G_elastic_local = NULL;
     50        IssmDouble* U_elastic = NULL;
     51        IssmDouble* U_elastic_local = NULL;
     52        IssmDouble* H_elastic = NULL;
     53        IssmDouble* H_elastic_local = NULL;
    4954        int         M,m,lower_row,upper_row;
    5055        IssmDouble  degacc=.01;
    … …  
    7277                iomodel->FetchData(&love_h,&nl,NULL,"md.slr.love_h");
    7378                iomodel->FetchData(&love_k,&nl,NULL,"md.slr.love_k");
     79                iomodel->FetchData(&love_l,&nl,NULL,"md.slr.love_l");
    7480
    7581                /*compute elastic green function for a range of angles*/
    … …  
    7783                M=reCast<int,IssmDouble>(180./degacc+1.);
    7884                G_elastic=xNew<IssmDouble>(M);
     85                U_elastic=xNew<IssmDouble>(M);
     86                H_elastic=xNew<IssmDouble>(M);
    7987               
    8088                /*compute combined legendre + love number (elastic green function:*/
    … …  
    8290                GetOwnershipBoundariesFromRange(&lower_row,&upper_row,m,IssmComm::GetComm());
    8391                G_elastic_local=xNew<IssmDouble>(m);
     92                U_elastic_local=xNew<IssmDouble>(m);
     93                H_elastic_local=xNew<IssmDouble>(m);
    8494
    8595                for(int i=lower_row;i<upper_row;i++){
    … …  
    8898
    8999                        G_elastic_local[i-lower_row]= (love_k[nl-1]-love_h[nl-1])/2.0/sin(alpha/2.0);
     100                        U_elastic_local[i-lower_row]= (love_h[nl-1])/2.0/sin(alpha/2.0);
     101                        H_elastic_local[i-lower_row]= 0;
    90102                        IssmDouble Pn,Pn1,Pn2;
     103                        IssmDouble Pn_p,Pn_p1,Pn_p2;
    91104                        for (int n=0;n<nl;n++) {
    92                                 IssmDouble deltalove;
    93 
    94                                 deltalove = (love_k[n]-love_k[nl-1]-love_h[n]+love_h[nl-1]);
    95 
    96                                 if(n==0)Pn=1;
    97                                 else if(n==1)Pn=cos(alpha);
    98                                 else Pn= ( (2*n-1)*cos(alpha)*Pn1 - (n-1)*Pn2 ) /n;
     105                                IssmDouble deltalove_G;
     106                                IssmDouble deltalove_U;
     107
     108                                deltalove_G = (love_k[n]-love_k[nl-1]-love_h[n]+love_h[nl-1]);
     109                                deltalove_U = (love_h[n]-love_h[nl-1]);
     110               
     111                                /*compute legendre polynomials: P_n(cos\theta) & d P_n(cos\theta)/ d\theta: */
     112                                if(n==0){
     113                                        Pn=1;
     114                                        Pn_p=0;
     115                                }
     116                                else if(n==1){
     117                                        Pn = cos(alpha);
     118                                        Pn_p = 1;
     119                                }
     120                                else{
     121                                        Pn = ( (2*n-1)*cos(alpha)*Pn1 - (n-1)*Pn2 ) /n;
     122                                        Pn_p = ( (2*n-1)*(Pn1+cos(alpha)*Pn_p1) - (n-1)*Pn_p2 ) /n;
     123                                }
    99124                                Pn2=Pn1; Pn1=Pn;
    100 
    101                                 G_elastic_local[i-lower_row] += deltalove*Pn;
     125                                Pn_p2=Pn_p1; Pn_p1=Pn_p;
     126
     127                                G_elastic_local[i-lower_row] += deltalove_G*Pn;         // gravitational potential
     128                                U_elastic_local[i-lower_row] += deltalove_U*Pn;         // vertical (up) displacement
     129                                H_elastic_local[i-lower_row] += sin(alpha)*love_l[n]*Pn_p;              // horizontal displacements
    102130                        }
    103131                }
    104132
    105                 /*merge G_elastic_local into G_elastic:{{{*/
     133                /*merge G_elastic_local into G_elastic; U_elastic_local into U_elastic; H_elastic_local to H_elastic:{{{*/
    106134                int* recvcounts=xNew<int>(IssmComm::GetSize());
    107135                int* displs=xNew<int>(IssmComm::GetSize());
    … …  
    115143                /*All gather:*/
    116144                ISSM_MPI_Allgatherv(G_elastic_local, m, ISSM_MPI_DOUBLE, G_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
     145                ISSM_MPI_Allgatherv(U_elastic_local, m, ISSM_MPI_DOUBLE, U_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
     146                ISSM_MPI_Allgatherv(H_elastic_local, m, ISSM_MPI_DOUBLE, H_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
    117147                /*free ressources: */
    118148                xDelete<int>(recvcounts);
    … …  
    124154                G_elastic[0]=G_elastic[1];
    125155                parameters->AddObject(new DoubleVecParam(SealevelriseGElasticEnum,G_elastic,M));
     156                U_elastic[0]=U_elastic[1];
     157                parameters->AddObject(new DoubleVecParam(SealevelriseUElasticEnum,U_elastic,M));
     158                H_elastic[0]=H_elastic[1];
     159                parameters->AddObject(new DoubleVecParam(SealevelriseHElasticEnum,H_elastic,M));
    126160
    127161                /*free ressources: */
    128162                xDelete<IssmDouble>(love_h);
    129163                xDelete<IssmDouble>(love_k);
     164                xDelete<IssmDouble>(love_l);
    130165                xDelete<IssmDouble>(G_elastic);
    131166                xDelete<IssmDouble>(G_elastic_local);
     167                xDelete<IssmDouble>(U_elastic);
     168                xDelete<IssmDouble>(U_elastic_local);
     169                xDelete<IssmDouble>(H_elastic);
     170                xDelete<IssmDouble>(H_elastic_local);
    132171        }
    133172       

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

    r20965 r20999  
    15011501                                name==MaterialsRheologyEsbarEnum ||
    15021502                                name==SealevelEnum ||
     1503                                name==SealevelUmotionEnum ||
     1504                                name==SealevelNmotionEnum ||
     1505                                name==SealevelEmotionEnum ||
     1506                                name==SealevelAbsoluteEnum ||
    15031507                                name==SealevelEustaticEnum ||
    15041508                                name==SealevelriseDeltathicknessEnum ||

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

    r20871 r20999  
    303303                virtual void          SealevelriseEustatic(Vector<IssmDouble>* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea)=0;
    304304                virtual void          SealevelriseNonEustatic(Vector<IssmDouble>* pSgo,IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea)=0;
     305                virtual void          SealevelriseGeodetic(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,IssmDouble* Sg,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble oceanarea,IssmDouble eartharea)=0;
    305306                #endif
    306307

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

    r20871 r20999  
    188188                void    SealevelriseEustatic(Vector<IssmDouble>* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){_error_("not implemented yet!");};
    189189                void    SealevelriseNonEustatic(Vector<IssmDouble>* pSgo,IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){_error_("not implemented yet!");};
     190                void    SealevelriseGeodetic(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,IssmDouble* Sg,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble oceanarea,IssmDouble eartharea){_error_("not implemented yet!");};
    190191                #endif
    191192

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

    r20871 r20999  
    172172                void    SealevelriseEustatic(Vector<IssmDouble>* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){_error_("not implemented yet!");};
    173173                void    SealevelriseNonEustatic(Vector<IssmDouble>* pSgo,IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){_error_("not implemented yet!");};
     174                void    SealevelriseGeodetic(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,IssmDouble* Sg,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble oceanarea,IssmDouble eartharea){_error_("not implemented yet!");};
    174175                IssmDouble    OceanArea(void){_error_("not implemented yet!");};
    175176                IssmDouble    OceanAverage(IssmDouble* Sg){_error_("not implemented yet!");};

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

    r20871 r20999  
    178178                void    SealevelriseEustatic(Vector<IssmDouble>* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){_error_("not implemented yet!");};
    179179                void    SealevelriseNonEustatic(Vector<IssmDouble>* pSgo,IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){_error_("not implemented yet!");};
     180                void    SealevelriseGeodetic(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,IssmDouble* Sg,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble oceanarea,IssmDouble eartharea){_error_("not implemented yet!");};
    180181                IssmDouble    OceanArea(void){_error_("not implemented yet!");};
    181182                IssmDouble    OceanAverage(IssmDouble* Sg){_error_("not implemented yet!");};

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

    r20993 r20999  
    38833883}
    38843884/*}}}*/
     3885void    Tria::SealevelriseGeodetic(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,IssmDouble* Sg,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble oceanarea,IssmDouble eartharea){ /*{{{*/
     3886
     3887        /*diverse:*/
     3888        int gsize;
     3889        bool spherical=true;
     3890        IssmDouble llr_list[NUMVERTICES][3];
     3891        IssmDouble xyz_list[NUMVERTICES][3];
     3892        IssmDouble area;
     3893        IssmDouble I;           //change in relative sea level or ice thickness
     3894        IssmDouble late,longe,re;
     3895        IssmDouble lati,longi,ri;
     3896        IssmDouble rho_ice,rho_water,rho_earth;
     3897        IssmDouble minlong=400;
     3898        IssmDouble maxlong=-20;
     3899
     3900        /*precomputed elastic green functions:*/
     3901        IssmDouble* U_elastic_precomputed = NULL;
     3902        IssmDouble* H_elastic_precomputed = NULL;
     3903        int         M;
     3904       
     3905        /*computation of Green functions:*/
     3906        IssmDouble* U_elastic= NULL;
     3907        IssmDouble* N_elastic= NULL;
     3908        IssmDouble* E_elastic= NULL;
     3909
     3910        /*optimization:*/
     3911        bool store_green_functions=false;
     3912
     3913        /*computational flags:*/
     3914        bool computerigid = true;
     3915        bool computeelastic= true;
     3916
     3917        /*early return if we are not on the ocean or on an ice cap:*/
     3918        if(!(this->inputs->Max(MaskIceLevelsetEnum)<0) && !IsWaterInElement()) return;
     3919
     3920        /*recover computational flags: */
     3921        this->parameters->FindParam(&computerigid,SealevelriseRigidEnum);
     3922        this->parameters->FindParam(&computeelastic,SealevelriseElasticEnum);
     3923       
     3924        /*early return if rigid or elastic not requested:*/
     3925        if(!computerigid && !computeelastic) return;
     3926
     3927        /*recover material parameters: */
     3928        rho_ice=matpar->GetMaterialParameter(MaterialsRhoIceEnum);
     3929        rho_water=matpar->GetMaterialParameter(MaterialsRhoFreshwaterEnum);
     3930        rho_earth=matpar->GetMaterialParameter(MaterialsEarthDensityEnum);
     3931
     3932        /*how many dofs are we working with here? */
     3933        this->parameters->FindParam(&gsize,MeshNumberofverticesEnum);
     3934
     3935        /*compute area of element:*/
     3936        area=GetAreaSpherical();
     3937
     3938        /*element centroid (spherical): */
     3939        /* Where is the centroid of this element?:{{{*/
     3940        ::GetVerticesCoordinates(&llr_list[0][0],this->vertices,NUMVERTICES,spherical);
     3941
     3942        minlong=400; maxlong=-20;
     3943        for (int i=0;i<NUMVERTICES;i++){
     3944                llr_list[i][0]=(90-llr_list[i][0]);
     3945                if(llr_list[i][1]<0)llr_list[i][1]=180+(180+llr_list[i][1]);
     3946                if(llr_list[i][1]>maxlong)maxlong=llr_list[i][1];
     3947                if(llr_list[i][1]<minlong)minlong=llr_list[i][1];
     3948        }
     3949        if(minlong==0 && maxlong>180){
     3950                if (llr_list[0][1]==0)llr_list[0][1]=360;
     3951                if (llr_list[1][1]==0)llr_list[1][1]=360;
     3952                if (llr_list[2][1]==0)llr_list[2][1]=360;
     3953        }
     3954
     3955        // correction at the north pole
     3956        if(llr_list[0][0]==0)llr_list[0][1]=(llr_list[1][1]+llr_list[2][1])/2.0;
     3957        if(llr_list[1][0]==0)llr_list[1][1]=(llr_list[0][1]+llr_list[2][1])/2.0;
     3958        if(llr_list[2][0]==0)llr_list[2][1]=(llr_list[0][1]+llr_list[1][1])/2.0;
     3959
     3960        //correction at the south pole
     3961        if(llr_list[0][0]==180)llr_list[0][1]=(llr_list[1][1]+llr_list[2][1])/2.0;
     3962        if(llr_list[1][0]==180)llr_list[1][1]=(llr_list[0][1]+llr_list[2][1])/2.0;
     3963        if(llr_list[2][0]==180)llr_list[2][1]=(llr_list[0][1]+llr_list[1][1])/2.0;
     3964
     3965        late=(llr_list[0][0]+llr_list[1][0]+llr_list[2][0])/3.0;
     3966        longe=(llr_list[0][1]+llr_list[1][1]+llr_list[2][1])/3.0;
     3967
     3968        late=90-late;
     3969        if(longe>180)longe=(longe-180)-180;
     3970
     3971        late=late/180*PI;
     3972        longe=longe/180*PI;
     3973        /*}}}*/
     3974
     3975        /*figure out gravity center of our element (Cartesian): */
     3976        IssmDouble x_element, y_element, z_element;
     3977        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     3978        x_element=(xyz_list[0][0]+xyz_list[1][0]+xyz_list[2][0])/3.0;
     3979        y_element=(xyz_list[0][1]+xyz_list[1][1]+xyz_list[2][1])/3.0;
     3980        z_element=(xyz_list[0][2]+xyz_list[1][2]+xyz_list[2][2])/3.0;
     3981
     3982        if(computeelastic){
     3983       
     3984                /*recover elastic Green's functions for displacement:*/
     3985                DoubleVecParam* U_parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelriseUElasticEnum)); _assert_(parameter);
     3986                DoubleVecParam* H_parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelriseHElasticEnum)); _assert_(parameter);
     3987                U_parameter->GetParameterValueByPointer(&U_elastic_precomputed,&M);
     3988                H_parameter->GetParameterValueByPointer(&H_elastic_precomputed,&M);
     3989
     3990                /*initialize: */
     3991                U_elastic=xNewZeroInit<IssmDouble>(gsize);
     3992                N_elastic=xNewZeroInit<IssmDouble>(gsize);
     3993                E_elastic=xNewZeroInit<IssmDouble>(gsize);
     3994        }
     3995
     3996        int* indices=xNew<int>(gsize);
     3997        IssmDouble* U_values=xNewZeroInit<IssmDouble>(gsize);
     3998        IssmDouble* N_values=xNewZeroInit<IssmDouble>(gsize);
     3999        IssmDouble* E_values=xNewZeroInit<IssmDouble>(gsize);
     4000
     4001        for(int i=0;i<gsize;i++){
     4002
     4003                indices[i]=i;
     4004                if(computeelastic){
     4005       
     4006                        IssmDouble alpha;
     4007                        IssmDouble delPhi,delLambda;
     4008
     4009                        /*Compute alpha angle between centroid and current vertex: */
     4010                        lati=latitude[i]/180*PI; longi=longitude[i]/180*PI;
     4011
     4012                        delPhi=fabs(lati-late); delLambda=fabs(longi-longe);
     4013                        alpha=2.*asin(sqrt(pow(sin(delPhi/2),2.0)+cos(lati)*cos(late)*pow(sin(delLambda/2),2)));
     4014
     4015                        /*Compute azimuths, both north and east components: */
     4016                        IssmDouble dx, dy, dz, x, y, z;
     4017                        IssmDouble N_azim, E_azim;
     4018                        x = xx[i]; y = yy[i]; z = zz[i];
     4019                        if(latitude[i]==90){
     4020                                x=1e-12; y=1e-12;
     4021                        }
     4022                        if(latitude[i]==-90){
     4023                                x=1e-12; y=1e-12;
     4024                        }
     4025                        dx = x_element-x; dy = y_element-y; dz = z_element-z;
     4026                        N_azim = (-z*x*dx-z*y*dy+(pow(x,2)+pow(y,2))*dz) /pow((pow(x,2)+pow(y,2))*(pow(x,2)+pow(y,2)+pow(z,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5);
     4027         E_azim = (-y*dx+x*dy) /pow((pow(x,2)+pow(y,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5);
     4028                       
     4029                        /*Elastic component  (from Eq 17 in Adhikari et al, GMD 2015): */
     4030                        int index=reCast<int,IssmDouble>(alpha/PI*(M-1));
     4031                        U_elastic[i] += U_elastic_precomputed[index];
     4032                        N_elastic[i] += H_elastic_precomputed[index]*N_azim;
     4033                        E_elastic[i] += H_elastic_precomputed[index]*E_azim;
     4034                }
     4035
     4036                /*Add all components to the pUp solution vectors:*/
     4037                if(computerigid){
     4038                        U_values[i]+=0; N_values[i]+=0; E_values[i]+=0;
     4039                }
     4040                if(computeelastic){
     4041                       
     4042                        if(this->inputs->Max(MaskIceLevelsetEnum)<0){
     4043                               
     4044                                /*Compute ice thickness change: */
     4045                                Input*  deltathickness_input=inputs->GetInput(SealevelriseDeltathicknessEnum);
     4046                                if (!deltathickness_input)_error_("delta thickness input needed to compute sea level rise!");
     4047                                deltathickness_input->GetInputAverage(&I);
     4048                       
     4049                                U_values[i]+=3*rho_ice/rho_earth*area/eartharea*I*U_elastic[i];
     4050                                N_values[i]+=3*rho_ice/rho_earth*area/eartharea*I*N_elastic[i];
     4051                                E_values[i]+=3*rho_ice/rho_earth*area/eartharea*I*E_elastic[i];
     4052                        }
     4053                        else if(IsWaterInElement()) {
     4054                       
     4055                                /*From Sg, recover water sea level rise:*/
     4056                                I=0; for(int i=0;i<NUMVERTICES;i++) I+=Sg[this->vertices[i]->Sid()]/NUMVERTICES;
     4057
     4058                                U_values[i]+=3*rho_water/rho_earth*area/eartharea*I*U_elastic[i];
     4059                                N_values[i]+=3*rho_water/rho_earth*area/eartharea*I*N_elastic[i];
     4060                                E_values[i]+=3*rho_water/rho_earth*area/eartharea*I*E_elastic[i];
     4061                        }
     4062                }
     4063        }
     4064        pUp->SetValues(gsize,indices,U_values,ADD_VAL);
     4065        pNorth->SetValues(gsize,indices,N_values,ADD_VAL);
     4066        pEast->SetValues(gsize,indices,E_values,ADD_VAL);
     4067
     4068        /*free ressources:*/
     4069        xDelete<int>(indices);
     4070        xDelete<IssmDouble>(U_values); xDelete<IssmDouble>(N_values); xDelete<IssmDouble>(E_values);
     4071
     4072        /*Free ressources:*/
     4073        if(computeelastic) {
     4074                xDelete<IssmDouble>(U_elastic); xDelete<IssmDouble>(N_elastic); xDelete<IssmDouble>(E_elastic);
     4075        }
     4076
     4077        return;
     4078}
     4079/*}}}*/
    38854080#endif
    3886 
    38874081
    38884082#ifdef _HAVE_DAKOTA_

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

    r20871 r20999  
    149149                void    SealevelriseEustatic(Vector<IssmDouble>* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea);
    150150                void    SealevelriseNonEustatic(Vector<IssmDouble>* pSgo,IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea);
     151                void    SealevelriseGeodetic(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,IssmDouble* Sg,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble oceanarea,IssmDouble eartharea);
    151152                #endif
    152153                /*}}}*/

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

    r20983 r20999  
    23992399}
    24002400/*}}}*/
     2401void FemModel::SealevelriseGeodetic(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, Vector<IssmDouble>* pSg, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz){/*{{{*/
     2402
     2403        /*serialized vectors:*/
     2404        IssmDouble* Sg=NULL;
     2405       
     2406        IssmDouble  oceanarea=0;
     2407        IssmDouble  oceanarea_cpu=0;
     2408        IssmDouble  eartharea=0;
     2409        IssmDouble  eartharea_cpu=0;
     2410
     2411        int         ns,nsmax;
     2412       
     2413        /*Serialize vectors from previous iteration:*/
     2414        Sg=pSg->ToMPISerial();
     2415
     2416        /*Go through elements, and add contribution from each element to the deflection vector wg:*/
     2417        ns = elements->Size();
     2418       
     2419        /*First, figure out the area of the ocean, which is needed to compute the eustatic component: */
     2420        for(int i=0;i<ns;i++){
     2421                Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
     2422                oceanarea_cpu += element->OceanArea();
     2423                eartharea_cpu += element->GetAreaSpherical();
     2424        }
     2425        ISSM_MPI_Reduce (&oceanarea_cpu,&oceanarea,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
     2426        ISSM_MPI_Bcast(&oceanarea,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
     2427       
     2428        ISSM_MPI_Reduce (&eartharea_cpu,&eartharea,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
     2429        ISSM_MPI_Bcast(&eartharea,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
     2430
     2431        /*Figure out max of ns: */
     2432        ISSM_MPI_Reduce(&ns,&nsmax,1,ISSM_MPI_INT,ISSM_MPI_MAX,0,IssmComm::GetComm());
     2433        ISSM_MPI_Bcast(&nsmax,1,ISSM_MPI_INT,0,IssmComm::GetComm());
     2434
     2435        /*Call the sea level rise core: */
     2436        for(int i=0;i<nsmax;i++){
     2437                if(i<ns){
     2438                        Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
     2439                        element->SealevelriseGeodetic(pUp,pNorth,pEast,Sg,latitude,longitude,radius,xx,yy,zz,oceanarea,eartharea);
     2440                }
     2441                if(i%100==0){
     2442                        pUp->Assemble();
     2443                        pNorth->Assemble();
     2444                        pEast->Assemble();
     2445                }
     2446        }
     2447       
     2448        /*One last time: */
     2449        pUp->Assemble();
     2450        pNorth->Assemble();
     2451        pEast->Assemble();
     2452
     2453        /*Free ressources:*/
     2454        xDelete<IssmDouble>(Sg);
     2455        xDelete<IssmDouble>(latitude);
     2456        xDelete<IssmDouble>(longitude);
     2457        xDelete<IssmDouble>(radius);
     2458        xDelete<IssmDouble>(xx);
     2459        xDelete<IssmDouble>(yy);
     2460        xDelete<IssmDouble>(zz);
     2461}
     2462/*}}}*/
    24012463IssmDouble FemModel::SealevelriseOceanAverage(Vector<IssmDouble>* Sg) { /*{{{*/
    24022464
    … …  
    24232485        ISSM_MPI_Reduce (&oceanvalue_cpu,&oceanvalue,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
    24242486        ISSM_MPI_Bcast(&oceanvalue,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
    2425 
    24262487
    24272488        /*Free ressources:*/

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

    r20935 r20999  
    116116                void SealevelriseEustatic(Vector<IssmDouble>* pSgi, IssmDouble* peustatic, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius);
    117117                void SealevelriseNonEustatic(Vector<IssmDouble>* pSgo, Vector<IssmDouble>* pSg_old, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius,bool verboseconvolution);
     118                void SealevelriseGeodetic(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, Vector<IssmDouble>* pSg_old, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz);
    118119                IssmDouble SealevelriseOceanAverage(Vector<IssmDouble>* Sg);
    119120                #endif

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

    r20348 r20999  
    1313
    1414        Vector<IssmDouble> *Sg    = NULL;
    15         Vector<IssmDouble> *Sg_eustatic    = NULL;
     15        Vector<IssmDouble> *Sg_absolute  = NULL;
     16        Vector<IssmDouble> *Sg_eustatic  = NULL;
     17        Vector<IssmDouble> *U_radial  = NULL;
     18        Vector<IssmDouble> *U_north   = NULL;
     19        Vector<IssmDouble> *U_east    = NULL;
    1620        bool save_results,isslr,iscoupler;
    1721        int configuration_type;
    … …  
    2024        char     **requested_outputs = NULL;
    2125       
     26        /*additional parameters: */
     27        int  gsize;
     28        bool spherical=true;
     29        IssmDouble          *latitude   = NULL;
     30        IssmDouble          *longitude  = NULL;
     31        IssmDouble          *radius     = NULL;
     32        IssmDouble          *xx     = NULL;
     33        IssmDouble          *yy     = NULL;
     34        IssmDouble          *zz     = NULL;
     35
    2236        /*Recover some parameters: */
    2337        femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
    … …  
    2640        femmodel->parameters->FindParam(&isslr,TransientIsslrEnum);
    2741        femmodel->parameters->FindParam(&iscoupler,TransientIscouplerEnum);
     42
     43        /*first, recover lat,long and radius vectors from vertices: */
     44        VertexCoordinatesx(&latitude,&longitude,&radius,femmodel->vertices,spherical);
     45
     46        /*recover x,y,z vectors from vertices: */
     47        VertexCoordinatesx(&xx,&yy,&zz,femmodel->vertices);
     48
     49        /*Figure out size of g-set deflection vector and allocate solution vector: */
     50        gsize      = femmodel->nodes->NumberOfDofs(configuration_type,GsetEnum);
    2851
    2952        /*several cases here, depending on value of iscoupler and isslr:
    … …  
    5780
    5881                Sg=sealevelrise_core_noneustatic(femmodel,Sg_eustatic); //ocean loading tems  (2nd and 5th terms on the RHS of Farrel and Clark)
    59 
     82               
    6083                /*get results into elements:*/
    61                 InputUpdateFromSolutionx(femmodel,Sg);
    62 
     84                //InputUpdateFromSolutionx(femmodel,Sg);                // from Eric
     85                InputUpdateFromVectorx(femmodel,Sg,SealevelEnum,VertexSIdEnum);
     86
     87                /*compute other geodetic signatures, such as absolute sea level chagne, components of 3-D crustal motion: */
     88                /*Initialize:*/
     89                U_radial = new Vector<IssmDouble>(gsize);
     90                U_north = new Vector<IssmDouble>(gsize);
     91                U_east = new Vector<IssmDouble>(gsize);
     92                Sg_absolute = new Vector<IssmDouble>(gsize);
     93               
     94                /*call the geodetic main modlule:*/
     95                femmodel->SealevelriseGeodetic(U_radial,U_north,U_east,Sg,latitude,longitude,radius,xx,yy,zz);
     96
     97                /*compute: absolute sea level change = relative sea level change + vertical motion*/
     98                Sg->Copy(Sg_absolute); Sg_absolute->AXPY(U_radial,1);
     99               
     100                /*get results into elements:*/
     101                InputUpdateFromVectorx(femmodel,U_radial,SealevelUmotionEnum,VertexSIdEnum);    // radial displacement
     102                InputUpdateFromVectorx(femmodel,U_north,SealevelNmotionEnum,VertexSIdEnum);     // north motion
     103                InputUpdateFromVectorx(femmodel,U_east,SealevelEmotionEnum,VertexSIdEnum);              // east motion
     104                InputUpdateFromVectorx(femmodel,Sg_absolute,SealevelAbsoluteEnum,VertexSIdEnum);
     105               
    63106                if(save_results){
    64107                        if(VerboseSolution()) _printf0_("   saving results\n");
    … …  
    66109                        femmodel->RequestedOutputsx(&femmodel->results,requested_outputs,numoutputs);
    67110                }
    68 
     111                       
    69112                if(solution_type==SealevelriseSolutionEnum)femmodel->RequestedDependentsx();
    70113
    … …  
    72115                delete Sg;
    73116                delete Sg_eustatic;
     117                delete U_radial;
     118                delete U_north;
     119                delete U_east;
     120                delete Sg_absolute;
    74121                if(numoutputs){for(int i=0;i<numoutputs;i++){xDelete<char>(requested_outputs[i]);} xDelete<char*>(requested_outputs);}
    75122        }

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

    r20367 r20999  
    5858        return Sgi;
    5959}
     60

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

    r20007 r20999  
    1212void slrconvergence(bool* pconverged, Vector<IssmDouble>* Sg,Vector<IssmDouble>* Sg_old,IssmDouble eps_rel,IssmDouble eps_abs);
    1313
    14 Vector<IssmDouble>* sealevelrise_core_noneustatic(FemModel* femmodel,Vector<IssmDouble>* Sg_eustatic){
     14Vector<IssmDouble>* sealevelrise_core_noneustatic(FemModel* femmodel,Vector<IssmDouble>* Sg_eustatic){ /*{{{*/
    1515
    1616        Vector<IssmDouble> *Sg    = NULL;
    … …  
    3535        IssmDouble          *radius    = NULL;
    3636        IssmDouble           eustatic;
    37 
    3837
    3938        /*Recover some parameters: */
    … …  
    112111
    113112        return Sg;
    114 }
     113} /*}}}*/
    115114
    116115void slrconvergence(bool* pconverged, Vector<IssmDouble>* Sg,Vector<IssmDouble>* Sg_old,IssmDouble eps_rel,IssmDouble eps_abs){ /*{{{*/
    … …  
    157156
    158157} /*}}}*/
     158

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

    r20983 r20999  
    10651065        SealevelriseAnalysisEnum,
    10661066        SealevelEnum,
     1067        SealevelUmotionEnum,
     1068        SealevelNmotionEnum,
     1069        SealevelEmotionEnum,
     1070        SealevelAbsoluteEnum,
    10671071        SealevelEustaticEnum,
    10681072        SealevelriseDeltathicknessEnum,
    … …  
    10781082        SealevelriseRotationEnum,
    10791083        SealevelriseGElasticEnum,
     1084        SealevelriseUElasticEnum,
     1085        SealevelriseHElasticEnum,
    10801086        SealevelriseDegaccEnum,
    10811087        SealevelriseTransitionsEnum,

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

    r20983 r20999  
    10171017                case SealevelriseAnalysisEnum : return "SealevelriseAnalysis";
    10181018                case SealevelEnum : return "Sealevel";
     1019                case SealevelUmotionEnum : return "SealevelUmotion";
     1020                case SealevelNmotionEnum : return "SealevelNmotion";
     1021                case SealevelEmotionEnum : return "SealevelEmotion";
     1022                case SealevelAbsoluteEnum : return "SealevelAbsolute";
    10191023                case SealevelEustaticEnum : return "SealevelEustatic";
    10201024                case SealevelriseDeltathicknessEnum : return "SealevelriseDeltathickness";
    … …  
    10301034                case SealevelriseRotationEnum : return "SealevelriseRotation";
    10311035                case SealevelriseGElasticEnum : return "SealevelriseGElastic";
     1036                case SealevelriseUElasticEnum : return "SealevelriseUElastic";
     1037                case SealevelriseHElasticEnum : return "SealevelriseHElastic";
    10321038                case SealevelriseDegaccEnum : return "SealevelriseDegacc";
    10331039                case SealevelriseTransitionsEnum : return "SealevelriseTransitions";

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

    r20983 r20999  
    10411041              else if (strcmp(name,"SealevelriseAnalysis")==0) return SealevelriseAnalysisEnum;
    10421042              else if (strcmp(name,"Sealevel")==0) return SealevelEnum;
     1043              else if (strcmp(name,"SealevelUmotion")==0) return SealevelUmotionEnum;
     1044              else if (strcmp(name,"SealevelNmotion")==0) return SealevelNmotionEnum;
     1045              else if (strcmp(name,"SealevelEmotion")==0) return SealevelEmotionEnum;
     1046              else if (strcmp(name,"SealevelAbsolute")==0) return SealevelAbsoluteEnum;
    10431047              else if (strcmp(name,"SealevelEustatic")==0) return SealevelEustaticEnum;
    10441048              else if (strcmp(name,"SealevelriseDeltathickness")==0) return SealevelriseDeltathicknessEnum;
    … …  
    10541058              else if (strcmp(name,"SealevelriseRotation")==0) return SealevelriseRotationEnum;
    10551059              else if (strcmp(name,"SealevelriseGElastic")==0) return SealevelriseGElasticEnum;
     1060              else if (strcmp(name,"SealevelriseUElastic")==0) return SealevelriseUElasticEnum;
     1061              else if (strcmp(name,"SealevelriseHElastic")==0) return SealevelriseHElasticEnum;
    10561062              else if (strcmp(name,"SealevelriseDegacc")==0) return SealevelriseDegaccEnum;
    10571063              else if (strcmp(name,"SealevelriseTransitions")==0) return SealevelriseTransitionsEnum;
Note: See TracChangeset for help on using the changeset viewer.