Changeset 24938


Ignore:
Timestamp:
05/31/20 00:21:13 (5 years ago)
Author:
Eric.Larour
Message:

CHG: optimizing sea level masks. Create new SealevelMasks classs
and new sealevelrise_core_masks that pre compute masks, and then
gets cores quick access to them.

Location:
issm/trunk-jpl/src/c
Files:
1 added
11 edited

Legend:

Unmodified
Added
Removed

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

    r24935 r24938  
    3030class DatasetInput2;
    3131class IoModel;
     32class SealevelMasks;
    3233class Gauss;
    3334class ElementVector;
    … …  
    373374                #endif
    374375                #ifdef _HAVE_SEALEVELRISE_
     376                virtual void          SetSealevelMasks(SealevelMasks* masks)=0;
    375377                virtual IssmDouble    GetArea3D(void)=0;
    376378                virtual IssmDouble    GetAreaSpherical(void)=0;
    377                 virtual IssmDouble    OceanAverage(IssmDouble* Sg)=0;
    378                 virtual IssmDouble    OceanArea(void)=0;
    379                 virtual void          SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old,IssmDouble eartharea)=0;
    380                 virtual void          SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea)=0;
     379                virtual IssmDouble    OceanAverage(IssmDouble* Sg, SealevelMasks* masks)=0;
     380                virtual void          SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks,IssmDouble eartharea)=0;
     381                virtual void          SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic,SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea)=0;
    381382                virtual void          SealevelriseGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius)=0;
    382                 virtual void          SealevelriseNonEustatic(IssmDouble* Sgo, IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea)=0;
    383                 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;
     383                virtual void          SealevelriseNonEustatic(IssmDouble* Sgo, IssmDouble* Sg_old,SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea)=0;
     384                virtual void          SealevelriseGeodetic(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg,SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea,int horiz)=0;
    384385                #endif
    385386

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

    r24924 r24938  
    212212                #endif
    213213                #ifdef _HAVE_SEALEVELRISE_
    214                 IssmDouble    OceanArea(void){_error_("not implemented yet!");};
    215                 IssmDouble    OceanAverage(IssmDouble* Sg){_error_("not implemented yet!");};
    216                 void    SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old,IssmDouble eartharea){_error_("not implemented yet!");};
     214                IssmDouble    OceanAverage(IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet!");};
     215                void    SetSealevelMasks(SealevelMasks* masks){_error_("not implemented yet!");};
     216                void    SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks,IssmDouble eartharea){_error_("not implemented yet!");};
    217217                void    SealevelriseGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius){_error_("not implemented yet!");};
    218                 void    SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){_error_("not implemented yet!");};
    219                 void    SealevelriseNonEustatic(IssmDouble* Sgo, IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea){_error_("not implemented yet!");};
    220                 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){_error_("not implemented yet!");};
     218                void    SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){_error_("not implemented yet!");};
     219                void    SealevelriseNonEustatic(IssmDouble* Sgo, IssmDouble* Sg_old, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea){_error_("not implemented yet!");};
     220                void    SealevelriseGeodetic(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea,int horiz){_error_("not implemented yet!");};
    221221                #endif
    222222

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

    r24924 r24938  
    172172#endif
    173173#ifdef _HAVE_SEALEVELRISE_
    174                 void    SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old,IssmDouble eartharea){_error_("not implemented yet!");};
     174                void    SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks,IssmDouble eartharea){_error_("not implemented yet!");};
     175                void    SetSealevelMasks(SealevelMasks* masks){_error_("not implemented yet!");};
    175176                void    SealevelriseGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius){_error_("not implemented yet!");};
    176                 void    SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){_error_("not implemented yet!");};
    177                 void    SealevelriseNonEustatic(IssmDouble* Sgo,IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea){_error_("not implemented yet!");};
    178                 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){_error_("not implemented yet!");};
    179                 IssmDouble    OceanArea(void){_error_("not implemented yet!");};
    180                 IssmDouble    OceanAverage(IssmDouble* Sg){_error_("not implemented yet!");};
     177                void    SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){_error_("not implemented yet!");};
     178                void    SealevelriseNonEustatic(IssmDouble* Sgo,IssmDouble* Sg_old, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea){_error_("not implemented yet!");};
     179                void    SealevelriseGeodetic(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea,int horiz){_error_("not implemented yet!");};
     180                IssmDouble    OceanAverage(IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet!");};
    181181#endif
    182182

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

    r24924 r24938  
    178178#endif
    179179#ifdef _HAVE_SEALEVELRISE_
    180                 void    SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old,IssmDouble eartharea){_error_("not implemented yet!");};
     180                void    SetSealevelMasks(SealevelMasks* masks){_error_("not implemented yet!");};
     181                void    SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks,IssmDouble eartharea){_error_("not implemented yet!");};
    181182                void    SealevelriseGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius){_error_("not implemented yet!");};
    182                 void    SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){_error_("not implemented yet!");};
    183                 void    SealevelriseNonEustatic(IssmDouble* Sgo, IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea){_error_("not implemented yet!");};
    184                 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){_error_("not implemented yet!");};
    185                 IssmDouble    OceanArea(void){_error_("not implemented yet!");};
    186                 IssmDouble    OceanAverage(IssmDouble* Sg){_error_("not implemented yet!");};
     183                void    SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){_error_("not implemented yet!");};
     184                void    SealevelriseNonEustatic(IssmDouble* Sgo, IssmDouble* Sg_old, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea){_error_("not implemented yet!");};
     185                void    SealevelriseGeodetic(IssmDouble* Up ,IssmDouble* North, IssmDouble* East, IssmDouble* Sg, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea,int horiz){_error_("not implemented yet!");};
     186                IssmDouble    OceanAverage(IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet!");};
    187187#endif
    188188

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

    r24935 r24938  
    30163016
    30173017        if(!IsIceInElement())return 0.;
    3018         //if(!IsIceOnlyInElement()) return 0;
    30193018
    30203019        int domaintype;
    … …  
    54475446#endif
    54485447#ifdef _HAVE_SEALEVELRISE_
    5449 IssmDouble    Tria::OceanArea(void){ /*{{{*/
    5450 
    5451         if(IsOceanInElement()) return GetAreaSpherical();
    5452         else return 0;
    5453 
    5454 }
    5455 /*}}}*/
    5456 IssmDouble Tria::OceanAverage(IssmDouble* Sg){ /*{{{*/
    5457 
    5458         if(IsOceanInElement()){
     5448IssmDouble Tria::OceanAverage(IssmDouble* Sg, SealevelMasks* masks){ /*{{{*/
     5449
     5450        if(masks->isoceanin[this->lid]){
    54595451
    54605452                IssmDouble area;
    … …  
    54735465}
    54745466/*}}}*/
    5475 void    Tria::SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old,IssmDouble eartharea){/*{{{*/
     5467void    Tria::SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks, IssmDouble eartharea){/*{{{*/
    54765468        /*early return if we are not on an ice cap OR ocean:*/
    5477         if(!IsIceOnlyInElement() && !IsOceanInElement()){
     5469        if(!masks->isiceonly[this->lid] && !masks->isoceanin[this->lid]){
    54785470                dI_list[0] = 0.0; // this is important!!!
    54795471                dI_list[1] = 0.0; // this is important!!!
    … …  
    55285520        re=(llr_list[0][2]+llr_list[1][2]+llr_list[2][2])/3.0;
    55295521
    5530         if(IsOceanInElement()){
     5522        if(masks->isoceanin[this->lid]){
    55315523                IssmDouble rho_water, S;
    55325524
    … …  
    55465538                dI_list[2] = +4*PI*(rho_water*S*area)*pow(re,4)*(1-pow(sin(late),2))/eartharea;
    55475539        }
    5548         else if(IsIceOnlyInElement()){
     5540        else if(masks->isiceonly[this->lid]){
    55495541                IssmDouble rho_ice, I;
    55505542
    … …  
    55645556        return;
    55655557}/*}}}*/
     5558void    Tria::SetSealevelMasks(SealevelMasks* masks){ /*{{{*/
     5559
     5560        masks->isiceonly[this->lid]=this->IsIceOnlyInElement();
     5561        masks->isoceanin[this->lid]=this->IsOceanInElement();
     5562       
     5563        /*are we fully floating:*/
     5564        Input2* gr_input=this->GetInput2(MaskOceanLevelsetEnum); _assert_(gr_input);
     5565        if (gr_input->GetInputMax()<=0)masks->isfullyfloating[this->lid]=true;
     5566        else masks->isfullyfloating[this->lid]=false;
     5567
     5568        /*are we not fully grounded: */
     5569        if ((gr_input->GetInputMin())<0) masks->notfullygrounded[this->lid]=true;
     5570        else masks->notfullygrounded[this->lid]=false;
     5571
     5572}
     5573/*}}}*/
    55665574void    Tria::SealevelriseGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius){ /*{{{*/
    55675575
    … …  
    56535661}
    56545662/*}}}*/
    5655 void    Tria::SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){ /*{{{*/
     5663void    Tria::SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic, SealevelMasks* masks, IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){ /*{{{*/
    56565664
    56575665        /*Computational flags:*/
    … …  
    56615669        this->parameters->FindParam(&bp_compute_fingerprints,DslComputeFingerprintsEnum);
    56625670
    5663         if(!IsOceanInElement()){
     5671        if(!masks->isoceanin[this->lid]){
    56645672                /*ok, there is ocean in this element, we should compute eustatic loads for the ocean if we have requested
    56655673                 *bottom pressure fingerprints:*/
    … …  
    56685676        //if(!IsIceInElement()){
    56695677                /*there is ice in this eleemnt, let's compute the eustatic response for ice changes:*/
    5670                 this->SealevelriseEustaticIce(Sgi,peustatic,latitude,longitude,radius,oceanarea,eartharea);
     5678                this->SealevelriseEustaticIce(Sgi,peustatic,masks, latitude,longitude,radius,oceanarea,eartharea);
    56715679        //}
    56725680
    56735681}
    56745682/*}}}*/
    5675 void    Tria::SealevelriseEustaticIce(IssmDouble* Sgi, IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){ /*{{{*/
     5683void    Tria::SealevelriseEustaticIce(IssmDouble* Sgi, IssmDouble* peustatic, SealevelMasks* masks, IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){ /*{{{*/
    56765684
    56775685        /*diverse:*/
    … …  
    57085716
    57095717        /*early return if we are not on an ice cap:*/
    5710         if(!IsIceOnlyInElement()){
     5718        if(!masks->isiceonly[this->lid]){
    57115719                constant=0; this->AddInput2(SealevelEustaticMaskEnum,&constant,P0Enum);
    57125720                *peustatic=0; //do not forget to assign this pointer, otherwise, global eustatic will be garbage!
    … …  
    57155723
    57165724        /*early return if we are fully floating:*/
    5717         Input2* gr_input=this->GetInput2(MaskOceanLevelsetEnum); _assert_(gr_input);
    5718         if (gr_input->GetInputMax()<=0){
     5725        if (masks->isfullyfloating[this->lid]){
    57195726                constant=0; this->AddInput2(SealevelEustaticMaskEnum,&constant,P0Enum);
    57205727                *peustatic=0; //do not forget to assign this pointer, otherwise, global eustatic will be garbage!
    … …  
    57235730
    57245731        /*If we are here, we are on ice that is fully grounded or half-way to floating:*/
    5725         if ((gr_input->GetInputMin())<0){
    5726                 notfullygrounded=true; //used later on.
    5727         }
     5732        if (masks->notfullygrounded[this->lid]) notfullygrounded=true; //used later on.
    57285733   
    57295734        /*Inform mask: */
    … …  
    58315836}
    58325837/*}}}*/
    5833 void    Tria::SealevelriseEustaticBottomPressure(IssmDouble* Sgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){ /*{{{*/
     5838void    Tria::SealevelriseEustaticBottomPressure(IssmDouble* Sgi,IssmDouble* peustatic, IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){ /*{{{*/
    58345839
    58355840        /*diverse:*/
    … …  
    59775982}
    59785983/*}}}*/
    5979 void    Tria::SealevelriseNonEustatic(IssmDouble* Sgo,IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea){ /*{{{*/
     5984void    Tria::SealevelriseNonEustatic(IssmDouble* Sgo,IssmDouble* Sg_old, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea){ /*{{{*/
    59805985
    59815986        /*diverse:*/
    … …  
    60116016
    60126017        /*early return if we are not on the ocean:*/
    6013         if (!IsOceanInElement()){
     6018        if (!masks->isoceanin[this->lid]){
    60146019                constant=0; this->AddInput2(SealevelEustaticOceanMaskEnum,&constant,P0Enum);
    60156020                return;
    … …  
    60766081}
    60776082/*}}}*/
    6078 void    Tria::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){ /*{{{*/
     6083void    Tria::SealevelriseGeodetic(IssmDouble* Up, IssmDouble* North ,IssmDouble* East,IssmDouble* Sg, SealevelMasks* masks, IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea,int horiz){ /*{{{*/
    60796084
    60806085        /*diverse:*/
    … …  
    61166121
    61176122        /*early return if we are not on the ocean or on an ice cap:*/
    6118         if(!IsIceOnlyInElement() && !IsOceanInElement()) return;
     6123        if(!masks->isiceonly[this->lid] && !masks->isoceanin[this->lid]) return;
    61196124
    61206125        /*early return if we are fully floating: */
    6121         Input2* gr_input=this->GetInput2(MaskOceanLevelsetEnum); _assert_(gr_input);
    6122         if(gr_input->GetInputMax()<=0)return;
     6126        if(masks->isfullyfloating[this->lid])return;
    61236127
    61246128        /*recover computational flags: */
    … …  
    61856189
    61866190        /*we are going to use the result of these masks a lot: */
    6187         bool isiceonlyinelement=IsIceOnlyInElement();
    6188         bool isoceaninelement=IsOceanInElement();
     6191        bool isiceonlyinelement=masks->isiceonly[this->lid];
     6192        bool isoceaninelement=masks->isoceanin[this->lid];
    61896193
    61906194        for(int i=0;i<gsize;i++){

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

    r24924 r24938  
    163163                #endif
    164164                #ifdef _HAVE_SEALEVELRISE_
    165                 IssmDouble OceanArea(void);
    166                 IssmDouble OceanAverage(IssmDouble* Sg);
    167                 void    SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old,IssmDouble eartharea);
     165                IssmDouble OceanAverage(IssmDouble* Sg, SealevelMasks* masks);
     166                void    SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old,SealevelMasks* masks, IssmDouble eartharea);
     167                void    SetSealevelMasks(SealevelMasks* masks);
    168168                void    SealevelriseGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius);
    169                 void    SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea);
    170                 void    SealevelriseEustaticIce(IssmDouble* Sgi, IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea);
     169                void    SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic,SealevelMasks* masks, IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea);
     170                void    SealevelriseEustaticIce(IssmDouble* Sgi, IssmDouble* peustatic,SealevelMasks* masks, IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea);
    171171                void    SealevelriseEustaticBottomPressure(IssmDouble* Sgi, IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea);
    172                 void    SealevelriseNonEustatic(IssmDouble* Sgo,IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea);
    173                 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);
     172                void    SealevelriseNonEustatic(IssmDouble* Sgo,IssmDouble* Sg_old,SealevelMasks* masks, IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea);
     173                void    SealevelriseGeodetic(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg,SealevelMasks* masks, IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea,int horiz);
    174174                #endif
    175175                /*}}}*/

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

    r24924 r24938  
    46944694}
    46954695/*}}}*/
    4696 void FemModel::SealevelriseEustatic(Vector<IssmDouble>* pRSLgi, IssmDouble* peartharea, IssmDouble* poceanarea, IssmDouble* peustatic, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius,int loop) { /*{{{*/
     4696void FemModel::SealevelriseEustatic(Vector<IssmDouble>* pRSLgi, IssmDouble* peartharea, IssmDouble* poceanarea, IssmDouble* peustatic, SealevelMasks* masks, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius,int loop) { /*{{{*/
    46974697
    46984698        /*serialized vectors:*/
    … …  
    47174717                IssmDouble area=element->GetAreaSpherical();
    47184718                eartharea_cpu += area;
    4719                 if(element->IsOceanInElement()) oceanarea_cpu += area;
     4719                if (masks->isoceanin[i]) oceanarea_cpu += area;
    47204720        }
    47214721        ISSM_MPI_Reduce (&oceanarea_cpu,&oceanarea,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
    … …  
    47294729        for(int i=0;i<elements->Size();i++){
    47304730                Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
    4731                 element->SealevelriseEustatic(RSLgi,&eustatic_cpu_e,latitude,longitude,radius,oceanarea,eartharea);
     4731                element->SealevelriseEustatic(RSLgi,&eustatic_cpu_e,masks, latitude,longitude,radius,oceanarea,eartharea);
    47324732                eustatic_cpu+=eustatic_cpu_e;
    47334733        }
    … …  
    47534753}
    47544754/*}}}*/
    4755 void FemModel::SealevelriseNonEustatic(Vector<IssmDouble>* pRSLgo, Vector<IssmDouble>* pRSLg_old, IssmDouble eartharea,IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, bool verboseconvolution,int loop){/*{{{*/
     4755void FemModel::SealevelriseNonEustatic(Vector<IssmDouble>* pRSLgo, Vector<IssmDouble>* pRSLg_old, IssmDouble eartharea, SealevelMasks* masks, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, bool verboseconvolution,int loop){/*{{{*/
    47564756
    47574757        /*serialized vectors:*/
    … …  
    47744774        for(int i=0;i<elements->Size();i++){
    47754775                Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
    4776                 element->SealevelriseNonEustatic(RSLgo,RSLg_old,latitude,longitude,radius,eartharea);
     4776                element->SealevelriseNonEustatic(RSLgo,RSLg_old,masks, latitude,longitude,radius,eartharea);
    47774777        }
    47784778        pRSLgo->SetValues(gsize,indices,RSLgo,ADD_VAL);
    … …  
    47864786}
    47874787/*}}}*/
    4788 void FemModel::SealevelriseRotationalFeedback(Vector<IssmDouble>* pRSLgo_rot, Vector<IssmDouble>* pRSLg_old, IssmDouble* pIxz, IssmDouble* pIyz, IssmDouble* pIzz, IssmDouble eartharea, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius){/*{{{*/
     4788void FemModel::SealevelriseRotationalFeedback(Vector<IssmDouble>* pRSLgo_rot, Vector<IssmDouble>* pRSLg_old, IssmDouble* pIxz, IssmDouble* pIyz, IssmDouble* pIzz, IssmDouble eartharea, SealevelMasks* masks, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius){/*{{{*/
    47894789
    47904790        /*serialized vectors:*/
    … …  
    48024802        for(int i=0;i<elements->Size();i++){
    48034803                Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
    4804                 element->SealevelriseMomentOfInertia(&moi_list[0],RSLg_old,eartharea);
     4804                element->SealevelriseMomentOfInertia(&moi_list[0],RSLg_old,masks, eartharea);
    48054805                moi_list_cpu[0] += moi_list[0];
    48064806                moi_list_cpu[1] += moi_list[1];
    … …  
    48624862}
    48634863/*}}}*/
    4864 void FemModel::SealevelriseElastic(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, Vector<IssmDouble>* pRSLg, IssmDouble eartharea, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz,int loop,int horiz){/*{{{*/
     4864void FemModel::SealevelriseElastic(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, Vector<IssmDouble>* pRSLg, IssmDouble eartharea, SealevelMasks* masks, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz,int loop,int horiz){/*{{{*/
    48654865
    48664866        /*serialized vectors:*/
    … …  
    48884888        for(int i=0;i<elements->Size();i++){
    48894889                Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
    4890                 element->SealevelriseGeodetic(Up,North,East,RSLg,latitude,longitude,radius,xx,yy,zz,eartharea,horiz);
     4890                element->SealevelriseGeodetic(Up,North,East,RSLg,masks, latitude,longitude,radius,xx,yy,zz,eartharea,horiz);
    48914891        }
    48924892
    … …  
    49084908}
    49094909/*}}}*/
    4910 IssmDouble FemModel::SealevelriseOceanAverage(Vector<IssmDouble>* RSLg,IssmDouble oceanarea) { /*{{{*/
     4910IssmDouble FemModel::SealevelriseOceanAverage(Vector<IssmDouble>* RSLg,SealevelMasks* masks, IssmDouble oceanarea) { /*{{{*/
    49114911
    49124912        IssmDouble* RSLg_serial=NULL;
    … …  
    49224922        for(int i=0;i<elements->Size();i++){
    49234923                Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
    4924                 oceanvalue_cpu += element->OceanAverage(RSLg_serial);
     4924                oceanvalue_cpu += element->OceanAverage(RSLg_serial,masks);
    49254925        }
    49264926       

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

    r24924 r24938  
    1818class Loads;
    1919class Materials;
     20class SealevelMasks;
    2021class Profiler;
    2122class Elements;
    … …  
    162163                #endif
    163164                #ifdef _HAVE_SEALEVELRISE_
    164                 void SealevelriseEustatic(Vector<IssmDouble>* pSgi, IssmDouble* peartharea, IssmDouble* poceanarea, IssmDouble* peustatic, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius,int loop);
     165                void SealevelriseEustatic(Vector<IssmDouble>* pSgi, IssmDouble* peartharea, IssmDouble* poceanarea, IssmDouble* peustatic, SealevelMasks* masks, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius,int loop);
    165166                void SealevelriseGeometry(IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius);
    166                 void SealevelriseNonEustatic(Vector<IssmDouble>* pSgo, Vector<IssmDouble>* pSg_old, IssmDouble eartharea, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius,bool verboseconvolution,int loop);
    167                 void SealevelriseRotationalFeedback(Vector<IssmDouble>* pRSLgo_rot, Vector<IssmDouble>* pRSLg_old, IssmDouble* pIxz, IssmDouble* pIyz, IssmDouble* pIzz, IssmDouble eartharea,IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius);
    168                 void SealevelriseElastic(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, Vector<IssmDouble>* pSg_old, IssmDouble eartharea,IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz,int loop,int horiz);
    169                 IssmDouble SealevelriseOceanAverage(Vector<IssmDouble>* Sg,IssmDouble oceanarea);
     167                void SealevelriseNonEustatic(Vector<IssmDouble>* pSgo, Vector<IssmDouble>* pSg_old, IssmDouble eartharea, SealevelMasks* masks, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius,bool verboseconvolution,int loop);
     168                void SealevelriseRotationalFeedback(Vector<IssmDouble>* pRSLgo_rot, Vector<IssmDouble>* pRSLg_old, IssmDouble* pIxz, IssmDouble* pIyz, IssmDouble* pIzz, IssmDouble eartharea,SealevelMasks* masks, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius);
     169                void SealevelriseElastic(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, Vector<IssmDouble>* pSg_old, IssmDouble eartharea,SealevelMasks* masks, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz,int loop,int horiz);
     170                IssmDouble SealevelriseOceanAverage(Vector<IssmDouble>* Sg,SealevelMasks* masks, IssmDouble oceanarea);
    170171                #endif
    171172                void HydrologyEPLupdateDomainx(IssmDouble* pEplcount);

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

    r24379 r24938  
    1818#include "./Massfluxatgate.h"
    1919#include "./Misfit.h"
     20#include "./SealevelMasks.h"
    2021#include "./Nodalvalue.h"
    2122#include "./Numberedcostfunction.h"

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

    r24924 r24938  
    99class FemModel;
    1010class Parameters;
     11class SealevelMasks;
    1112template <class doubletype> class Matrix;
    1213template <class doubletype> class Vector;
    … …  
    5657void steric_core(FemModel* femmodel);
    5758void sealevelrise_core_geometry(FemModel* femmodel);
    58 Vector<IssmDouble>* sealevelrise_core_eustatic(FemModel* femmodel,IssmDouble* peartharea,IssmDouble* poceanarea);
    59 Vector<IssmDouble>* sealevelrise_core_noneustatic(FemModel* femmodel,Vector<IssmDouble>* RSLg_eustatic,IssmDouble eartharea,IssmDouble oceanarea);
    60 void sealevelrise_core_elastic(Vector<IssmDouble>** pU_radial, Vector<IssmDouble>** pU_north,Vector<IssmDouble>** pU_east,FemModel* femmodel,Vector<IssmDouble>* RSLg,IssmDouble eartharea);
     59SealevelMasks* sealevelrise_core_masks(FemModel* femmodel);
     60Vector<IssmDouble>* sealevelrise_core_eustatic(FemModel* femmodel,SealevelMasks* mask, IssmDouble* peartharea,IssmDouble* poceanarea);
     61Vector<IssmDouble>* sealevelrise_core_noneustatic(FemModel* femmodel,SealevelMasks* masks, Vector<IssmDouble>* RSLg_eustatic,IssmDouble eartharea,IssmDouble oceanarea);
     62void sealevelrise_core_elastic(Vector<IssmDouble>** pU_radial, Vector<IssmDouble>** pU_north,Vector<IssmDouble>** pU_east,FemModel* femmodel,Vector<IssmDouble>* RSLg,IssmDouble eartharea, SealevelMasks* masks);
    6163void sealevelrise_core_viscous(Vector<IssmDouble>** pU_gia,Vector<IssmDouble>** pN_gia,FemModel*  femmodel,Vector<IssmDouble>* RSLg);
    6264void sealevelrise_diagnostics(FemModel* femmodel,Vector<IssmDouble>* RSLg);

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

    r24934 r24938  
    8585        Vector<IssmDouble> *N_gia_rate= NULL;
    8686        Vector<IssmDouble> *U_gia_rate= NULL;
     87        SealevelMasks* masks=NULL;
    8788
    8889        /*parameters:*/
    … …  
    148149        if(modelid==earthid){
    149150
     151                /*call masks core: */
     152                masks=sealevelrise_core_masks(femmodel);
     153
    150154                /*call eustatic core  (generalized eustatic - Farrel and Clark, Eq 4, 1st, 3rd and 4rd terms on the RHS) */
    151                 RSLg_eustatic=sealevelrise_core_eustatic(femmodel,&eartharea,&oceanarea);
     155                RSLg_eustatic=sealevelrise_core_eustatic(femmodel,masks,&eartharea,&oceanarea);
    152156
    153157                /*call non-eustatic core (ocean loading tems  - 2nd and 5th terms on the RHS of Farrel and Clark) */
    154                 RSLg=sealevelrise_core_noneustatic(femmodel,RSLg_eustatic,eartharea,oceanarea);
     158                RSLg=sealevelrise_core_noneustatic(femmodel,masks,RSLg_eustatic,eartharea,oceanarea);
    155159
    156160                /*compute other elastic geodetic signatures, such as components of 3-D crustal motion: */
    157                 sealevelrise_core_elastic(&U_esa,&U_north_esa,&U_east_esa,femmodel,RSLg,eartharea);
     161                sealevelrise_core_elastic(&U_esa,&U_north_esa,&U_east_esa,femmodel,RSLg,eartharea,masks);
    158162
    159163                /*compute viscosus (GIA) geodetic signatures:*/
    … …  
    223227        delete N_gia_rate;
    224228        delete U_gia_rate;
     229        //delete masks;
    225230        /*}}}*/
    226231
    … …  
    307312}
    308313/*}}}*/
     314SealevelMasks* sealevelrise_core_masks(FemModel* femmodel) {  /*{{{*/
     315
     316        if(VerboseSolution()) _printf0_("         computing sea level masks\n");
     317       
     318        /*initialize SealevelMasks structure: */
     319        SealevelMasks* masks=new SealevelMasks(femmodel->elements->Size());
     320
     321        /*go through elements and fill the masks: */
     322        for (int i=0;i<femmodel->elements->Size();i++){
     323
     324                Element*   element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
     325                element->SetSealevelMasks(masks);
     326        }
     327
     328        return masks;
     329}/*}}}*/
    309330void sealevelrise_core_geometry(FemModel* femmodel) {  /*{{{*/
    310331
    … …  
    333354
    334355}/*}}}*/
    335 Vector<IssmDouble>* sealevelrise_core_eustatic(FemModel* femmodel,IssmDouble* peartharea,IssmDouble* poceanarea){ /*{{{*/
     356Vector<IssmDouble>* sealevelrise_core_eustatic(FemModel* femmodel,SealevelMasks* masks, IssmDouble* peartharea,IssmDouble* poceanarea){ /*{{{*/
    336357
    337358        /*Eustatic core of the SLR solution (terms that are constant with respect to sea-level)*/
    … …  
    367388
    368389        /*call the eustatic main module: */
    369         femmodel->SealevelriseEustatic(RSLgi,&eartharea, &oceanarea,&eustatic, latitude, longitude, radius,loop); //this computes
     390        femmodel->SealevelriseEustatic(RSLgi,&eartharea, &oceanarea,&eustatic, masks, latitude, longitude, radius,loop); //this computes
    370391
    371392        /*we need to average RSLgi over the ocean: RHS term  4 in Eq.4 of Farrel and clarke. Only the elements can do that: */
    372         RSLgi_oceanaverage=femmodel->SealevelriseOceanAverage(RSLgi,oceanarea);
     393        RSLgi_oceanaverage=femmodel->SealevelriseOceanAverage(RSLgi,masks, oceanarea);
    373394
    374395        /*RSLg is the sum of the pure eustatic component (term 3) and the contribution from the perturbation to the graviation potential due to the
    … …  
    389410        return RSLgi;
    390411}/*}}}*/
    391 Vector<IssmDouble>* sealevelrise_core_noneustatic(FemModel* femmodel,Vector<IssmDouble>* RSLg_eustatic,IssmDouble eartharea,IssmDouble oceanarea){ /*{{{*/
     412Vector<IssmDouble>* sealevelrise_core_noneustatic(FemModel* femmodel, SealevelMasks* masks, Vector<IssmDouble>* RSLg_eustatic,IssmDouble eartharea,IssmDouble oceanarea){ /*{{{*/
    392413
    393414        /*sealevelrise_core_noneustatic.cpp //this computes the contributions from Eq.4 of Farrel and Clarke, rhs terms 2 and 5.
    … …  
    458479
    459480                /*call the non eustatic module: */
    460                 femmodel->SealevelriseNonEustatic(RSLgo, RSLg_old,  eartharea,latitude,longitude, radius,verboseconvolution,loop);
     481                femmodel->SealevelriseNonEustatic(RSLgo, RSLg_old,  eartharea, masks, latitude,longitude, radius,verboseconvolution,loop);
    461482
    462483                /*assemble solution vector: */
    … …  
    467488                        /*call rotational feedback  module: */
    468489                        RSLgo_rot = new Vector<IssmDouble>(gsize); RSLgo_rot->Assemble();
    469                         femmodel->SealevelriseRotationalFeedback(RSLgo_rot,RSLg_old,&Ixz,&Iyz,&Izz,eartharea,latitude,longitude,radius);
     490                        femmodel->SealevelriseRotationalFeedback(RSLgo_rot,RSLg_old,&Ixz,&Iyz,&Izz, eartharea, masks, latitude,longitude,radius);
    470491                        RSLgo_rot->Assemble();
    471492
    … …  
    479500
    480501                /*we need to average RSLgo over the ocean: RHS term  5 in Eq.4 of Farrel and clarke. Only the elements can do that: */
    481                 RSLgo_oceanaverage=femmodel->SealevelriseOceanAverage(RSLgo,oceanarea);
     502                RSLgo_oceanaverage=femmodel->SealevelriseOceanAverage(RSLgo,masks, oceanarea);
    482503
    483504                /*RSLg is the sum of the eustatic term, and the ocean terms: */
    … …  
    516537        return RSLg;
    517538} /*}}}*/
    518 void sealevelrise_core_elastic(Vector<IssmDouble>** pU_esa, Vector<IssmDouble>** pU_north_esa,Vector<IssmDouble>** pU_east_esa,FemModel* femmodel,Vector<IssmDouble>* RSLg,IssmDouble eartharea){ /*{{{*/
     539void sealevelrise_core_elastic(Vector<IssmDouble>** pU_esa, Vector<IssmDouble>** pU_north_esa,Vector<IssmDouble>** pU_east_esa,FemModel* femmodel,Vector<IssmDouble>* RSLg,IssmDouble eartharea, SealevelMasks* masks){ /*{{{*/
    519540
    520541        Vector<IssmDouble> *U_esa  = NULL;
    … …  
    556577
    557578        /*call the elastic main modlule:*/
    558         femmodel->SealevelriseElastic(U_esa,U_north_esa,U_east_esa,RSLg,eartharea,latitude,longitude,radius,xx,yy,zz,loop,horiz);
     579        femmodel->SealevelriseElastic(U_esa,U_north_esa,U_east_esa,RSLg,eartharea, masks, latitude,longitude,radius,xx,yy,zz,loop,horiz);
    559580
    560581        /*Assign output pointers:*/
Note: See TracChangeset for help on using the changeset viewer.