Changeset 21272


Ignore:
Timestamp:
10/14/16 12:30:14 (8 years ago)
Author:
adhikari
Message:

CHG: Elastostatic adjustment on a 2D plan geometry

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

Legend:

Unmodified
Added
Removed

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

    r21261 r21272  
    3232        /*Create inputs: */
    3333        iomodel->FetchDataToInput(elements,"md.mask.ice_levelset",MaskIceLevelsetEnum);
    34         iomodel->FetchDataToInput(elements,"md.mask.ocean_levelset",MaskOceanLevelsetEnum);
    35         iomodel->FetchDataToInput(elements,"md.mask.land_levelset",MaskLandLevelsetEnum);
    3634        iomodel->FetchDataToInput(elements,"md.esa.deltathickness",EsaDeltathicknessEnum);
    3735

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

    r21261 r21272  
    297297                #endif
    298298                #ifdef _HAVE_ESA_
    299                 virtual void          EsaGeodetic(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea)=0;
     299                virtual void          EsaGeodetic2D(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,IssmDouble* xx,IssmDouble* yy)=0;
     300                virtual void          EsaGeodetic3D(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea)=0;
    300301                #endif
    301302                #ifdef _HAVE_SEALEVELRISE_

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

    r21261 r21272  
    184184                #endif
    185185                #ifdef _HAVE_ESA_
    186                 void    EsaGeodetic(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea){_error_("not implemented yet!");};
     186                void    EsaGeodetic2D(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,IssmDouble* xx,IssmDouble* yy){_error_("not implemented yet!");};
     187                void    EsaGeodetic3D(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea){_error_("not implemented yet!");};
    187188                #endif
    188189                #ifdef _HAVE_SEALEVELRISE_

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

    r21261 r21272  
    169169#endif
    170170#ifdef _HAVE_ESA_
    171                 void    EsaGeodetic(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea){_error_("not implemented yet!");};
     171                void    EsaGeodetic2D(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,IssmDouble* xx,IssmDouble* yy){_error_("not implemented yet!");};
     172                void    EsaGeodetic3D(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea){_error_("not implemented yet!");};
    172173#endif
    173174#ifdef _HAVE_SEALEVELRISE_

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

    r21261 r21272  
    176176#endif
    177177#ifdef _HAVE_ESA_
    178                 void    EsaGeodetic(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea){_error_("not implemented yet!");};
     178                void    EsaGeodetic2D(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,IssmDouble* xx,IssmDouble* yy){_error_("not implemented yet!");};
     179                void    EsaGeodetic3D(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea){_error_("not implemented yet!");};
    179180#endif
    180181#ifdef _HAVE_SEALEVELRISE_

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

    r21261 r21272  
    35813581#endif
    35823582#ifdef _HAVE_ESA_
    3583 void    Tria::EsaGeodetic(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea){ /*{{{*/
     3583void    Tria::EsaGeodetic2D(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,IssmDouble* xx,IssmDouble* yy){ /*{{{*/
     3584
     3585        /*diverse:*/
     3586        int gsize;
     3587        IssmDouble xyz_list[NUMVERTICES][3];
     3588        IssmDouble area;
     3589        IssmDouble earth_radius = 6371012.0;    // Earth's radius [m]
     3590        IssmDouble I;           //ice/water loading
     3591        IssmDouble rho_ice, rho_earth;
     3592
     3593        /*precomputed elastic green functions:*/
     3594        IssmDouble* U_elastic_precomputed = NULL;
     3595        IssmDouble* H_elastic_precomputed = NULL;
     3596        int         M;
     3597       
     3598        /*computation of Green functions:*/
     3599        IssmDouble* U_elastic= NULL;
     3600        IssmDouble* N_elastic= NULL;
     3601        IssmDouble* E_elastic= NULL;
     3602       
     3603        /*optimization:*/
     3604        bool store_green_functions=false;
     3605
     3606        /*Compute ice thickness change: */
     3607        Input*  deltathickness_input=inputs->GetInput(EsaDeltathicknessEnum);
     3608        if (!deltathickness_input)_error_("delta thickness input needed to compute elastic adjustment!");
     3609        deltathickness_input->GetInputAverage(&I);
     3610               
     3611        /*early return if we are not on the (ice) loading point: */
     3612        if(I==0) return;
     3613
     3614        /*recover material parameters: */
     3615        rho_ice=matpar->GetMaterialParameter(MaterialsRhoIceEnum);
     3616        rho_earth=matpar->GetMaterialParameter(MaterialsEarthDensityEnum);
     3617
     3618        /*how many dofs are we working with here? */
     3619        this->parameters->FindParam(&gsize,MeshNumberofverticesEnum);
     3620
     3621        /*compute area of element:*/
     3622        area=GetArea();
     3623
     3624        /*figure out gravity center of our element (Cartesian): */
     3625        IssmDouble x_element, y_element;
     3626        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     3627        x_element=(xyz_list[0][0]+xyz_list[1][0]+xyz_list[2][0])/3.0;
     3628        y_element=(xyz_list[0][1]+xyz_list[1][1]+xyz_list[2][1])/3.0;
     3629
     3630        /*recover elastic Green's functions for displacement:*/
     3631        DoubleVecParam* U_parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(EsaUElasticEnum)); _assert_(U_parameter);
     3632        DoubleVecParam* H_parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(EsaHElasticEnum)); _assert_(H_parameter);
     3633        U_parameter->GetParameterValueByPointer(&U_elastic_precomputed,&M);
     3634        H_parameter->GetParameterValueByPointer(&H_elastic_precomputed,&M);
     3635
     3636        /*initialize: */
     3637        U_elastic=xNewZeroInit<IssmDouble>(gsize);
     3638        N_elastic=xNewZeroInit<IssmDouble>(gsize);
     3639        E_elastic=xNewZeroInit<IssmDouble>(gsize);
     3640
     3641        int* indices=xNew<int>(gsize);
     3642        IssmDouble* U_values=xNewZeroInit<IssmDouble>(gsize);
     3643        IssmDouble* N_values=xNewZeroInit<IssmDouble>(gsize);
     3644        IssmDouble* E_values=xNewZeroInit<IssmDouble>(gsize);
     3645        IssmDouble dx, dy;
     3646        IssmDouble dist, alpha, ang;
     3647        IssmDouble N_azim, E_azim;
     3648
     3649        for(int i=0;i<gsize;i++){
     3650
     3651                indices[i]=i;
     3652                IssmDouble N_azim=0;
     3653                IssmDouble E_azim=0;
     3654
     3655                /*Compute alpha angle between centroid and current vertex: */
     3656                dx = x_element - xx[i];         dy = y_element - yy[i];
     3657                dist = sqrt(pow(dx,2)+pow(dy,2));                                               // distance between vertex and elemental centroid [m]
     3658                alpha = dist*360.0/(2*PI*earth_radius) * PI/180.0;      // [in radians] 360 degree = 2*pi*earth_radius
     3659
     3660                /*Compute azimuths, both north and east components: */
     3661                ang = PI/2 - atan2(dy,dx);              // this is bearing angle!
     3662                N_azim = cos(ang);
     3663                E_azim = sin(ang);
     3664
     3665                /*Elastic component  (from Eq 17 in Adhikari et al, GMD 2015): */
     3666                int index=reCast<int,IssmDouble>(alpha/PI*(M-1));
     3667                U_elastic[i] += U_elastic_precomputed[index];
     3668                N_elastic[i] += H_elastic_precomputed[index]*N_azim;
     3669                E_elastic[i] += H_elastic_precomputed[index]*E_azim;
     3670
     3671                /*Add all components to the pUp solution vectors:*/
     3672                U_values[i]+=3*rho_ice/rho_earth*area/(4*PI*pow(earth_radius,2))*I*U_elastic[i];
     3673                N_values[i]+=3*rho_ice/rho_earth*area/(4*PI*pow(earth_radius,2))*I*N_elastic[i];
     3674                E_values[i]+=3*rho_ice/rho_earth*area/(4*PI*pow(earth_radius,2))*I*E_elastic[i];
     3675        }
     3676
     3677        pUp->SetValues(gsize,indices,U_values,ADD_VAL);
     3678        pNorth->SetValues(gsize,indices,N_values,ADD_VAL);
     3679        pEast->SetValues(gsize,indices,E_values,ADD_VAL);
     3680       
     3681        /*free ressources:*/
     3682        xDelete<int>(indices);
     3683        xDelete<IssmDouble>(U_values); xDelete<IssmDouble>(N_values); xDelete<IssmDouble>(E_values);
     3684        xDelete<IssmDouble>(U_elastic); xDelete<IssmDouble>(N_elastic); xDelete<IssmDouble>(E_elastic);
     3685
     3686        return;
     3687}
     3688/*}}}*/
     3689void    Tria::EsaGeodetic3D(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea){ /*{{{*/
    35843690
    35853691        /*diverse:*/

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

    r21261 r21272  
    145145                #endif
    146146                #ifdef _HAVE_ESA_
    147                 void    EsaGeodetic(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea);
     147                void    EsaGeodetic2D(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,IssmDouble* xx,IssmDouble* yy);
     148                void    EsaGeodetic3D(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea);
    148149                #endif
    149150                #ifdef _HAVE_SEALEVELRISE_

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

    r21261 r21272  
    23002300#endif
    23012301#ifdef _HAVE_ESA_
    2302 void FemModel::EsaGeodetic(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz){/*{{{*/
     2302void FemModel::EsaGeodetic2D(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, IssmDouble* xx, IssmDouble* yy){/*{{{*/
     2303
     2304        int         ns,nsmax;
     2305       
     2306        /*Go through elements, and add contribution from each element to the deflection vector wg:*/
     2307        ns = elements->Size();
     2308       
     2309        /*Figure out max of ns: */
     2310        ISSM_MPI_Reduce(&ns,&nsmax,1,ISSM_MPI_INT,ISSM_MPI_MAX,0,IssmComm::GetComm());
     2311        ISSM_MPI_Bcast(&nsmax,1,ISSM_MPI_INT,0,IssmComm::GetComm());
     2312
     2313        /*Call the esa geodetic core: */
     2314        for(int i=0;i<nsmax;i++){
     2315                if(i<ns){
     2316                        Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
     2317                        element->EsaGeodetic2D(pUp,pNorth,pEast,xx,yy);
     2318                }
     2319                if(i%100==0){
     2320                        pUp->Assemble();
     2321                        pNorth->Assemble();
     2322                        pEast->Assemble();
     2323                }
     2324        }
     2325       
     2326        /*One last time: */
     2327        pUp->Assemble();
     2328        pNorth->Assemble();
     2329        pEast->Assemble();
     2330
     2331        /*Free ressources:*/
     2332        xDelete<IssmDouble>(xx);
     2333        xDelete<IssmDouble>(yy);
     2334}
     2335/*}}}*/
     2336void FemModel::EsaGeodetic3D(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz){/*{{{*/
    23032337
    23042338        IssmDouble  eartharea=0;
    … …  
    23262360                if(i<ns){
    23272361                        Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
    2328                         element->EsaGeodetic(pUp,pNorth,pEast,latitude,longitude,radius,xx,yy,zz,eartharea);
     2362                        element->EsaGeodetic3D(pUp,pNorth,pEast,latitude,longitude,radius,xx,yy,zz,eartharea);
    23292363                }
    23302364                if(i%100==0){

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

    r21261 r21272  
    114114                #endif
    115115                #ifdef _HAVE_ESA_
    116                 void EsaGeodetic(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz);
     116                void EsaGeodetic2D(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, IssmDouble* xx, IssmDouble* yy);
     117                void EsaGeodetic3D(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz);
    117118                #endif
    118119                #ifdef _HAVE_SEALEVELRISE_

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

    r21261 r21272  
    1717        bool save_results,isesa,iscoupler;
    1818        int configuration_type;
     19        int domaintype;
    1920        int solution_type;
    2021        int        numoutputs        = 0;
    … …  
    3233
    3334        /*Recover some parameters: */
     35        femmodel->parameters->FindParam(&domaintype,DomainTypeEnum);
    3436        femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
    3537        femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
    … …  
    3840        femmodel->parameters->FindParam(&iscoupler,TransientIscouplerEnum);
    3941
    40         /*first, recover lat,long and radius vectors from vertices: */
     42        /* recover coordinates of vertices: */
    4143        VertexCoordinatesx(&latitude,&longitude,&radius,femmodel->vertices,spherical);
    42 
    43         /*recover x,y,z vectors from vertices: */
    4444        VertexCoordinatesx(&xx,&yy,&zz,femmodel->vertices);
    4545
    … …  
    7878               
    7979                /*call the geodetic main modlule:*/
    80                 femmodel->EsaGeodetic(U_radial,U_north,U_east,latitude,longitude,radius,xx,yy,zz);
     80                if(domaintype==Domain3DsurfaceEnum){
     81                        femmodel->EsaGeodetic3D(U_radial,U_north,U_east,latitude,longitude,radius,xx,yy,zz);
     82                }
     83                if(domaintype==Domain2DhorizontalEnum){
     84                        femmodel->EsaGeodetic2D(U_radial,U_north,U_east,xx,yy);
     85                }
    8186
    8287                /*get results into elements:*/
Note: See TracChangeset for help on using the changeset viewer.