Changeset 24924


Ignore:
Timestamp:
05/29/20 19:08:11 (5 years ago)
Author:
Eric.Larour
Message:

CHG: further optimization of sea level core by adding a geometry core
that precomputes indices for each centroid and every vertex to index into
the gridded Green function kernels. These indices are the most computational
aspects of the sea level solver.

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

Legend:

Unmodified
Added
Removed

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

    r24915 r24924  
    378378                virtual void          SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old,IssmDouble eartharea)=0;
    379379                virtual void          SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea)=0;
     380                virtual void          SealevelriseGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius)=0;
    380381                virtual void          SealevelriseNonEustatic(IssmDouble* Sgo, IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea)=0;
    381382                virtual void          SealevelriseGeodetic(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea,int horiz)=0;

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

    r24915 r24924  
    215215                IssmDouble    OceanAverage(IssmDouble* Sg){_error_("not implemented yet!");};
    216216                void    SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old,IssmDouble eartharea){_error_("not implemented yet!");};
     217                void    SealevelriseGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius){_error_("not implemented yet!");};
    217218                void    SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){_error_("not implemented yet!");};
    218219                void    SealevelriseNonEustatic(IssmDouble* Sgo, IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea){_error_("not implemented yet!");};

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

    r24915 r24924  
    173173#ifdef _HAVE_SEALEVELRISE_
    174174                void    SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old,IssmDouble eartharea){_error_("not implemented yet!");};
     175                void    SealevelriseGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius){_error_("not implemented yet!");};
    175176                void    SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){_error_("not implemented yet!");};
    176177                void    SealevelriseNonEustatic(IssmDouble* Sgo,IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea){_error_("not implemented yet!");};

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

    r24915 r24924  
    179179#ifdef _HAVE_SEALEVELRISE_
    180180                void    SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old,IssmDouble eartharea){_error_("not implemented yet!");};
     181                void    SealevelriseGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius){_error_("not implemented yet!");};
    181182                void    SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){_error_("not implemented yet!");};
    182183                void    SealevelriseNonEustatic(IssmDouble* Sgo, IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea){_error_("not implemented yet!");};

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

    r24921 r24924  
    55645564        return;
    55655565}/*}}}*/
    5566 void    Tria::SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){ /*{{{*/
    5567 
    5568         /*Computational flags:*/
    5569         int bp_compute_fingerprints= 0;
    5570 
    5571         /*some paramters first: */
    5572         this->parameters->FindParam(&bp_compute_fingerprints,DslComputeFingerprintsEnum);
    5573 
    5574         if(!IsOceanInElement()){
    5575                 /*ok, there is ocean in this element, we should compute eustatic loads for the ocean if we have requested
    5576                  *bottom pressure fingerprints:*/
    5577                 if(bp_compute_fingerprints)this->SealevelriseEustaticBottomPressure(Sgi,peustatic,latitude,longitude,radius,oceanarea,eartharea);
    5578         }
    5579         //if(!IsIceInElement()){
    5580                 /*there is ice in this eleemnt, let's compute the eustatic response for ice changes:*/
    5581                 this->SealevelriseEustaticIce(Sgi,peustatic,latitude,longitude,radius,oceanarea,eartharea);
    5582         //}
    5583 
    5584 }
    5585 /*}}}*/
    5586 void    Tria::SealevelriseEustaticIce(IssmDouble* Sgi, IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){ /*{{{*/
     5566void    Tria::SealevelriseGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius){ /*{{{*/
    55875567
    55885568        /*diverse:*/
    … …  
    55955575        IssmDouble late,longe,re;
    55965576        IssmDouble lati,longi,ri;
    5597         bool notfullygrounded=false;
    55985577
    55995578        /*elastic green function:*/
    5600         int         index;
    5601         IssmDouble* G_elastic_precomputed=NULL;
    5602         IssmDouble* G_rigid_precomputed=NULL;
     5579        IssmDouble* indices=NULL;
    56035580        int         M;
    5604 
    5605         /*ice properties: */
    5606         IssmDouble rho_ice,rho_water,rho_earth;
    5607 
    5608         /*constants:*/
    5609         IssmDouble constant=0;
    5610 
    5611         /*Initialize eustatic component: do not skip this step :):*/
    5612         IssmDouble eustatic = 0.;
     5581        IssmDouble* dummypointer=NULL;
    56135582
    56145583        /*Computational flags:*/
    56155584        bool computerigid = true;
    5616         bool computeelastic= true;
    5617         bool scaleoceanarea= false;
    5618 
    5619         /*early return if we are not on an ice cap:*/
    5620         if(!IsIceOnlyInElement()){
    5621                 constant=0; this->AddInput2(SealevelEustaticMaskEnum,&constant,P0Enum);
    5622                 *peustatic=0; //do not forget to assign this pointer, otherwise, global eustatic will be garbage!
    5623                 return;
    5624         }
    5625 
    5626         /*early return if we are fully floating:*/
    5627         Input2* gr_input=this->GetInput2(MaskOceanLevelsetEnum); _assert_(gr_input);
    5628         if (gr_input->GetInputMax()<=0){
    5629                 constant=0; this->AddInput2(SealevelEustaticMaskEnum,&constant,P0Enum);
    5630                 *peustatic=0; //do not forget to assign this pointer, otherwise, global eustatic will be garbage!
    5631                 return;
    5632         }
    5633 
    5634         /*If we are here, we are on ice that is fully grounded or half-way to floating:*/
    5635         if ((gr_input->GetInputMin())<0){
    5636                 notfullygrounded=true; //used later on.
    5637         }
    5638 
    5639         /*Inform mask: */
    5640         constant=1; this->AddInput2(SealevelEustaticMaskEnum,&constant,P0Enum);
    5641 
    5642         /*recover material parameters: */
    5643         rho_ice=FindParam(MaterialsRhoIceEnum);
    5644         rho_water=FindParam(MaterialsRhoSeawaterEnum);
    5645         rho_earth=FindParam(MaterialsEarthDensityEnum);
    5646 
    5647         /*recover love numbers and computational flags: */
     5585
     5586        /*recover computational flags: */
    56485587        this->parameters->FindParam(&computerigid,SealevelriseRigidEnum);
    5649         this->parameters->FindParam(&computeelastic,SealevelriseElasticEnum);
    5650         this->parameters->FindParam(&scaleoceanarea,SealevelriseOceanAreaScalingEnum);
    5651 
    5652         /*recover elastic green function:*/
    5653         DoubleVecParam* parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelriseGRigidEnum)); _assert_(parameter);
    5654         parameter->GetParameterValueByPointer(&G_rigid_precomputed,&M);
    5655 
    5656         if(computeelastic){
    5657                 parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelriseGElasticEnum)); _assert_(parameter);
    5658                 parameter->GetParameterValueByPointer(&G_elastic_precomputed,&M);
    5659         }
     5588        if(!computerigid) return; //we are running eustatic solution only.
    56605589
    56615590        /*how many dofs are we working with here? */
    56625591        this->parameters->FindParam(&gsize,MeshNumberofverticesEnum);
     5592
     5593        /*recover size of indexed tables:*/
     5594        DoubleVecParam* parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelriseGRigidEnum)); _assert_(parameter);
     5595        parameter->GetParameterValueByPointer(&dummypointer,&M);
     5596
     5597        /*allocate indices:*/
     5598        indices=xNew<IssmDouble>(gsize);
    56635599
    56645600        /* Where is the centroid of this element?:{{{*/
    … …  
    57005636        longe=longe/180*PI;
    57015637        /*}}}*/
     5638
     5639        for(int i=0;i<gsize;i++){
     5640                IssmDouble alpha;
     5641                IssmDouble delPhi,delLambda;
     5642
     5643                /*Compute alpha angle between centroid and current vertex and index into precomputed tables: */
     5644                lati=latitude[i]/180*PI; longi=longitude[i]/180*PI;
     5645                delPhi=fabs(lati-late); delLambda=fabs(longi-longe);
     5646                alpha=2.*asin(sqrt(pow(sin(delPhi/2),2.0)+cos(lati)*cos(late)*pow(sin(delLambda/2),2)));
     5647                indices[i]=reCast<int,IssmDouble>(alpha/PI*reCast<IssmDouble,int>(M-1));
     5648        }
     5649
     5650        /*Add in inputs:*/
     5651    this->inputs2->SetArrayInput(SealevelriseIndicesEnum,this->lid,indices,gsize);
     5652
     5653        return;
     5654}
     5655/*}}}*/
     5656void    Tria::SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){ /*{{{*/
     5657
     5658        /*Computational flags:*/
     5659        int bp_compute_fingerprints= 0;
     5660
     5661        /*some paramters first: */
     5662        this->parameters->FindParam(&bp_compute_fingerprints,DslComputeFingerprintsEnum);
     5663
     5664        if(!IsOceanInElement()){
     5665                /*ok, there is ocean in this element, we should compute eustatic loads for the ocean if we have requested
     5666                 *bottom pressure fingerprints:*/
     5667                if(bp_compute_fingerprints)this->SealevelriseEustaticBottomPressure(Sgi,peustatic,latitude,longitude,radius,oceanarea,eartharea);
     5668        }
     5669        //if(!IsIceInElement()){
     5670                /*there is ice in this eleemnt, let's compute the eustatic response for ice changes:*/
     5671                this->SealevelriseEustaticIce(Sgi,peustatic,latitude,longitude,radius,oceanarea,eartharea);
     5672        //}
     5673
     5674}
     5675/*}}}*/
     5676void    Tria::SealevelriseEustaticIce(IssmDouble* Sgi, IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){ /*{{{*/
     5677
     5678        /*diverse:*/
     5679        int gsize,dummy;
     5680        bool spherical=true;
     5681        IssmDouble llr_list[NUMVERTICES][3];
     5682        IssmDouble area;
     5683        IssmDouble I;  //change in ice thickness or water level(Farrel and Clarke, Equ. 4)
     5684        IssmDouble rho;
     5685        IssmDouble late,longe,re;
     5686        IssmDouble lati,longi,ri;
     5687        bool notfullygrounded=false;
     5688
     5689        /*elastic green function:*/
     5690        int         index;
     5691        IssmDouble* indices=NULL;
     5692        IssmDouble* G_elastic_precomputed=NULL;
     5693        IssmDouble* G_rigid_precomputed=NULL;
     5694        int         M;
     5695
     5696        /*ice properties: */
     5697        IssmDouble rho_ice,rho_water,rho_earth;
     5698
     5699        /*constants:*/
     5700        IssmDouble constant=0;
     5701
     5702        /*Initialize eustatic component: do not skip this step :):*/
     5703        IssmDouble eustatic = 0.;
     5704
     5705        /*Computational flags:*/
     5706        bool computerigid = true;
     5707        bool computeelastic= true;
     5708        bool scaleoceanarea= false;
     5709
     5710        /*early return if we are not on an ice cap:*/
     5711        if(!IsIceOnlyInElement()){
     5712                constant=0; this->AddInput2(SealevelEustaticMaskEnum,&constant,P0Enum);
     5713                *peustatic=0; //do not forget to assign this pointer, otherwise, global eustatic will be garbage!
     5714                return;
     5715        }
     5716
     5717        /*early return if we are fully floating:*/
     5718        Input2* gr_input=this->GetInput2(MaskOceanLevelsetEnum); _assert_(gr_input);
     5719        if (gr_input->GetInputMax()<=0){
     5720                constant=0; this->AddInput2(SealevelEustaticMaskEnum,&constant,P0Enum);
     5721                *peustatic=0; //do not forget to assign this pointer, otherwise, global eustatic will be garbage!
     5722                return;
     5723        }
     5724
     5725        /*If we are here, we are on ice that is fully grounded or half-way to floating:*/
     5726        if ((gr_input->GetInputMin())<0){
     5727                notfullygrounded=true; //used later on.
     5728        }
     5729   
     5730        /*Inform mask: */
     5731        constant=1; this->AddInput2(SealevelEustaticMaskEnum,&constant,P0Enum);
     5732
     5733        /*recover material parameters: */
     5734        rho_ice=FindParam(MaterialsRhoIceEnum);
     5735        rho_water=FindParam(MaterialsRhoSeawaterEnum);
     5736        rho_earth=FindParam(MaterialsEarthDensityEnum);
     5737
     5738        /*recover love numbers and computational flags: */
     5739        this->parameters->FindParam(&computerigid,SealevelriseRigidEnum);
     5740        this->parameters->FindParam(&computeelastic,SealevelriseElasticEnum);
     5741        this->parameters->FindParam(&scaleoceanarea,SealevelriseOceanAreaScalingEnum);
     5742
     5743        /*recover elastic green function:*/
     5744        DoubleVecParam* parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelriseGRigidEnum)); _assert_(parameter);
     5745        parameter->GetParameterValueByPointer(&G_rigid_precomputed,&M);
     5746
     5747        if(computeelastic){
     5748                parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelriseGElasticEnum)); _assert_(parameter);
     5749                parameter->GetParameterValueByPointer(&G_elastic_precomputed,&M);
     5750        }
     5751
     5752        /*how many dofs are we working with here? */
     5753        this->parameters->FindParam(&gsize,MeshNumberofverticesEnum);
     5754
     5755        /*retrieve indices:*/
     5756        if(computerigid){this->inputs2->GetArray(SealevelriseIndicesEnum,this->lid,&indices,&dummy); _assert_(dummy==gsize);}
    57025757
    57035758        /*Compute area of element. Scale it by grounded fraction if not fully grounded: */
    … …  
    57555810
    57565811                        /*Compute alpha angle between centroid and current vertex and index into precomputed tables: */
    5757                         lati=latitude[i]/180*PI; longi=longitude[i]/180*PI;
     5812                        /*lati=latitude[i]/180*PI; longi=longitude[i]/180*PI;
    57585813                        delPhi=fabs(lati-late); delLambda=fabs(longi-longe);
    57595814                        alpha=2.*asin(sqrt(pow(sin(delPhi/2),2.0)+cos(lati)*cos(late)*pow(sin(delLambda/2),2)));
    5760                         index=reCast<int,IssmDouble>(alpha/PI*reCast<IssmDouble,int>(M-1));
     5815                        index=reCast<int,IssmDouble>(alpha/PI*reCast<IssmDouble,int>(M-1));*/
     5816                        index=reCast<int,IssmDouble>(indices[i]);
    57615817
    57625818                        //Elastic component  (from Eq 17 in Adhikari et al, GMD 2015)
    … …  
    57675823                }
    57685824        }
     5825
     5826        /*Free ressources:*/
     5827        if(computerigid)xDelete<IssmDouble>(indices);
    57695828
    57705829        /*Assign output pointer:*/
    … …  
    59245983
    59255984        /*diverse:*/
    5926         int gsize;
     5985        int gsize,dummy;
    59275986        bool spherical=true;
    59285987        IssmDouble llr_list[NUMVERTICES][3];
    … …  
    59405999        int         M;
    59416000        int         index;
     6001        IssmDouble* indices=NULL;
    59426002        IssmDouble alpha;
    59436003        IssmDouble delPhi,delLambda;
    … …  
    59746034        this->parameters->FindParam(&gsize,MeshNumberofverticesEnum);
    59756035
     6036        /*retrieve indices:*/
     6037        if(computerigid){this->inputs2->GetArray(SealevelriseIndicesEnum,this->lid,&indices,&dummy); _assert_(dummy==gsize);}
     6038
    59766039        /*From Sg_old, recover water sea level rise:*/
    59776040        S=0; for(int i=0;i<NUMVERTICES;i++) S+=Sg_old[this->vertices[i]->Sid()]/NUMVERTICES;
    … …  
    59796042        /*Compute area of element:*/
    59806043        area=GetAreaSpherical();
    5981 
    5982         /* Where is the centroid of this element?:{{{*/
    5983         ::GetVerticesCoordinates(&llr_list[0][0],this->vertices,NUMVERTICES,spherical);
    5984 
    5985         minlong=400; maxlong=-20;
    5986         for (int i=0;i<NUMVERTICES;i++){
    5987                 llr_list[i][0]=(90-llr_list[i][0]);
    5988                 if(llr_list[i][1]<0)llr_list[i][1]=180+(180+llr_list[i][1]);
    5989                 if(llr_list[i][1]>maxlong)maxlong=llr_list[i][1];
    5990                 if(llr_list[i][1]<minlong)minlong=llr_list[i][1];
    5991         }
    5992         if(minlong==0 && maxlong>180){
    5993                 if (llr_list[0][1]==0)llr_list[0][1]=360;
    5994                 if (llr_list[1][1]==0)llr_list[1][1]=360;
    5995                 if (llr_list[2][1]==0)llr_list[2][1]=360;
    5996         }
    5997 
    5998         // correction at the north pole
    5999         if(llr_list[0][0]==0)llr_list[0][1]=(llr_list[1][1]+llr_list[2][1])/2.0;
    6000         if(llr_list[1][0]==0)llr_list[1][1]=(llr_list[0][1]+llr_list[2][1])/2.0;
    6001         if(llr_list[2][0]==0)llr_list[2][1]=(llr_list[0][1]+llr_list[1][1])/2.0;
    6002 
    6003         //correction at the south pole
    6004         if(llr_list[0][0]==180)llr_list[0][1]=(llr_list[1][1]+llr_list[2][1])/2.0;
    6005         if(llr_list[1][0]==180)llr_list[1][1]=(llr_list[0][1]+llr_list[2][1])/2.0;
    6006         if(llr_list[2][0]==180)llr_list[2][1]=(llr_list[0][1]+llr_list[1][1])/2.0;
    6007 
    6008         late=(llr_list[0][0]+llr_list[1][0]+llr_list[2][0])/3.0;
    6009         longe=(llr_list[0][1]+llr_list[1][1]+llr_list[2][1])/3.0;
    6010 
    6011         late=90-late;
    6012         if(longe>180)longe=(longe-180)-180;
    6013 
    6014         late=late/180*PI;
    6015         longe=longe/180*PI;
    6016         /*}}}*/
    60176044
    60186045        /*recover rigid and elastic green functions:*/
    … …  
    60306057
    60316058                /*Compute alpha angle between centroid and current vertex : */
    6032                 lati=latitude[i]/180*PI; longi=longitude[i]/180*PI;
     6059                /*lati=latitude[i]/180*PI; longi=longitude[i]/180*PI;
    60336060                delPhi=fabs(lati-late); delLambda=fabs(longi-longe);
    60346061                alpha=2.*asin(sqrt(pow(sin(delPhi/2),2.0)+cos(lati)*cos(late)*pow(sin(delLambda/2),2)));
    6035                 index=reCast<int,IssmDouble>(alpha/PI*(M-1));
     6062                index=reCast<int,IssmDouble>(alpha/PI*(M-1));*/
     6063                       
     6064                index=reCast<int,IssmDouble>(indices[i]);
    60366065
    60376066                /*Rigid earth gravitational perturbation: */
    … …  
    60466075        }
    60476076
     6077        /*Free ressources:*/
     6078        if(computerigid)xDelete<IssmDouble>(indices);
     6079
    60486080
    60496081        return;
    … …  
    60536085
    60546086        /*diverse:*/
    6055         int gsize;
     6087        int gsize,dummy;
    60566088        bool spherical=true;
    60576089        IssmDouble llr_list[NUMVERTICES][3];
    … …  
    60696101        IssmDouble* H_elastic_precomputed = NULL;
    60706102        int         M;
     6103        IssmDouble* indices=NULL;
     6104        int         index;
    60716105
    60726106        /*computation of Green functions:*/
    … …  
    61096143        this->parameters->FindParam(&gsize,MeshNumberofverticesEnum);
    61106144
     6145        /*retrieve indices:*/
     6146        if(computerigid){this->inputs2->GetArray(SealevelriseIndicesEnum,this->lid,&indices,&dummy); _assert_(dummy==gsize);}
     6147
    61116148        /*compute area of element:*/
    61126149        area=GetAreaSpherical();
    6113 
    6114         /*element centroid (spherical): */
    6115         /* Where is the centroid of this element?:{{{*/
    6116         ::GetVerticesCoordinates(&llr_list[0][0],this->vertices,NUMVERTICES,spherical);
    6117 
    6118         minlong=400; maxlong=-20;
    6119         for (int i=0;i<NUMVERTICES;i++){
    6120                 llr_list[i][0]=(90-llr_list[i][0]);
    6121                 if(llr_list[i][1]<0)llr_list[i][1]=180+(180+llr_list[i][1]);
    6122                 if(llr_list[i][1]>maxlong)maxlong=llr_list[i][1];
    6123                 if(llr_list[i][1]<minlong)minlong=llr_list[i][1];
    6124         }
    6125         if(minlong==0 && maxlong>180){
    6126                 if (llr_list[0][1]==0)llr_list[0][1]=360;
    6127                 if (llr_list[1][1]==0)llr_list[1][1]=360;
    6128                 if (llr_list[2][1]==0)llr_list[2][1]=360;
    6129         }
    6130 
    6131         // correction at the north pole
    6132         if(llr_list[0][0]==0)llr_list[0][1]=(llr_list[1][1]+llr_list[2][1])/2.0;
    6133         if(llr_list[1][0]==0)llr_list[1][1]=(llr_list[0][1]+llr_list[2][1])/2.0;
    6134         if(llr_list[2][0]==0)llr_list[2][1]=(llr_list[0][1]+llr_list[1][1])/2.0;
    6135 
    6136         //correction at the south pole
    6137         if(llr_list[0][0]==180)llr_list[0][1]=(llr_list[1][1]+llr_list[2][1])/2.0;
    6138         if(llr_list[1][0]==180)llr_list[1][1]=(llr_list[0][1]+llr_list[2][1])/2.0;
    6139         if(llr_list[2][0]==180)llr_list[2][1]=(llr_list[0][1]+llr_list[1][1])/2.0;
    6140 
    6141         late=(llr_list[0][0]+llr_list[1][0]+llr_list[2][0])/3.0;
    6142         longe=(llr_list[0][1]+llr_list[1][1]+llr_list[2][1])/3.0;
    6143 
    6144         late=90-late;
    6145         if(longe>180)longe=(longe-180)-180;
    6146 
    6147         late=late/180*PI;
    6148         longe=longe/180*PI;
    6149         /*}}}*/
    61506150
    61516151        /*figure out gravity center of our element (Cartesian): */
    … …  
    61926192        for(int i=0;i<gsize;i++){
    61936193
    6194                 //indices[i]=i;
    6195 
    61966194                /*Compute alpha angle between centroid and current vertex: */
    6197                 lati=latitude[i]/180*PI; longi=longitude[i]/180*PI;
    6198 
     6195                /*lati=latitude[i]/180*PI; longi=longitude[i]/180*PI;
    61996196                delPhi=fabs(lati-late); delLambda=fabs(longi-longe);
    6200                 alpha=2.*asin(sqrt(pow(sin(delPhi/2),2.0)+cos(lati)*cos(late)*pow(sin(delLambda/2),2)));
     6197                alpha=2.*asin(sqrt(pow(sin(delPhi/2),2.0)+cos(lati)*cos(late)*pow(sin(delLambda/2),2)));*/
     6198
     6199                index=reCast<int,IssmDouble>(indices[i]);
    62016200
    62026201                /*Compute azimuths, both north and east components: */
    … …  
    62156214
    62166215                /*Elastic component  (from Eq 17 in Adhikari et al, GMD 2015): */
    6217                 int index=reCast<int,IssmDouble>(alpha/PI*(M-1));
    62186216                U_elastic[i] += U_elastic_precomputed[index];
    62196217                if(horiz){
    … …  
    62406238
    62416239        /*free ressources:*/
     6240        if(computerigid)xDelete<IssmDouble>(indices);
    62426241        xDelete<IssmDouble>(U_values);
    62436242        xDelete<IssmDouble>(U_elastic);

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

    r24915 r24924  
    166166                IssmDouble OceanAverage(IssmDouble* Sg);
    167167                void    SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old,IssmDouble eartharea);
     168                void    SealevelriseGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius);
    168169                void    SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea);
    169170                void    SealevelriseEustaticIce(IssmDouble* Sgi, IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea);

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

    r24917 r24924  
    46844684#endif
    46854685#ifdef _HAVE_SEALEVELRISE_
     4686void FemModel::SealevelriseGeometry(IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius) { /*{{{*/
     4687
     4688        /*Run sealevelrie geometry routine in elements:*/
     4689        for(int i=0;i<elements->Size();i++){
     4690                Element*   element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
     4691                element->SealevelriseGeometry(latitude,longitude,radius);
     4692        }
     4693
     4694}
     4695/*}}}*/
    46864696void FemModel::SealevelriseEustatic(Vector<IssmDouble>* pRSLgi, IssmDouble* peartharea, IssmDouble* poceanarea, IssmDouble* peustatic, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius,int loop) { /*{{{*/
    46874697

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

    r24916 r24924  
    163163                #ifdef _HAVE_SEALEVELRISE_
    164164                void SealevelriseEustatic(Vector<IssmDouble>* pSgi, IssmDouble* peartharea, IssmDouble* poceanarea, IssmDouble* peustatic, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius,int loop);
     165                void SealevelriseGeometry(IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius);
    165166                void SealevelriseNonEustatic(Vector<IssmDouble>* pSgo, Vector<IssmDouble>* pSg_old, IssmDouble eartharea, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius,bool verboseconvolution,int loop);
    166167                void SealevelriseRotationalFeedback(Vector<IssmDouble>* pRSLgo_rot, Vector<IssmDouble>* pRSLg_old, IssmDouble* pIxz, IssmDouble* pIyz, IssmDouble* pIzz, IssmDouble eartharea,IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius);

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

    r24916 r24924  
    5555void geodetic_core(FemModel* femmodel);
    5656void steric_core(FemModel* femmodel);
     57void sealevelrise_core_geometry(FemModel* femmodel);
    5758Vector<IssmDouble>* sealevelrise_core_eustatic(FemModel* femmodel,IssmDouble* peartharea,IssmDouble* poceanarea);
    5859Vector<IssmDouble>* sealevelrise_core_noneustatic(FemModel* femmodel,Vector<IssmDouble>* RSLg_eustatic,IssmDouble eartharea,IssmDouble oceanarea);

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

    r24916 r24924  
    4141        /*set SLR configuration: */
    4242        femmodel->SetCurrentConfiguration(SealevelriseAnalysisEnum);
     43
     44        /*run geometry core:*/
     45        sealevelrise_core_geometry(femmodel);
    4346
    4447        /*Run geodetic:*/
    … …  
    304307}
    305308/*}}}*/
     309void sealevelrise_core_geometry(FemModel* femmodel) {  /*{{{*/
     310
     311        /*Geometry core where we compute indices into tables pre computed in the SealevelRiseAnalysis: */
     312
     313        /*parameters: */
     314        bool spherical=true;
     315        IssmDouble *latitude  = NULL;
     316        IssmDouble *longitude = NULL;
     317        IssmDouble *radius    = NULL;
     318
     319        /*Verbose: */
     320        if(VerboseSolution()) _printf0_("         computing geometrical offsets into precomputed Green tables \n");
     321
     322        /*first, recover lat,long and radius vectors from vertices: */
     323        VertexCoordinatesx(&latitude,&longitude,&radius,femmodel->vertices,spherical);
     324
     325        /*call the FemModel geometry module: */
     326        femmodel->SealevelriseGeometry(latitude, longitude, radius);
     327
     328        /*Free ressources:*/
     329        xDelete<IssmDouble>(latitude);
     330        xDelete<IssmDouble>(longitude);
     331        xDelete<IssmDouble>(radius);
     332
     333
     334}/*}}}*/
    306335Vector<IssmDouble>* sealevelrise_core_eustatic(FemModel* femmodel,IssmDouble* peartharea,IssmDouble* poceanarea){ /*{{{*/
    307336

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

    r24920 r24924  
    9999        }
    100100
    101         //if(isslr) sealevelrise_geometrycore(femmodel);
     101        if(isslr) sealevelrise_core_geometry(femmodel);
    102102
    103103        while(time < finaltime - (yts*DBL_EPSILON)){ //make sure we run up to finaltime.

  • issm/trunk-jpl/src/c/shared/Enum/Enum.vim

    r24921 r24924  
    683683syn keyword cConstant SealevelriseSpcthicknessEnum
    684684syn keyword cConstant SealevelriseHydroRateEnum
     685syn keyword cConstant SealevelriseIndicesEnum
    685686syn keyword cConstant SedimentHeadEnum
    686687syn keyword cConstant SedimentHeadOldEnum

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

    r24921 r24924  
    680680        SealevelriseSpcthicknessEnum,
    681681        SealevelriseHydroRateEnum,
     682        SealevelriseIndicesEnum,
    682683        SedimentHeadEnum,
    683684        SedimentHeadOldEnum,

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

    r24921 r24924  
    685685                case SealevelriseSpcthicknessEnum : return "SealevelriseSpcthickness";
    686686                case SealevelriseHydroRateEnum : return "SealevelriseHydroRate";
     687                case SealevelriseIndicesEnum : return "SealevelriseIndices";
    687688                case SedimentHeadEnum : return "SedimentHead";
    688689                case SedimentHeadOldEnum : return "SedimentHeadOld";

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

    r24921 r24924  
    700700              else if (strcmp(name,"SealevelriseSpcthickness")==0) return SealevelriseSpcthicknessEnum;
    701701              else if (strcmp(name,"SealevelriseHydroRate")==0) return SealevelriseHydroRateEnum;
     702              else if (strcmp(name,"SealevelriseIndices")==0) return SealevelriseIndicesEnum;
    702703              else if (strcmp(name,"SedimentHead")==0) return SedimentHeadEnum;
    703704              else if (strcmp(name,"SedimentHeadOld")==0) return SedimentHeadOldEnum;
    … …  
    751752              else if (strcmp(name,"SmbMassBalanceClimate")==0) return SmbMassBalanceClimateEnum;
    752753              else if (strcmp(name,"SmbMassBalance")==0) return SmbMassBalanceEnum;
    753               else if (strcmp(name,"SmbMassBalanceSubstep")==0) return SmbMassBalanceSubstepEnum;
    754754         else stage=7;
    755755   }
    756756   if(stage==7){
    757               if (strcmp(name,"SmbMassBalanceTransient")==0) return SmbMassBalanceTransientEnum;
     757              if (strcmp(name,"SmbMassBalanceSubstep")==0) return SmbMassBalanceSubstepEnum;
     758              else if (strcmp(name,"SmbMassBalanceTransient")==0) return SmbMassBalanceTransientEnum;
    758759              else if (strcmp(name,"SmbMeanLHF")==0) return SmbMeanLHFEnum;
    759760              else if (strcmp(name,"SmbMeanSHF")==0) return SmbMeanSHFEnum;
    … …  
    874875              else if (strcmp(name,"Outputdefinition13")==0) return Outputdefinition13Enum;
    875876              else if (strcmp(name,"Outputdefinition14")==0) return Outputdefinition14Enum;
    876               else if (strcmp(name,"Outputdefinition15")==0) return Outputdefinition15Enum;
    877877         else stage=8;
    878878   }
    879879   if(stage==8){
    880               if (strcmp(name,"Outputdefinition16")==0) return Outputdefinition16Enum;
     880              if (strcmp(name,"Outputdefinition15")==0) return Outputdefinition15Enum;
     881              else if (strcmp(name,"Outputdefinition16")==0) return Outputdefinition16Enum;
    881882              else if (strcmp(name,"Outputdefinition17")==0) return Outputdefinition17Enum;
    882883              else if (strcmp(name,"Outputdefinition18")==0) return Outputdefinition18Enum;
    … …  
    997998              else if (strcmp(name,"BeckmannGoosseFloatingMeltRate")==0) return BeckmannGoosseFloatingMeltRateEnum;
    998999              else if (strcmp(name,"BedSlopeSolution")==0) return BedSlopeSolutionEnum;
    999               else if (strcmp(name,"BoolExternalResult")==0) return BoolExternalResultEnum;
    10001000         else stage=9;
    10011001   }
    10021002   if(stage==9){
    1003               if (strcmp(name,"BoolInput")==0) return BoolInputEnum;
     1003              if (strcmp(name,"BoolExternalResult")==0) return BoolExternalResultEnum;
     1004              else if (strcmp(name,"BoolInput")==0) return BoolInputEnum;
    10041005              else if (strcmp(name,"BoolInput2")==0) return BoolInput2Enum;
    10051006              else if (strcmp(name,"IntInput2")==0) return IntInput2Enum;
    … …  
    11201121              else if (strcmp(name,"IceVolumeAboveFloatation")==0) return IceVolumeAboveFloatationEnum;
    11211122              else if (strcmp(name,"IceVolumeAboveFloatationScaled")==0) return IceVolumeAboveFloatationScaledEnum;
    1122               else if (strcmp(name,"IceVolume")==0) return IceVolumeEnum;
    11231123         else stage=10;
    11241124   }
    11251125   if(stage==10){
    1126               if (strcmp(name,"IceVolumeScaled")==0) return IceVolumeScaledEnum;
     1126              if (strcmp(name,"IceVolume")==0) return IceVolumeEnum;
     1127              else if (strcmp(name,"IceVolumeScaled")==0) return IceVolumeScaledEnum;
    11271128              else if (strcmp(name,"IcefrontMassFlux")==0) return IcefrontMassFluxEnum;
    11281129              else if (strcmp(name,"IcefrontMassFluxLevelset")==0) return IcefrontMassFluxLevelsetEnum;
    … …  
    12431244              else if (strcmp(name,"PentaInput")==0) return PentaInputEnum;
    12441245              else if (strcmp(name,"Profiler")==0) return ProfilerEnum;
    1245               else if (strcmp(name,"ProfilingCurrentFlops")==0) return ProfilingCurrentFlopsEnum;
    12461246         else stage=11;
    12471247   }
    12481248   if(stage==11){
    1249               if (strcmp(name,"ProfilingCurrentMem")==0) return ProfilingCurrentMemEnum;
     1249              if (strcmp(name,"ProfilingCurrentFlops")==0) return ProfilingCurrentFlopsEnum;
     1250              else if (strcmp(name,"ProfilingCurrentMem")==0) return ProfilingCurrentMemEnum;
    12501251              else if (strcmp(name,"ProfilingSolutionTime")==0) return ProfilingSolutionTimeEnum;
    12511252              else if (strcmp(name,"Regionaloutput")==0) return RegionaloutputEnum;
Note: See TracChangeset for help on using the changeset viewer.