Changeset 22349


Ignore:
Timestamp:
01/14/18 12:54:39 (7 years ago)
Author:
adhikari
Message:

CHG: ESA2D can now compute north east components of 2D crustal disp vector

Location:
issm/trunk-jpl/src
Files:
15 edited

Legend:

Unmodified
Added
Removed

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

    r21908 r22349  
    5656        int         *transitions_N    = NULL;
    5757        int          ntransitions;
     58
     59        /*some constant parameters: */
     60        parameters->AddObject(iomodel->CopyConstantObject("md.esa.hemisphere",EsaHemisphereEnum));
    5861
    5962        /*love numbers: */

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

    r22326 r22349  
    17501750                                name==EsaNmotionEnum ||
    17511751                                name==EsaEmotionEnum ||
     1752                                name==EsaXmotionEnum ||
     1753                                name==EsaYmotionEnum ||
    17521754                                name==EsaStrainratexxEnum ||
    17531755                                name==EsaStrainratexyEnum ||

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

    r22326 r22349  
    310310                #endif
    311311                #ifdef _HAVE_ESA_
    312                 virtual void          EsaGeodetic2D(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,IssmDouble* xx,IssmDouble* yy)=0;
     312                virtual void          EsaGeodetic2D(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast, Vector<IssmDouble>* pX, Vector<IssmDouble>* pY,IssmDouble* xx,IssmDouble* yy)=0;
    313313                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;
    314314                #endif

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

    r22326 r22349  
    184184                #endif
    185185                #ifdef _HAVE_ESA_
    186                 void    EsaGeodetic2D(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,IssmDouble* xx,IssmDouble* yy){_error_("not implemented yet!");};
     186                void    EsaGeodetic2D(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,Vector<IssmDouble>* pX,Vector<IssmDouble>* pY,IssmDouble* xx,IssmDouble* yy){_error_("not implemented yet!");};
    187187                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!");};
    188188                #endif

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

    r22326 r22349  
    167167#endif
    168168#ifdef _HAVE_ESA_
    169                 void    EsaGeodetic2D(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,IssmDouble* xx,IssmDouble* yy){_error_("not implemented yet!");};
     169                void    EsaGeodetic2D(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,Vector<IssmDouble>* pX,Vector<IssmDouble>* pY,IssmDouble* xx,IssmDouble* yy){_error_("not implemented yet!");};
    170170                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!");};
    171171#endif

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

    r22326 r22349  
    174174#endif
    175175#ifdef _HAVE_ESA_
    176                 void    EsaGeodetic2D(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,IssmDouble* xx,IssmDouble* yy){_error_("not implemented yet!");};
     176                void    EsaGeodetic2D(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,Vector<IssmDouble>* pX,Vector<IssmDouble>* pY,IssmDouble* xx,IssmDouble* yy){_error_("not implemented yet!");};
    177177                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!");};
    178178#endif

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

    r22330 r22349  
    539539
    540540        /*Retrieve all inputs we will be needing: */
    541         Input* vx_input=this->GetInput(EsaEmotionEnum); _assert_(vx_input);
    542         Input* vy_input=this->GetInput(EsaNmotionEnum); _assert_(vy_input);
     541        Input* vx_input=this->GetInput(EsaXmotionEnum); _assert_(vx_input);
     542        Input* vy_input=this->GetInput(EsaYmotionEnum); _assert_(vy_input);
    543543       
    544544        /* Start looping on the number of vertices: */
    … …  
    37503750#endif
    37513751#ifdef _HAVE_ESA_
    3752 void    Tria::EsaGeodetic2D(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,IssmDouble* xx,IssmDouble* yy){ /*{{{*/
     3752void    Tria::EsaGeodetic2D(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,Vector<IssmDouble>* pX,Vector<IssmDouble>* pY,IssmDouble* xx,IssmDouble* yy){ /*{{{*/
    37533753
    37543754        /*diverse:*/
    … …  
    37633763        IssmDouble* U_elastic_precomputed = NULL;
    37643764        IssmDouble* H_elastic_precomputed = NULL;
    3765         int         M;
     3765        int         M, hemi;
    37663766       
    37673767        /*computation of Green functions:*/
    … …  
    37693769        IssmDouble* N_elastic= NULL;
    37703770        IssmDouble* E_elastic= NULL;
     3771        IssmDouble* X_elastic= NULL;
     3772        IssmDouble* Y_elastic= NULL;
    37713773       
    37723774        /*optimization:*/
    … …  
    37883790        this->parameters->FindParam(&gsize,MeshNumberofverticesEnum);
    37893791
     3792        /*which hemisphere? for north-south, east-west components*/
     3793        this->parameters->FindParam(&hemi,EsaHemisphereEnum);
     3794       
    37903795        /*compute area of element:*/
    37913796        area=GetArea();
    … …  
    38073812        N_elastic=xNewZeroInit<IssmDouble>(gsize);
    38083813        E_elastic=xNewZeroInit<IssmDouble>(gsize);
     3814        X_elastic=xNewZeroInit<IssmDouble>(gsize);
     3815        Y_elastic=xNewZeroInit<IssmDouble>(gsize);
    38093816
    38103817        int* indices=xNew<int>(gsize);
    … …  
    38123819        IssmDouble* N_values=xNewZeroInit<IssmDouble>(gsize);
    38133820        IssmDouble* E_values=xNewZeroInit<IssmDouble>(gsize);
    3814         IssmDouble dx, dy;
    3815         IssmDouble dist, alpha, ang;
    3816         IssmDouble N_azim, E_azim;
     3821        IssmDouble* X_values=xNewZeroInit<IssmDouble>(gsize);
     3822        IssmDouble* Y_values=xNewZeroInit<IssmDouble>(gsize);
     3823        IssmDouble dx, dy, dist, alpha, ang, ang2;
     3824        IssmDouble N_azim, E_azim, X_azim, Y_azim;
    38173825
    38183826        for(int i=0;i<gsize;i++){
    38193827
    38203828                indices[i]=i;
     3829
    38213830                IssmDouble N_azim=0;
    38223831                IssmDouble E_azim=0;
    … …  
    38293838                /*Compute azimuths, both north and east components: */
    38303839                ang = PI/2 - atan2(dy,dx);              // this is bearing angle!
    3831                 N_azim = cos(ang);
    3832                 E_azim = sin(ang);
     3840                Y_azim = cos(ang);
     3841                X_azim = sin(ang);
    38333842
    38343843                /*Elastic component  (from Eq 17 in Adhikari et al, GMD 2015): */
    38353844                int index=reCast<int,IssmDouble>(alpha/PI*(M-1));
    38363845                U_elastic[i] += U_elastic_precomputed[index];
    3837                 N_elastic[i] += H_elastic_precomputed[index]*N_azim;
    3838                 E_elastic[i] += H_elastic_precomputed[index]*E_azim;
     3846                Y_elastic[i] += H_elastic_precomputed[index]*Y_azim;
     3847                X_elastic[i] += H_elastic_precomputed[index]*X_azim;
    38393848
    38403849                /*Add all components to the pUp solution vectors:*/
    38413850                U_values[i]+=3*rho_ice/rho_earth*area/(4*PI*pow(earth_radius,2))*I*U_elastic[i];
    3842                 N_values[i]+=3*rho_ice/rho_earth*area/(4*PI*pow(earth_radius,2))*I*N_elastic[i];
    3843                 E_values[i]+=3*rho_ice/rho_earth*area/(4*PI*pow(earth_radius,2))*I*E_elastic[i];
     3851                Y_values[i]+=3*rho_ice/rho_earth*area/(4*PI*pow(earth_radius,2))*I*Y_elastic[i];
     3852                X_values[i]+=3*rho_ice/rho_earth*area/(4*PI*pow(earth_radius,2))*I*X_elastic[i];
     3853               
     3854                /*North-south, East-west components */
     3855                if (hemi == -1) {
     3856                        ang2 = PI/2 - atan2(yy[i],xx[i]);
     3857                }
     3858                else if (hemi == 1) {
     3859                        ang2 = PI/2 - atan2(-yy[i],-xx[i]);
     3860                }
     3861                if (hemi != 0){
     3862                        N_azim = Y_azim*cos(ang2) + X_azim*sin(ang2);
     3863                        E_azim = X_azim*cos(ang2) - Y_azim*sin(ang2);
     3864                        N_elastic[i] += H_elastic_precomputed[index]*N_azim;
     3865                        E_elastic[i] += H_elastic_precomputed[index]*E_azim;
     3866                        N_values[i]+=3*rho_ice/rho_earth*area/(4*PI*pow(earth_radius,2))*I*N_elastic[i];
     3867                        E_values[i]+=3*rho_ice/rho_earth*area/(4*PI*pow(earth_radius,2))*I*E_elastic[i];
     3868                }
    38443869        }
    38453870
    … …  
    38473872        pNorth->SetValues(gsize,indices,N_values,ADD_VAL);
    38483873        pEast->SetValues(gsize,indices,E_values,ADD_VAL);
     3874        pX->SetValues(gsize,indices,X_values,ADD_VAL);
     3875        pY->SetValues(gsize,indices,Y_values,ADD_VAL);
    38493876       
    38503877        /*free ressources:*/
    … …  
    38523879        xDelete<IssmDouble>(U_values); xDelete<IssmDouble>(N_values); xDelete<IssmDouble>(E_values);
    38533880        xDelete<IssmDouble>(U_elastic); xDelete<IssmDouble>(N_elastic); xDelete<IssmDouble>(E_elastic);
     3881        xDelete<IssmDouble>(X_values); xDelete<IssmDouble>(Y_values);
     3882        xDelete<IssmDouble>(X_elastic); xDelete<IssmDouble>(Y_elastic);
    38543883
    38553884        return;

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

    r22326 r22349  
    143143                #endif
    144144                #ifdef _HAVE_ESA_
    145                 void    EsaGeodetic2D(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,IssmDouble* xx,IssmDouble* yy);
     145                void    EsaGeodetic2D(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,  Vector<IssmDouble>* pX,Vector<IssmDouble>* pY,IssmDouble* xx,IssmDouble* yy);
    146146                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);
    147147                #endif

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

    r22326 r22349  
    39763976#endif
    39773977#ifdef _HAVE_ESA_
    3978 void FemModel::EsaGeodetic2D(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, IssmDouble* xx, IssmDouble* yy){/*{{{*/
     3978void FemModel::EsaGeodetic2D(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, Vector<IssmDouble>* pX, Vector<IssmDouble>* pY, IssmDouble* xx, IssmDouble* yy){/*{{{*/
    39793979
    39803980        int         ns,nsmax;
    … …  
    39913991                if(i<ns){
    39923992                        Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
    3993                         element->EsaGeodetic2D(pUp,pNorth,pEast,xx,yy);
     3993                        element->EsaGeodetic2D(pUp,pNorth,pEast,pX,pY,xx,yy);
    39943994                }
    39953995                if(i%100==0){
    … …  
    39973997                        pNorth->Assemble();
    39983998                        pEast->Assemble();
     3999                        pX->Assemble();
     4000                        pY->Assemble();
    39994001                }
    40004002        }
    … …  
    40044006        pNorth->Assemble();
    40054007        pEast->Assemble();
     4008        pX->Assemble();
     4009        pY->Assemble();
    40064010
    40074011        /*Free ressources:*/

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

    r22326 r22349  
    134134                #endif
    135135                #ifdef _HAVE_ESA_
    136                 void EsaGeodetic2D(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, IssmDouble* xx, IssmDouble* yy);
     136                void EsaGeodetic2D(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, Vector<IssmDouble>* pX, Vector<IssmDouble>* pY, IssmDouble* xx, IssmDouble* yy);
    137137                void EsaGeodetic3D(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz);
    138138                #endif

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

    r21272 r22349  
    1515        Vector<IssmDouble> *U_north   = NULL;
    1616        Vector<IssmDouble> *U_east    = NULL;
     17        Vector<IssmDouble> *U_x   = NULL;
     18        Vector<IssmDouble> *U_y    = NULL;
    1719        bool save_results,isesa,iscoupler;
    1820        int configuration_type;
    … …  
    7678                U_north = new Vector<IssmDouble>(gsize);
    7779                U_east = new Vector<IssmDouble>(gsize);
     80                U_x = new Vector<IssmDouble>(gsize);
     81                U_y = new Vector<IssmDouble>(gsize);
    7882               
    7983                /*call the geodetic main modlule:*/
    … …  
    8286                }
    8387                if(domaintype==Domain2DhorizontalEnum){
    84                         femmodel->EsaGeodetic2D(U_radial,U_north,U_east,xx,yy);
     88                        femmodel->EsaGeodetic2D(U_radial,U_north,U_east,U_x,U_y,xx,yy);
     89                        InputUpdateFromVectorx(femmodel,U_x,EsaXmotionEnum,VertexSIdEnum);
     90                        InputUpdateFromVectorx(femmodel,U_y,EsaYmotionEnum,VertexSIdEnum);
    8591                }
    8692
    … …  
    102108                delete U_north;
    103109                delete U_east;
     110                delete U_x;
     111                delete U_y;
    104112                if(numoutputs){for(int i=0;i<numoutputs;i++){xDelete<char>(requested_outputs[i]);} xDelete<char*>(requested_outputs);}
    105113        }

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

    r22326 r22349  
    872872        EsaNmotionEnum,
    873873        EsaEmotionEnum,
     874        EsaXmotionEnum,
     875        EsaYmotionEnum,
     876        EsaHemisphereEnum,
    874877        EsaStrainratexxEnum,
    875878   EsaStrainratexyEnum,

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

    r22326 r22349  
    847847                case EsaNmotionEnum : return "EsaNmotion";
    848848                case EsaEmotionEnum : return "EsaEmotion";
     849                case EsaXmotionEnum : return "EsaXmotion";
     850                case EsaYmotionEnum : return "EsaYmotion";
     851                case EsaHemisphereEnum : return "EsaHemisphere";
    849852                case EsaStrainratexxEnum : return "EsaStrainratexx";
    850853                case EsaStrainratexyEnum : return "EsaStrainratexy";

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

    r22326 r22349  
    865865              else if (strcmp(name,"EsaNmotion")==0) return EsaNmotionEnum;
    866866              else if (strcmp(name,"EsaEmotion")==0) return EsaEmotionEnum;
     867              else if (strcmp(name,"EsaXmotion")==0) return EsaXmotionEnum;
     868              else if (strcmp(name,"EsaYmotion")==0) return EsaYmotionEnum;
     869              else if (strcmp(name,"EsaHemisphere")==0) return EsaHemisphereEnum;
    867870              else if (strcmp(name,"EsaStrainratexx")==0) return EsaStrainratexxEnum;
    868871              else if (strcmp(name,"EsaStrainratexy")==0) return EsaStrainratexyEnum;
    … …  
    872875              else if (strcmp(name,"EsaUElastic")==0) return EsaUElasticEnum;
    873876              else if (strcmp(name,"EsaHElastic")==0) return EsaHElasticEnum;
    874               else if (strcmp(name,"EsaTransitions")==0) return EsaTransitionsEnum;
    875               else if (strcmp(name,"EsaRequestedOutputs")==0) return EsaRequestedOutputsEnum;
    876               else if (strcmp(name,"EsaNumRequestedOutputs")==0) return EsaNumRequestedOutputsEnum;
    877877         else stage=8;
    878878   }
    879879   if(stage==8){
    880               if (strcmp(name,"TransientAmrFrequency")==0) return TransientAmrFrequencyEnum;
     880              if (strcmp(name,"EsaTransitions")==0) return EsaTransitionsEnum;
     881              else if (strcmp(name,"EsaRequestedOutputs")==0) return EsaRequestedOutputsEnum;
     882              else if (strcmp(name,"EsaNumRequestedOutputs")==0) return EsaNumRequestedOutputsEnum;
     883              else if (strcmp(name,"TransientAmrFrequency")==0) return TransientAmrFrequencyEnum;
    881884              else if (strcmp(name,"AmrType")==0) return AmrTypeEnum;
    882885              else if (strcmp(name,"AmrRestart")==0) return AmrRestartEnum;
    … …  
    995998              else if (strcmp(name,"EnthalpyAnalysis")==0) return EnthalpyAnalysisEnum;
    996999              else if (strcmp(name,"HydrologyShreveAnalysis")==0) return HydrologyShreveAnalysisEnum;
    997               else if (strcmp(name,"HydrologyDCInefficientAnalysis")==0) return HydrologyDCInefficientAnalysisEnum;
    998               else if (strcmp(name,"HydrologyDCEfficientAnalysis")==0) return HydrologyDCEfficientAnalysisEnum;
    999               else if (strcmp(name,"HydrologySommersAnalysis")==0) return HydrologySommersAnalysisEnum;
    10001000         else stage=9;
    10011001   }
    10021002   if(stage==9){
    1003               if (strcmp(name,"HydrologySolution")==0) return HydrologySolutionEnum;
     1003              if (strcmp(name,"HydrologyDCInefficientAnalysis")==0) return HydrologyDCInefficientAnalysisEnum;
     1004              else if (strcmp(name,"HydrologyDCEfficientAnalysis")==0) return HydrologyDCEfficientAnalysisEnum;
     1005              else if (strcmp(name,"HydrologySommersAnalysis")==0) return HydrologySommersAnalysisEnum;
     1006              else if (strcmp(name,"HydrologySolution")==0) return HydrologySolutionEnum;
    10041007              else if (strcmp(name,"MeltingAnalysis")==0) return MeltingAnalysisEnum;
    10051008              else if (strcmp(name,"MasstransportAnalysis")==0) return MasstransportAnalysisEnum;

  • issm/trunk-jpl/src/m/classes/esa.m

    r21620 r22349  
    1010                love_l         = 0; %ideam
    1111                degacc         = 0;
     12                hemisphere              = 0;
    1213                requested_outputs      = {};
    1314                transitions    = {};
    … …  
    2627                %numerical discretization accuracy
    2728                self.degacc=.01;
    28                
     29       
     30                %computational flags:
     31                self.hemisphere=0;
     32
    2933                %output default:
    3034                self.requested_outputs={'default'};
    … …  
    6670                        fielddisplay(self,'love_h','load Love number for radial displacement');
    6771                        fielddisplay(self,'love_l','load Love number for horizontal displacements');
     72                        fielddisplay(self,'hemisphere','North-south, East-west components of 2-D horiz displacement vector: -1 south, 1 north');
    6873                        fielddisplay(self,'degacc','accuracy (default .01 deg) for numerical discretization of the Green''s functions');
    6974                        fielddisplay(self,'transitions','indices into parts of the mesh that will be icecaps');
    … …  
    7580                        WriteData(fid,prefix,'object',self,'fieldname','love_h','format','DoubleMat','mattype',1);
    7681                        WriteData(fid,prefix,'object',self,'fieldname','love_l','format','DoubleMat','mattype',1);
     82                        WriteData(fid,prefix,'object',self,'fieldname','hemisphere','format','Integer');
    7783                        WriteData(fid,prefix,'object',self,'fieldname','degacc','format','Double');
    7884                        WriteData(fid,prefix,'object',self,'fieldname','transitions','format','MatArray');
    … …  
    9399                        writejs1Darray(fid,[modelname '.esa.love_h'],self.love_h);
    94100                        writejs1Darray(fid,[modelname '.esa.love_l'],self.love_l);
     101                        writejsdouble(fid,[modelname '.esa.hemisphere'],self.hemisphere);
    95102                        writejsdouble(fid,[modelname '.esa.degacc'],self.degacc);
    96103                        writejscellstring(fid,[modelname '.esa.requested_outputs'],self.requested_outputs);
Note: See TracChangeset for help on using the changeset viewer.