Changeset 26110


Ignore:
Timestamp:
03/18/21 20:17:37 (4 years ago)
Author:
Eric.Larour
Message:

CHG: new synthesis of the GRD module in sealevelchange_core using the new approach
for convolutions.

Location:
issm/trunk-jpl
Files:
2 added
26 edited

Legend:

Unmodified
Added
Removed

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

    r26099 r26110  
    9797        ./classes/Hook.cpp \
    9898        ./classes/Radar.cpp \
     99        ./classes/BarystaticContributions.cpp \
    99100        ./classes/ExternalResults/Results.cpp \
    100101        ./classes/Elements/Element.cpp \

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

    r26099 r26110  
    183183                xDelete<IssmDouble>(partitionocean);
    184184        }
    185 
     185        /*New optimized code:*/
     186        BarystaticContributions* barystaticcontributions=new BarystaticContributions(iomodel);
     187        parameters->AddObject(new GenericParam<BarystaticContributions*>(barystaticcontributions,BarystaticContributionsEnum));
    186188       
    187189        /*Deal with external multi-model ensembles: {{{*/

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

    r26099 r26110  
    3838class ElementMatrix;
    3939class ElementVector;
     40class BarystaticContributions;
    4041/*}}}*/
    4142
    … …  
    389390                virtual void          DeformationFromSurfaceLoads(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg,SealevelMasks* masks)=0;
    390391                virtual void       GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt,Matlitho* litho, IssmDouble* x,IssmDouble* y)=0;
     392
     393                virtual void       SealevelchangeGeometryOptim(IssmDouble* lat,IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze)=0;
     394                virtual void       SealevelchangeBarystaticLoads(Vector<IssmDouble>* loads, BarystaticContributions* barycontrib, SealevelMasks* masks)=0;
     395                virtual void       SealevelchangeInitialConvolution(Vector<IssmDouble>* sealevelloads,Vector<IssmDouble>* oceanareas,IssmDouble* loads, SealevelMasks* masks)=0;
     396                virtual void       SealevelchangeOceanConvolution(Vector<IssmDouble>* newsealevelloads, IssmDouble* sealevelloads, IssmDouble* loads, SealevelMasks* masks)=0;
     397                virtual void       OceanAverageOptim(IssmDouble* poceanaverage, IssmDouble* poceanarea, IssmDouble* Sg, SealevelMasks* masks)=0;
    391398                #endif
    392399

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

    r26099 r26110  
    227227                void    DeformationFromSurfaceLoads(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet!");};
    228228                void           GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt,Matlitho* litho, IssmDouble* x,IssmDouble* y){_error_("not implemented yet");};
     229
     230                void       SealevelchangeGeometryOptim(IssmDouble* lat,IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze){_error_("not implemented yet");};
     231                void       SealevelchangeBarystaticLoads(Vector<IssmDouble>* loads, BarystaticContributions* barycontrib, SealevelMasks* masks){_error_("not implemented yet");};
     232                void       SealevelchangeInitialConvolution(Vector<IssmDouble>* sealevelloads,Vector<IssmDouble>* oceanareas,IssmDouble* loads, SealevelMasks* masks){_error_("not implemented yet");};
     233                void       SealevelchangeOceanConvolution(Vector<IssmDouble>* newsealevelloads, IssmDouble* sealevelloads, IssmDouble* loads, SealevelMasks* masks){_error_("not implemented yet");};
     234                void       OceanAverageOptim(IssmDouble* poceanaverage, IssmDouble* poceanarea, IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet");};
    229235                #endif
    230236

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

    r26099 r26110  
    182182                IssmDouble    OceanAverage(IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet!");};
    183183                void        GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt,Matlitho* litho, IssmDouble* x,IssmDouble* y){_error_("not implemented yet");};
     184
     185                void       SealevelchangeGeometryOptim(IssmDouble* lat,IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze){_error_("not implemented yet");};
     186                void       SealevelchangeBarystaticLoads(Vector<IssmDouble>* loads, BarystaticContributions* barycontrib, SealevelMasks* masks){_error_("not implemented yet");};
     187                void       SealevelchangeInitialConvolution(Vector<IssmDouble>* sealevelloads,Vector<IssmDouble>* oceanareas,IssmDouble* loads, SealevelMasks* masks){_error_("not implemented yet");};
     188                void       SealevelchangeOceanConvolution(Vector<IssmDouble>* newsealevelloads, IssmDouble* sealevelloads, IssmDouble* loads, SealevelMasks* masks){_error_("not implemented yet");};
     189                void       OceanAverageOptim(IssmDouble* poceanaverage, IssmDouble* poceanarea, IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet");};
     190
    184191#endif
    185192

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

    r26099 r26110  
    189189                IssmDouble    OceanAverage(IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet!");};
    190190                void        GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt, Matlitho* litho, IssmDouble* x,IssmDouble* y){_error_("not implemented yet");};
     191                void       SealevelchangeGeometryOptim(IssmDouble* lat,IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze){_error_("not implemented yet");};
     192                void       SealevelchangeBarystaticLoads(Vector<IssmDouble>* loads, BarystaticContributions* barycontrib, SealevelMasks* masks){_error_("not implemented yet");};
     193                void       SealevelchangeInitialConvolution(Vector<IssmDouble>* sealevelloads,Vector<IssmDouble>* oceanareas,IssmDouble* loads, SealevelMasks* masks){_error_("not implemented yet");};
     194                void       SealevelchangeOceanConvolution(Vector<IssmDouble>* newsealevelloads, IssmDouble* sealevelloads, IssmDouble* loads, SealevelMasks* masks){_error_("not implemented yet");};
     195                void       OceanAverageOptim(IssmDouble* poceanaverage, IssmDouble* poceanarea, IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet");};
     196
    191197#endif
    192198

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

    r26107 r26110  
    54325432#endif
    54335433#ifdef _HAVE_SEALEVELCHANGE_
     5434//old code
    54345435IssmDouble Tria::OceanAverage(IssmDouble* Sg, SealevelMasks* masks){ /*{{{*/
    54355436
    … …  
    54515452}
    54525453/*}}}*/
    5453 void    Tria::SealevelchangeMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks){/*{{{*/
     5454void          Tria::SealevelchangeMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks){/*{{{*/
    54545455        /*early return if we are not on an ice cap OR ocean:*/
    54555456        if(!masks->isiceonly[this->lid] && !masks->isoceanin[this->lid]){
    … …  
    55455546        return;
    55465547}/*}}}*/
    5547 void    Tria::SetSealevelMasks(SealevelMasks* masks){ /*{{{*/
     5548void       Tria::SetSealevelMasks(SealevelMasks* masks){ /*{{{*/
    55485549
    55495550        masks->isiceonly[this->lid]=this->IsIceOnlyInElement();
    … …  
    55585559        if ((gr_input->GetInputMin())<0) masks->notfullygrounded[this->lid]=true;
    55595560        else masks->notfullygrounded[this->lid]=false;
    5560 
    5561 }
    5562 /*}}}*/
    5563 void    Tria::SealevelchangeGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze){ /*{{{*/
     5561}
     5562/*}}}*/
     5563void       Tria::SealevelchangeGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze){ /*{{{*/
    55645564
    55655565        /*diverse:*/
    … …  
    59515951}
    59525952/*}}}*/
    5953 void    Tria::SealevelchangeBarystaticBottomPressure(IssmDouble* Sgi,SealevelMasks* masks){ /*{{{*/
     5953void       Tria::SealevelchangeBarystaticBottomPressure(IssmDouble* Sgi,SealevelMasks* masks){ /*{{{*/
    59545954
    59555955        /*diverse:*/
    … …  
    60016001}
    60026002/*}}}*/
    6003 void    Tria::SealevelchangeSal(IssmDouble* Sgo,IssmDouble* Sg_old, SealevelMasks* masks){ /*{{{*/
     6003void       Tria::SealevelchangeSal(IssmDouble* Sgo,IssmDouble* Sg_old, SealevelMasks* masks){ /*{{{*/
    60046004
    60056005        /*diverse:*/
    … …  
    60436043}
    60446044/*}}}*/
    6045 void    Tria::DeformationFromSurfaceLoads(IssmDouble* Up, IssmDouble* North ,IssmDouble* East,IssmDouble* Sg, SealevelMasks* masks){ /*{{{*/
     6045void       Tria::DeformationFromSurfaceLoads(IssmDouble* Up, IssmDouble* North ,IssmDouble* East,IssmDouble* Sg, SealevelMasks* masks){ /*{{{*/
    60466046
    60476047        /*diverse:*/
    … …  
    61326132}
    61336133/*}}}*/
    6134 void    Tria::GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt, Matlitho* litho, IssmDouble* x, IssmDouble* y){/*{{{*/
     6134void       Tria::GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt, Matlitho* litho, IssmDouble* x, IssmDouble* y){/*{{{*/
    61356135
    61366136        IssmDouble xyz_list[NUMVERTICES][3];
    … …  
    62386238/*}}}*/
    62396239
    6240 //new logic
    6241 void    Tria::SealevelchangeGeometryOptim(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze){ /*{{{*/
     6240//new code
     6241void       Tria::SealevelchangeGeometryOptim(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze){ /*{{{*/
    62426242
    62436243        /*diverse:*/
    … …  
    64286428}
    64296429/*}}}*/
    6430 void    Tria::SealevelchangeBarystaticLoads(IssmDouble* barystatic_contribution,IssmDouble* localloads, SealevelMasks* masks, Matrix<IssmDouble>* barystatic_contribution_onpartition, IssmDouble* partition, IssmDouble oceanarea){ /*{{{*/
     6430void       Tria::SealevelchangeBarystaticLoads(Vector<IssmDouble>* loads, BarystaticContributions* barycontrib, SealevelMasks* masks){ /*{{{*/
    64316431
    64326432        /*diverse:*/
    … …  
    64366436        IssmDouble I,W,BP;  //change in ice thickness, water column or bottom pressure (Farrel and Clarke, Equ. 4)
    64376437        bool notfullygrounded=false;
    6438         bool scaleoceanarea= false;
    64396438        bool computerigid= false;
    64406439        int  glfraction=1;
    … …  
    64696468                        constant=0; this->AddInput(SealevelBarystaticMaskEnum,&constant,P0Enum);
    64706469                #endif
    6471                 /*barystatic_contribution[0]+=0;
    6472                 barystatic_contribution[1]+=0;
    6473                 barystatic_contribution[2]+=0;*/
    64746470                return;
    64756471                }
    … …  
    64836479                        this->AddInput(SealevelBarystaticMaskEnum,&constant,P0Enum);
    64846480                        #endif
    6485                         /*barystatic_contribution[0]+=0;
    6486                         barystatic_contribution[1]+=0;
    6487                         barystatic_contribution[2]+=0;*/
    64886481                        return;
    64896482                }
    … …  
    64946487        if(!isice  && !ishydro){
    64956488                if(!masks->isoceanin[this->lid]){
    6496                         /*barystatic_contribution[0]+=0;
    6497                         barystatic_contribution[1]+=0;
    6498                         barystatic_contribution[2]+=0;*/
    64996489                        return;
    65006490                }
    … …  
    65096499        rho_water=FindParam(MaterialsRhoSeawaterEnum);
    65106500        this->parameters->FindParam(&computerigid,SolidearthSettingsRigidEnum);
    6511         this->parameters->FindParam(&scaleoceanarea,SolidearthSettingsOceanAreaScalingEnum);
    65126501        this->parameters->FindParam(&glfraction,SolidearthSettingsGlfractionEnum);
    65136502        this->parameters->FindParam(&npartice,SolidearthNpartIceEnum);
    … …  
    65186507
    65196508        /*Deal with ice loads if we are on grounded ice:*/
    6520         if(masks->isiceonly[this->lid] && !masks->isfullyfloating[this->lid]){ /*{{{*/
    6521 
    6522                 /*Compute fraction of the element that is grounded: */
     6509        if(masks->isiceonly[this->lid] && !masks->isfullyfloating[this->lid]){
     6510
     6511                /*Compute fraction of the element that is grounded: {{{*/
    65236512                if(notfullygrounded){
    65246513                        IssmDouble xyz_list[NUMVERTICES][3];
    … …  
    65326521                }
    65336522                else phi_ice=1.0;
     6523                /*}}}*/
    65346524
    65356525                /*Inform mask: */
    … …  
    65696559                /*}}}*/
    65706560
    6571                 /*Compute barystatic contribution:*/
    6572                 _assert_(oceanarea>0.);
    6573                 if(scaleoceanarea) oceanarea=3.619e+14; // use true ocean area, m^2
    6574                 bslcice = rho_ice*area*I/(oceanarea*rho_water);
     6561                /*Compute barystatic contribution in kg:*/
     6562                bslcice = rho_ice*area*I;
    65756563                _assert_(!xIsNan<IssmDouble>(bslcice));
    65766564
    65776565                /*Transfer thickness change into kg/m^2:*/
    65786566                I=I*rho_ice*phi_ice;
    6579         } /*}}}*/
     6567        }
    65806568
    65816569        /*Deal with water loads if we are on ground:*/
    … …  
    65876575                deltathickness_input->GetInputAverage(&W);
    65886576
    6589                 /*Compute barystatic component:*/
    6590                 _assert_(oceanarea>0.);
    6591                 if(scaleoceanarea) oceanarea=3.619e+14; // use true ocean area, m^2
    6592                 bslchydro = rho_freshwater*area*phi_water*W/(oceanarea*rho_water);
     6577                /*Compute barystatic component in kg:*/
     6578                bslchydro = rho_freshwater*area*phi_water*W;
    65936579                _assert_(!xIsNan<IssmDouble>(bslchydro));
    65946580
    … …  
    66056591                bottompressure_change_input->GetInputAverage(&BP);
    66066592               
    6607                 /*Compute barystatic component:*/
    6608                 bslcbp = rho_water*area*BP/(oceanarea*rho_water);
     6593                /*Compute barystatic component in kg:*/
     6594                bslcbp = rho_water*area*BP;
    66096595
    66106596                /*convert from m to kg/m^2:*/
    … …  
    66136599
    66146600        /*Plug all loads into total load vector:*/
    6615         localloads[this->lid]+=I+W+BP;
     6601        loads->SetValue(this->sid,I+W+BP,INS_VAL);
    66166602       
    6617         /*Plug bslcice into barystatic contribution vector:*/
    6618         if(barystatic_contribution_onpartition){
    6619                 int idi=reCast<int>(partition[this->Sid()])+1;
    6620                 int idj=0;
    6621                 idj=0;barystatic_contribution_onpartition->SetValues(1,&idi,1,&idj,&bslcice,ADD_VAL);
    6622                 idj=1;barystatic_contribution_onpartition->SetValues(1,&idi,1,&idj,&bslchydro,ADD_VAL);
    6623                 idj=2;barystatic_contribution_onpartition->SetValues(1,&idi,1,&idj,&bslcbp,ADD_VAL);
    6624         }
    6625 
    6626         barystatic_contribution[0]+=bslcice;
    6627         barystatic_contribution[1]+=bslchydro;
    6628         barystatic_contribution[2]+=bslcbp;
    6629 
    6630 }
    6631 /*}}}*/
    6632 IssmDouble Tria::SealevelchangeConvolution(IssmDouble* loads){ /*{{{*/
     6603        /*Keep track of barystatic contributions:*/
     6604        barycontrib->Set(this->Sid(),bslcice,bslchydro,bslcbp);
     6605
     6606}
     6607/*}}}*/
     6608void       Tria::SealevelchangeInitialConvolution(Vector<IssmDouble>* sealevelloads,Vector<IssmDouble>* oceanareas,IssmDouble* loads, SealevelMasks* masks){ /*{{{*/
    66336609
    66346610        /*sal green function:*/
    66356611        IssmDouble* G=NULL;
    6636         IssmDouble Sealevel[NUMVERTICES]={0,0,0};
     6612        IssmDouble SealevelGRD[NUMVERTICES]={0,0,0};
     6613        IssmDouble oceanaverage,oceanarea=0;
    66376614
    66386615        bool sal = false;
    … …  
    66476624                for(int i=0;i<NUMVERTICES;i++) {
    66486625                        for (int e=0;e<nel;e++){
    6649                                 Sealevel[i]+=G[i*nel+e]*loads[e];
    6650                         }
    6651                 }
    6652 
    6653                 this->AddInput(SealevelRSLEnum,Sealevel,P1Enum);
    6654         }
    6655 
     6626                                SealevelGRD[i]+=G[i*nel+e]*loads[e];
     6627                        }
     6628                }
     6629        }
     6630
     6631        /*compute ocean average over element:*/
     6632        OceanAverageOptim(&oceanaverage,&oceanarea,SealevelGRD,masks);
     6633       
     6634        /*add ocean average in the global sealevelloads vector:*/
     6635        sealevelloads->SetValue(this->sid,oceanaverage,INS_VAL);
     6636        oceanareas->SetValue(this->sid,oceanarea,INS_VAL);
     6637       
     6638        return;
    66566639} /*}}}*/
    6657 
     6640void       Tria::SealevelchangeOceanConvolution(Vector<IssmDouble>* newsealevelloads, IssmDouble* sealevelloads, IssmDouble* loads, SealevelMasks* masks){ /*{{{*/
     6641
     6642        bool converged=false;
     6643        IssmDouble SealevelGRD[3]={0,0,0};
     6644        IssmDouble oceanaverage,oceanarea=0;
     6645        int nel;
     6646        bool sal = false;
     6647        IssmDouble* G=NULL;
     6648        int size;
     6649       
     6650        this->parameters->FindParam(&nel,MeshNumberofelementsEnum);
     6651        this->parameters->FindParam(&sal,SolidearthSettingsRigidEnum);
     6652               
     6653        if(sal){
     6654                this->inputs->GetArrayPtr(SealevelchangeGEnum,this->lid,&G,&size);
     6655
     6656                for(int i=0;i<NUMVERTICES;i++) {
     6657                        for (int e=0;e<nel;e++){
     6658                                SealevelGRD[i]+=G[i*nel+e]*(sealevelloads[e]+loads[e]);
     6659                        }
     6660                }
     6661        }
     6662        OceanAverageOptim(&oceanaverage,&oceanarea,SealevelGRD,masks);
     6663        newsealevelloads->SetValue(this->sid,oceanaverage,INS_VAL);
     6664
     6665} /*}}}*/
     6666void       Tria::OceanAverageOptim(IssmDouble* poceanaverage, IssmDouble* poceanarea, IssmDouble* Sg, SealevelMasks* masks){ /*{{{*/
     6667
     6668        IssmDouble phi=1.0;
     6669        bool iscoastline=false;
     6670        IssmDouble area;
     6671        IssmDouble Sg_avg=0; //output
     6672       
     6673        /*Do we have an ocean?:*/
     6674        if(!masks->isoceanin[this->lid]){
     6675                *poceanarea=0;
     6676                *poceanaverage=0;
     6677        }
     6678
     6679        /*Do we have a coastline?:*/
     6680        if(!masks->isfullyfloating[this->lid])iscoastline=true;
     6681
     6682        if(iscoastline){
     6683                IssmDouble xyz_list[NUMVERTICES][3];
     6684                ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     6685                phi=1.0-this->GetGroundedPortion(&xyz_list[0][0]);
     6686        }
     6687
     6688        /*Get area of element:*/
     6689        this->Element::GetInputValue(&area,AreaEnum);
     6690
     6691        /*Average over ocean if there is no coastline, area of the ocean
     6692         *is the areaa of the element:*/
     6693        if(!iscoastline){
     6694
     6695                /*Average Sg over vertices:*/
     6696                for(int i=0;i<NUMVERTICES;i++) Sg_avg+=Sg[i]/NUMVERTICES;
     6697
     6698                *poceanaverage=Sg_avg;
     6699                *poceanarea=area;
     6700                return;
     6701        }
     6702
     6703        /*Average over  the ocean only if there is a coastline. Area of the
     6704         * ocean will be the fraction times the area of the element:*/
     6705        area=phi*area;
     6706   
     6707        IssmDouble total_weight=0;
     6708        bool mainlyfloating = true;
     6709        int         point1;
     6710        IssmDouble  fraction1,fraction2;
     6711
     6712        /*Recover portion of element that is grounded*/
     6713        this->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating);
     6714        //!mainlyfloating so that the integration is done on the ocean (and not the grounded) part
     6715        Gauss* gauss = this->NewGauss(point1,fraction1,fraction2,!mainlyfloating,2);
     6716
     6717        /* Start  looping on the number of gaussian points and average over these gaussian points: */
     6718        total_weight=0;
     6719        Sg_avg=0;
     6720        while(gauss->next()){
     6721                IssmDouble Sg_gauss=0;
     6722                TriaRef::GetInputValue(&Sg_gauss, Sg, gauss,P1Enum);
     6723                Sg_avg+=Sg_gauss*gauss->weight;
     6724                total_weight+=gauss->weight;
     6725        }
     6726        Sg_avg=Sg_avg/total_weight;
     6727        delete gauss;
     6728
     6729        *poceanaverage=Sg_avg;
     6730        *poceanarea=area;
     6731        return;
     6732
     6733}
     6734/*}}}*/
    66586735#endif
    66596736

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

    r26099 r26110  
    162162                #endif
    163163                #ifdef _HAVE_SEALEVELCHANGE_
     164                void       SealevelchangeMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old,SealevelMasks* masks);
     165                void       SealevelchangeGeometry(IssmDouble* lat, IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze);
     166                IssmDouble SealevelchangeBarystaticIce(IssmDouble* Sgi, SealevelMasks* masks, Vector<IssmDouble>* barystatic_contribution, IssmDouble* partition, IssmDouble oceanarea);
     167                IssmDouble SealevelchangeBarystaticHydro(IssmDouble* Sgi, SealevelMasks* masks, Vector<IssmDouble>* barystatic_contribution, IssmDouble* partition,IssmDouble oceanarea);
     168                void       SealevelchangeBarystaticBottomPressure(IssmDouble* Sgi,SealevelMasks* masks);
     169                void       SealevelchangeSal(IssmDouble* Sgo,IssmDouble* Sg_old,SealevelMasks* masks);
     170                void       DeformationFromSurfaceLoads(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg,SealevelMasks* masks);
     171                void       GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt,Matlitho* litho, IssmDouble* x,IssmDouble* y);
    164172                IssmDouble OceanAverage(IssmDouble* Sg, SealevelMasks* masks);
    165                 void    SealevelchangeMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old,SealevelMasks* masks);
    166                 void    SetSealevelMasks(SealevelMasks* masks);
    167                 void    SealevelchangeGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze);
    168                 void    SealevelchangeGeometryOptim(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze);
    169 
    170                 IssmDouble    SealevelchangeBarystaticIce(IssmDouble* Sgi, SealevelMasks* masks, Vector<IssmDouble>* barystatic_contribution, IssmDouble* partition, IssmDouble oceanarea);
    171                 IssmDouble    SealevelchangeBarystaticHydro(IssmDouble* Sgi, SealevelMasks* masks, Vector<IssmDouble>* barystatic_contribution, IssmDouble* partition,IssmDouble oceanarea);
    172                 void    SealevelchangeBarystaticBottomPressure(IssmDouble* Sgi,SealevelMasks* masks);
    173                 void    SealevelchangeSal(IssmDouble* Sgo,IssmDouble* Sg_old,SealevelMasks* masks);
    174                 void    DeformationFromSurfaceLoads(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg,SealevelMasks* masks);
    175                 void    GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt,Matlitho* litho, IssmDouble* x,IssmDouble* y);
    176                 void    SealevelchangeBarystaticLoads(IssmDouble* barystatic_contribution,IssmDouble* localloads, SealevelMasks* masks, Matrix<IssmDouble>* barystatic_contribution_onpartition, IssmDouble* partition, IssmDouble oceanarea);
    177                 IssmDouble SealevelchangeConvolution(IssmDouble* loads);
     173                void       SetSealevelMasks(SealevelMasks* masks);
     174               
     175                void       SealevelchangeGeometryOptim(IssmDouble* lat,IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze);
     176                void       SealevelchangeBarystaticLoads(Vector<IssmDouble>* loads, BarystaticContributions* barycontrib, SealevelMasks* masks);
     177                void       SealevelchangeInitialConvolution(Vector<IssmDouble>* sealevelloads,Vector<IssmDouble>* oceanareas,IssmDouble* loads, SealevelMasks* masks);
     178                void       SealevelchangeOceanConvolution(Vector<IssmDouble>* newsealevelloads, IssmDouble* sealevelloads, IssmDouble* loads, SealevelMasks* masks);
     179                void       OceanAverageOptim(IssmDouble* poceanaverage, IssmDouble* poceanarea, IssmDouble* Sg, SealevelMasks* masks);
    178180                #endif
    179181                /*}}}*/

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

    r26099 r26110  
    47494749#endif
    47504750#ifdef _HAVE_SEALEVELCHANGE_
    4751 void FemModel::SealevelchangeBarystatic(Vector<IssmDouble>* pRSLgi, IssmDouble* poceanarea, IssmDouble* pbslc,IssmDouble* pbslcice, IssmDouble* pbslchydro, IssmDouble** pbslcice_partition,IssmDouble** pbslchydro_partition,SealevelMasks* masks) { /*{{{*/
    4752 
    4753         /*serialized vectors:*/
    4754         IssmDouble  bslcice       = 0.;
    4755         IssmDouble  bslcice_cpu   = 0.;
    4756         IssmDouble  bslchydro       = 0.;
    4757         IssmDouble  bslchydro_cpu   = 0.;
    4758         IssmDouble  area      = 0.;
    4759         IssmDouble  oceanarea      = 0.;
    4760         IssmDouble  oceanarea_cpu  = 0.;
    4761         int bp_compute_fingerprints= 0;
    4762         bool isoceantransport=false;
    4763 
    4764         Vector<IssmDouble>* bslcice_partition=NULL;
    4765         IssmDouble* bslcice_partition_serial=NULL;
    4766         IssmDouble* partitionice=NULL;
    4767         int npartice,nel;
    4768 
    4769         Vector<IssmDouble>* bslchydro_partition=NULL;
    4770         IssmDouble* bslchydro_partition_serial=NULL;
    4771         IssmDouble* partitionhydro=NULL;
    4772         bool istws=0;
    4773         int nparthydro;
    4774                
    4775         int npartocean;
    4776         Vector<IssmDouble>* bslcocean_partition=NULL;
    4777         IssmDouble* bslcocean_partition_serial=NULL;
    4778         IssmDouble* partitionocean=NULL;
    4779 
    4780    /*Initialize temporary vector that will be used to sum barystatic components
    4781     * on all local elements, prior to assembly:*/
    4782         int gsize = this->nodes->NumberOfDofs(GsetEnum);
    4783         IssmDouble* RSLgi=xNewZeroInit<IssmDouble>(gsize);
    4784         int* indices=xNew<int>(gsize);
    4785    for(int i=0;i<gsize;i++) indices[i]=i;
    4786 
    4787         /*First, figure out the area of the ocean, which is needed to compute the barystatic component: */
    4788         int i = -1;
    4789         for(Object* & object : this->elements->objects){
    4790                 i +=1;
    4791                 Element* element = xDynamicCast<Element*>(object);
    4792                 element->GetInputValue(&area,AreaEnum);
    4793                 if (masks->isoceanin[i]) oceanarea_cpu += area;
    4794         }
    4795         ISSM_MPI_Reduce (&oceanarea_cpu,&oceanarea,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
    4796         ISSM_MPI_Bcast(&oceanarea,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
    4797         _assert_(oceanarea>0.);
    4798 
    4799         /*Initialize partition vectors to retrieve barystatic contributions: */
    4800         this->parameters->FindParam(&npartice,SolidearthNpartIceEnum);
    4801         if(npartice){
    4802                 this->parameters->FindParam(&partitionice,&nel,NULL,SolidearthPartitionIceEnum);
    4803                 bslcice_partition= new Vector<IssmDouble>(npartice);
    4804         }
    4805 
    4806         this->parameters->FindParam(&nparthydro,SolidearthNpartHydroEnum);
    4807         if(nparthydro){
    4808                 this->parameters->FindParam(&partitionhydro,&nel,NULL,SolidearthPartitionHydroEnum);
    4809                 bslchydro_partition= new Vector<IssmDouble>(nparthydro);
    4810         }
    4811 
    4812         this->parameters->FindParam(&npartocean,SolidearthNpartOceanEnum);
    4813         if(npartocean){
    4814                 this->parameters->FindParam(&partitionocean,&nel,NULL,SolidearthPartitionOceanEnum);
    4815                 bslchydro_partition= new Vector<IssmDouble>(npartocean);
    4816         }
    4817         /*For later:
    4818         npartbarystatic=npartice;
    4819         if(nparthydro>npartbarystatic)npartbarystatic=nparthydro;
    4820         if(npartocean>npartbarystatic)npartbarystatic=npartocean;
    4821         bslc_partition=new Matrix(IssmDouble>(npartbarystatic,3);
    4822 
    4823         bslc_cpu[0]=0; bslc_cpu[1]=0; bslc_cpu[2]=0;
    4824         for(Object* & object : this->elements->objects){
    4825                 Element* element = xDynamicCast<Element*>(object);
    4826                 element->SealevelchangeBarystaticLoads(&bslc_cpu[0], localloads,masks, bslcice_partition,partitionice,oceanarea);
    4827         }
    4828         MPI Bcast localloads -> loads
    4829 
    4830         for(Object* & object : this->elements->objects){
    4831                 Element* element = xDynamicCast<Element*>(object);
    4832                 element->SealevelchangeConvolution(loads);
    4833         }
    4834         */
    4835 
    4836 
    4837 
    4838 
    4839 
    4840         /*Call the barystatic sea level change core for ice : */
    4841         bslcice_cpu=0;
    4842         for(Object* & object : this->elements->objects){
    4843                 Element* element = xDynamicCast<Element*>(object);
    4844                 bslcice_cpu+=element->SealevelchangeBarystaticIce(RSLgi,masks, bslcice_partition,partitionice,oceanarea);
    4845         }
    4846 
    4847         /*Call the barystatic sea level change core for hydro: */
    4848         bslchydro_cpu=0; //make sure to initialize this, so we have a total barystatic contribution computed at 0.
    4849         this->parameters->FindParam(&istws,TransientIshydrologyEnum);
    4850         if(istws){
    4851                 for(int i=0;i<elements->Size();i++){
    4852                         Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
    4853                         bslchydro_cpu+=element->SealevelchangeBarystaticHydro(RSLgi,masks, bslchydro_partition,partitionhydro,oceanarea);
    4854                 }
    4855         }
    4856 
    4857         /*Call the barystatic sea level change core for bottom pressures: */
    4858         this->parameters->FindParam(&bp_compute_fingerprints,SolidearthSettingsComputeBpGrdEnum);
    4859         this->parameters->FindParam(&isoceantransport,TransientIsoceantransportEnum);
    4860         if(bp_compute_fingerprints && isoceantransport){
    4861                 for(int i=0;i<elements->Size();i++){
    4862                         Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
    4863                         element->SealevelchangeBarystaticBottomPressure(RSLgi,masks);
    4864                 }
    4865         }
    4866 
    4867         /*Plug values once and assemble: */
    4868         pRSLgi->SetValues(gsize,indices,RSLgi,ADD_VAL);
    4869         pRSLgi->Assemble();
    4870 
    4871         /*Sum all barystatic components from all cpus:*/
    4872         ISSM_MPI_Reduce (&bslcice_cpu,&bslcice,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
    4873         ISSM_MPI_Bcast(&bslcice,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
    4874         _assert_(!xIsNan<IssmDouble>(bslcice));
    4875 
    4876         ISSM_MPI_Reduce (&bslchydro_cpu,&bslchydro,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
    4877         ISSM_MPI_Bcast(&bslchydro,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
    4878         _assert_(!xIsNan<IssmDouble>(bslchydro));
    4879 
    4880         /*Take care of partition vectors:*/
    4881         if(bslcice_partition){
    4882                 bslcice_partition->Assemble();
    4883                 bslcice_partition_serial=bslcice_partition->ToMPISerial();
    4884         }
    4885         if(bslchydro_partition){
    4886                 bslchydro_partition->Assemble();
    4887                 bslchydro_partition_serial=bslchydro_partition->ToMPISerial();
    4888         }
    4889 
    4890 
    4891         /*Free ressources:*/
    4892         xDelete<int>(indices);
    4893         xDelete<IssmDouble>(RSLgi);
    4894         if(bslchydro_partition)delete bslchydro_partition;
    4895         if(bslcice_partition)delete bslcice_partition;
    4896         if(partitionhydro)xDelete<IssmDouble>(partitionhydro);
    4897         if(partitionice)xDelete<IssmDouble>(partitionice);
    4898 
    4899         /*Assign output pointers:*/
    4900         *poceanarea = oceanarea;
    4901         *pbslcice  = bslcice;
    4902         *pbslchydro  = bslchydro;
    4903         *pbslc=bslchydro+bslcice;
    4904         *pbslcice_partition=bslcice_partition_serial;
    4905         *pbslchydro_partition=bslchydro_partition_serial;
    4906 
    4907 }
    4908 /*}}}*/
    49094751void FemModel::SealevelchangeSal(Vector<IssmDouble>* pRSLgo, Vector<IssmDouble>* pRSLg_old,  SealevelMasks* masks, bool verboseconvolution){/*{{{*/
    49104752

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

    r25379 r26110  
    1919#include "./Misfit.h"
    2020#include "./SealevelMasks.h"
     21#include "./BarystaticContributions.h"
    2122#include "./Nodalvalue.h"
    2223#include "./Numberedcostfunction.h"

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

    r26104 r26110  
    8181void ResetBoundaryConditions(FemModel* femmodel, int analysis_type);
    8282void PrintBanner(void);
    83 void TransferForcing(FemModel* femmodel,int forcingenum);
    84 void TransferSealevel(FemModel* femmodel,int forcingenum);
    8583void EarthMassTransport(FemModel* femmodel);
    86 void slcconvergence(bool* pconverged, Vector<IssmDouble>* RSLg,Vector<IssmDouble>* RSLg_old,IssmDouble eps_rel,IssmDouble eps_abs);
    8784
    8885//solution configuration

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

    r26109 r26110  
    1212#include "../modules/modules.h"
    1313#include "../solutionsequences/solutionsequences.h"
     14/*support routines local definitions:{{{*/
     15void TransferForcing(FemModel* femmodel,int forcingenum);
     16void TransferSealevel(FemModel* femmodel,int forcingenum);
     17void slcconvergence(bool* pconverged, Vector<IssmDouble>* RSLg,Vector<IssmDouble>* RSLg_old,IssmDouble eps_rel,IssmDouble eps_abs);
     18IssmDouble SealevelloadsOceanAverage(Vector<IssmDouble>* sealevelloads,Vector<IssmDouble>* oceanareas, IssmDouble oceanarea);
     19/*}}}*/
    1420
    1521/*main cores:*/
    … …  
    383389        }
    384390}; /*}}}*/
     391
     392void grd_core_optim(FemModel* femmodel,SealevelMasks* masks) { /*{{{*/
     393
     394        /*variables:{{{*/
     395        int nel;
     396        BarystaticContributions* barycontrib=NULL;
     397        GenericParam<BarystaticContributions*>* barycontribparam=NULL;
     398        IssmDouble               barystatic;
     399       
     400        Vector<IssmDouble>*    loads=NULL;
     401        IssmDouble*            allloads=NULL;
     402        Vector<IssmDouble>*    sealevelloads=NULL;
     403        Vector<IssmDouble>*    oldsealevelloads=NULL;
     404        IssmDouble             sealevelloadsaverage;
     405        IssmDouble*            allsealevelloads=NULL;
     406        Vector<IssmDouble>*    oceanareas=NULL;
     407        IssmDouble             oceanarea;
     408        bool                   scaleoceanarea=false;
     409        IssmDouble             rho_water;
     410
     411        IssmDouble           eps_rel;
     412        IssmDouble           eps_abs;
     413        int                  step;
     414        IssmDouble           time;
     415       
     416        IssmDouble cumbslc;
     417        IssmDouble cumbslcice;
     418        IssmDouble cumbslchydro;
     419        /*}}}*/
     420
     421        /*retrieve parameters: */
     422        femmodel->parameters->FindParam(&scaleoceanarea,SolidearthSettingsOceanAreaScalingEnum);
     423        barycontribparam = xDynamicCast<GenericParam<BarystaticContributions*>*>(femmodel->parameters->FindParamObject(BarystaticContributionsEnum));
     424        barycontrib=barycontribparam->GetParameterValue();
     425        femmodel->parameters->FindParam(&rho_water,MaterialsRhoSeawaterEnum);
     426        femmodel->parameters->FindParam(&eps_rel,SolidearthSettingsReltolEnum);
     427        femmodel->parameters->FindParam(&eps_abs,SolidearthSettingsAbstolEnum);
     428       
     429        /*initialize matrices and vectors:*/
     430        femmodel->parameters->FindParam(&nel,MeshNumberofelementsEnum);
     431        loads=new Vector<IssmDouble>(nel);
     432        sealevelloads=new Vector<IssmDouble>(nel);
     433        oceanareas=new Vector<IssmDouble>(nel);
     434       
     435        /*buildup loads: */
     436        for(Object* & object : femmodel->elements->objects){
     437                Element* element = xDynamicCast<Element*>(object);
     438                element->SealevelchangeBarystaticLoads(loads, barycontrib,masks);
     439        }
     440
     441        //Communicate loads from local to global:
     442        loads->Assemble(); allloads=loads->ToMPISerial();
     443
     444        /*convolve loads:*/
     445        for(Object* & object : femmodel->elements->objects){
     446                Element* element = xDynamicCast<Element*>(object);
     447                element->SealevelchangeInitialConvolution(sealevelloads,oceanareas,allloads,masks);
     448        }
     449
     450        //Get ocean area:
     451        oceanareas->Assemble(); oceanareas->Sum(&oceanarea); _assert_(oceanarea>0.);
     452        if(scaleoceanarea) oceanarea=3.619e+14; // use true ocean area, m^2
     453
     454        //Get sea level loads ocean average:
     455        sealevelloadsaverage=SealevelloadsOceanAverage(sealevelloads,oceanareas,oceanarea);
     456       
     457        //substract ocean average and barystatic contributionfrom sea level loads:
     458        barystatic=barycontrib->Total()/oceanarea/rho_water;
     459        sealevelloads->Shift(-sealevelloadsaverage+barystatic);
     460        allsealevelloads=sealevelloads->ToMPISerial();
     461
     462        bool converged=false;
     463        for(;;){
     464                       
     465                oldsealevelloads=sealevelloads->Duplicate();
     466
     467                /*convolve load and sealevel loads on oceans:*/
     468                for(Object* & object : femmodel->elements->objects){
     469                        Element* element = xDynamicCast<Element*>(object);
     470                        element->SealevelchangeOceanConvolution(sealevelloads, allsealevelloads, allloads,masks);
     471                }
     472                sealevelloads->Assemble();
     473               
     474                //substract ocean average and barystatic contribution
     475                sealevelloadsaverage=SealevelloadsOceanAverage(sealevelloads,oceanareas,oceanarea);
     476                sealevelloads->Shift(-sealevelloadsaverage+barystatic);
     477
     478                //broadcast loads
     479                allsealevelloads=sealevelloads->ToMPISerial();
     480
     481                //convergence?
     482                slcconvergence(&converged,sealevelloads,oldsealevelloads,eps_rel,eps_abs);
     483                if (converged)break;
     484        }
     485
     486        /*cumulate barystatic contributions and save to results: */
     487        barycontrib->Cumulate(femmodel->parameters);
     488        barycontrib->Save(femmodel->results,femmodel->parameters,oceanarea);
     489
     490}
     491/*}}}*/
    385492
    386493//Geometry:
    … …  
    10621169
    10631170} /*}}}*/
     1171IssmDouble SealevelloadsOceanAverage(Vector<IssmDouble>* sealevelloads,Vector<IssmDouble>* oceanareas, IssmDouble oceanarea){ /*{{{*/
     1172        IssmDouble sealevelloadsaverage;       
     1173        Vector<IssmDouble>* sealevelloadsvolume=sealevelloads->Duplicate();
     1174        sealevelloadsvolume->PointwiseMult(sealevelloads,oceanareas);
     1175        sealevelloadsvolume->Sum(&sealevelloadsaverage);
     1176        delete sealevelloadsvolume;
     1177        return sealevelloadsaverage/oceanarea;
     1178} /*}}}*/

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

    r26099 r26110  
    747747syn keyword cConstant BslcIceEnum
    748748syn keyword cConstant BslcHydroEnum
     749syn keyword cConstant BslcOceanEnum
    749750syn keyword cConstant BslcRateEnum
    750751syn keyword cConstant GmtslcEnum
    … …  
    10711072syn keyword cConstant BasalforcingsIsmip6Enum
    10721073syn keyword cConstant BasalforcingsPicoEnum
     1074syn keyword cConstant BarystaticContributionsEnum
    10731075syn keyword cConstant BeckmannGoosseFloatingMeltRateEnum
    10741076syn keyword cConstant BedSlopeSolutionEnum
    … …  
    14331435syn keyword cType AmrNeopz
    14341436syn keyword cType ArrayInput
     1437syn keyword cType BarystaticContributions
    14351438syn keyword cType BoolInput
    14361439syn keyword cType BoolParam

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

    r26099 r26110  
    743743        BslcIceEnum,
    744744        BslcHydroEnum,
     745        BslcOceanEnum,
    745746        BslcRateEnum,
    746747        GmtslcEnum,
    … …  
    10701071        BasalforcingsIsmip6Enum,
    10711072        BasalforcingsPicoEnum,
     1073        BarystaticContributionsEnum,
    10721074        BeckmannGoosseFloatingMeltRateEnum,
    10731075        BedSlopeSolutionEnum,

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

    r26099 r26110  
    749749                case BslcIceEnum : return "BslcIce";
    750750                case BslcHydroEnum : return "BslcHydro";
     751                case BslcOceanEnum : return "BslcOcean";
    751752                case BslcRateEnum : return "BslcRate";
    752753                case GmtslcEnum : return "Gmtslc";
    … …  
    10731074                case BasalforcingsIsmip6Enum : return "BasalforcingsIsmip6";
    10741075                case BasalforcingsPicoEnum : return "BasalforcingsPico";
     1076                case BarystaticContributionsEnum : return "BarystaticContributions";
    10751077                case BeckmannGoosseFloatingMeltRateEnum : return "BeckmannGoosseFloatingMeltRate";
    10761078                case BedSlopeSolutionEnum : return "BedSlopeSolution";

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

    r26099 r26110  
    767767              else if (strcmp(name,"BslcIce")==0) return BslcIceEnum;
    768768              else if (strcmp(name,"BslcHydro")==0) return BslcHydroEnum;
     769              else if (strcmp(name,"BslcOcean")==0) return BslcOceanEnum;
    769770              else if (strcmp(name,"BslcRate")==0) return BslcRateEnum;
    770771              else if (strcmp(name,"Gmtslc")==0) return GmtslcEnum;
    … …  
    874875              else if (strcmp(name,"SmbSmbref")==0) return SmbSmbrefEnum;
    875876              else if (strcmp(name,"SmbSzaValue")==0) return SmbSzaValueEnum;
    876               else if (strcmp(name,"SmbT")==0) return SmbTEnum;
    877877         else stage=8;
    878878   }
    879879   if(stage==8){
    880               if (strcmp(name,"SmbTa")==0) return SmbTaEnum;
     880              if (strcmp(name,"SmbT")==0) return SmbTEnum;
     881              else if (strcmp(name,"SmbTa")==0) return SmbTaEnum;
    881882              else if (strcmp(name,"SmbTeValue")==0) return SmbTeValueEnum;
    882883              else if (strcmp(name,"SmbTemperaturesAnomaly")==0) return SmbTemperaturesAnomalyEnum;
    … …  
    997998              else if (strcmp(name,"Outputdefinition32")==0) return Outputdefinition32Enum;
    998999              else if (strcmp(name,"Outputdefinition33")==0) return Outputdefinition33Enum;
    999               else if (strcmp(name,"Outputdefinition34")==0) return Outputdefinition34Enum;
    10001000         else stage=9;
    10011001   }
    10021002   if(stage==9){
    1003               if (strcmp(name,"Outputdefinition35")==0) return Outputdefinition35Enum;
     1003              if (strcmp(name,"Outputdefinition34")==0) return Outputdefinition34Enum;
     1004              else if (strcmp(name,"Outputdefinition35")==0) return Outputdefinition35Enum;
    10041005              else if (strcmp(name,"Outputdefinition36")==0) return Outputdefinition36Enum;
    10051006              else if (strcmp(name,"Outputdefinition37")==0) return Outputdefinition37Enum;
    … …  
    10971098              else if (strcmp(name,"BasalforcingsIsmip6")==0) return BasalforcingsIsmip6Enum;
    10981099              else if (strcmp(name,"BasalforcingsPico")==0) return BasalforcingsPicoEnum;
     1100              else if (strcmp(name,"BarystaticContributions")==0) return BarystaticContributionsEnum;
    10991101              else if (strcmp(name,"BeckmannGoosseFloatingMeltRate")==0) return BeckmannGoosseFloatingMeltRateEnum;
    11001102              else if (strcmp(name,"BedSlopeSolution")==0) return BedSlopeSolutionEnum;
    … …  
    11191121              else if (strcmp(name,"Closed")==0) return ClosedEnum;
    11201122              else if (strcmp(name,"Colinear")==0) return ColinearEnum;
    1121               else if (strcmp(name,"Constraints")==0) return ConstraintsEnum;
    1122               else if (strcmp(name,"Contact")==0) return ContactEnum;
    11231123         else stage=10;
    11241124   }
    11251125   if(stage==10){
    1126               if (strcmp(name,"Contour")==0) return ContourEnum;
     1126              if (strcmp(name,"Constraints")==0) return ConstraintsEnum;
     1127              else if (strcmp(name,"Contact")==0) return ContactEnum;
     1128              else if (strcmp(name,"Contour")==0) return ContourEnum;
    11271129              else if (strcmp(name,"Contours")==0) return ContoursEnum;
    11281130              else if (strcmp(name,"ControlInput")==0) return ControlInputEnum;
    … …  
    12421244              else if (strcmp(name,"L2ProjectionEPLAnalysis")==0) return L2ProjectionEPLAnalysisEnum;
    12431245              else if (strcmp(name,"LACrouzeixRaviart")==0) return LACrouzeixRaviartEnum;
    1244               else if (strcmp(name,"LATaylorHood")==0) return LATaylorHoodEnum;
    1245               else if (strcmp(name,"LambdaS")==0) return LambdaSEnum;
    12461246         else stage=11;
    12471247   }
    12481248   if(stage==11){
    1249               if (strcmp(name,"LevelsetAnalysis")==0) return LevelsetAnalysisEnum;
     1249              if (strcmp(name,"LATaylorHood")==0) return LATaylorHoodEnum;
     1250              else if (strcmp(name,"LambdaS")==0) return LambdaSEnum;
     1251              else if (strcmp(name,"LevelsetAnalysis")==0) return LevelsetAnalysisEnum;
    12501252              else if (strcmp(name,"LevelsetfunctionPicard")==0) return LevelsetfunctionPicardEnum;
    12511253              else if (strcmp(name,"LinearFloatingMeltRate")==0) return LinearFloatingMeltRateEnum;
    … …  
    13651367              else if (strcmp(name,"SMBpddSicopolis")==0) return SMBpddSicopolisEnum;
    13661368              else if (strcmp(name,"SMBsemic")==0) return SMBsemicEnum;
    1367               else if (strcmp(name,"SSAApproximation")==0) return SSAApproximationEnum;
    1368               else if (strcmp(name,"SSAFSApproximation")==0) return SSAFSApproximationEnum;
    13691369         else stage=12;
    13701370   }
    13711371   if(stage==12){
    1372               if (strcmp(name,"SSAHOApproximation")==0) return SSAHOApproximationEnum;
     1372              if (strcmp(name,"SSAApproximation")==0) return SSAApproximationEnum;
     1373              else if (strcmp(name,"SSAFSApproximation")==0) return SSAFSApproximationEnum;
     1374              else if (strcmp(name,"SSAHOApproximation")==0) return SSAHOApproximationEnum;
    13731375              else if (strcmp(name,"Scaled")==0) return ScaledEnum;
    13741376              else if (strcmp(name,"SealevelAbsolute")==0) return SealevelAbsoluteEnum;

  • issm/trunk-jpl/src/c/toolkits/issm/IssmAbsVec.h

    r25995 r26110  
    5252                virtual void PointwiseMult(IssmAbsVec* x,IssmAbsVec* y)=0;
    5353                virtual void Pow(doubletype scale_factor)=0;
     54                virtual void Sum(doubletype* pvalue)=0;
    5455};
    5556

  • issm/trunk-jpl/src/c/toolkits/issm/IssmMpiVec.h

    r25995 r26110  
    594594                }
    595595                /*}}}*/
     596                void Sum(doubletype* pvalue){/*{{{*/
     597                        _error_("not support yet!");
     598                }
     599                /*}}}*/
    596600                void BucketsBuildScatterBuffers(int** pnumvalues_forcpu,int** prow_indices_forcpu,doubletype** pvalues_forcpu,int** pmodes_forcpu,DataSet** bucketsforcpu,int num_procs){/*{{{*/
    597601

  • issm/trunk-jpl/src/c/toolkits/issm/IssmSeqVec.h

    r25995 r26110  
    324324                }
    325325                /*}}}*/
     326                void Sum(doubletype* pvalue){/*{{{*/
     327
     328                        doubletype value=0;
     329                        int i;
     330                        for(i=0;i<this->M;i++)value+=this->vector[i];
     331                        *pvalue=value;
     332                }
     333                /*}}}*/
    326334               
    327335};

  • issm/trunk-jpl/src/c/toolkits/issm/IssmVec.h

    r25995 r26110  
    220220                }
    221221                /*}}}*/
     222                void Sum(doubletype*  pvalue){/*{{{*/
     223                        vector->Sum(pvalue);
     224                }
     225                /*}}}*/
    222226};
    223227

  • issm/trunk-jpl/src/c/toolkits/objects/Matrix.h

    r24679 r26110  
    282282                }
    283283                /*}}}*/
     284                doubletype* ToMPISerial(void){/*{{{*/
     285
     286                        doubletype* output=NULL;
     287
     288                        if(type==PetscMatType){
     289                                #ifdef _HAVE_PETSC_
     290                                output=this->pmatrix->ToMPISerial();
     291                                #endif
     292                        }
     293                        else{
     294                                _error_("not implemented yet!");
     295                        }
     296
     297                        return output;
     298                }
     299                /*}}}*/
    284300                void SetValues(int m,int* idxm,int n,int* idxn,IssmDouble* values,InsMode mode){/*{{{*/
    285301
    … …  
    307323                }
    308324                /*}}}*/
    309                 /*
    310                 * sets all values to 0 but keeps the structure of a sparse matrix
    311                 */
    312325                void SetZero(void) {/*{{{*/
     326                        // sets all values to 0 but keeps the structure of a sparse matrix
    313327                        if(type==PetscMatType){
    314328                                #ifdef _HAVE_PETSC_

  • issm/trunk-jpl/src/c/toolkits/objects/Vector.h

    r25995 r26110  
    406406                }
    407407                /*}}}*/
    408 };
     408void Sum(doubletype* pvalue){ /*{{{*/
     409        _assert_(this);/*{{{*/
     410
     411        if(type==PetscVecType){
     412                #ifdef _HAVE_PETSC_
     413                this->pvector->Sum(pvalue);
     414                #endif
     415        }
     416        else this->ivector->Sum(pvalue);
     417}
     418/*}}}*/
     419}; /*}}}*/
    409420#endif //#ifndef _VECTOR_H_

  • issm/trunk-jpl/src/c/toolkits/petsc/objects/PetscVec.cpp

    r24240 r26110  
    249249}
    250250/*}}}*/
     251void PetscVec::Sum(IssmDouble* pvalue){/*{{{*/
     252
     253        _assert_(this->vector);
     254        VecSum(this->vector,pvalue);
     255
     256}
     257/*}}}*/
    251258IssmDouble PetscVec::Dot(PetscVec* input){/*{{{*/
    252259

  • issm/trunk-jpl/src/c/toolkits/petsc/objects/PetscVec.h

    r24240 r26110  
    5858                void        Scale(IssmDouble scale_factor);
    5959                void        Pow(IssmDouble scale_factor);
     60                void        Sum(IssmDouble* pvalue);
    6061                void        PointwiseDivide(PetscVec* x,PetscVec* y);
    6162                void        PointwiseMult(PetscVec* x,PetscVec* y);

  • issm/trunk-jpl/src/c/toolkits/petsc/patches/petscpatches.h

    r23643 r26110  
    3232void MatMultPatch(Mat A,Vec X, Vec AX,ISSM_MPI_Comm comm);
    3333void MatToSerial(double** poutmatrix,Mat matrix,ISSM_MPI_Comm comm);
     34void MatToMPISerial(double** poutmatrix,Mat matrix,ISSM_MPI_Comm comm);
    3435Vec  SerialToVec(double* vector,int vector_size);
    3536InsertMode ISSMToPetscInsertMode(InsMode mode);

  • issm/trunk-jpl/test/NightlyRun/test2004.m

    r26077 r26110  
    134134                disp('      reading bedrock');
    135135                md.geometry.bed=-ones(md.mesh.numberofvertices,1);
     136                md.geometry.base=md.geometry.bed;
     137                md.geometry.thickness=1000*ones(md.mesh.numberofvertices,1);
     138                md.geometry.surface=md.geometry.bed+md.geometry.thickness;
     139
    136140        end % }}}
    137141        %Slc: {{{
    … …  
    167171                end
    168172
    169                 md.solidearth.initialsealevel=zeros(md.mesh.numberofvertices,1);
    170 
    171                 md.dsl.global_average_thermosteric_sea_level_change=[0;0];
    172                 md.dsl.sea_surface_height_change_above_geoid=zeros(md.mesh.numberofvertices+1,1);
    173                 md.dsl.sea_water_pressure_change_at_sea_floor=zeros(md.mesh.numberofvertices+1,1);
     173                md.initialization.sealevel=zeros(md.mesh.numberofvertices,1);
     174
     175                md.dsl.global_average_thermosteric_sea_level=[0;0];
     176                md.dsl.sea_surface_height_above_geoid=zeros(md.mesh.numberofvertices+1,1);
     177                md.dsl.sea_water_pressure_at_sea_floor=zeros(md.mesh.numberofvertices+1,1);
    174178
    175179        end %}}}
    … …  
    313317        di=md.materials.rho_ice/md.materials.rho_water;
    314318        md.geometry.bed=-ones(md.mesh.numberofvertices,1);
     319        md.geometry.base=md.geometry.bed;
     320        md.geometry.thickness=1000*ones(md.mesh.numberofvertices,1);
     321        md.geometry.surface=md.geometry.bed+md.geometry.thickness;
    315322        % }}}
    316323        %materials:  {{{
    … …  
    346353sl.transfer('mask.ocean_levelset');
    347354sl.transfer('geometry.bed');
     355sl.transfer('geometry.surface');
     356sl.transfer('geometry.thickness');
     357sl.transfer('geometry.base');
    348358sl.transfer('mesh.lat');
    349359sl.transfer('mesh.long');
    … …  
    400410md.transient.isslc=1;
    401411
     412%Initializations:
     413md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1);
     414md.basalforcings.floatingice_melting_rate=zeros(md.mesh.numberofvertices,1);
     415md.initialization.vx=zeros(md.mesh.numberofvertices,1);
     416md.initialization.vy=zeros(md.mesh.numberofvertices,1);
     417md.initialization.sealevel=zeros(md.mesh.numberofvertices,1);
     418md.initialization.bottompressure=zeros(md.mesh.numberofvertices,1);
     419md.initialization.dsl=zeros(md.mesh.numberofvertices,1);
     420md.initialization.str=0;
     421md.smb.mass_balance=zeros(md.mesh.numberofvertices,1);
    402422
    403423%max number of iterations reverted back to 10 (i.e. the original default value)
Note: See TracChangeset for help on using the changeset viewer.