Changeset 19984


Ignore:
Timestamp:
01/22/16 09:56:50 (9 years ago)
Author:
Eric.Larour
Message:

CHG: new sea level rise solution. Not valdated yet!

This solution requires one Sealevelrise analysis (so added the corresponding EnumToAnalysis and analysis.h files + SealevelriseAnalysis.* files).
In terms of solution core: we have a new sealevelrise_core.cpp files + corresponding hook up in CorePointerFromSolutionEnum.
This core calls the new FemModel Sealevelrise module, which loops through the elements, hence the mods to Element.* along with all
derivatives classes, in particular Tria.
In terms of ModelProcessorx, we have a modified CreateElementsVerticesAndMaterials to include lat,long and radius, which also translsates
into modifications for the Vertex object. The VertexCoordinatesx module is also modified, to be able to retrieve x,y,z but also lat,long,r
from vertices.
Of course, this is all driven from matlab, where we have a new field in the model, the slr class.

Location:
issm/trunk-jpl/src
Files:
4 added
25 edited

Legend:

Unmodified
Added
Removed

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

    r19749 r19984  
    111111                                        ./shared/Numerics/NewtonSolveDnorm.cpp\
    112112                                        ./shared/Numerics/extrema.cpp\
     113                                        ./shared/Numerics/legendre.cpp\
    113114                                        ./shared/Numerics/XZvectorsToCoordinateSystem.cpp\
    114115                                        ./shared/Exceptions/Exceptions.cpp\
    … …  
    464465endif
    465466#}}}
     467#Slr sources  {{{
     468if SEALEVELRISE
     469issm_sources +=  ./cores/sealevelrise_core.cpp\
     470                                        ./analyses/SealevelriseAnalysis.cpp
     471endif
     472#}}}
    466473#Metis sources  {{{
    467474if METIS

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

    r19720 r19984  
    116116                case LsfReinitializationAnalysisEnum : return new LsfReinitializationAnalysis();
    117117                #endif
     118                #ifdef _HAVE_SEALEVELRISE_
     119                case SealevelriseAnalysisEnum : return new SealevelriseAnalysis();
     120                #endif
    118121                default : _error_("enum provided not supported ("<<EnumToStringx(analysis_enum)<<")");
    119122        }

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

    r19720 r19984  
    3232#include "./MasstransportAnalysis.h"
    3333#include "./SmbAnalysis.h"
     34#include "./SealevelriseAnalysis.h"
    3435#include "./MeltingAnalysis.h"
    3536#include "./MeshdeformationAnalysis.h"

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

    r19825 r19984  
    14641464bool       Element::IsIceInElement(){/*{{{*/
    14651465        return (this->inputs->Min(MaskIceLevelsetEnum)<0.);
     1466}
     1467/*}}}*/
     1468bool       Element::IsWaterInElement(){/*{{{*/
     1469        return (this->inputs->Max(MaskOceanLevelsetEnum)>0.);
    14661470}
    14671471/*}}}*/
    … …  
    15161520                                name==MaterialsRheologyBbarEnum ||
    15171521                                name==MaterialsRheologyNEnum ||
     1522                                name==SealevelriseSEnum ||
    15181523                                name==GiaWEnum ||
    15191524                                name==GiadWdtEnum ||

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

    r19764 r19984  
    119119                bool               IsFloating();
    120120                bool               IsIceInElement();
     121                bool               IsWaterInElement();
    121122                bool                 IsInput(int name);
    122123                void               LinearFloatingiceMeltingRate();
    … …  
    299300                virtual void       GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt,IssmDouble* x,IssmDouble* y)=0;
    300301                #endif
     302                #ifdef _HAVE_SEALEVELRISE_
     303                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 IssmDouble    OceanArea(void)=0;
     305                #endif
     306
    301307};
    302308#endif

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

    r19764 r19984  
    180180                void           GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt,IssmDouble* x,IssmDouble* y);
    181181                #endif
     182                #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                IssmDouble    OceanArea(void){_error_("not implemented yet!");};
     185                #endif
     186
    182187                /*}}}*/
    183188};

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

    r19764 r19984  
    166166#endif
    167167
     168#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                IssmDouble    OceanArea(void){_error_("not implemented yet!");};
     171#endif
     172
    168173#ifdef _HAVE_DAKOTA_
    169174                void        InputUpdateFromMatrixDakota(IssmDouble* matrix, int nows, int ncols, int name, int type){_error_("not implemented yet");};

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

    r19764 r19984  
    172172                void        GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt,IssmDouble* x,IssmDouble* y){_error_("not implemented yet");};
    173173#endif
     174#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                IssmDouble    OceanArea(void){_error_("not implemented yet!");};
     177#endif
    174178
    175179#ifdef _HAVE_DAKOTA_

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

    r19825 r19984  
    10341034        _assert_(x2*y3 - y2*x3 + x1*y2 - y1*x2 + x3*y1 - y3*x1>0);
    10351035        return (x2*y3 - y2*x3 + x1*y2 - y1*x2 + x3*y1 - y3*x1)/2;
     1036}
     1037/*}}}*/
     1038IssmDouble Tria::GetArea3D(void){/*{{{*/
     1039
     1040        IssmDouble xyz_list[NUMVERTICES][3];
     1041        IssmDouble x1,y1,z1,x2,y2,z2,x3,y3,z3;
     1042        IssmDouble detm1,detm2,detm3;
     1043
     1044        /*Get xyz list: */
     1045        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     1046        x1=xyz_list[0][0]; y1=xyz_list[0][1]; z1=xyz_list[0][2];
     1047        x2=xyz_list[1][0]; y2=xyz_list[1][1]; z2=xyz_list[1][2];
     1048        x3=xyz_list[2][0]; y3=xyz_list[2][1]; z3=xyz_list[2][2];
     1049
     1050        detm1=x1*y2 - x2*y1 - x1*y3 + x3*y1 + x2*y3 - x3*y2;
     1051        detm2=y1*z2 - y2*z1 - y1*z3 + y3*z1 + y2*z3 - y3*z2;
     1052        detm3=x2*z1 - x1*z2 + x1*z3 - x3*z1 - x2*z3 + x3*z2;
     1053
     1054        return sqrt(pow(detm1,2) + pow(detm2,2) + pow(detm3,2))/2;
    10361055}
    10371056/*}}}*/
    … …  
    34833502#endif
    34843503
     3504#ifdef _HAVE_SEALEVELRISE_
     3505void    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){ /*{{{*/
     3506
     3507        /*diverse:*/
     3508        int gsize;
     3509        bool spherical=true;
     3510        IssmDouble x0,y0,z0,a;
     3511        IssmDouble xyz_list[NUMVERTICES][3];
     3512        IssmDouble area;
     3513        IssmDouble I;  //change in ice thickness or water level(Farrel and Clarke, Equ. 4)
     3514        IssmDouble Me;
     3515        IssmDouble rho;
     3516       
     3517        /*love numbers:*/
     3518        IssmDouble* love_h=NULL;
     3519        IssmDouble* love_k=NULL;
     3520        IssmDouble* deltalove=NULL;
     3521        int nl;
     3522
     3523        /*ice properties: */
     3524        IssmDouble rho_ice,rho_water,rho_earth;
     3525
     3526        /*other constants: */
     3527        //IssmDouble g;
     3528
     3529        /*Solution outputs: */
     3530        Vector<IssmDouble>* pSolution=NULL;
     3531       
     3532        /*Initialize eustatic component: do not skip this step :):*/
     3533        IssmDouble eustatic = 0;
     3534
     3535        /*Computational flags:*/
     3536        bool computerigid = true;
     3537        bool computeelastic= true;
     3538        bool computeeustatic = true;
     3539       
     3540        /*recover material parameters: */
     3541        rho_ice=matpar->GetMaterialParameter(MaterialsRhoIceEnum);
     3542        rho_water=matpar->GetMaterialParameter(MaterialsRhoSeawaterEnum);
     3543        rho_earth=matpar->GetMaterialParameter(MaterialsEarthDensityEnum);
     3544
     3545        /*recover love numbers and computational flags: */
     3546        this->parameters->FindParam(&love_h,&nl,SealevelriseLoveHEnum);
     3547        this->parameters->FindParam(&love_k,&nl,SealevelriseLoveKEnum);
     3548        this->parameters->FindParam(&computerigid,SealevelriseRigidEnum);
     3549        this->parameters->FindParam(&computeelastic,SealevelriseElasticEnum);
     3550        this->parameters->FindParam(&computeeustatic,SealevelriseEustaticEnum);
     3551
     3552        /*how many dofs are we working with here? */
     3553        this->parameters->FindParam(&gsize,MeshNumberofverticesEnum);
     3554
     3555        /*retrieve coordinates: */
     3556        ::GetVerticesCoordinates(&xyz_list[0][0],this->vertices,NUMVERTICES);
     3557
     3558        /* Where is the centroid of this element?. To avoid issues with lat,long
     3559         * being between [0,180] and [0 360], and issues with jumps at the central
     3560         * meridian and poles, we do everything in cartesian coordinate system: */
     3561        x0=(xyz_list[0][0]+xyz_list[1][0]+xyz_list[2][0])/3;
     3562        y0=(xyz_list[0][1]+xyz_list[1][1]+xyz_list[2][1])/3;
     3563        z0=(xyz_list[0][2]+xyz_list[1][2]+xyz_list[2][2])/3;
     3564
     3565        /*radius at this location?:*/
     3566        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.
     3567
     3568        /*Compute mass of the earth:*/
     3569        Me= rho_earth*4/3*PI*pow(a,3);
     3570
     3571        /*Compute area of element:*/
     3572        area=GetArea3D();
     3573
     3574        /*Compute ice thickness or sea level change: */
     3575        if(IsWaterInElement()){
     3576
     3577                IssmDouble phi_I_I_O,phi_O_O_O;
     3578               
     3579                /*From Sg_old, recover water sea level rise:*/
     3580                I=0; for(int i=0;i<NUMVERTICES;i++) I+=Sg_old[this->vertices[i]->Sid()]/NUMVERTICES;
     3581                rho=rho_water;
     3582                pSolution=pSgo;
     3583
     3584                /*Compute eustatic component: */
     3585                if(computeeustatic){
     3586                        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;
     3587                        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;
     3588                        eustatic -= (phi_I_I_O + phi_O_O_O); //Watch out the sign "-"!
     3589                }
     3590
     3591        }
     3592        else{
     3593                Input*  deltathickness_input=inputs->GetInput(SealevelriseDeltathicknessEnum);
     3594                if (!deltathickness_input)_error_("delta thickness input needed to compute sea level rise!");
     3595                deltathickness_input->GetInputAverage(&I);
     3596                rho=rho_ice;
     3597                pSolution=pSgi;
     3598
     3599                /*Compute eustatic compoent:*/
     3600                if(computeeustatic)eustatic -= rho_ice*area*I/(oceanarea*rho_water); //Watch out the sign "-"!
     3601        }
     3602
     3603                               
     3604        /*Speed up of love number comutation: */
     3605        if(computeelastic){
     3606                deltalove=xNew<IssmDouble>(nl);
     3607                for (int n=0;n<nl;n++) deltalove[n] = (love_k[n]-love_k[nl-1]-love_h[n]+love_h[nl-1]);
     3608        }
     3609
     3610        if(computeelastic | computerigid){
     3611                for(int i=0;i<gsize;i++){
     3612
     3613                        IssmDouble alpha;
     3614                        IssmDouble G_rigid=0;  //do not remove =0!
     3615                        IssmDouble G_elastic=0;  //do not remove =0!
     3616                        IssmDouble Pn1,Pn2; //first two legendre polynomials
     3617                        IssmDouble cosalpha;
     3618
     3619                        /*Compute alpha angle between centroid and current vertex and cosine of this angle: */
     3620                        alpha = acos(
     3621                                        (x[i]*x0+y[i]*y0+z[i]*z0)  //scalar product  of two position vectors
     3622                                        / sqrt(pow(x0,2.0)+pow(y0,2.0)+pow(z0,2.0))  //norm of first position vector
     3623                                        / sqrt(pow(x[i],2.0)+pow(y[i],2.0)+pow(z[i],2.0))  //norm of second position vector
     3624                                        );
     3625                        cosalpha=cos(alpha); //compute this to speed up the elastic component.
     3626
     3627                        //Rigid earth gravitational perturbation:
     3628                        if(computerigid)G_rigid=1.0/2.0/sin(alpha/2.0);
     3629
     3630                        //Elastic component  (from Eq 17 in Adhikari et al, GMD 2015)
     3631                        if(computeelastic){
     3632                                G_elastic = (love_k[nl-1]-love_h[nl-1])/2.0/sin(alpha/2.0);
     3633                                for(int n=0;n<nl;n++) G_elastic += deltalove[n]*legendre(Pn1,Pn2,n,cosalpha);
     3634                        }
     3635
     3636                        /*Add all components to the pSgi or pSgo solution vectors:*/
     3637                        pSolution->SetValue(i,rho*a/Me*I*area*(G_rigid+G_elastic),ADD_VAL);
     3638                }
     3639        }
     3640       
     3641        /*Assign output pointer:*/
     3642        *peustatic=eustatic;
     3643
     3644        /*Free ressources:*/
     3645        xDelete<IssmDouble>(love_h);
     3646        xDelete<IssmDouble>(love_k);
     3647        xDelete<IssmDouble>(deltalove);
     3648        return;
     3649}
     3650/*}}}*/
     3651IssmDouble    Tria::OceanArea(void){ /*{{{*/
     3652
     3653        if(IsWaterInElement()) return GetArea3D();
     3654        else return 0;
     3655
     3656}
     3657/*}}}*/
     3658#endif
     3659
     3660
    34853661#ifdef _HAVE_DAKOTA_
    34863662void       Tria::InputUpdateFromMatrixDakota(IssmDouble* matrix, int nrows, int ncols, int name, int type){/*{{{*/

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

    r19764 r19984  
    143143                void   GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt,IssmDouble* x,IssmDouble* y);
    144144                #endif
     145                #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);
     147                IssmDouble OceanArea(void);
     148                #endif
    145149                /*}}}*/
    146150                /*Tria specific routines:{{{*/
    … …  
    148152                void           AddInput(int input_enum, IssmDouble* values, int interpolation_enum);
    149153                IssmDouble     GetArea(void);
    150                 IssmDouble      GetAreaIce(void);
     154                IssmDouble         GetArea3D(void);
     155                IssmDouble         GetAreaIce(void);
    151156                void           GetAreaCoordinates(IssmDouble *area_coordinates,IssmDouble* xyz_zero,IssmDouble* xyz_list,int numpoints);
    152157                void            GetLevelsetIntersection(int** pindices, int* pnumiceverts, IssmDouble* fraction, int levelset_enum, IssmDouble level);

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

    r19803 r19984  
    534534                        break;
    535535               
     536                case SealevelriseSolutionEnum:
     537                        analyses_temp[numanalyses++]=SealevelriseAnalysisEnum;
     538                        break;
     539
    536540                case SmbSolutionEnum:
    537541                        analyses_temp[numanalyses++]=SmbAnalysisEnum;
    … …  
    22012205/*}}}*/
    22022206#endif
     2207#ifdef _HAVE_SEALEVELRISE_
     2208void 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) { /*{{{*/
     2209
     2210        /*serialized vectors:*/
     2211        IssmDouble* Sg_old=NULL;
     2212        IssmDouble* Sgi_old=NULL;
     2213        IssmDouble* Sgo_old=NULL;
     2214        IssmDouble  eustatic=0;
     2215        IssmDouble  eustatic_cpu=0;
     2216        IssmDouble  eustatic_cpu_e=0;
     2217        IssmDouble  oceanarea=0;
     2218        IssmDouble  oceanarea_cpu=0;
     2219        int         ns;
     2220       
     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
     2227        /*Go through elements, and add contribution from each element to the deflection vector wg:*/
     2228        ns = elements->Size();
     2229
     2230        /*First, figure out the area of the ocean, which is needed to compute the eustatic component: */
     2231        for(int i=0;i<ns;i++){
     2232                Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
     2233                oceanarea_cpu += element->OceanArea();
     2234        }
     2235        ISSM_MPI_Reduce (&oceanarea_cpu,&oceanarea,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
     2236        ISSM_MPI_Bcast(&oceanarea,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
     2237
     2238        /*Call the sea level rise core: */
     2239        for(int i=0;i<ns;i++){
     2240               
     2241                if(IssmComm::GetRank()==0)if(VerboseConvergence())if(i%1000)_printf_("\r" << "      convolution progress: " << (float)i/(float)ns*100 << "\%");
     2242       
     2243                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);
     2245                eustatic_cpu+=eustatic_cpu_e;
     2246        }
     2247        if(IssmComm::GetRank()==0)if(VerboseConvergence())_printf_("\n");
     2248
     2249        /*Sum all eustatic components from all cpus:*/
     2250        ISSM_MPI_Reduce (&eustatic_cpu,&eustatic,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
     2251        ISSM_MPI_Bcast(&eustatic,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
     2252
     2253        /*Assign output pointers:*/
     2254        *peustatic=eustatic;
     2255
     2256        /*Free ressources:*/
     2257        xDelete<IssmDouble>(Sg_old);
     2258        xDelete<IssmDouble>(Sgi_old);
     2259        xDelete<IssmDouble>(Sgo_old);
     2260}
     2261/*}}}*/
     2262#endif
     2263
    22032264void FemModel::HydrologyEPLupdateDomainx(IssmDouble* pEplcount){ /*{{{*/
    22042265

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

    r19803 r19984  
    112112                void Deflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt, IssmDouble* x, IssmDouble* y);
    113113                #endif
     114                #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);
     116                #endif
    114117                void TimeAdaptx(IssmDouble* pdt);
    115118                void UpdateConstraintsx(void);

  • issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp

    r19720 r19984  
    5858        mantle_shear_modulus=0;
    5959        mantle_density=0;
     60       
     61        earth_density=0;
    6062
    6163        poisson=0;
    … …  
    160162                        iomodel->Constant(&this->mantle_shear_modulus,MaterialsMantleShearModulusEnum);
    161163                        iomodel->Constant(&this->mantle_density,MaterialsMantleDensityEnum);
     164
     165                        /*slr:*/
     166                        iomodel->Constant(&this->earth_density,MaterialsEarthDensityEnum);
     167
    162168                        break;
    163169                default:
    … …  
    255261        matpar->mantle_shear_modulus=this->mantle_shear_modulus;
    256262        matpar->mantle_density=this->mantle_density;
     263       
     264        matpar->earth_density=this->earth_density;
    257265
    258266        return matpar;
    … …  
    302310        MARSHALLING(mantle_shear_modulus);
    303311        MARSHALLING(mantle_density);
     312       
     313        //slr:
     314        MARSHALLING(earth_density);
    304315
    305316        //Sea ice:
    … …  
    516527                case MaterialsMantleDensityEnum:             return this->mantle_density;
    517528                case MaterialsMantleShearModulusEnum:        return this->mantle_shear_modulus;
     529                case MaterialsEarthDensityEnum:              return this->earth_density;
    518530                default: _error_("Enum "<<EnumToStringx(enum_in)<<" not supported yet");
    519531        }

  • issm/trunk-jpl/src/c/classes/Materials/Matpar.h

    r19554 r19984  
    5858                IssmDouble mantle_shear_modulus;
    5959                IssmDouble mantle_density;
     60
     61                /*slr:*/
     62                IssmDouble earth_density;
    6063
    6164                /*Sea ice*/

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

    r19199 r19984  
    3636                        _assert_(iomodel->Data(BaseEnum) && iomodel->Data(ThicknessEnum));
    3737                        this->sigma = (iomodel->Data(MeshZEnum)[i]-iomodel->Data(BaseEnum)[i])/(iomodel->Data(ThicknessEnum)[i]);
     38                        break;
     39                case Domain3DsurfaceEnum:
     40                        this->latitute     = iomodel->Data(MeshLatEnum)[i];
     41                        this->longitude    = iomodel->Data(MeshLongEnum)[i];
     42                        this->R            = iomodel->Data(MeshREnum)[i];
    3843                        break;
    3944                case Domain2DhorizontalEnum:
    … …  
    121126IssmDouble Vertex::GetZ(){/*{{{*/
    122127        return this->z;
     128}
     129/*}}}*/
     130IssmDouble Vertex::GetLatitude(){/*{{{*/
     131        return this->latitute;
     132}
     133/*}}}*/
     134IssmDouble Vertex::GetLongitude(){/*{{{*/
     135        return this->longitude;
     136}
     137/*}}}*/
     138IssmDouble Vertex::GetRadius(){/*{{{*/
     139        return this->R;
    123140}
    124141/*}}}*/
    … …  
    239256}
    240257/*}}}*/
    241 void       Vertex::VertexCoordinates(Vector<IssmDouble>* vx,Vector<IssmDouble>* vy,Vector<IssmDouble>* vz){/*{{{*/
     258void       Vertex::VertexCoordinates(Vector<IssmDouble>* vx,Vector<IssmDouble>* vy,Vector<IssmDouble>* vz, bool spherical){/*{{{*/
    242259
    243260        if (this->clone==true) return;
    244261
    245         vx->SetValue(this->sid,this->x,INS_VAL);
    246         vy->SetValue(this->sid,this->y,INS_VAL);
    247         vz->SetValue(this->sid,this->z,INS_VAL);
     262        if(!spherical){
     263                vx->SetValue(this->sid,this->x,INS_VAL);
     264                vy->SetValue(this->sid,this->y,INS_VAL);
     265                vz->SetValue(this->sid,this->z,INS_VAL);
     266        }
     267        else{
     268                vx->SetValue(this->sid,this->latitute,INS_VAL);
     269                vy->SetValue(this->sid,this->longitude,INS_VAL);
     270                vz->SetValue(this->sid,this->R,INS_VAL);
     271        }
    248272
    249273        return;
    … …  
    252276
    253277/*Methods relating to Vertex, but not internal methods: */
    254 void GetVerticesCoordinates(IssmDouble* xyz,Vertex** vertices, int numvertices){ /*{{{*/
     278void GetVerticesCoordinates(IssmDouble* xyz,Vertex** vertices, int numvertices,bool spherical){ /*{{{*/
    255279
    256280        _assert_(vertices);
    257281        _assert_(xyz);
    258282
    259         for(int i=0;i<numvertices;i++) {
    260                 xyz[i*3+0]=vertices[i]->GetX();
    261                 xyz[i*3+1]=vertices[i]->GetY();
    262                 xyz[i*3+2]=vertices[i]->GetZ();
     283        if(!spherical){
     284                for(int i=0;i<numvertices;i++) {
     285                        xyz[i*3+0]=vertices[i]->GetX();
     286                        xyz[i*3+1]=vertices[i]->GetY();
     287                        xyz[i*3+2]=vertices[i]->GetZ();
     288                }
     289        }
     290        else{
     291                for(int i=0;i<numvertices;i++) {
     292                        xyz[i*3+0]=vertices[i]->GetLatitude();
     293                        xyz[i*3+1]=vertices[i]->GetLongitude();
     294                        xyz[i*3+2]=vertices[i]->GetRadius();
     295                }
    263296        }
    264297}/*}}}*/

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

    r19199 r19984  
    2828                IssmDouble y;
    2929                IssmDouble z;
     30                IssmDouble latitute;
     31                IssmDouble longitude;
     32                IssmDouble R;
    3033                IssmDouble sigma;        //sigma coordinate: (z-bed)/thickness
    3134                int        connectivity; //number of vertices connected to this vertex
    … …  
    5255                IssmDouble GetY(void);
    5356                IssmDouble GetZ(void);
     57                IssmDouble GetLatitude(void);
     58                IssmDouble GetLongitude(void);
     59                IssmDouble GetRadius(void);
    5460                void       UpdatePosition(Vector<IssmDouble>* vx,Vector<IssmDouble>* vy,Vector<IssmDouble>* vz,Parameters* parameters,IssmDouble* thickness,IssmDouble* bed);
    5561                void       DistributePids(int* ppidcount);
    … …  
    5965                void       SetClone(int* minranks);
    6066                void       ToXYZ(Matrix<IssmDouble>* matrix);
    61                 void       VertexCoordinates(Vector<IssmDouble>* vx,Vector<IssmDouble>* vy,Vector<IssmDouble>* vz);
     67                void       VertexCoordinates(Vector<IssmDouble>* vx,Vector<IssmDouble>* vy,Vector<IssmDouble>* vz,bool spherical=false);
    6268};
    6369
    6470/*Methods relating to Vertex object: */
    65 void GetVerticesCoordinates(IssmDouble* xyz,Vertex** vertices, int numvertices);
     71void GetVerticesCoordinates(IssmDouble* xyz,Vertex** vertices, int numvertices,bool spherical=false);
    6672
    6773#endif  /* _VERTEX_H */

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

    r19087 r19984  
    5959                        solutioncore=&masstransport_core;
    6060                        break;
     61                case SealevelriseSolutionEnum:
     62                        solutioncore=&sealevelrise_core;
     63                        break;
    6164                case GiaSolutionEnum:
    6265                        #if _HAVE_GIA_

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

    r19527 r19984  
    4646void dummy_core(FemModel* femmodel);
    4747void gia_core(FemModel* femmodel);
     48void sealevelrise_core(FemModel* femmodel);
    4849void smb_core(FemModel* femmodel);
    4950void damage_core(FemModel* femmodel);

  • issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp

    r19087 r19984  
    9595        /*Fetch data:*/
    9696        iomodel->FetchData(6,MeshXEnum,MeshYEnum,MeshZEnum,BaseEnum,ThicknessEnum,MaskIceLevelsetEnum);
     97        if (iomodel->domaintype == Domain3DsurfaceEnum) iomodel->FetchData(3,MeshLatEnum,MeshLongEnum,MeshREnum);
     98       
    9799        CreateNumberNodeToElementConnectivity(iomodel);
    98100
    … …  
    103105        /*Free data: */
    104106        iomodel->DeleteData(6,MeshXEnum,MeshYEnum,MeshZEnum,BaseEnum,ThicknessEnum,MaskIceLevelsetEnum);
     107        if (iomodel->domaintype == Domain3DsurfaceEnum) iomodel->DeleteData(3,MeshLatEnum,MeshLongEnum,MeshREnum);
    105108}

  • issm/trunk-jpl/src/c/modules/VertexCoordinatesx/VertexCoordinatesx.cpp

    r14999 r19984  
    99#include "../../toolkits/toolkits.h"
    1010
    11 void VertexCoordinatesx( IssmDouble** px, IssmDouble** py, IssmDouble** pz, Vertices* vertices) {
     11void VertexCoordinatesx( IssmDouble** px, IssmDouble** py, IssmDouble** pz, Vertices* vertices,bool spherical) {
    1212
    1313        /*output: */
    … …  
    3434        for(i=0;i<vertices->Size();i++){
    3535                Vertex* vertex=(Vertex*)vertices->GetObjectByOffset(i);
    36                 vertex->VertexCoordinates(vx,vy,vz);
     36                vertex->VertexCoordinates(vx,vy,vz,spherical);
    3737        }
    3838

  • issm/trunk-jpl/src/c/modules/VertexCoordinatesx/VertexCoordinatesx.h

    r15000 r19984  
    88
    99/* local prototypes: */
    10 void VertexCoordinatesx( IssmDouble** px, IssmDouble** py, IssmDouble** pz, Vertices* vertices);
     10void VertexCoordinatesx( IssmDouble** px, IssmDouble** py, IssmDouble** pz, Vertices* vertices,bool spherical=false);
    1111
    1212#endif  /* _VERTEX_COORDINATESX_H */

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

    r19958 r19984  
    6464                end % }}}
    6565                function marshall(self,md,fid) % {{{
    66                         WriteData(fid,'object',self,'fieldname','groundedice_levelset','format','DoubleMat','mattype',1);
    67                         WriteData(fid,'object',self,'fieldname','ice_levelset','format','DoubleMat','mattype',1);
    68                         WriteData(fid,'object',self,'fieldname','ocean_levelset','format','DoubleMat','mattype',1);
     66                        WriteData(fid,'object',self,'class','mask','fieldname','groundedice_levelset','format','DoubleMat','mattype',1);
     67                        WriteData(fid,'object',self,'class','mask','fieldname','ice_levelset','format','DoubleMat','mattype',1);
     68                        WriteData(fid,'object',self,'class','mask','fieldname','ocean_levelset','format','DoubleMat','mattype',1);
    6969
    7070                        % get mask of vertices of elements with ice

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

    r19958 r19984  
    2828                mantle_density             = 0.;
    2929
     30                %slr
     31                earth_density              = 0;
     32
    3033        end
    3134        methods
    … …  
    100103                        self.mantle_density             = 3.34;       % (g/cm^-3)
    101104
     105                        %SLR
     106                        self.earth_density= 5512;  % average density of the Earth, (kg/m^3)
     107
    102108                end % }}}
    103109                function md = checkconsistency(self,md,solution,analyses) % {{{
    … …  
    115121                                md = checkfield(md,'fieldname','materials.mantle_shear_modulus','>',0,'numel',1);
    116122                                md = checkfield(md,'fieldname','materials.mantle_density','>',0,'numel',1);
     123                        end
     124                        if ismember(SealevelriseAnalysisEnum(),analyses),
     125                                md = checkfield(md,'fieldname','materials.earth_density','>',0,'numel',1);
    117126                        end
    118127
    … …  
    140149                        fielddisplay(self,'mantle_shear_modulus','Mantle shear modulus [Pa]');
    141150                        fielddisplay(self,'mantle_density','Mantle density [g/cm^-3]');
     151                        fielddisplay(self,'earth_density','Mantle density [kg/m^-3]');
    142152                end % }}}
    143153                function marshall(self,md,fid) % {{{
    … …  
    163173                        WriteData(fid,'object',self,'class','materials','fieldname','mantle_shear_modulus','format','Double');
    164174                        WriteData(fid,'object',self,'class','materials','fieldname','mantle_density','format','Double','scale',10^3);
     175                        WriteData(fid,'object',self,'class','materials','fieldname','earth_density','format','Double');
    165176                end % }}}
    166177                function savemodeljs(self,fid,modelname) % {{{
    … …  
    186197                        writejsdouble(fid,[modelname '.materials.mantle_shear_modulus'],self.mantle_shear_modulus);
    187198                        writejsdouble(fid,[modelname '.materials.mantle_density'],self.mantle_density);
     199                        writejsdouble(fid,[modelname '.materials.earth_density'],self.earth_density);
    188200
    189201                end % }}}

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

    r19957 r19984  
    2222                initialization   = 0;
    2323                rifts            = 0;
     24                slr              = 0;
    2425
    2526                debug            = 0;
    … …  
    10641065                        md.friction         = friction();
    10651066                        md.rifts            = rifts();
     1067                        md.slr              = slr();
    10661068                        md.timestepping     = timestepping();
    10671069                        md.groundingline    = groundingline();
    … …  
    12381240                        disp(sprintf('%19s: %-22s -- %s','initialization'  ,['[1x1 ' class(self.initialization) ']'],'initial guess/state'));
    12391241                        disp(sprintf('%19s: %-22s -- %s','rifts'           ,['[1x1 ' class(self.rifts) ']'],'rifts properties'));
     1242                        disp(sprintf('%19s: %-22s -- %s','slr'             ,['[1x1 ' class(self.slr) ']'],'slr forcings'));
    12401243                        disp(sprintf('%19s: %-22s -- %s','debug'           ,['[1x1 ' class(self.debug) ']'],'debugging tools (valgrind, gprof)'));
    12411244                        disp(sprintf('%19s: %-22s -- %s','verbose'         ,['[1x1 ' class(self.verbose) ']'],'verbosity level in solve'));

  • issm/trunk-jpl/src/m/consistency/ismodelselfconsistent.m

    r19744 r19984  
    8585                analyses=[FlaimAnalysisEnum()];
    8686
     87        case SealevelriseSolutionEnum(),
     88                analyses=[SealevelriseAnalysisEnum()];
     89
    8790        case HydrologySolutionEnum(),
    8891                analyses=[L2ProjectionBaseAnalysisEnum();HydrologyShreveAnalysisEnum();HydrologyDCInefficientAnalysisEnum();HydrologyDCEfficientAnalysisEnum()];
Note: See TracChangeset for help on using the changeset viewer.