Changeset 26165


Ignore:
Timestamp:
03/30/21 19:32:13 (4 years ago)
Author:
Eric.Larour
Message:

CHG: hooking the new optimized code and removing the old one. Moving everything
away from FemModel into the sea level core (including Ivins deformation model).

Location:
issm/trunk-jpl
Files:
18 edited

Legend:

Unmodified
Added
Removed

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

    r26126 r26165  
    145145        parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.glfraction",SolidearthSettingsGlfractionEnum));
    146146        parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.cross_section_shape",SolidearthSettingsCrossSectionShapeEnum));
    147         parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.optim",SolidearthSettingsOptimEnum));
    148147        parameters->AddObject(new DoubleParam(CumBslcEnum,0.0));
    149148        parameters->AddObject(new DoubleParam(CumBslcIceEnum,0.0));

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

    r26164 r26165  
    383383                virtual IssmDouble    GetArea3D(void)=0;
    384384                virtual IssmDouble    GetAreaSpherical(void)=0;
    385                 virtual IssmDouble    OceanAverage(IssmDouble* Sg, SealevelMasks* masks)=0;
    386                 virtual void          SealevelchangeMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks)=0;
    387                 virtual void          SealevelchangeMomentOfInertiaElement(IssmDouble* dI_list, IssmDouble* loads, IssmDouble* sealevelloads)=0;
    388                 virtual IssmDouble    SealevelchangeBarystaticIce(IssmDouble* Sgi, SealevelMasks* masks,Vector<IssmDouble>* barystatic_contribution,IssmDouble* partition,IssmDouble oceanarea)=0;
    389                 virtual IssmDouble    SealevelchangeBarystaticHydro(IssmDouble* Sgi, SealevelMasks* masks,Vector<IssmDouble>* barystatic_contribution, IssmDouble* partition, IssmDouble oceanarea)=0;
    390                 virtual void          SealevelchangeBarystaticBottomPressure(IssmDouble* Sgi, SealevelMasks* masks)=0;
    391                 virtual void          SealevelchangeGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx, IssmDouble* yy, IssmDouble* zz,
    392                                 IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze)=0;
    393                 virtual void          SealevelchangeSal(IssmDouble* Sgo, IssmDouble* Sg_old,SealevelMasks* mask)=0;
    394                 virtual void          DeformationFromSurfaceLoads(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg,SealevelMasks* masks)=0;
    395                 virtual void       GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt,Matlitho* litho, IssmDouble* x,IssmDouble* y)=0;
    396 
    397                 virtual void       SealevelchangeGeometryOptim(IssmDouble* lat,IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae)=0;
    398                 virtual void       SealevelchangeBarystaticLoads(Vector<IssmDouble>* loads, BarystaticContributions* barycontrib, SealevelMasks* masks)=0;
    399                 virtual void       SealevelchangeConvolution(Vector<IssmDouble>* sealevelloads, Vector<IssmDouble>* oceanareas, IssmDouble* allsealevelloads, IssmDouble* allloads,IssmDouble* rotationaxismotionvector)=0;
    400                 virtual void       SealevelchangeDeformationConvolution(IssmDouble* sealevelloads, IssmDouble* loads, IssmDouble* rotationvector)=0;
    401                 virtual void       SealevelchangeShift(Vector<IssmDouble>* loads, IssmDouble offset, SealevelMasks* masks)=0;
     385                virtual void          SealevelchangeMomentOfInertia(IssmDouble* dI_list, IssmDouble* loads, IssmDouble* sealevelloads)=0;
     386                virtual void          GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt,Matlitho* litho, IssmDouble* x,IssmDouble* y)=0;
     387
     388                virtual void          SealevelchangeGeometry(IssmDouble* lat,IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae)=0;
     389                virtual void          SealevelchangeBarystaticLoads(Vector<IssmDouble>* loads, BarystaticContributions* barycontrib, SealevelMasks* masks)=0;
     390                virtual void          SealevelchangeConvolution(Vector<IssmDouble>* sealevelloads, Vector<IssmDouble>* oceanareas, IssmDouble* allsealevelloads, IssmDouble* allloads,IssmDouble* rotationaxismotionvector)=0;
     391                virtual void          SealevelchangeDeformationConvolution(IssmDouble* sealevelloads, IssmDouble* loads, IssmDouble* rotationvector)=0;
     392                virtual void          SealevelchangeShift(Vector<IssmDouble>* loads, IssmDouble offset, SealevelMasks* masks)=0;
    402393                #endif
    403394

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

    r26164 r26165  
    219219                #endif
    220220                #ifdef _HAVE_SEALEVELCHANGE_
    221                 IssmDouble    OceanAverage(IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet!");};
    222                 void    SetSealevelMasks(SealevelMasks* masks){_error_("not implemented yet!");};
    223                 void    SealevelchangeMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks){_error_("not implemented yet!");};
    224                 void    SealevelchangeMomentOfInertiaElement(IssmDouble* dI_list, IssmDouble* loads, IssmDouble* sealevelloads){_error_("not implemented yet!");};
    225                 void    SealevelchangeGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx, IssmDouble* yy, IssmDouble* zz,
    226                                         IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze){_error_("not implemented yet!");};
    227                 IssmDouble    SealevelchangeBarystaticIce(IssmDouble* Sgi, SealevelMasks* masks, Vector<IssmDouble>* barystatic_contribution, IssmDouble* partition, IssmDouble oceanarea){_error_("not implemented yet!");};
    228                 IssmDouble    SealevelchangeBarystaticHydro(IssmDouble* Sgi, SealevelMasks* masks, Vector<IssmDouble>* barystatic_contribution, IssmDouble* partition,IssmDouble oceanarea){_error_("not implemented yet!");};
    229                 void    SealevelchangeBarystaticBottomPressure(IssmDouble* Sgi,SealevelMasks* masks){_error_("not implemented yet!");};
    230                 void    SealevelchangeSal(IssmDouble* Sgo, IssmDouble* Sg_old, SealevelMasks* masks){_error_("not implemented yet!");};
    231                 void    DeformationFromSurfaceLoads(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet!");};
    232                 void           GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt,Matlitho* litho, IssmDouble* x,IssmDouble* y){_error_("not implemented yet");};
    233 
    234                 void       SealevelchangeGeometryOptim(IssmDouble* lat,IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae){_error_("not implemented yet");};
     221                void       SetSealevelMasks(SealevelMasks* masks){_error_("not implemented yet!");};
     222                void       SealevelchangeMomentOfInertia(IssmDouble* dI_list, IssmDouble* loads, IssmDouble* sealevelloads){_error_("not implemented yet!");};
     223                void       GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt,Matlitho* litho, IssmDouble* x,IssmDouble* y){_error_("not implemented yet");};
     224
     225                void       SealevelchangeGeometry(IssmDouble* lat,IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae){_error_("not implemented yet");};
    235226                void       SealevelchangeBarystaticLoads(Vector<IssmDouble>* loads, BarystaticContributions* barycontrib, SealevelMasks* masks){_error_("not implemented yet");};
    236227                void       SealevelchangeConvolution(Vector<IssmDouble>* sealevelloads, Vector<IssmDouble>* oceanareas, IssmDouble* allsealevelloads, IssmDouble* allloads,IssmDouble* rotationaxismotionvector){_error_("not implemented yet");};

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

    r26164 r26165  
    174174#endif
    175175#ifdef _HAVE_SEALEVELCHANGE_
    176                 void    SealevelchangeMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks){_error_("not implemented yet!");};
    177                 void    SealevelchangeMomentOfInertiaElement(IssmDouble* dI_list, IssmDouble* loads, IssmDouble* sealevelloads){_error_("not implemented yet!");};
    178                 void    SetSealevelMasks(SealevelMasks* masks){_error_("not implemented yet!");};
    179                 void    SealevelchangeGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx, IssmDouble* yy, IssmDouble* zz,
    180                                         IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze){_error_("not implemented yet!");};
    181                 IssmDouble    SealevelchangeBarystaticIce(IssmDouble* Sgi, SealevelMasks* masks, Vector<IssmDouble>* barystatic_contribution, IssmDouble* partition, IssmDouble oceanarea){_error_("not implemented yet!");};
    182                 IssmDouble    SealevelchangeBarystaticHydro(IssmDouble* Sgi, SealevelMasks* masks, Vector<IssmDouble>* barystatic_contribution, IssmDouble* partition,IssmDouble oceanarea){_error_("not implemented yet!");};
    183                 void    SealevelchangeBarystaticBottomPressure(IssmDouble* Sgi,SealevelMasks* masks){_error_("not implemented yet!");};
    184                 void    SealevelchangeSal(IssmDouble* Sgo,IssmDouble* Sg_old, SealevelMasks* masks){_error_("not implemented yet!");};
    185                 void    DeformationFromSurfaceLoads(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet!");};
    186                 IssmDouble    OceanAverage(IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet!");};
    187                 void        GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt,Matlitho* litho, IssmDouble* x,IssmDouble* y){_error_("not implemented yet");};
    188 
    189                 void       SealevelchangeGeometryOptim(IssmDouble* lat,IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae){_error_("not implemented yet");};
     176                void       SealevelchangeMomentOfInertia(IssmDouble* dI_list, IssmDouble* loads, IssmDouble* sealevelloads){_error_("not implemented yet!");};
     177                void       SetSealevelMasks(SealevelMasks* masks){_error_("not implemented yet!");};
     178                void       GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt,Matlitho* litho, IssmDouble* x,IssmDouble* y){_error_("not implemented yet");};
     179                void       SealevelchangeGeometry(IssmDouble* lat,IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae){_error_("not implemented yet");};
    190180                void       SealevelchangeBarystaticLoads(Vector<IssmDouble>* loads, BarystaticContributions* barycontrib, SealevelMasks* masks){_error_("not implemented yet");};
    191181                void       SealevelchangeConvolution(Vector<IssmDouble>* sealevelloads, Vector<IssmDouble>* oceanareas, IssmDouble* allsealevelloads, IssmDouble* allloads,IssmDouble* rotationaxismotionvector){_error_("not implemented yet");};

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

    r26164 r26165  
    181181#endif
    182182#ifdef _HAVE_SEALEVELCHANGE_
    183                 void    SetSealevelMasks(SealevelMasks* masks){_error_("not implemented yet!");};
    184                 void    SealevelchangeMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks){_error_("not implemented yet!");};
    185                 void    SealevelchangeMomentOfInertiaElement(IssmDouble* dI_list, IssmDouble* loads, IssmDouble* sealevelloads){_error_("not implemented yet!");};
    186                 void    SealevelchangeGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx, IssmDouble* yy, IssmDouble* zz,
    187                                         IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze){_error_("not implemented yet!");};
    188                 IssmDouble    SealevelchangeBarystaticIce(IssmDouble* Sgi, SealevelMasks* masks, Vector<IssmDouble>* barystatic_contribution, IssmDouble* partition, IssmDouble oceanarea){_error_("not implemented yet!");};
    189                 IssmDouble    SealevelchangeBarystaticHydro(IssmDouble* Sgi, SealevelMasks* masks, Vector<IssmDouble>* barystatic_contribution, IssmDouble* partition,IssmDouble oceanarea){_error_("not implemented yet!");};
    190                 void    SealevelchangeBarystaticBottomPressure(IssmDouble* Sgi,SealevelMasks* masks){_error_("not implemented yet!");};
    191                 void    SealevelchangeSal(IssmDouble* Sgo, IssmDouble* Sg_old, SealevelMasks* masks){_error_("not implemented yet!");};
    192                 void    DeformationFromSurfaceLoads(IssmDouble* Up ,IssmDouble* North, IssmDouble* East, IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet!");};
    193                 IssmDouble    OceanAverage(IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet!");};
    194                 void        GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt, Matlitho* litho, IssmDouble* x,IssmDouble* y){_error_("not implemented yet");};
    195                 void       SealevelchangeGeometryOptim(IssmDouble* lat,IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae){_error_("not implemented yet");};
     183                void       SetSealevelMasks(SealevelMasks* masks){_error_("not implemented yet!");};
     184                void       SealevelchangeMomentOfInertia(IssmDouble* dI_list, IssmDouble* loads, IssmDouble* sealevelloads){_error_("not implemented yet!");};
     185                void       GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt, Matlitho* litho, IssmDouble* x,IssmDouble* y){_error_("not implemented yet");};
     186                void       SealevelchangeGeometry(IssmDouble* lat,IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae){_error_("not implemented yet");};
    196187                void       SealevelchangeBarystaticLoads(Vector<IssmDouble>* loads, BarystaticContributions* barycontrib, SealevelMasks* masks){_error_("not implemented yet");};
    197188                void       SealevelchangeConvolution(Vector<IssmDouble>* sealevelloads, Vector<IssmDouble>* oceanareas, IssmDouble* allsealevelloads, IssmDouble* allloads,IssmDouble* rotationaxismotionvector){_error_("not implemented yet");};

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

    r26164 r26165  
    55565556#ifdef _HAVE_SEALEVELCHANGE_
    55575557//old code
    5558 IssmDouble Tria::OceanAverage(IssmDouble* Sg, SealevelMasks* masks){ /*{{{*/
    5559 
    5560         if(masks->isoceanin[this->lid]){
    5561 
    5562                 IssmDouble area;
    5563 
    5564                 /*Get area of element:*/
    5565                 this->Element::GetInputValue(&area,AreaEnum);
    5566 
    5567                 /*Average Sg over vertices:*/
    5568                 IssmDouble Sg_avg=0; for(int i=0;i<NUMVERTICES;i++) Sg_avg+=Sg[this->vertices[i]->Sid()]/NUMVERTICES;
    5569 
    5570                 /*return: */
    5571                 return area*Sg_avg;
    5572         }
    5573         else return 0;
    5574 
    5575 }
    5576 /*}}}*/
    5577 void          Tria::SealevelchangeMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks){/*{{{*/
    5578         /*early return if we are not on an ice cap OR ocean:*/
    5579         if(!masks->isiceonly[this->lid] && !masks->isoceanin[this->lid]){
    5580                 dI_list[0] = 0.0; // this is important!!!
    5581                 dI_list[1] = 0.0; // this is important!!!
    5582                 dI_list[2] = 0.0; // this is important!!!
    5583                 return;
    5584         }
    5585 
    5586         /*Compute area of element:*/
    5587         IssmDouble area,planetarea;
    5588         this->Element::GetInputValue(&area,AreaEnum);
    5589 
    5590         /*recover earth area: */
    5591         this->parameters->FindParam(&planetarea,SolidearthPlanetAreaEnum);
    5592 
    5593         /*Compute lat,long,radius of elemental centroid: */
    5594         bool spherical=true;
    5595         IssmDouble llr_list[NUMVERTICES][3];
    5596         IssmDouble late,longe,re;
    5597         /* Where is the centroid of this element?:{{{*/
    5598         ::GetVerticesCoordinates(&llr_list[0][0],this->vertices,NUMVERTICES,spherical);
    5599 
    5600         IssmDouble minlong=400;
    5601         IssmDouble maxlong=-20;
    5602         for (int i=0;i<NUMVERTICES;i++){
    5603                 llr_list[i][0]=(90-llr_list[i][0]);
    5604                 if(llr_list[i][1]<0)llr_list[i][1]=180+(180+llr_list[i][1]);
    5605                 if(llr_list[i][1]>maxlong)maxlong=llr_list[i][1];
    5606                 if(llr_list[i][1]<minlong)minlong=llr_list[i][1];
    5607         }
    5608         if(minlong==0 && maxlong>180){
    5609                 if (llr_list[0][1]==0)llr_list[0][1]=360;
    5610                 if (llr_list[1][1]==0)llr_list[1][1]=360;
    5611                 if (llr_list[2][1]==0)llr_list[2][1]=360;
    5612         }
    5613 
    5614         // correction at the north pole
    5615         if(llr_list[0][0]==0)llr_list[0][1]=(llr_list[1][1]+llr_list[2][1])/2.0;
    5616         if(llr_list[1][0]==0)llr_list[1][1]=(llr_list[0][1]+llr_list[2][1])/2.0;
    5617         if(llr_list[2][0]==0)llr_list[2][1]=(llr_list[0][1]+llr_list[1][1])/2.0;
    5618 
    5619         //correction at the south pole
    5620         if(llr_list[0][0]==180)llr_list[0][1]=(llr_list[1][1]+llr_list[2][1])/2.0;
    5621         if(llr_list[1][0]==180)llr_list[1][1]=(llr_list[0][1]+llr_list[2][1])/2.0;
    5622         if(llr_list[2][0]==180)llr_list[2][1]=(llr_list[0][1]+llr_list[1][1])/2.0;
    5623 
    5624         late=(llr_list[0][0]+llr_list[1][0]+llr_list[2][0])/3.0;
    5625         longe=(llr_list[0][1]+llr_list[1][1]+llr_list[2][1])/3.0;
    5626 
    5627         late=90.-late;
    5628         if(longe>180.)longe=(longe-180.)-180.;
    5629 
    5630         late=late/180.*M_PI;
    5631         longe=longe/180.*M_PI;
    5632         /*}}}*/
    5633         re=(llr_list[0][2]+llr_list[1][2]+llr_list[2][2])/3.0;
    5634 
    5635         if(masks->isoceanin[this->lid]){
    5636                 IssmDouble rho_water, S;
    5637 
    5638                 /*recover material parameters: */
    5639                 rho_water=FindParam(MaterialsRhoSeawaterEnum);
    5640 
    5641                 /*From Sg_old, recover water sea level change:*/
    5642                 S=0; for(int i=0;i<NUMVERTICES;i++) S+=Sg_old[this->vertices[i]->Sid()]/NUMVERTICES;
    5643 
    5644                 /* Perturbation terms for moment of inertia (moi_list):
    5645                  * computed analytically (see Wu & Peltier, eqs 10 & 32)
    5646                  * also consistent with my GMD formulation!
    5647                  * ALL in geographic coordinates
    5648                  * */
    5649                 dI_list[0] = -4*M_PI*(rho_water*S*area)*pow(re,4)*(sin(late)*cos(late)*cos(longe))/planetarea;
    5650                 dI_list[1] = -4*M_PI*(rho_water*S*area)*pow(re,4)*(sin(late)*cos(late)*sin(longe))/planetarea;
    5651                 dI_list[2] = +4*M_PI*(rho_water*S*area)*pow(re,4)*(1-pow(sin(late),2))/planetarea;
    5652         }
    5653         if(masks->isiceonly[this->lid]){
    5654                
    5655                 IssmDouble rho_ice,I;
    5656                 bool computerigid= false;
    5657        
    5658                 bool notfullygrounded=false;
    5659        
    5660                 bool scaleoceanarea= false;
    5661                 int  glfraction=1;
    5662                 IssmDouble phi=1.0;
    5663                 /*{{{*/
    5664                 if (masks->notfullygrounded[this->lid]) notfullygrounded=true; //used later on.
    5665 
    5666                 /*recover some parameters:*/
    5667                 rho_ice=FindParam(MaterialsRhoIceEnum);
    5668                 this->parameters->FindParam(&computerigid,SolidearthSettingsRigidEnum);
    5669                 this->parameters->FindParam(&scaleoceanarea,SolidearthSettingsOceanAreaScalingEnum);
    5670                 this->parameters->FindParam(&glfraction,SolidearthSettingsGlfractionEnum);
    5671 
    5672                 /*Get area of element: precomputed in the sealevelchange_geometry:*/
    5673                 this->Element::GetInputValue(&area,AreaEnum);
    5674 
    5675                 /*Compute fraction of the element that is grounded: */
    5676                 if(notfullygrounded){
    5677                         IssmDouble xyz_list[NUMVERTICES][3];
    5678                         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    5679 
    5680                         phi=this->GetGroundedPortion(&xyz_list[0][0]); //watch out, this only works because of the Thales theorem! We are in 3D, but this routine is inherently for 2D trias
    5681                         if(glfraction==0)phi=1;
    5682                         #ifdef _ISSM_DEBUG_
    5683                         this->AddInput(SealevelBarystaticMaskEnum,&phi,P0Enum);
    5684                         #endif
    5685                 }
    5686                 else phi=1.0;
    5687 
    5688                 /*Retrieve surface load for ice: */
    5689                 Input* deltathickness_input=this->GetInput(DeltaIceThicknessEnum);
    5690                 if (!deltathickness_input)_error_("delta thickness input needed to compute sea level change!");
    5691 
    5692                 /*/Average ice thickness over grounded area of the element only: {{{*/
    5693                 if(!notfullygrounded)deltathickness_input->GetInputAverage(&I);
    5694                 else{
    5695                         IssmDouble total_weight=0;
    5696                         bool mainlyfloating = true;
    5697                         int         point1;
    5698                         IssmDouble  fraction1,fraction2;
    5699 
    5700                         /*Recover portion of element that is grounded*/
    5701                         this->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating);
    5702                         Gauss* gauss = this->NewGauss(point1,fraction1,fraction2,mainlyfloating,2);
    5703 
    5704                         /* Start  looping on the number of gaussian points and average over these gaussian points: */
    5705                         total_weight=0;
    5706                         I=0;
    5707                         while(gauss->next()){
    5708                                 IssmDouble Ig=0;
    5709                                 deltathickness_input->GetInputValue(&Ig,gauss);
    5710                                 I+=Ig*gauss->weight;
    5711                                 total_weight+=gauss->weight;
    5712                         }
    5713                         if(total_weight) I=I/total_weight;
    5714                         delete gauss;
    5715                 }
    5716                 /*}}}*/
    5717 
    5718                 /*}}}*/
    5719                 /*convert to kg/m^2*/
    5720                 I=I*phi;
    5721 
    5722                 dI_list[0] = -4*M_PI*(rho_ice*I*area)*pow(re,4)*(sin(late)*cos(late)*cos(longe))/planetarea;
    5723                 dI_list[1] = -4*M_PI*(rho_ice*I*area)*pow(re,4)*(sin(late)*cos(late)*sin(longe))/planetarea;
    5724                 dI_list[2] = +4*M_PI*(rho_ice*I*area)*pow(re,4)*(1-pow(sin(late),2))/planetarea;
    5725         }
    5726 
    5727         return;
    5728 }/*}}}*/
    57295558void       Tria::SetSealevelMasks(SealevelMasks* masks){ /*{{{*/
    57305559
    … …  
    57405569        if ((gr_input->GetInputMin())<0) masks->notfullygrounded[this->lid]=true;
    57415570        else masks->notfullygrounded[this->lid]=false;
    5742 }
    5743 /*}}}*/
    5744 void       Tria::SealevelchangeGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze){ /*{{{*/
    5745 
    5746         /*diverse:*/
    5747         int gsize;
    5748         bool spherical=true;
    5749         IssmDouble llr_list[NUMVERTICES][3];
    5750         IssmDouble xyz_list[NUMVERTICES][3];
    5751         IssmDouble area,planetarea,planetradius;
    5752         IssmDouble I;  //change in ice thickness or water level(Farrel and Clarke, Equ. 4)
    5753         IssmDouble rho_earth;
    5754         IssmDouble late,longe,re;
    5755         IssmDouble lati,longi,ri;
    5756         IssmDouble constant;
    5757         IssmDouble x_element,y_element,z_element,x,y,z,dx,dy,dz,N_azim,E_azim;
    5758 
    5759         #ifdef _HAVE_RESTRICT_
    5760         IssmDouble* __restrict__ G=NULL;
    5761         IssmDouble* __restrict__ GU=NULL;
    5762         IssmDouble* __restrict__ GN=NULL;
    5763         IssmDouble* __restrict__ GE=NULL;
    5764         IssmDouble* __restrict__ G_elastic_precomputed=NULL;
    5765         IssmDouble* __restrict__ G_rigid_precomputed=NULL;
    5766         IssmDouble* __restrict__ U_elastic_precomputed=NULL;
    5767         IssmDouble* __restrict__ H_elastic_precomputed=NULL;
    5768         IssmDouble* __restrict__ indices=NULL;
    5769         #else
    5770         IssmDouble* G=NULL;
    5771         IssmDouble* GU=NULL;
    5772         IssmDouble* GN=NULL;
    5773         IssmDouble* GE=NULL;
    5774         IssmDouble* G_elastic_precomputed=NULL;
    5775         IssmDouble* G_rigid_precomputed=NULL;
    5776         IssmDouble* U_elastic_precomputed=NULL;
    5777         IssmDouble* H_elastic_precomputed=NULL;
    5778         IssmDouble* indices=NULL;
    5779         #endif
    5780 
    5781         /*elastic green function:*/
    5782         int index;
    5783         int M,Mtest;
    5784 
    5785         /*Computational flags:*/
    5786         bool computerigid = false;
    5787         bool computeelastic = false;
    5788         int  horiz;
    5789 
    5790         /*recover parameters: */
    5791         rho_earth=FindParam(MaterialsEarthDensityEnum);
    5792         this->parameters->FindParam(&computerigid,SolidearthSettingsRigidEnum);
    5793         this->parameters->FindParam(&computeelastic,SolidearthSettingsElasticEnum);
    5794         this->parameters->FindParam(&gsize,MeshNumberofverticesEnum);
    5795         this->parameters->FindParam(&planetarea,SolidearthPlanetAreaEnum);
    5796         this->parameters->FindParam(&planetradius,SolidearthPlanetRadiusEnum);
    5797         this->parameters->FindParam(&horiz,SolidearthSettingsHorizEnum);
    5798 
    5799         /*compute area and add to inputs:*/
    5800         area=GetAreaSpherical();
    5801         this->inputs->SetDoubleInput(AreaEnum,this->lid,area);
    5802         this->AddInput(SealevelAreaEnum,&area,P0Enum);
    5803         if(!computerigid)return;
    5804 
    5805         /*recover precomputed green function kernels:*/
    5806         DoubleVecParam* parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelchangeGRigidEnum)); _assert_(parameter);
    5807         parameter->GetParameterValueByPointer((IssmDouble**)&G_rigid_precomputed,&M);
    5808         _assert_(M>0);
    5809 
    5810         /*allocate indices:*/
    5811         indices=xNew<IssmDouble>(gsize);
    5812         G=xNewZeroInit<IssmDouble>(gsize);
    5813 
    5814         if(computeelastic){
    5815                 parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelchangeGElasticEnum)); _assert_(parameter);
    5816                 parameter->GetParameterValueByPointer((IssmDouble**)&G_elastic_precomputed,&Mtest); _assert_(Mtest==M);
    5817 
    5818                 parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelchangeHElasticEnum)); _assert_(parameter);
    5819                 parameter->GetParameterValueByPointer((IssmDouble**)&H_elastic_precomputed,&Mtest); _assert_(Mtest==M);
    5820 
    5821                 parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelchangeUElasticEnum)); _assert_(parameter);
    5822                 parameter->GetParameterValueByPointer((IssmDouble**)&U_elastic_precomputed,&Mtest); _assert_(Mtest==M);
    5823 
    5824                 GU=xNewZeroInit<IssmDouble>(gsize);
    5825                 if(horiz){
    5826                         GN=xNewZeroInit<IssmDouble>(gsize);
    5827                         GE=xNewZeroInit<IssmDouble>(gsize);
    5828                 }
    5829         }
    5830 
    5831         /* Where is the centroid of this element:*/
    5832 
    5833         /*retrieve coordinates: */
    5834         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    5835 
    5836         /*figure out gravity center of our element (Cartesian): */
    5837         x_element=(xyz_list[0][0]+xyz_list[1][0]+xyz_list[2][0])/3.0;
    5838         y_element=(xyz_list[0][1]+xyz_list[1][1]+xyz_list[2][1])/3.0;
    5839         z_element=(xyz_list[0][2]+xyz_list[1][2]+xyz_list[2][2])/3.0;
    5840 
    5841         /*compute gravity center in lat,long: */
    5842         late= asin(z_element/planetradius);
    5843         longe = atan2(y_element,x_element);
    5844 
    5845         constant=3/rho_earth*area/planetarea;
    5846 
    5847         for(int i=0;i<gsize;i++){
    5848                 IssmDouble alpha;
    5849                 IssmDouble delPhi,delLambda;
    5850 
    5851                 /*Compute alpha angle between centroid and current vertex and index into precomputed tables: */
    5852                 lati=latitude[i]/180.*M_PI; longi=longitude[i]/180.*M_PI;
    5853                 delPhi=fabs(lati-late); delLambda=fabs(longi-longe); if (delLambda>M_PI)delLambda=2*M_PI-delLambda;
    5854                 alpha=2.*asin(sqrt(pow(sin(delPhi/2),2)+cos(lati)*cos(late)*pow(sin(delLambda/2),2)));
    5855                 indices[i]=alpha/M_PI*reCast<IssmDouble,int>(M-1);
    5856                 index=reCast<int,IssmDouble>(indices[i]);
    5857                 _assert_(index>=0); _assert_(index<M);
    5858 
    5859                 /*Rigid earth gravitational perturbation: */
    5860                 G[i]+=G_rigid_precomputed[index];
    5861                
    5862                 if(computeelastic){
    5863                         G[i]+=G_elastic_precomputed[index];
    5864                 }
    5865                 G[i]=G[i]*constant;
    5866 
    5867                 /*Elastic components:*/
    5868                 if(computeelastic){
    5869                         GU[i] =  constant * U_elastic_precomputed[index];
    5870                         if(horiz){
    5871                                 /*Compute azimuths, both north and east components: */
    5872                                 x = xx[i]; y = yy[i]; z = zz[i];
    5873                                 if(latitude[i]==90){
    5874                                         x=1e-12; y=1e-12;
    5875                                 }
    5876                                 if(latitude[i]==-90){
    5877                                         x=1e-12; y=1e-12;
    5878                                 }
    5879                                 dx = x_element-x; dy = y_element-y; dz = z_element-z;
    5880                                 N_azim = (-z*x*dx-z*y*dy+(pow(x,2)+pow(y,2))*dz) /pow((pow(x,2)+pow(y,2))*(pow(x,2)+pow(y,2)+pow(z,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5);
    5881                                 E_azim = (-y*dx+x*dy) /pow((pow(x,2)+pow(y,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5);
    5882 
    5883                                 GN[i] = constant*H_elastic_precomputed[index]*N_azim;
    5884                                 GE[i] = constant*H_elastic_precomputed[index]*E_azim;
    5885                         }
    5886                 }
    5887         }
    5888 
    5889         /*Add in inputs:*/
    5890         this->inputs->SetArrayInput(SealevelchangeGEnum,this->lid,G,gsize);
    5891         if(computeelastic){
    5892                 this->inputs->SetArrayInput(SealevelchangeGUEnum,this->lid,GU,gsize);
    5893                 if(horiz){
    5894                         this->inputs->SetArrayInput(SealevelchangeGNEnum,this->lid,GN,gsize);
    5895                         this->inputs->SetArrayInput(SealevelchangeGEEnum,this->lid,GE,gsize);
    5896                 }
    5897         }
    5898         this->inputs->SetDoubleInput(AreaEnum,this->lid,area);
    5899         this->AddInput(SealevelAreaEnum,&area,P0Enum);
    5900 
    5901         /*Free allocations:*/
    5902         #ifdef _HAVE_RESTRICT_
    5903         delete indices;
    5904         delete G;
    5905         if(computeelastic){
    5906                 delete GU;
    5907                 if(horiz){
    5908                         delete GN;
    5909                         delete GE;
    5910                 }
    5911         }
    5912         #else
    5913         xDelete(indices);
    5914         xDelete(G);
    5915         if(computeelastic){
    5916                 xDelete(GU);
    5917                 if(horiz){
    5918                         xDelete(GN);
    5919                         xDelete(GE);
    5920                 }
    5921         }
    5922         #endif
    5923 
    5924         return;
    5925 }
    5926 /*}}}*/
    5927 IssmDouble Tria::SealevelchangeBarystaticIce(IssmDouble* Sgi, SealevelMasks* masks, Vector<IssmDouble>* barystatic_contribution, IssmDouble* partition, IssmDouble oceanarea){ /*{{{*/
    5928 
    5929         /*diverse:*/
    5930         int gsize;
    5931         IssmDouble area;
    5932         IssmDouble phi=1.0; //WARNING: do not touch this, default is entire elemnt contributes barystatic
    5933         IssmDouble I;  //change in ice thickness or water level(Farrel and Clarke, Equ. 4)
    5934         bool notfullygrounded=false;
    5935         bool scaleoceanarea= false;
    5936         bool computerigid= false;
    5937         int  glfraction=1;
    5938 
    5939         /*output: */
    5940         IssmDouble bslcice=0;
    5941 
    5942         /*elastic green function:*/
    5943         IssmDouble* G=NULL;
    5944 
    5945         /*ice properties: */
    5946         IssmDouble rho_ice,rho_water;
    5947 
    5948         /*constants:*/
    5949         IssmDouble constant=0;
    5950 
    5951         /*early return if we are not on an ice cap:*/
    5952         if(!masks->isiceonly[this->lid]){
    5953                 #ifdef _ISSM_DEBUG_
    5954                 constant=0; this->AddInput(SealevelBarystaticMaskEnum,&constant,P0Enum);
    5955                 #endif
    5956                 bslcice=0;
    5957                 return bslcice;
    5958         }
    5959 
    5960         /*early return if we are fully floating:*/
    5961         if (masks->isfullyfloating[this->lid]){
    5962                 constant=0;
    5963                 #ifdef _ISSM_DEBUG_
    5964                 this->AddInput(SealevelBarystaticMaskEnum,&constant,P0Enum);
    5965                 #endif
    5966                 bslcice=0;
    5967                 return bslcice;
    5968         }
    5969 
    5970         /*If we are here, we are on ice that is fully grounded or half-way to floating:*/
    5971         if (masks->notfullygrounded[this->lid]) notfullygrounded=true; //used later on.
    5972 
    5973         /*Inform mask: */
    5974         constant=1;
    5975         #ifdef _ISSM_DEBUG_
    5976         this->AddInput(SealevelBarystaticMaskEnum,&constant,P0Enum);
    5977         #endif
    5978 
    5979         /*recover some parameters:*/
    5980         rho_ice=FindParam(MaterialsRhoIceEnum);
    5981         rho_water=FindParam(MaterialsRhoSeawaterEnum);
    5982         this->parameters->FindParam(&computerigid,SolidearthSettingsRigidEnum);
    5983         this->parameters->FindParam(&scaleoceanarea,SolidearthSettingsOceanAreaScalingEnum);
    5984         this->parameters->FindParam(&glfraction,SolidearthSettingsGlfractionEnum);
    5985 
    5986         /*retrieve precomputed G:*/
    5987         if(computerigid)this->inputs->GetArrayPtr(SealevelchangeGEnum,this->lid,&G,&gsize);
    5988 
    5989         /*Get area of element: precomputed in the sealevelchange_geometry:*/
    5990         this->Element::GetInputValue(&area,AreaEnum);
    5991 
    5992         /*Compute fraction of the element that is grounded: */
    5993         if(notfullygrounded){
    5994                 IssmDouble xyz_list[NUMVERTICES][3];
    5995                 ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    5996 
    5997                 phi=this->GetGroundedPortion(&xyz_list[0][0]); //watch out, this only works because of the Thales theorem! We are in 3D, but this routine is inherently for 2D trias
    5998                 if(glfraction==0)phi=1;
    5999                 #ifdef _ISSM_DEBUG_
    6000                 this->AddInput(SealevelBarystaticMaskEnum,&phi,P0Enum);
    6001                 #endif
    6002         }
    6003         else phi=1.0;
    6004 
    6005         /*Retrieve surface load for ice: */
    6006         Input* deltathickness_input=this->GetInput(DeltaIceThicknessEnum);
    6007         if (!deltathickness_input)_error_("delta thickness input needed to compute sea level change!");
    6008 
    6009         /*/Average ice thickness over grounded area of the element only: {{{*/
    6010         if(!notfullygrounded)deltathickness_input->GetInputAverage(&I);
    6011         else{
    6012                 IssmDouble total_weight=0;
    6013                 bool mainlyfloating = true;
    6014                 int         point1;
    6015                 IssmDouble  fraction1,fraction2;
    6016 
    6017                 /*Recover portion of element that is grounded*/
    6018                 this->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating);
    6019                 Gauss* gauss = this->NewGauss(point1,fraction1,fraction2,mainlyfloating,2);
    6020 
    6021                 /* Start  looping on the number of gaussian points and average over these gaussian points: */
    6022                 total_weight=0;
    6023                 I=0;
    6024                 while(gauss->next()){
    6025                         IssmDouble Ig=0;
    6026                         deltathickness_input->GetInputValue(&Ig,gauss);
    6027                         I+=Ig*gauss->weight;
    6028                         total_weight+=gauss->weight;
    6029                 }
    6030                 I=I/total_weight;
    6031                 delete gauss;
    6032         }
    6033         /*}}}*/
    6034 
    6035         /*Compute barystatic contribution:*/
    6036         _assert_(oceanarea>0.);
    6037         if(scaleoceanarea) oceanarea=3.619e+14; // use true ocean area, m^2
    6038         bslcice = rho_ice*area*phi*I/(oceanarea*rho_water);
    6039         _assert_(!xIsNan<IssmDouble>(bslcice));
    6040 
    6041         if(computerigid){
    6042                 /*convert from m to kg/m^2:*/
    6043                 I=I*rho_ice*phi;
    6044 
    6045                 /*convolve:*/
    6046                 for(int i=0;i<gsize;i++) Sgi[i]+=G[i]*I;
    6047         }
    6048 
    6049         /*Plug bslcice into barystatic contribution vector:*/
    6050         if(barystatic_contribution){
    6051                 id=reCast<int>(partition[this->Sid()])+1;
    6052                 barystatic_contribution->SetValue(id,bslcice,ADD_VAL);
    6053         }
    6054         /*return :*/
    6055         return bslcice;
    6056 }
    6057 /*}}}*/
    6058 IssmDouble Tria::SealevelchangeBarystaticHydro(IssmDouble* Sgi, SealevelMasks* masks, Vector<IssmDouble>* barystatic_contribution, IssmDouble* partition, IssmDouble oceanarea){ /*{{{*/
    6059 
    6060         /*diverse:*/
    6061         int gsize;
    6062         IssmDouble area;
    6063         IssmDouble phi=1.0; //WARNING: do not touch this, default is entire elemnt contributes barystatic
    6064         IssmDouble W;  //change in water height thickness (Farrel and Clarke, Equ. 4)
    6065         bool notfullygrounded=false;
    6066         bool scaleoceanarea= false;
    6067         bool computeelastic= false;
    6068 
    6069         /*elastic green function:*/
    6070         IssmDouble* G=NULL;
    6071 
    6072         /*ice properties: */
    6073         IssmDouble rho_water;
    6074         IssmDouble rho_freshwater;
    6075 
    6076         /*output:*/
    6077         IssmDouble bslchydro = 0;
    6078 
    6079         /*early return if we are on an ice cap. Nop, we may have hydro
    6080          * and ice on the same masscon:*/
    6081         /*if(masks->isiceonly[this->lid]){
    6082                 bslchydro=0;
    6083                 return bslchydro;
    6084         }*/
    6085 
    6086         /*early return if we are fully floating:*/
    6087         if (masks->isfullyfloating[this->lid]){
    6088                 bslchydro=0;
    6089                 return bslchydro;
    6090         }
    6091 
    6092         /*If we are here, we are on earth, not on ice: */
    6093 
    6094         /*recover parameters: */
    6095         rho_water=FindParam(MaterialsRhoSeawaterEnum);
    6096         rho_freshwater=FindParam(MaterialsRhoFreshwaterEnum);
    6097         this->parameters->FindParam(&computeelastic,SolidearthSettingsElasticEnum);
    6098         this->parameters->FindParam(&scaleoceanarea,SolidearthSettingsOceanAreaScalingEnum);
    6099 
    6100         /*retrieve precomputed G:*/
    6101         if(computeelastic)this->inputs->GetArrayPtr(SealevelchangeGEnum,this->lid,&G,&gsize);
    6102 
    6103         /*Get area of element: precomputed in the sealevelchange_geometry:*/
    6104         this->Element::GetInputValue(&area,AreaEnum);
    6105 
    6106         /*Retrieve water height at vertices: */
    6107         Input* deltathickness_input=this->GetInput(DeltaTwsEnum);
    6108         if (!deltathickness_input)_error_("DeltaTwsEnum input needed to compute sea level change!");
    6109         deltathickness_input->GetInputAverage(&W);
    6110 
    6111         /*Compute barystatic component:*/
    6112         _assert_(oceanarea>0.);
    6113         if(scaleoceanarea) oceanarea=3.619e+14; // use true ocean area, m^2
    6114         bslchydro = rho_freshwater*area*phi*W/(oceanarea*rho_water);
    6115         _assert_(!xIsNan<IssmDouble>(bslchydro));
    6116 
    6117         if(computeelastic){
    6118                 /*convert from m to kg/m^2:*/
    6119                 W=W*rho_freshwater*phi;
    6120 
    6121                 /*convolve:*/
    6122                 for(int i=0;i<gsize;i++) Sgi[i]+=G[i]*W;
    6123         }
    6124 
    6125         /*Plug bslcice into barystatic contribution vector:*/
    6126         if(barystatic_contribution){
    6127                 id=reCast<int>(partition[this->Sid()])+1;
    6128                 barystatic_contribution->SetValue(id,bslchydro,ADD_VAL);
    6129         }
    6130         /*output:*/
    6131         return bslchydro;
    6132 }
    6133 /*}}}*/
    6134 void       Tria::SealevelchangeBarystaticBottomPressure(IssmDouble* Sgi,SealevelMasks* masks){ /*{{{*/
    6135 
    6136         /*diverse:*/
    6137         int gsize;
    6138         IssmDouble area;
    6139         IssmDouble BP;  //change in bottom pressure (Farrel and Clarke, Equ. 4)
    6140         IssmDouble constant;
    6141         bool computeelastic= false;
    6142 
    6143         /*elastic green function:*/
    6144         IssmDouble* G=NULL;
    6145 
    6146         /*water properties: */
    6147         IssmDouble rho_water;
    6148 
    6149         /*exit now?:*/
    6150         this->parameters->FindParam(&computeelastic,SolidearthSettingsElasticEnum);
    6151         if(!computeelastic)return;
    6152 
    6153         /*we are here to compute fingerprints originating fromn bottom pressure changes:*/
    6154         if(!masks->isoceanin[this->lid])return;
    6155 
    6156         /*Inform mask: */
    6157         #ifdef _ISSM_DEBUG_
    6158         constant=1; this->AddInput(SealevelBarystaticMaskEnum,&constant,P0Enum);
    6159         #endif
    6160 
    6161         /*recover material parameters: */
    6162         rho_water=FindParam(MaterialsRhoSeawaterEnum);
    6163 
    6164         /*retrieve precomputed G:*/
    6165         this->inputs->GetArrayPtr(SealevelchangeGEnum,this->lid,&G,&gsize);
    6166 
    6167         /*Get area of element: precomputed in the sealevelchange_geometry:*/
    6168         this->Element::GetInputValue(&area,AreaEnum);
    6169 
    6170         /*Retrieve bottom pressure change and average over the element: */
    6171         Input* bottompressure_change_input=this->GetInput(DeltaBottomPressureEnum);
    6172         if (!bottompressure_change_input)_error_("DeltaBottomPressureEnum pressure input needed to compute sea level change fingerprint!");
    6173         bottompressure_change_input->GetInputAverage(&BP);
    6174 
    6175         /*convert from m to kg/m^2:*/
    6176         BP=BP*rho_water;
    6177 
    6178         /*convolve:*/
    6179         for(int i=0;i<gsize;i++) Sgi[i]+=G[i]*BP;
    6180 
    6181         return;
    6182 }
    6183 /*}}}*/
    6184 void       Tria::SealevelchangeSal(IssmDouble* Sgo,IssmDouble* Sg_old, SealevelMasks* masks){ /*{{{*/
    6185 
    6186         /*diverse:*/
    6187         int gsize,dummy;
    6188         IssmDouble S;  //change in water water level(Farrel and Clarke, Equ. 4)
    6189         IssmDouble constant=0;
    6190         IssmDouble rho_water;
    6191         IssmDouble* G=NULL;
    6192 
    6193         /*retrieve parameters:*/
    6194         rho_water=FindParam(MaterialsRhoSeawaterEnum);
    6195 
    6196         /*early return if we are not on the ocean:*/
    6197         if (!masks->isoceanin[this->lid]){
    6198                 constant=0;
    6199                 #ifdef _ISSM_DEBUG_
    6200                 this->AddInput(SealevelBarystaticOceanMaskEnum,&constant,P0Enum);
    6201                 #endif
    6202                 return;
    6203         }
    6204         constant=1;
    6205         #ifdef _ISSM_DEBUG_
    6206         this->AddInput(SealevelBarystaticOceanMaskEnum,&constant,P0Enum);
    6207         #endif
    6208 
    6209         /*how many dofs are we working with here? */
    6210         this->parameters->FindParam(&gsize,MeshNumberofverticesEnum);
    6211 
    6212         /*retrieve precomputed G:*/
    6213         this->inputs->GetArrayPtr(SealevelchangeGEnum,this->lid,&G,&dummy); _assert_(dummy==gsize);
    6214 
    6215         /*From Sg_old, recover water sea level change:*/
    6216         S=0; for(int i=0;i<NUMVERTICES;i++) S+=Sg_old[this->vertices[i]->Sid()]/NUMVERTICES;
    6217 
    6218         /*convert to kg/m^2: */
    6219         S=S*rho_water;
    6220 
    6221         for(int i=0;i<gsize;i++) Sgo[i]+=G[i]*S;
    6222 
    6223         return;
    6224 }
    6225 /*}}}*/
    6226 void       Tria::DeformationFromSurfaceLoads(IssmDouble* Up, IssmDouble* North ,IssmDouble* East,IssmDouble* Sg, SealevelMasks* masks){ /*{{{*/
    6227 
    6228         /*diverse:*/
    6229         int gsize;
    6230         IssmDouble I=0;
    6231         IssmDouble S=0;
    6232         IssmDouble BP=0;
    6233         IssmDouble rho_ice,rho_water;
    6234         int horiz;
    6235         int  bp_compute_fingerprints= 0;
    6236 
    6237         /*precomputed elastic green functions:*/
    6238         IssmDouble* GU=NULL;
    6239         IssmDouble* GN=NULL;
    6240         IssmDouble* GE=NULL;
    6241 
    6242         /*computational flags:*/
    6243         bool computeelastic= false;
    6244         bool computerigid= false;
    6245         bool notfullygrounded=false;
    6246         bool scaleoceanarea= false;
    6247         int  glfraction=1;
    6248         IssmDouble area;
    6249         IssmDouble phi=1.0;
    6250 
    6251         /*early return if we are not on the ocean or on an ice cap:*/
    6252         if(!masks->isiceonly[this->lid] && !masks->isoceanin[this->lid]) return;
    6253 
    6254         /*recover parameters:*/
    6255         this->parameters->FindParam(&computeelastic,SolidearthSettingsElasticEnum);
    6256         this->parameters->FindParam(&horiz,SolidearthSettingsHorizEnum);
    6257         this->parameters->FindParam(&bp_compute_fingerprints,SolidearthSettingsComputeBpGrdEnum);
    6258 
    6259         /*early return if elastic not requested:*/
    6260         if(!computeelastic) return;
    6261 
    6262         /*recover material parameters: */
    6263         rho_ice=FindParam(MaterialsRhoIceEnum);
    6264         rho_water=FindParam(MaterialsRhoSeawaterEnum);
    6265 
    6266         /*recover elastic Green's functions for displacement:*/
    6267         this->inputs->GetArrayPtr(SealevelchangeGUEnum,this->lid,&GU,&gsize);
    6268         if(horiz){
    6269                 this->inputs->GetArrayPtr(SealevelchangeGEEnum,this->lid,&GE,&gsize);
    6270                 this->inputs->GetArrayPtr(SealevelchangeGNEnum,this->lid,&GN,&gsize);
    6271         }
    6272 
    6273 
    6274         if(masks->isoceanin[this->lid]){
    6275                 /*From Sg, recover water sea level change:*/
    6276                 S=0; for(int i=0;i<NUMVERTICES;i++) S+=Sg[this->vertices[i]->Sid()]/NUMVERTICES;
    6277 
    6278                 /*convert to kg/m^2:*/
    6279                 S=rho_water*S;
    6280                
    6281                 /*If bottom pressures are available, retrieve them to load the bedrock:*/
    6282                 if(bp_compute_fingerprints){
    6283                         Input* bottompressure_change_input=this->GetInput(DeltaBottomPressureEnum);
    6284                         if (!bottompressure_change_input)_error_("bottom pressure input needed to compute sea level change fingerprint!");
    6285                         bottompressure_change_input->GetInputAverage(&BP);
    6286 
    6287                         /*convert from m to kg/m^2:*/
    6288                         BP=BP*rho_water;
    6289 
    6290                         S+=BP;
    6291                 }
    6292 
    6293                 for(int i=0;i<gsize;i++){
    6294                         Up[i]+=S*GU[i];
    6295                         if(horiz){
    6296                                 North[i]+=S*GN[i];
    6297                                 East[i]+=S*GE[i];
    6298                         }
    6299                 }
    6300         }
    6301         if (masks->isiceonly[this->lid]){
    6302 
    6303                 if (masks->isfullyfloating[this->lid]) return;
    6304 
    6305        
    6306                 if (masks->notfullygrounded[this->lid]) notfullygrounded=true; //used later on.
    6307 
    6308                 /*recover some parameters:*/
    6309                 rho_ice=FindParam(MaterialsRhoIceEnum);
    6310                 rho_water=FindParam(MaterialsRhoSeawaterEnum);
    6311                 this->parameters->FindParam(&computerigid,SolidearthSettingsRigidEnum);
    6312                 this->parameters->FindParam(&scaleoceanarea,SolidearthSettingsOceanAreaScalingEnum);
    6313                 this->parameters->FindParam(&glfraction,SolidearthSettingsGlfractionEnum);
    6314 
    6315                 /*Get area of element: precomputed in the sealevelchange_geometry:*/
    6316                 this->Element::GetInputValue(&area,AreaEnum);
    6317 
    6318                 /*Compute fraction of the element that is grounded: */
    6319                 if(notfullygrounded){
    6320                         IssmDouble xyz_list[NUMVERTICES][3];
    6321                         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    6322 
    6323                         phi=this->GetGroundedPortion(&xyz_list[0][0]); //watch out, this only works because of the Thales theorem! We are in 3D, but this routine is inherently for 2D trias
    6324                         if(glfraction==0)phi=1;
    6325                         #ifdef _ISSM_DEBUG_
    6326                         this->AddInput(SealevelBarystaticMaskEnum,&phi,P0Enum);
    6327                         #endif
    6328                 }
    6329                 else phi=1.0;
    6330 
    6331                 /*Retrieve surface load for ice: */
    6332                 Input* deltathickness_input=this->GetInput(DeltaIceThicknessEnum);
    6333                 if (!deltathickness_input)_error_("delta thickness input needed to compute sea level change!");
    6334 
    6335                 /*/Average ice thickness over grounded area of the element only: {{{*/
    6336                 if(!notfullygrounded)deltathickness_input->GetInputAverage(&I);
    6337                 else{
    6338                         IssmDouble total_weight=0;
    6339                         bool mainlyfloating = true;
    6340                         int         point1;
    6341                         IssmDouble  fraction1,fraction2;
    6342 
    6343                         /*Recover portion of element that is grounded*/
    6344                         this->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating);
    6345                         Gauss* gauss = this->NewGauss(point1,fraction1,fraction2,mainlyfloating,2);
    6346 
    6347                         /* Start  looping on the number of gaussian points and average over these gaussian points: */
    6348                         total_weight=0;
    6349                         I=0;
    6350                         while(gauss->next()){
    6351                                 IssmDouble Ig=0;
    6352                                 deltathickness_input->GetInputValue(&Ig,gauss);
    6353                                 I+=Ig*gauss->weight;
    6354                                 total_weight+=gauss->weight;
    6355                         }
    6356                         I=I/total_weight;
    6357                         delete gauss;
    6358                 }
    6359                 /*}}}*/
    6360 
    6361                 /*convert to kg/m^2*/
    6362                 I=I*rho_ice*phi;
    6363                
    6364                 for(int i=0;i<gsize;i++){
    6365                         Up[i]+=I*GU[i];
    6366                         if(horiz){
    6367                                 North[i]+=I*GN[i];
    6368                                 East[i]+=I*GE[i];
    6369                         }
    6370                 }
    6371         }
    6372 
    6373         return;
    63745571}
    63755572/*}}}*/
    … …  
    64795676}
    64805677/*}}}*/
    6481 
    6482 //new code
    6483 void       Tria::SealevelchangeGeometryOptim(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae){ /*{{{*/
     5678void       Tria::SealevelchangeGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae){ /*{{{*/
    64845679
    64855680        /*diverse:*/
    … …  
    70986293}
    70996294/*}}}*/
    7100 void       Tria::SealevelchangeMomentOfInertiaElement(IssmDouble* dI_list, IssmDouble* loads, IssmDouble* sealevelloads){/*{{{*/
     6295void       Tria::SealevelchangeMomentOfInertia(IssmDouble* dI_list, IssmDouble* loads, IssmDouble* sealevelloads){/*{{{*/
    71016296               
    71026297        IssmDouble S=0;

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

    r26164 r26165  
    165165                #endif
    166166                #ifdef _HAVE_SEALEVELCHANGE_
    167                 void       SealevelchangeMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old,SealevelMasks* masks);
    168                 void       SealevelchangeMomentOfInertiaElement(IssmDouble* dI_list, IssmDouble* loads, IssmDouble* sealevelloads);
     167                void       SealevelchangeMomentOfInertia(IssmDouble* dI_list, IssmDouble* loads, IssmDouble* sealevelloads);
    169168                void       SealevelchangeGeometry(IssmDouble* lat, IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze);
    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);
    175169                void       GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt,Matlitho* litho, IssmDouble* x,IssmDouble* y);
    176                 IssmDouble OceanAverage(IssmDouble* Sg, SealevelMasks* masks);
    177170                void       SetSealevelMasks(SealevelMasks* masks);
    178                
    179                 void       SealevelchangeGeometryOptim(IssmDouble* lat,IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae);
     171                void       SealevelchangeGeometry(IssmDouble* lat,IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae);
    180172                void       SealevelchangeBarystaticLoads(Vector<IssmDouble>* loads, BarystaticContributions* barycontrib, SealevelMasks* masks);
    181173                void       SealevelchangeConvolution(Vector<IssmDouble>* sealevelloads, Vector<IssmDouble>* oceanareas, IssmDouble* allsealevelloads, IssmDouble* allloads,IssmDouble* rotationaxismotionvector);

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

    r26114 r26165  
    47484748/*}}}*/
    47494749#endif
    4750 #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 /*}}}*/
    4909 void FemModel::SealevelchangeSal(Vector<IssmDouble>* pRSLgo, Vector<IssmDouble>* pRSLg_old,  SealevelMasks* masks, bool verboseconvolution){/*{{{*/
    4910 
    4911         /*serialized vectors:*/
    4912         IssmDouble* RSLg_old=NULL;
    4913 
    4914         IssmDouble* RSLgo  = NULL;
    4915         int* indices = NULL;
    4916         int         gsize;
    4917 
    4918         bool computerigid = true;
    4919         bool computeelastic= true;
    4920 
    4921         /*recover computational flags: */
    4922         this->parameters->FindParam(&computerigid,SolidearthSettingsRigidEnum);
    4923 
    4924         /*Initialize temporary vector that will be used to sum barystatic components on all local elements, prior
    4925          * to assembly:*/
    4926         gsize = this->nodes->NumberOfDofs(GsetEnum);
    4927         RSLgo=xNewZeroInit<IssmDouble>(gsize);
    4928         indices=xNew<int>(gsize); for (int i=0;i<gsize;i++)indices[i]=i;
    4929 
    4930         /*Serialize vectors from previous iteration:*/
    4931         RSLg_old=pRSLg_old->ToMPISerial();
    4932 
    4933         /*Call the sal sea level change core only if required: */
    4934         if(computerigid){
    4935                 for(Object* & object : this->elements->objects){
    4936                         Element* element = xDynamicCast<Element*>(object);
    4937                         element->SealevelchangeSal(RSLgo,RSLg_old,masks);
    4938                 }
    4939         }
    4940         pRSLgo->SetValues(gsize,indices,RSLgo,ADD_VAL);
    4941         pRSLgo->Assemble();
    4942 
    4943         /*Free ressources:*/
    4944         xDelete<int>(indices);
    4945         xDelete<IssmDouble>(RSLgo);
    4946         xDelete<IssmDouble>(RSLg_old);
    4947 
    4948 }
    4949 /*}}}*/
    4950 void FemModel::SealevelchangeRotationalFeedback(Vector<IssmDouble>* pRSLgo_rot, Vector<IssmDouble>* pRSLg_old, IssmDouble* pIxz, IssmDouble* pIyz, IssmDouble* pIzz,  SealevelMasks* masks){/*{{{*/
    4951 
    4952         /*serialized vectors:*/
    4953         bool spherical=true;
    4954         IssmDouble* RSLg_old=NULL;
    4955         IssmDouble*     tide_love_h  = NULL;
    4956         IssmDouble*     tide_love_k  = NULL;
    4957         IssmDouble*     load_love_k  = NULL;
    4958         IssmDouble  tide_love_k2secular;
    4959         IssmDouble  moi_e, moi_p, omega, g;
    4960         IssmDouble      m1, m2, m3;
    4961         IssmDouble      lati, longi, radi, value;
    4962         IssmDouble          *latitude    = NULL;
    4963         IssmDouble          *longitude    = NULL;
    4964         IssmDouble          *radius    = NULL;
    4965 
    4966         /*Serialize vectors from previous iteration:*/
    4967         RSLg_old=pRSLg_old->ToMPISerial();
    4968 
    4969         IssmDouble moi_list[3]={0,0,0};
    4970         IssmDouble moi_list_cpu[3]={0,0,0};
    4971         for(Object* & object : this->elements->objects){
    4972                 Element* element = xDynamicCast<Element*>(object);
    4973                 element->SealevelchangeMomentOfInertia(&moi_list[0],RSLg_old,masks );
    4974                 moi_list_cpu[0] += moi_list[0];
    4975                 moi_list_cpu[1] += moi_list[1];
    4976                 moi_list_cpu[2] += moi_list[2];
    4977         }
    4978         ISSM_MPI_Reduce (&moi_list_cpu[0],&moi_list[0],1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
    4979         ISSM_MPI_Bcast(&moi_list[0],1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
    4980         //
    4981         ISSM_MPI_Reduce (&moi_list_cpu[1],&moi_list[1],1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
    4982         ISSM_MPI_Bcast(&moi_list[1],1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
    4983         //
    4984         ISSM_MPI_Reduce (&moi_list_cpu[2],&moi_list[2],1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
    4985         ISSM_MPI_Bcast(&moi_list[2],1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
    4986 
    4987         /*pull out some useful parameters: */
    4988         parameters->FindParam(&load_love_k,NULL,NULL,LoadLoveKEnum);
    4989         parameters->FindParam(&tide_love_h,NULL,NULL,TidalLoveHEnum);
    4990         parameters->FindParam(&tide_love_k,NULL,NULL,TidalLoveKEnum);
    4991         parameters->FindParam(&tide_love_k2secular,TidalLoveK2SecularEnum);
    4992         parameters->FindParam(&moi_e,RotationalEquatorialMoiEnum);
    4993         parameters->FindParam(&moi_p,RotationalPolarMoiEnum);
    4994         parameters->FindParam(&omega,RotationalAngularVelocityEnum);
    4995 
    4996         /*compute perturbation terms for angular velocity vector: */
    4997         m1 = 1/(1-tide_love_k[2]/tide_love_k2secular) * (1+load_love_k[2])/(moi_p-moi_e) * moi_list[0];
    4998         m2 = 1/(1-tide_love_k[2]/tide_love_k2secular) * (1+load_love_k[2])/(moi_p-moi_e) * moi_list[1];
    4999         m3 = -(1+load_love_k[2])/moi_p * moi_list[2];   // term associated with fluid number (3-order-of-magnitude smaller) is negelected
    5000 
    5001         /*recover lat,long and radius vectors from vertices: */
    5002         VertexCoordinatesx(&latitude,&longitude,&radius,this->vertices,spherical);
    5003 
    5004         /* Green's function (1+k_2-h_2/g): checked against Glenn Milne's thesis Chapter 3 (eqs: 3.3-4, 3.10-11)
    5005          * Perturbation terms for angular velocity vector (m1, m2, m3): checked against Mitrovica (2005 Appendix) & Jensen et al (2013 Appendix A3)
    5006          * Sea level rotational feedback: checked against GMD eqs 8-9 (only first order terms, i.e., degree 2 order 0 & 1 considered)
    5007          * all DONE in Geographic coordinates: theta \in [-90,90], lambda \in [-180 180]
    5008          */
    5009         for(Object* & object : vertices->objects){
    5010                 Vertex* vertex=xDynamicCast<Vertex*>(object);
    5011                 int sid=vertex->Sid();
    5012 
    5013                 lati=latitude[sid]/180*PI;
    5014                 longi=longitude[sid]/180*PI;
    5015                 radi=radius[sid];
    5016 
    5017                 /*only first order terms are considered now: */
    5018                 value=((1.0+tide_love_k[2]-tide_love_h[2])/9.81)*pow(omega*radi,2.0)*
    5019                                                 (-m3/6.0 + 0.5*m3*cos(2.0*lati) - 0.5*sin(2.*lati)*(m1*cos(longi)+m2*sin(longi)));
    5020 
    5021                 pRSLgo_rot->SetValue(sid,value,INS_VAL); //INS_VAL ensures that you don't add several times
    5022         }
    5023 
    5024         /*Assemble mesh velocity*/
    5025         pRSLgo_rot->Assemble();
    5026 
    5027         /*Assign output pointers:*/
    5028         if(pIxz)*pIxz=moi_list[0];
    5029         if(pIyz)*pIyz=moi_list[1];
    5030         if(pIzz)*pIzz=moi_list[2];
    5031         xDelete<IssmDouble>(latitude);
    5032         xDelete<IssmDouble>(longitude);
    5033         xDelete<IssmDouble>(tide_love_h);
    5034         xDelete<IssmDouble>(tide_love_k);
    5035         xDelete<IssmDouble>(load_love_k);
    5036 
    5037         xDelete<IssmDouble>(radius);
    5038 
    5039         /*Free ressources:*/
    5040         xDelete<IssmDouble>(RSLg_old);
    5041 
    5042 }
    5043 /*}}}*/
    5044 void FemModel::SealevelchangeDeformation(Vector<IssmDouble>* pgeoid, Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, Vector<IssmDouble>* pRSLg, SealevelMasks* masks){/*{{{*/
    5045 
    5046         /*serialized vectors:*/
    5047         IssmDouble* RSLg=NULL;
    5048 
    5049         IssmDouble* Up  = NULL;
    5050         IssmDouble* North  = NULL;
    5051         IssmDouble* East  = NULL;
    5052         int* indices = NULL;
    5053         int  gsize;
    5054         int  horiz;
    5055 
    5056         /*retrieve parameters:*/
    5057         this->parameters->FindParam(&horiz,SolidearthSettingsHorizEnum);
    5058 
    5059         /*Serialize vectors from previous iteration:*/
    5060         RSLg=pRSLg->ToMPISerial();
    5061 
    5062         /*Initialize temporary vector that will be used to sum barystatic components on all local elements, prior
    5063          * to assembly:*/
    5064         gsize = this->nodes->NumberOfDofs(GsetEnum);
    5065         Up=xNewZeroInit<IssmDouble>(gsize);
    5066         if(horiz){
    5067                 North=xNewZeroInit<IssmDouble>(gsize);
    5068                 East=xNewZeroInit<IssmDouble>(gsize);
    5069         }
    5070         indices=xNew<int>(gsize); for (int i=0;i<gsize;i++)indices[i]=i;
    5071 
    5072         /*Call the deformation from loading routines:*/
    5073         for(Object* & object : this->elements->objects){
    5074                 Element* element = xDynamicCast<Element*>(object);
    5075                 element->DeformationFromSurfaceLoads(Up,North,East,RSLg,masks);
    5076         }
    5077 
    5078         pUp->SetValues(gsize,indices,Up,ADD_VAL);
    5079         pUp->Assemble();
    5080 
    5081 
    5082         /*Add RSL to Up to find the geoid: */
    5083         pUp->Copy(pgeoid); pgeoid->AXPY(pRSLg,1);
    5084 
    5085         if (horiz){
    5086                 pNorth->SetValues(gsize,indices,North,ADD_VAL);
    5087                 pNorth->Assemble();
    5088                 pEast->SetValues(gsize,indices,East,ADD_VAL);
    5089                 pEast->Assemble();
    5090         }
    5091 
    5092         /*Free ressources:*/
    5093         xDelete<IssmDouble>(Up);
    5094         if(horiz){
    5095                 xDelete<IssmDouble>(North);
    5096                 xDelete<IssmDouble>(East);
    5097         }
    5098         xDelete<int>(indices);
    5099         xDelete<IssmDouble>(RSLg);
    5100 }
    5101 /*}}}*/
    5102 IssmDouble FemModel::SealevelchangeOceanAverage(Vector<IssmDouble>* RSLg,SealevelMasks* masks, IssmDouble oceanarea) { /*{{{*/
    5103 
    5104         IssmDouble* RSLg_serial=NULL;
    5105         IssmDouble  oceanvalue,oceanvalue_cpu;
    5106 
    5107         /*Serialize vectors from previous iteration:*/
    5108         RSLg_serial=RSLg->ToMPISerial();
    5109 
    5110         /*Initialize:*/
    5111         oceanvalue_cpu=0;
    5112 
    5113         /*Go through elements, and add contribution from each element and divide by overall ocean area:*/
    5114         for(Object* & object : this->elements->objects){
    5115                 Element* element = xDynamicCast<Element*>(object);
    5116                 oceanvalue_cpu += element->OceanAverage(RSLg_serial,masks);
    5117         }
    5118 
    5119         ISSM_MPI_Reduce (&oceanvalue_cpu,&oceanvalue,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
    5120         ISSM_MPI_Bcast(&oceanvalue,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
    5121 
    5122         /*Free ressources:*/
    5123         xDelete<IssmDouble>(RSLg_serial);
    5124 
    5125         return oceanvalue/oceanarea;
    5126 }
    5127 /*}}}*/
    5128 void FemModel::IvinsDeformation(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt, IssmDouble* x, IssmDouble* y){ /*{{{*/
    5129 
    5130         /*Find the litho material to be used by all the elements:*/
    5131         Matlitho* matlitho=NULL;
    5132         for (Object* & object: this->materials->objects){
    5133                 Material* material=xDynamicCast<Material*>(object);
    5134                 if(material->ObjectEnum()==MatlithoEnum){
    5135                         matlitho=xDynamicCast<Matlitho*>(material);
    5136                         break;
    5137                 }
    5138         }
    5139         _assert_(matlitho);
    5140 
    5141 
    5142         /*Go through elements, and add contribution from each element to the deflection vector wg:*/
    5143         for(Object* & object : this->elements->objects){
    5144                 Element* element = xDynamicCast<Element*>(object);
    5145                 element->GiaDeflection(wg,dwgdt, matlitho, x,y);
    5146         }
    5147 
    5148         /*Assemble parallel vector:*/
    5149         dwgdt->Assemble();
    5150         wg->Assemble();
    5151 }
    5152 /*}}}*/
    5153 #endif
    51544750void FemModel::HydrologyEPLupdateDomainx(IssmDouble* pEplcount){ /*{{{*/
    51554751

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

    r26114 r26165  
    163163                void EsaGeodetic3D(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz);
    164164                #endif
    165                 #ifdef _HAVE_SEALEVELCHANGE_
    166                 void SealevelchangeBarystatic(Vector<IssmDouble>* pRSLgi, IssmDouble* poceanarea, IssmDouble* pbslc,IssmDouble* pbslcice, IssmDouble* pbslchydro, IssmDouble** pbslcice_partition,IssmDouble** pbslchydro_partition, SealevelMasks* masks);
    167                 void SealevelchangeSal(Vector<IssmDouble>* pSgo, Vector<IssmDouble>* pSg_old,  SealevelMasks* masks,bool verboseconvolution);
    168                 void SealevelchangeRotationalFeedback(Vector<IssmDouble>* pRSLgo_rot, Vector<IssmDouble>* pRSLg_old, IssmDouble* pIxz, IssmDouble* pIyz, IssmDouble* pIzz, SealevelMasks* masks);
    169                 void SealevelchangeDeformation(Vector<IssmDouble>* pgeoid,Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, Vector<IssmDouble>* pSg_old, SealevelMasks* masks);
    170                 IssmDouble SealevelchangeOceanAverage(Vector<IssmDouble>* Sg,SealevelMasks* masks, IssmDouble oceanarea);
    171                 void IvinsDeformation(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt, IssmDouble* x, IssmDouble* y);
    172                 #endif
    173165                void HydrologyEPLupdateDomainx(IssmDouble* pEplcount);
    174166                void HydrologyIDSupdateDomainx(IssmDouble* pIDScount);

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

    r26119 r26165  
    6363#endif
    6464void grd_core(FemModel* femmodel);
    65 void grd_core_optim(FemModel* femmodel);
    6665void solidearthexternal_core(FemModel* femmodel);
    6766void dynstr_core(FemModel* femmodel);
    68 Vector<IssmDouble>* sealevelchange_core_barystatic(FemModel* femmodel,SealevelMasks* mask, IssmDouble* poceanarea);
    69 Vector<IssmDouble>* sealevelchange_core_sal(FemModel* femmodel,SealevelMasks* masks, Vector<IssmDouble>* RSLg_barystatic,IssmDouble oceanarea);
    70 void sealevelchange_core_deformation(Vector<IssmDouble>** pN_radial, Vector<IssmDouble>** pU_radial, Vector<IssmDouble>** pU_north,Vector<IssmDouble>** pU_east,FemModel* femmodel,Vector<IssmDouble>* RSLg, SealevelMasks* masks);
    7167void couplerinput_core(FemModel* femmodel);
    7268void coupleroutput_core(FemModel* femmodel);

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

    r26164 r26165  
    2121void RotationAxisMotion(IssmDouble* m, FemModel* femmodel,IssmDouble* loads, IssmDouble* sealeveloads);
    2222void ConserveOceanMass(FemModel* femmodel,Vector<IssmDouble>* sealevelloads, IssmDouble offset, SealevelMasks* masks);
     23void ivins_deformation_core(FemModel* femmodel);
    2324/*}}}*/
    2425
    … …  
    3132        /*Parameters, variables:*/
    3233        bool save_results;
    33         int  optim=0;
    3434
    3535        /*Retrieve parameters:*/
    3636        femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
    37         femmodel->parameters->FindParam(&optim,SolidearthSettingsOptimEnum);
    3837       
    3938        /*Verbose: */
    … …  
    5352
    5453        /*Run geodetic:*/
    55         if(!optim) grd_core(femmodel);
    56         else grd_core_optim(femmodel);
     54        grd_core(femmodel);
    5755
    5856        /*Run steric core for sure:*/
    … …  
    208206       
    209207}; /*}}}*/
    210 void grd_core(FemModel* femmodel){ /*{{{*/
    211 
    212         /*Gravity rotation deformation core GRD: */
    213 
    214         /*variables:*/
    215         Vector<IssmDouble> *RSLg    = NULL;
    216         Vector<IssmDouble> *RSLg_barystatic  = NULL;
    217         Vector<IssmDouble> *U_grd  = NULL;
    218         Vector<IssmDouble> *N_grd  = NULL;
    219         Vector<IssmDouble> *U_north_grd   = NULL;
    220         Vector<IssmDouble> *U_east_grd    = NULL;
    221         Vector<IssmDouble> *bedrock  = NULL;
    222         Vector<IssmDouble> *bedrockeast  = NULL;
    223         Vector<IssmDouble> *bedrocknorth  = NULL;
    224         Vector<IssmDouble> *geoid= NULL;
     208void grd_core(FemModel* femmodel) { /*{{{*/
     209
     210        /*variables:{{{*/
     211        int nel;
     212        BarystaticContributions* barycontrib=NULL;
    225213        SealevelMasks* masks=NULL;
    226 
    227         /*parameters:*/
     214        GenericParam<BarystaticContributions*>* barycontribparam=NULL;
     215        IssmDouble rotationaxismotionvector[3]={0};
     216       
     217        Vector<IssmDouble>*    loads=NULL;
     218        IssmDouble*            allloads=NULL;
     219        Vector<IssmDouble>*    sealevelloads=NULL;
     220        Vector<IssmDouble>*    oldsealevelloads=NULL;
     221        IssmDouble*            allsealevelloads=NULL;
     222        Vector<IssmDouble>*    oceanareas=NULL;
     223        IssmDouble             oceanarea;
     224        IssmDouble             oceanaverage;
     225        bool                   scaleoceanarea=false;
     226        IssmDouble             rho_water;
     227
     228        IssmDouble           eps_rel;
     229        IssmDouble           eps_abs;
     230        int                  max_nonlinear_iterations;
     231        int                  iterations=0;
     232        int                  step;
     233        IssmDouble           time;
     234        bool converged=false;
     235
    228236        int  modelid,earthid;
    229237        int  horiz;
    230         IssmDouble oceanarea;
    231238        int  count,frequency,iscoupling;
    232239        int  grd=0;
     240        int  grdmodel;
     241        int  computesealevel=0;
     242
     243        /*}}}*/
    233244
    234245        /*Verbose: */
    235246        if(VerboseSolution()) _printf0_("         computing GRD patterns\n");
    236        
    237         /*retrieve parameters:*/
     247
     248        /*retrieve parameters:{{{*/
    238249        femmodel->parameters->FindParam(&grd,SolidearthSettingsGRDEnum);
    239250        femmodel->parameters->FindParam(&frequency,SolidearthSettingsRunFrequencyEnum);
    240251        femmodel->parameters->FindParam(&count,SealevelchangeRunCountEnum);
     252        femmodel->parameters->FindParam(&computesealevel,SolidearthSettingsComputesealevelchangeEnum);
     253        femmodel->parameters->FindParam(&max_nonlinear_iterations,SolidearthSettingsMaxiterEnum);
     254        femmodel->parameters->FindParam(&grdmodel,GrdModelEnum);
     255        /*}}}*/
    241256
    242257        /*only run if grd was requested, if we are the earth, and we have reached
    … …  
    251266        }
    252267
     268        /*branch directly to Ivins deformation core if requested:*/
     269        if(grdmodel==IvinsEnum){
     270                ivins_deformation_core(femmodel);
     271                return;
     272        }
     273
     274        /*retrieve parameters: {{{*/
     275        femmodel->parameters->FindParam(&scaleoceanarea,SolidearthSettingsOceanAreaScalingEnum);
     276        barycontribparam = xDynamicCast<GenericParam<BarystaticContributions*>*>(femmodel->parameters->FindParamObject(BarystaticContributionsEnum));
     277        barycontrib=barycontribparam->GetParameterValue();
     278        femmodel->parameters->FindParam(&rho_water,MaterialsRhoSeawaterEnum);
     279        femmodel->parameters->FindParam(&eps_rel,SolidearthSettingsReltolEnum);
     280        femmodel->parameters->FindParam(&eps_abs,SolidearthSettingsAbstolEnum);
     281        femmodel->parameters->FindParam(&horiz,SolidearthSettingsHorizEnum);
     282        /*}}}*/
     283
     284        /*initialize matrices and vectors:*/
     285        femmodel->parameters->FindParam(&nel,MeshNumberofelementsEnum);
     286        loads=new Vector<IssmDouble>(nel);
     287        oceanareas=new Vector<IssmDouble>(nel);
     288        sealevelloads=new Vector<IssmDouble>(nel);
     289        sealevelloads->Set(0);sealevelloads->Assemble();
     290
    253291        /*call masks core: */
    254292        masks=sealevel_masks(femmodel);
    255 
    256         /*call barystatic core  (generalized eustatic - Farrel and Clark, Eq 4, 1st, 3rd and 4rd terms on the RHS) */
    257         RSLg_barystatic=sealevelchange_core_barystatic(femmodel,masks,&oceanarea);
    258 
    259         /*call self attraction and loading module (ocean loading tems  - 2nd and 5th terms on the RHS of Farrel and Clark) */
    260         RSLg=sealevelchange_core_sal(femmodel,masks,RSLg_barystatic,oceanarea);
    261 
    262         /*compute bedrock motion and derive geoid: */
    263         sealevelchange_core_deformation(&N_grd,&U_grd,&U_north_grd,&U_east_grd,femmodel,RSLg,masks);
     293        if(VerboseSolution()) _printf0_("         starting  GRD convolutions\n");
     294       
     295        /*buildup loads: */
     296        for(Object* & object : femmodel->elements->objects){
     297                Element* element = xDynamicCast<Element*>(object);
     298                element->SealevelchangeBarystaticLoads(loads, barycontrib,masks);
     299        }
     300        loads->Assemble();
     301
     302        //broadcast loads
     303        allloads=loads->ToMPISerial();
     304
     305        //compute rotation axis motion:
     306        RotationAxisMotion(&rotationaxismotionvector[0],femmodel,allloads,NULL);
     307
     308        /*skip computation of sea level if requested, which means sea level loads should be zeroed */
     309        if(!computesealevel){
     310                allsealevelloads=xNewZeroInit<IssmDouble>(nel);
     311                goto deformation;
     312        }
     313
     314        if(VerboseSolution()) _printf0_("         converging GRD convolutions\n");
     315        for(;;){
     316                       
     317                oldsealevelloads=sealevelloads->Duplicate(); sealevelloads->Copy(oldsealevelloads);
     318
     319                /*convolve load and sealevel loads on oceans:*/
     320                for(Object* & object : femmodel->elements->objects){
     321                        Element* element = xDynamicCast<Element*>(object);
     322                        element->SealevelchangeConvolution(sealevelloads, oceanareas, allsealevelloads, allloads,rotationaxismotionvector);
     323                }
     324                sealevelloads->Assemble();
     325
     326                /*compute ocean areas:*/
     327                if(!allsealevelloads){ //first time in the loop
     328                        oceanareas->Assemble(); oceanareas->Sum(&oceanarea); _assert_(oceanarea>0.);
     329                        if(scaleoceanarea) oceanarea=3.619e+14; // use true ocean area, m^2
     330                }
     331               
     332                //Conserve ocean mass:
     333                oceanaverage=SealevelloadsOceanAverage(sealevelloads,oceanareas,oceanarea);
     334                ConserveOceanMass(femmodel,sealevelloads,barycontrib->Total()/oceanarea - oceanaverage,masks);
     335
     336                //broadcast sea level loads
     337                allsealevelloads=sealevelloads->ToMPISerial();
     338
     339                //compute rotation axis motion:
     340                RotationAxisMotion(&rotationaxismotionvector[0],femmodel,allloads,allsealevelloads);
     341
     342                //convergence?
     343                slcconvergence(&converged,sealevelloads,oldsealevelloads,eps_rel,eps_abs);
     344                if (converged)break;
     345
     346                //early return?
     347                if(iterations>=max_nonlinear_iterations)break;
     348                else iterations++;
     349        }
     350
     351        deformation:
     352       
     353        if(VerboseSolution()) _printf0_("         deformation GRD convolutions\n");
     354
     355        /*convolve loads and sea level loads to get the deformation:*/
     356        for(Object* & object : femmodel->elements->objects){
     357                Element* element = xDynamicCast<Element*>(object);
     358                element->SealevelchangeDeformationConvolution(allsealevelloads, allloads, rotationaxismotionvector);
     359        }
     360       
     361        if(VerboseSolution()) _printf0_("         updating GRD fields\n");
    264362
    265363        /*Update bedrock motion and geoid:*/
    266         GetVectorFromInputsx(&geoid,femmodel,SealevelEnum,VertexSIdEnum);
    267         GetVectorFromInputsx(&bedrock,femmodel,BedEnum,VertexSIdEnum);
     364        if(computesealevel){
     365                femmodel->inputs->Shift(SealevelGRDEnum,barycontrib->Total()/rho_water/oceanarea- oceanaverage/rho_water);
     366
     367                /*cumulate barystatic contributions and save to results: */
     368                barycontrib->Cumulate(femmodel->parameters);
     369                barycontrib->Save(femmodel->results,femmodel->parameters,oceanarea);
     370        }
     371
     372        femmodel->inputs->AXPY(1,SealevelGRDEnum,SealevelEnum);
     373        femmodel->inputs->AXPY(1,BedGRDEnum,BedEnum);
    268374        if(horiz){
    269                 GetVectorFromInputsx(&bedrockeast,femmodel,BedEastEnum,VertexSIdEnum);
    270                 GetVectorFromInputsx(&bedrocknorth,femmodel,BedNorthEnum,VertexSIdEnum);
    271         }
    272 
    273         geoid->AXPY(N_grd,1);
    274         bedrock->AXPY(U_grd,1);
    275         if(horiz){
    276                 bedrockeast->AXPY(U_east_grd,1);
    277                 bedrocknorth->AXPY(U_north_grd,1);
    278         }
    279 
    280         /*get some of the updates into elements:*/
    281         InputUpdateFromVectorx(femmodel,U_grd,SealevelUGrdEnum,VertexSIdEnum);
    282         InputUpdateFromVectorx(femmodel,N_grd,SealevelEnum,VertexSIdEnum);
    283         InputUpdateFromVectorx(femmodel,N_grd,SealevelNGrdEnum,VertexSIdEnum);
    284         InputUpdateFromVectorx(femmodel,bedrock,BedEnum,VertexSIdEnum);
    285         if(RSLg)InputUpdateFromVectorx(femmodel,RSLg,SealevelRSLEnum,VertexSIdEnum);
    286         if(RSLg_barystatic)InputUpdateFromVectorx(femmodel,RSLg_barystatic,SealevelRSLBarystaticEnum,VertexSIdEnum);
    287         if (horiz){
    288                 InputUpdateFromVectorx(femmodel,U_north_grd,SealevelUNorthEsaEnum,VertexSIdEnum);       
    289                 InputUpdateFromVectorx(femmodel,U_east_grd,SealevelUEastEsaEnum,VertexSIdEnum);
    290         }
    291 
    292         /*reset counter to 1:*/
    293         femmodel->parameters->SetParam(1,SealevelchangeRunCountEnum); //reset counter.
    294 
    295         /*free ressources:{{{*/
    296         delete RSLg;
    297         delete RSLg_barystatic;
    298         delete U_grd;
    299         delete N_grd;
    300         delete bedrock;
    301         delete geoid;
    302         if(horiz){
    303                 delete U_north_grd;
    304                 delete U_east_grd;
    305                 delete bedrockeast;
    306                 delete bedrocknorth;
    307         }
    308         delete masks;
    309         /*}}}*/
    310 
    311 }
     375                femmodel->inputs->AXPY(1,BedEastGRDEnum,BedEastEnum);
     376                femmodel->inputs->AXPY(1,BedNorthGRDEnum, BedNorthEnum);
     377        }
     378
     379}
    312380/*}}}*/
    313381void dynstr_core(FemModel* femmodel){ /*{{{*/
    … …  
    397465        }
    398466}; /*}}}*/
    399 
    400 void grd_core_optim(FemModel* femmodel) { /*{{{*/
    401 
    402         /*variables:{{{*/
    403         int nel;
    404         BarystaticContributions* barycontrib=NULL;
    405         SealevelMasks* masks=NULL;
    406         GenericParam<BarystaticContributions*>* barycontribparam=NULL;
    407         IssmDouble rotationaxismotionvector[3]={0};
    408        
    409         Vector<IssmDouble>*    loads=NULL;
    410         IssmDouble*            allloads=NULL;
    411         Vector<IssmDouble>*    sealevelloads=NULL;
    412         Vector<IssmDouble>*    oldsealevelloads=NULL;
    413         IssmDouble*            allsealevelloads=NULL;
    414         Vector<IssmDouble>*    oceanareas=NULL;
    415         IssmDouble             oceanarea;
    416         IssmDouble             oceanaverage;
    417         bool                   scaleoceanarea=false;
    418         IssmDouble             rho_water;
    419 
    420         IssmDouble           eps_rel;
    421         IssmDouble           eps_abs;
    422         int                  max_nonlinear_iterations;
    423         int                  iterations=0;
    424         int                  step;
    425         IssmDouble           time;
    426         bool converged=false;
    427 
    428         int  modelid,earthid;
    429         int  horiz;
    430         int  count,frequency,iscoupling;
    431         int  grd=0;
    432         int  computesealevel=0;
    433 
    434         /*}}}*/
    435 
    436         /*Verbose: */
    437         if(VerboseSolution()) _printf0_("         computing GRD patterns\n");
    438 
    439         /*retrieve parameters:{{{*/
    440         femmodel->parameters->FindParam(&grd,SolidearthSettingsGRDEnum);
    441         femmodel->parameters->FindParam(&frequency,SolidearthSettingsRunFrequencyEnum);
    442         femmodel->parameters->FindParam(&count,SealevelchangeRunCountEnum);
    443         femmodel->parameters->FindParam(&computesealevel,SolidearthSettingsComputesealevelchangeEnum);
    444         femmodel->parameters->FindParam(&max_nonlinear_iterations,SolidearthSettingsMaxiterEnum);
    445         /*}}}*/
    446 
    447         /*only run if grd was requested, if we are the earth, and we have reached
    448          * the necessary number of time steps dictated by :*/
    449         if(!grd)            return;
    450         if(count!=frequency)return;
    451         femmodel->parameters->FindParam(&iscoupling,IsSlcCouplingEnum);
    452         if(iscoupling){
    453                 femmodel->parameters->FindParam(&modelid,ModelIdEnum);
    454                 femmodel->parameters->FindParam(&earthid,EarthIdEnum);
    455                 if(modelid!=earthid)return;
    456         }
    457         /*retrieve parameters: {{{*/
    458         femmodel->parameters->FindParam(&scaleoceanarea,SolidearthSettingsOceanAreaScalingEnum);
    459         barycontribparam = xDynamicCast<GenericParam<BarystaticContributions*>*>(femmodel->parameters->FindParamObject(BarystaticContributionsEnum));
    460         barycontrib=barycontribparam->GetParameterValue();
    461         femmodel->parameters->FindParam(&rho_water,MaterialsRhoSeawaterEnum);
    462         femmodel->parameters->FindParam(&eps_rel,SolidearthSettingsReltolEnum);
    463         femmodel->parameters->FindParam(&eps_abs,SolidearthSettingsAbstolEnum);
    464         femmodel->parameters->FindParam(&horiz,SolidearthSettingsHorizEnum);
    465         /*}}}*/
    466 
    467         /*initialize matrices and vectors:*/
    468         femmodel->parameters->FindParam(&nel,MeshNumberofelementsEnum);
    469         loads=new Vector<IssmDouble>(nel);
    470         oceanareas=new Vector<IssmDouble>(nel);
    471         sealevelloads=new Vector<IssmDouble>(nel);
    472         sealevelloads->Set(0);sealevelloads->Assemble();
    473 
    474         /*call masks core: */
    475         masks=sealevel_masks(femmodel);
    476         if(VerboseSolution()) _printf0_("         starting  GRD convolutions\n");
    477        
    478         /*buildup loads: */
     467void ivins_deformation_core(FemModel* femmodel){ /*{{{*/
     468
     469        int  gsize;
     470        Vector<IssmDouble> *bedup  = NULL;
     471        Vector<IssmDouble> *beduprate= NULL;
     472        IssmDouble          *xx     = NULL;
     473        IssmDouble          *yy     = NULL;
     474       
     475        if(VerboseSolution()) _printf0_("         computing vertical deformation using Ivins model. \n");
     476
     477        /*find size of vectors:*/
     478        gsize      = femmodel->nodes->NumberOfDofs(GsetEnum);
     479
     480        /*Find the litho material to be used by all the elements:*/
     481        Matlitho* matlitho=NULL;
     482        for (Object* & object: femmodel->materials->objects){
     483                Material* material=xDynamicCast<Material*>(object);
     484                if(material->ObjectEnum()==MatlithoEnum){
     485                        matlitho=xDynamicCast<Matlitho*>(material);
     486                        break;
     487                }
     488        }
     489
     490        /*initialize vectors:*/
     491        bedup = new Vector<IssmDouble>(gsize);
     492        beduprate = new Vector<IssmDouble>(gsize);
     493       
     494        /*retrieve geometric information: */
     495        VertexCoordinatesx(&xx,&yy,NULL,femmodel->vertices);
     496
     497        /*Go through elements, and add contribution from each element to the deflection vector wg:*/
    479498        for(Object* & object : femmodel->elements->objects){
    480499                Element* element = xDynamicCast<Element*>(object);
    481                 element->SealevelchangeBarystaticLoads(loads, barycontrib,masks);
    482         }
    483         loads->Assemble();
    484 
    485         //broadcast loads
    486         allloads=loads->ToMPISerial();
    487 
    488         //compute rotation axis motion:
    489         RotationAxisMotion(&rotationaxismotionvector[0],femmodel,allloads,NULL);
    490 
    491         /*skip computation of sea level if requested, which means sea level loads should be zeroed */
    492         if(!computesealevel){
    493                 allsealevelloads=xNewZeroInit<IssmDouble>(nel);
    494                 goto deformation;
    495         }
    496 
    497         if(VerboseSolution()) _printf0_("         converging GRD convolutions\n");
    498         for(;;){
    499                        
    500                 oldsealevelloads=sealevelloads->Duplicate(); sealevelloads->Copy(oldsealevelloads);
    501 
    502                 /*convolve load and sealevel loads on oceans:*/
    503                 for(Object* & object : femmodel->elements->objects){
    504                         Element* element = xDynamicCast<Element*>(object);
    505                         element->SealevelchangeConvolution(sealevelloads, oceanareas, allsealevelloads, allloads,rotationaxismotionvector);
    506                 }
    507                 sealevelloads->Assemble();
    508 
    509                 /*compute ocean areas:*/
    510                 if(!allsealevelloads){ //first time in the loop
    511                         oceanareas->Assemble(); oceanareas->Sum(&oceanarea); _assert_(oceanarea>0.);
    512                         if(scaleoceanarea) oceanarea=3.619e+14; // use true ocean area, m^2
    513                 }
    514                
    515                 //Conserve ocean mass:
    516                 oceanaverage=SealevelloadsOceanAverage(sealevelloads,oceanareas,oceanarea);
    517                 ConserveOceanMass(femmodel,sealevelloads,barycontrib->Total()/oceanarea - oceanaverage,masks);
    518 
    519                 //broadcast sea level loads
    520                 allsealevelloads=sealevelloads->ToMPISerial();
    521 
    522                 //compute rotation axis motion:
    523                 RotationAxisMotion(&rotationaxismotionvector[0],femmodel,allloads,allsealevelloads);
    524 
    525                 //convergence?
    526                 slcconvergence(&converged,sealevelloads,oldsealevelloads,eps_rel,eps_abs);
    527                 if (converged)break;
    528 
    529                 //early return?
    530                 if(iterations>=max_nonlinear_iterations)break;
    531                 else iterations++;
    532         }
    533 
    534         deformation:
    535        
    536         if(VerboseSolution()) _printf0_("         deformation GRD convolutions\n");
    537 
    538         /*convolve loads and sea level loads to get the deformation:*/
    539         for(Object* & object : femmodel->elements->objects){
    540                 Element* element = xDynamicCast<Element*>(object);
    541                 element->SealevelchangeDeformationConvolution(allsealevelloads, allloads, rotationaxismotionvector);
    542         }
    543        
    544         if(VerboseSolution()) _printf0_("         updating GRD fields\n");
    545 
    546         /*Update bedrock motion and geoid:*/
    547         if(computesealevel){
    548                 femmodel->inputs->Shift(SealevelGRDEnum,barycontrib->Total()/rho_water/oceanarea- oceanaverage/rho_water);
    549 
    550                 /*cumulate barystatic contributions and save to results: */
    551                 barycontrib->Cumulate(femmodel->parameters);
    552                 barycontrib->Save(femmodel->results,femmodel->parameters,oceanarea);
    553         }
    554 
    555         femmodel->inputs->AXPY(1,SealevelGRDEnum,SealevelEnum);
     500                element->GiaDeflection(bedup,beduprate, matlitho, xx,yy);
     501        }
     502
     503        /*Assemble parallel vector:*/
     504        beduprate->Assemble();
     505        bedup->Assemble();
     506
     507        /*Save results:*/
     508        InputUpdateFromVectorx(femmodel,bedup,BedGRDEnum,VertexSIdEnum);
    556509        femmodel->inputs->AXPY(1,BedGRDEnum,BedEnum);
    557         if(horiz){
    558                 femmodel->inputs->AXPY(1,BedEastGRDEnum,BedEastEnum);
    559                 femmodel->inputs->AXPY(1,BedNorthGRDEnum, BedNorthEnum);
    560         }
    561 
     510
     511        /*Free ressources: */
     512        xDelete<IssmDouble>(xx);
     513        xDelete<IssmDouble>(yy);
     514        delete beduprate;
     515        delete bedup;
    562516}
    563517/*}}}*/
    … …  
    612566        femmodel->parameters->FindParam(&geometrydone,SealevelchangeGeometryDoneEnum);
    613567        femmodel->parameters->FindParam(&grdmodel,GrdModelEnum);
    614         femmodel->parameters->FindParam(&optim,SolidearthSettingsOptimEnum);
    615568       
    616569        /*early return?:*/
    … …  
    635588        for(Object* & object : femmodel->elements->objects){
    636589                Element*   element=xDynamicCast<Element*>(object);
    637                 if(optim) element->SealevelchangeGeometryOptim(latitude,longitude,radius,xx,yy,zz,xxe,yye,zze,areae);
    638                 else      element->SealevelchangeGeometry(latitude,longitude,radius,xx,yy,zz,xxe,yye,zze);
     590                element->SealevelchangeGeometry(latitude,longitude,radius,xx,yy,zz,xxe,yye,zze,areae);
    639591        }
    640592
    … …  
    658610
    659611}/*}}}*/
    660 
    661 //GRD:
    662 Vector<IssmDouble>* sealevelchange_core_barystatic(FemModel* femmodel,SealevelMasks* masks, IssmDouble* poceanarea){ /*{{{*/
    663 
    664         /*Barystatic core of the SLR solution (terms that are constant with respect to sea-level)*/
    665 
    666         Vector<IssmDouble> *RSLgi    = NULL;
    667         IssmDouble          RSLgi_oceanaverage   = 0;
    668 
    669         /*parameters: */
    670         int  gsize;
    671         IssmDouble oceanarea;
    672         int        step;
    673         IssmDouble time;
    674 
    675         /*barystatic contribution:*/
    676         IssmDouble bslc;
    677         IssmDouble bslcice;
    678         IssmDouble* bslcice_partition=NULL;
    679         IssmDouble bslchydro;
    680         IssmDouble* bslchydro_partition=NULL;
    681         IssmDouble cumbslc;
    682         IssmDouble cumbslcice;
    683         IssmDouble cumbslchydro;
    684         IssmDouble* cumbslcice_partition=NULL;
    685         int npartice;
    686         IssmDouble* cumbslchydro_partition=NULL;
    687         int nparthydro;
    688         int computesealevel=0;
    689 
    690         /*early return if we are not computing sea level, but rather deformation: */
    691         femmodel->parameters->FindParam(&computesealevel,SolidearthSettingsComputesealevelchangeEnum);
    692         if (!computesealevel)return NULL;
    693        
    694         if(VerboseSolution()) _printf0_("         computing bslc components on ice\n");
    695 
    696         /*Figure out size of g-set deflection vector and allocate solution vector: */
    697         gsize = femmodel->nodes->NumberOfDofs(GsetEnum);
    698 
    699         /*some parameters:*/
    700         femmodel->parameters->FindParam(&step,StepEnum);
    701         femmodel->parameters->FindParam(&time,TimeEnum);
    702         femmodel->parameters->FindParam(&cumbslc,CumBslcEnum);
    703         femmodel->parameters->FindParam(&cumbslcice,CumBslcIceEnum);
    704         femmodel->parameters->FindParam(&cumbslchydro,CumBslcHydroEnum);
    705         femmodel->parameters->FindParam(&npartice,SolidearthNpartIceEnum);
    706         femmodel->parameters->FindParam(&nparthydro,SolidearthNpartHydroEnum);
    707         if(npartice) femmodel->parameters->FindParam(&cumbslcice_partition,&npartice,NULL,CumBslcIcePartitionEnum);
    708         if(nparthydro) femmodel->parameters->FindParam(&cumbslchydro_partition,&nparthydro,NULL,CumBslcHydroPartitionEnum);
    709 
    710         /*Initialize:*/
    711         RSLgi = new Vector<IssmDouble>(gsize);
    712 
    713         /*call the bslc main module: */
    714         femmodel->SealevelchangeBarystatic(RSLgi,&oceanarea,&bslc, &bslcice, &bslchydro, &bslcice_partition, &bslchydro_partition,masks); //this computes
    715 
    716         /*we need to average RSLgi over the ocean: RHS term  4 in Eq.4 of Farrel and clarke. Only the elements can do that: */
    717         RSLgi_oceanaverage=femmodel->SealevelchangeOceanAverage(RSLgi,masks, oceanarea);
    718 
    719         /*RSLg is the sum of the pure bslc component (term 3) and the contribution from the perturbation to the graviation potential due to the
    720          * presence of ice (terms 1 and 4 in Eq.4 of Farrel and Clarke):*/
    721         RSLgi->Shift(-bslc-RSLgi_oceanaverage);
    722 
    723         /*save bslc and cumulated bslc value for results:{{{ */
    724         femmodel->results->AddResult(new GenericExternalResult<IssmDouble>(femmodel->results->Size()+1,BslcEnum,-bslc,step,time));
    725         femmodel->results->AddResult(new GenericExternalResult<IssmDouble>(femmodel->results->Size()+1,BslcIceEnum,-bslcice,step,time));
    726         femmodel->results->AddResult(new GenericExternalResult<IssmDouble>(femmodel->results->Size()+1,BslcHydroEnum,-bslchydro,step,time));
    727 
    728         //cumulative barystatic contribution:
    729         cumbslc=cumbslc-bslc;
    730         cumbslcice=cumbslcice-bslcice;
    731         cumbslchydro=cumbslchydro-bslchydro;
    732 
    733         femmodel->parameters->SetParam(cumbslc,CumBslcEnum);
    734         femmodel->parameters->SetParam(cumbslcice,CumBslcIceEnum);
    735         femmodel->parameters->SetParam(cumbslchydro,CumBslcHydroEnum);
    736 
    737         femmodel->results->AddResult(new GenericExternalResult<IssmDouble>(femmodel->results->Size()+1,CumBslcEnum,cumbslc,step,time));
    738         femmodel->results->AddResult(new GenericExternalResult<IssmDouble>(femmodel->results->Size()+1,CumBslcIceEnum,cumbslcice,step,time));
    739         femmodel->results->AddResult(new GenericExternalResult<IssmDouble>(femmodel->results->Size()+1,CumBslcHydroEnum,cumbslchydro,step,time));
    740 
    741         //cumulative barystatic contributions by partition:
    742         if(npartice){
    743                 for(int i=0;i<npartice;i++) cumbslcice_partition[i] -= bslcice_partition[i];
    744                 femmodel->parameters->SetParam(cumbslcice_partition,npartice,1,CumBslcIcePartitionEnum);
    745                 femmodel->results->AddResult(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,CumBslcIcePartitionEnum,cumbslcice_partition,npartice,1,step,time));
    746         }
    747 
    748         if(nparthydro){
    749                 for(int i=0;i<nparthydro;i++) cumbslchydro_partition[i] -= bslchydro_partition[i];
    750                 femmodel->parameters->SetParam(cumbslchydro_partition,nparthydro,1,CumBslcHydroPartitionEnum);
    751                 femmodel->results->AddResult(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,CumBslcHydroPartitionEnum,cumbslchydro_partition,nparthydro,1,step,time));
    752         }
    753         /*}}}*/
    754        
    755         /*Assign output pointers and return: */
    756         *poceanarea=oceanarea;
    757         return RSLgi;
    758 }/*}}}*/
    759 Vector<IssmDouble>* sealevelchange_core_sal(FemModel* femmodel, SealevelMasks* masks, Vector<IssmDouble>* RSLg_barystatic,IssmDouble oceanarea){ /*{{{*/
    760 
    761         /*this core computes the contributions from Eq.4 of Farrel and Clarke, rhs terms 2 and 5.
    762           sal core of the SLR solution */
    763 
    764         Vector<IssmDouble> *RSLg    = NULL;
    765         Vector<IssmDouble> *RSLg_old    = NULL;
    766 
    767         Vector<IssmDouble> *RSLgo    = NULL; //ocean convolution of the perturbation to gravity potential.
    768         Vector<IssmDouble> *RSLgo_rot= NULL; // rotational feedback
    769         IssmDouble          RSLgo_oceanaverage = 0;  //average of RSLgo over the ocean.
    770 
    771         /*parameters: */
    772         int count;
    773         bool save_results;
    774         int  gsize;
    775         bool converged=true;
    776         bool rotation=true;
    777         bool verboseconvolution=true;
    778         int max_nonlinear_iterations;
    779         IssmDouble           eps_rel;
    780         IssmDouble           eps_abs;
    781         IssmDouble                      Ixz, Iyz, Izz;
    782         int computesealevel=0;
    783 
    784         /*early return if we are not computing sea level, but rather deformation: */
    785         femmodel->parameters->FindParam(&computesealevel,SolidearthSettingsComputesealevelchangeEnum);
    786         if (!computesealevel)return NULL;
    787        
    788         if(VerboseSolution()) _printf0_("         converging on ocean components\n");
    789 
    790         /*Recover some parameters: */
    791         femmodel->parameters->FindParam(&max_nonlinear_iterations,SolidearthSettingsMaxiterEnum);
    792         femmodel->parameters->FindParam(&eps_rel,SolidearthSettingsReltolEnum);
    793         femmodel->parameters->FindParam(&eps_abs,SolidearthSettingsAbstolEnum);
    794 
    795         /*computational flag: */
    796         femmodel->parameters->FindParam(&rotation,SolidearthSettingsRotationEnum);
    797 
    798         /*Figure out size of g-set deflection vector and allocate solution vector: */
    799         gsize = femmodel->nodes->NumberOfDofs(GsetEnum);
    800 
    801         /*Initialize:*/
    802         RSLg = new Vector<IssmDouble>(gsize);
    803         RSLg->Assemble();
    804         RSLg_barystatic->Copy(RSLg);  //first initialize RSLg with the barystatic component computed in sealevelchange_core_barystatic.
    805 
    806         RSLg_old = new Vector<IssmDouble>(gsize);
    807         RSLg_old->Assemble();
    808 
    809         count=1;
    810         converged=false;
    811 
    812         /*Start loop: */
    813         for(;;){
    814 
    815                 //save pointer to old sea level
    816                 delete RSLg_old; RSLg_old=RSLg;
    817 
    818                 /*Initialize solution vector: */
    819                 RSLg  = new Vector<IssmDouble>(gsize); RSLg->Assemble();
    820                 RSLgo = new Vector<IssmDouble>(gsize); RSLgo->Assemble();
    821 
    822                 /*call the sal module: */
    823                 femmodel->SealevelchangeSal(RSLgo, RSLg_old,  masks, verboseconvolution);
    824 
    825                 /*assemble solution vector: */
    826                 RSLgo->Assemble();
    827 
    828                 if(rotation){
    829 
    830                         /*call rotational feedback  module: */
    831                         RSLgo_rot = new Vector<IssmDouble>(gsize); RSLgo_rot->Assemble();
    832                         femmodel->SealevelchangeRotationalFeedback(RSLgo_rot,RSLg_old,&Ixz,&Iyz,&Izz, masks);
    833                         RSLgo_rot->Assemble();
    834 
    835                         /*save changes in inertia tensor as results: */
    836                         femmodel->results->AddResult(new GenericExternalResult<IssmDouble>(femmodel->results->Size()+1,SealevelInertiaTensorXZEnum,Ixz));
    837                         femmodel->results->AddResult(new GenericExternalResult<IssmDouble>(femmodel->results->Size()+1,SealevelInertiaTensorYZEnum,Iyz));
    838                         femmodel->results->AddResult(new GenericExternalResult<IssmDouble>(femmodel->results->Size()+1,SealevelInertiaTensorZZEnum,Izz));
    839 
    840                         RSLgo->AXPY(RSLgo_rot,1);
    841                 }
    842 
    843                 /*we need to average RSLgo over the ocean: RHS term  5 in Eq.4 of Farrel and clarke. Only the elements can do that: */
    844                 RSLgo_oceanaverage=femmodel->SealevelchangeOceanAverage(RSLgo,masks, oceanarea);
    845 
    846                 /*RSLg is the sum of the barystatic term, and the ocean terms: */
    847                 RSLg_barystatic->Copy(RSLg); RSLg->AXPY(RSLgo,1);
    848                 RSLg->Shift(-RSLgo_oceanaverage);
    849 
    850                 /*convergence criterion:*/
    851                 slcconvergence(&converged,RSLg,RSLg_old,eps_rel,eps_abs);
    852                
    853        
    854                 /*free ressources: */
    855                 delete RSLgo;
    856                 delete RSLgo_rot;
    857 
    858                 /*Increase count: */
    859                 count++;
    860                 if(converged==true){
    861                         break;
    862                 }
    863                 if(count>=max_nonlinear_iterations){
    864                         _printf0_("   maximum number of nonlinear iterations (" << max_nonlinear_iterations << ") exceeded\n");
    865                         converged=true;
    866                         break;
    867                 }       
    868 
    869                 /*some minor verbosing adjustment:*/
    870                 if(count>1)verboseconvolution=false;
    871 
    872         }
    873         if(VerboseConvergence()) _printf0_("\n              total number of iterations: " << count-1 << "\n");
    874        
    875        
    876         delete RSLg_old;
    877 
    878         return RSLg;
    879 } /*}}}*/
    880 void sealevelchange_core_deformation(Vector<IssmDouble>** pN_grd, Vector<IssmDouble>** pU_grd, Vector<IssmDouble>** pU_north_grd,Vector<IssmDouble>** pU_east_grd,FemModel* femmodel,Vector<IssmDouble>* RSLg, SealevelMasks* masks){ /*{{{*/
    881 
    882         Vector<IssmDouble> *U_grd  = NULL;
    883         Vector<IssmDouble> *dUdt_grd  = NULL;
    884         Vector<IssmDouble> *N_grd  = NULL;
    885         Vector<IssmDouble> *U_north_grd   = NULL;
    886         Vector<IssmDouble> *U_east_grd    = NULL;
    887 
    888         /*parameters: */
    889         int  gsize;
    890         bool spherical=true;
    891 
    892         IssmDouble          *latitude   = NULL;
    893         IssmDouble          *longitude  = NULL;
    894         IssmDouble          *radius     = NULL;
    895         IssmDouble          *xx     = NULL;
    896         IssmDouble          *yy     = NULL;
    897         IssmDouble          *zz     = NULL;
    898         int  horiz;
    899         int  grdmodel;
    900        
    901         if(VerboseSolution()) _printf0_("         computing vertical and horizontal geodetic signatures\n");
    902 
    903         /*retrieve some parameters:*/
    904         femmodel->parameters->FindParam(&horiz,SolidearthSettingsHorizEnum);
    905         femmodel->parameters->FindParam(&grdmodel,GrdModelEnum);
    906 
    907         /*find size of vectors:*/
    908         gsize      = femmodel->nodes->NumberOfDofs(GsetEnum);
    909 
    910         /*intialize vectors:*/
    911         U_grd = new Vector<IssmDouble>(gsize);
    912         if(grdmodel==IvinsEnum) dUdt_grd = new Vector<IssmDouble>(gsize);
    913         N_grd = new Vector<IssmDouble>(gsize);
    914         if (horiz){
    915                 U_north_grd = new Vector<IssmDouble>(gsize);
    916                 U_east_grd = new Vector<IssmDouble>(gsize);
    917         }
    918        
    919         /*retrieve geometric information: */
    920         VertexCoordinatesx(&latitude,&longitude,&radius,femmodel->vertices,spherical);
    921         VertexCoordinatesx(&xx,&yy,&zz,femmodel->vertices);
    922 
    923         /*call the deformation module: */
    924         switch(grdmodel){
    925                 case NoneEnum:
    926                         //do nothing:
    927                         break;
    928                 case IvinsEnum:
    929                         femmodel->IvinsDeformation(U_grd,dUdt_grd,xx,yy);
    930                         break;
    931                 case ElasticEnum:
    932                         femmodel->SealevelchangeDeformation(N_grd, U_grd,U_north_grd,U_east_grd,RSLg, masks);
    933                         break;
    934                 default:
    935                         _error_("Grd model " << EnumToStringx(grdmodel) << " not supported yet");
    936         }
    937 
    938         /*Assign output pointers:*/
    939         *pU_grd=U_grd;
    940         *pN_grd=N_grd;
    941         if(horiz){
    942                 *pU_east_grd=U_east_grd;
    943                 *pU_north_grd=U_north_grd;
    944         }
    945 
    946         /*Free ressources: */
    947         xDelete<IssmDouble>(longitude);
    948         xDelete<IssmDouble>(latitude);
    949         xDelete<IssmDouble>(xx);
    950         xDelete<IssmDouble>(yy);
    951         xDelete<IssmDouble>(zz);
    952         xDelete<IssmDouble>(radius);
    953         if(grdmodel==IvinsEnum)delete dUdt_grd;
    954 }
    955 /*}}}*/
    956612
    957613/*Support routines:*/
    … …  
    1284940        for(Object* & object : femmodel->elements->objects){
    1285941                Element* element = xDynamicCast<Element*>(object);
    1286                 element->SealevelchangeMomentOfInertiaElement(&moi_list[0],loads,sealevelloads);
     942                element->SealevelchangeMomentOfInertia(&moi_list[0],loads,sealevelloads);
    1287943                moi_list_cpu[0] += moi_list[0];
    1288944                moi_list_cpu[1] += moi_list[1];

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

    r26155 r26165  
    350350syn keyword cConstant SolidearthSettingsAbstolEnum
    351351syn keyword cConstant SolidearthSettingsCrossSectionShapeEnum
    352 syn keyword cConstant SolidearthSettingsOptimEnum
    353352syn keyword cConstant RotationalAngularVelocityEnum
    354353syn keyword cConstant SolidearthSettingsElasticEnum
    … …  
    14521451syn keyword cType Cfsurfacesquare
    14531452syn keyword cType Channel
    1454 syn keyword cType classes
    14551453syn keyword cType Constraint
    14561454syn keyword cType Constraints
    … …  
    14591457syn keyword cType ControlInput
    14601458syn keyword cType Covertree
     1459syn keyword cType DataSetParam
    14611460syn keyword cType DatasetInput
    1462 syn keyword cType DataSetParam
    14631461syn keyword cType Definition
    14641462syn keyword cType DependentObject
    … …  
    14731471syn keyword cType ElementInput
    14741472syn keyword cType ElementMatrix
     1473syn keyword cType ElementVector
    14751474syn keyword cType Elements
    1476 syn keyword cType ElementVector
    14771475syn keyword cType ExponentialVariogram
    14781476syn keyword cType ExternalResult
    … …  
    14811479syn keyword cType Friction
    14821480syn keyword cType Gauss
    1483 syn keyword cType GaussianVariogram
    1484 syn keyword cType gaussobjects
    14851481syn keyword cType GaussPenta
    14861482syn keyword cType GaussSeg
    14871483syn keyword cType GaussTetra
    14881484syn keyword cType GaussTria
     1485syn keyword cType GaussianVariogram
    14891486syn keyword cType GenericExternalResult
    14901487syn keyword cType GenericOption
    … …  
    15011498syn keyword cType IssmDirectApplicInterface
    15021499syn keyword cType IssmParallelDirectApplicInterface
    1503 syn keyword cType krigingobjects
    15041500syn keyword cType Load
    15051501syn keyword cType Loads
    … …  
    15121508syn keyword cType Matice
    15131509syn keyword cType Matlitho
    1514 syn keyword cType matrixobjects
    15151510syn keyword cType MatrixParam
    15161511syn keyword cType Misfit
    … …  
    15251520syn keyword cType Observations
    15261521syn keyword cType Option
     1522syn keyword cType OptionUtilities
    15271523syn keyword cType Options
    1528 syn keyword cType OptionUtilities
    15291524syn keyword cType Param
    15301525syn keyword cType Parameters
    … …  
    15401535syn keyword cType Regionaloutput
    15411536syn keyword cType Results
     1537syn keyword cType RiftStruct
    15421538syn keyword cType Riftfront
    1543 syn keyword cType RiftStruct
    15441539syn keyword cType SealevelMasks
    15451540syn keyword cType Seg
    15461541syn keyword cType SegInput
     1542syn keyword cType SegRef
    15471543syn keyword cType Segment
    1548 syn keyword cType SegRef
    15491544syn keyword cType SpcDynamic
    15501545syn keyword cType SpcStatic
    … …  
    15651560syn keyword cType Vertex
    15661561syn keyword cType Vertices
     1562syn keyword cType classes
     1563syn keyword cType gaussobjects
     1564syn keyword cType krigingobjects
     1565syn keyword cType matrixobjects
    15671566syn keyword cType AdjointBalancethickness2Analysis
    15681567syn keyword cType AdjointBalancethicknessAnalysis

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

    r26155 r26165  
    344344        SolidearthSettingsAbstolEnum,
    345345        SolidearthSettingsCrossSectionShapeEnum,
    346         SolidearthSettingsOptimEnum,
    347346        RotationalAngularVelocityEnum,
    348347        SolidearthSettingsElasticEnum,

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

    r26155 r26165  
    352352                case SolidearthSettingsAbstolEnum : return "SolidearthSettingsAbstol";
    353353                case SolidearthSettingsCrossSectionShapeEnum : return "SolidearthSettingsCrossSectionShape";
    354                 case SolidearthSettingsOptimEnum : return "SolidearthSettingsOptim";
    355354                case RotationalAngularVelocityEnum : return "RotationalAngularVelocity";
    356355                case SolidearthSettingsElasticEnum : return "SolidearthSettingsElastic";

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

    r26155 r26165  
    358358              else if (strcmp(name,"SolidearthSettingsAbstol")==0) return SolidearthSettingsAbstolEnum;
    359359              else if (strcmp(name,"SolidearthSettingsCrossSectionShape")==0) return SolidearthSettingsCrossSectionShapeEnum;
    360               else if (strcmp(name,"SolidearthSettingsOptim")==0) return SolidearthSettingsOptimEnum;
    361360              else if (strcmp(name,"RotationalAngularVelocity")==0) return RotationalAngularVelocityEnum;
    362361              else if (strcmp(name,"SolidearthSettingsElastic")==0) return SolidearthSettingsElasticEnum;
    … …  
    383382              else if (strcmp(name,"SolidearthSettingsReltol")==0) return SolidearthSettingsReltolEnum;
    384383              else if (strcmp(name,"SealevelchangeRequestedOutputs")==0) return SealevelchangeRequestedOutputsEnum;
     384              else if (strcmp(name,"SolidearthSettingsRigid")==0) return SolidearthSettingsRigidEnum;
    385385         else stage=4;
    386386   }
    387387   if(stage==4){
    388               if (strcmp(name,"SolidearthSettingsRigid")==0) return SolidearthSettingsRigidEnum;
    389               else if (strcmp(name,"SolidearthSettingsRotation")==0) return SolidearthSettingsRotationEnum;
     388              if (strcmp(name,"SolidearthSettingsRotation")==0) return SolidearthSettingsRotationEnum;
    390389              else if (strcmp(name,"SealevelchangeRunCount")==0) return SealevelchangeRunCountEnum;
    391390              else if (strcmp(name,"SealevelchangeTransitions")==0) return SealevelchangeTransitionsEnum;
    … …  
    506505              else if (strcmp(name,"TransientIsgroundingline")==0) return TransientIsgroundinglineEnum;
    507506              else if (strcmp(name,"TransientIshydrology")==0) return TransientIshydrologyEnum;
     507              else if (strcmp(name,"TransientIsmasstransport")==0) return TransientIsmasstransportEnum;
    508508         else stage=5;
    509509   }
    510510   if(stage==5){
    511               if (strcmp(name,"TransientIsmasstransport")==0) return TransientIsmasstransportEnum;
    512               else if (strcmp(name,"TransientIsoceantransport")==0) return TransientIsoceantransportEnum;
     511              if (strcmp(name,"TransientIsoceantransport")==0) return TransientIsoceantransportEnum;
    513512              else if (strcmp(name,"TransientIsmovingfront")==0) return TransientIsmovingfrontEnum;
    514513              else if (strcmp(name,"TransientIsoceancoupling")==0) return TransientIsoceancouplingEnum;
    … …  
    629628              else if (strcmp(name,"EplHead")==0) return EplHeadEnum;
    630629              else if (strcmp(name,"EplHeadOld")==0) return EplHeadOldEnum;
     630              else if (strcmp(name,"EplHeadSlopeX")==0) return EplHeadSlopeXEnum;
    631631         else stage=6;
    632632   }
    633633   if(stage==6){
    634               if (strcmp(name,"EplHeadSlopeX")==0) return EplHeadSlopeXEnum;
    635               else if (strcmp(name,"EplHeadSlopeY")==0) return EplHeadSlopeYEnum;
     634              if (strcmp(name,"EplHeadSlopeY")==0) return EplHeadSlopeYEnum;
    636635              else if (strcmp(name,"EplHeadSubstep")==0) return EplHeadSubstepEnum;
    637636              else if (strcmp(name,"EplHeadTransient")==0) return EplHeadTransientEnum;
    … …  
    752751              else if (strcmp(name,"RadarIcePeriod")==0) return RadarIcePeriodEnum;
    753752              else if (strcmp(name,"RadarPowerMacGregor")==0) return RadarPowerMacGregorEnum;
     753              else if (strcmp(name,"RadarPowerWolff")==0) return RadarPowerWolffEnum;
    754754         else stage=7;
    755755   }
    756756   if(stage==7){
    757               if (strcmp(name,"RadarPowerWolff")==0) return RadarPowerWolffEnum;
    758               else if (strcmp(name,"RheologyBAbsGradient")==0) return RheologyBAbsGradientEnum;
     757              if (strcmp(name,"RheologyBAbsGradient")==0) return RheologyBAbsGradientEnum;
    759758              else if (strcmp(name,"RheologyBInitialguess")==0) return RheologyBInitialguessEnum;
    760759              else if (strcmp(name,"RheologyBInitialguessMisfit")==0) return RheologyBInitialguessMisfitEnum;
    … …  
    875874              else if (strcmp(name,"SmbReini")==0) return SmbReiniEnum;
    876875              else if (strcmp(name,"SmbRunoff")==0) return SmbRunoffEnum;
     876              else if (strcmp(name,"SmbRunoffSubstep")==0) return SmbRunoffSubstepEnum;
    877877         else stage=8;
    878878   }
    879879   if(stage==8){
    880               if (strcmp(name,"SmbRunoffSubstep")==0) return SmbRunoffSubstepEnum;
    881               else if (strcmp(name,"SmbRunoffTransient")==0) return SmbRunoffTransientEnum;
     880              if (strcmp(name,"SmbRunoffTransient")==0) return SmbRunoffTransientEnum;
    882881              else if (strcmp(name,"SmbS0gcm")==0) return SmbS0gcmEnum;
    883882              else if (strcmp(name,"SmbS0p")==0) return SmbS0pEnum;
    … …  
    998997              else if (strcmp(name,"Outputdefinition22")==0) return Outputdefinition22Enum;
    999998              else if (strcmp(name,"Outputdefinition23")==0) return Outputdefinition23Enum;
     999              else if (strcmp(name,"Outputdefinition24")==0) return Outputdefinition24Enum;
    10001000         else stage=9;
    10011001   }
    10021002   if(stage==9){
    1003               if (strcmp(name,"Outputdefinition24")==0) return Outputdefinition24Enum;
    1004               else if (strcmp(name,"Outputdefinition25")==0) return Outputdefinition25Enum;
     1003              if (strcmp(name,"Outputdefinition25")==0) return Outputdefinition25Enum;
    10051004              else if (strcmp(name,"Outputdefinition26")==0) return Outputdefinition26Enum;
    10061005              else if (strcmp(name,"Outputdefinition27")==0) return Outputdefinition27Enum;
    … …  
    11211120              else if (strcmp(name,"CalvingHab")==0) return CalvingHabEnum;
    11221121              else if (strcmp(name,"CalvingLevermann")==0) return CalvingLevermannEnum;
     1122              else if (strcmp(name,"CalvingVonmises")==0) return CalvingVonmisesEnum;
    11231123         else stage=10;
    11241124   }
    11251125   if(stage==10){
    1126               if (strcmp(name,"CalvingVonmises")==0) return CalvingVonmisesEnum;
    1127               else if (strcmp(name,"Cfdragcoeffabsgrad")==0) return CfdragcoeffabsgradEnum;
     1126              if (strcmp(name,"Cfdragcoeffabsgrad")==0) return CfdragcoeffabsgradEnum;
    11281127              else if (strcmp(name,"Cfsurfacelogvel")==0) return CfsurfacelogvelEnum;
    11291128              else if (strcmp(name,"Cfsurfacesquare")==0) return CfsurfacesquareEnum;
    … …  
    12441243              else if (strcmp(name,"IntParam")==0) return IntParamEnum;
    12451244              else if (strcmp(name,"IntVecParam")==0) return IntVecParamEnum;
     1245              else if (strcmp(name,"Inputs")==0) return InputsEnum;
    12461246         else stage=11;
    12471247   }
    12481248   if(stage==11){
    1249               if (strcmp(name,"Inputs")==0) return InputsEnum;
    1250               else if (strcmp(name,"Internal")==0) return InternalEnum;
     1249              if (strcmp(name,"Internal")==0) return InternalEnum;
    12511250              else if (strcmp(name,"Intersect")==0) return IntersectEnum;
    12521251              else if (strcmp(name,"InversionVzObs")==0) return InversionVzObsEnum;
    … …  
    13671366              else if (strcmp(name,"SMBd18opdd")==0) return SMBd18opddEnum;
    13681367              else if (strcmp(name,"SMBforcing")==0) return SMBforcingEnum;
     1368              else if (strcmp(name,"SMBgcm")==0) return SMBgcmEnum;
    13691369         else stage=12;
    13701370   }
    13711371   if(stage==12){
    1372               if (strcmp(name,"SMBgcm")==0) return SMBgcmEnum;
    1373               else if (strcmp(name,"SMBgemb")==0) return SMBgembEnum;
     1372              if (strcmp(name,"SMBgemb")==0) return SMBgembEnum;
    13741373              else if (strcmp(name,"SMBgradients")==0) return SMBgradientsEnum;
    13751374              else if (strcmp(name,"SMBgradientscomponents")==0) return SMBgradientscomponentsEnum;

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

    r26126 r26165  
    2222                grdmodel               = 0; %grd model (0 by default, 1 for elastic, 2 for Ivins)
    2323                cross_section_shape    = 0; %cross section only used when grd model is Ivins
    24                 optim                  = 0; %new optimized version of the GRD code.
    2524        end
    2625        methods
    … …  
    6665                self.cross_section_shape=1; %square as default (see iedge in GiaDeflectionCorex)
    6766
    68                 %optim?
    69                 self.optim=0;
    70 
    7167                %no grd model by default:
    7268                self.grdmodel=0;
    … …  
    8783                        md = checkfield(md,'fieldname','solidearth.settings.grdmodel','values',[1 2]);
    8884                        md = checkfield(md,'fieldname','solidearth.settings.cross_section_shape','numel',[1],'values',[1,2]);
    89                         md = checkfield(md,'fieldname','solidearth.settings.optim','values',[0,1]);
    9085
    9186                        %checks on computational flags
    … …  
    131126                        fielddisplay(self,'grdmodel','type of deformation model, 1 for elastic, 2 for visco-elastic from Ivins');
    132127                        fielddisplay(self,'cross_section_shape','1: square-edged (default). 2: elliptical. See iedge in GiaDeflectionCore');
    133                         fielddisplay(self,'optim','use optimized version of the GRD code? (default 0)');
    134128                end % }}}
    135129                function marshall(self,prefix,md,fid) % {{{
    … …  
    150144                        WriteData(fid,prefix,'object',self,'fieldname','grdmodel','name','md.solidearth.settings.grdmodel','format','Integer');
    151145                        WriteData(fid,prefix,'object',self,'fieldname','cross_section_shape','name','md.solidearth.settings.cross_section_shape','format','Integer');
    152                         WriteData(fid,prefix,'object',self,'fieldname','optim','name','md.solidearth.settings.optim','format','Integer');
    153146                end % }}}
    154147                function savemodeljs(self,fid,modelname) % {{{
    … …  
    165158                        writejsdouble(fid,[modelname '.solidearth.settings.glfraction'],self.glfraction);
    166159                        writejsdouble(fid,[modelname '.solidearth.settings.cross_section_shape'],self.cross_section_shape);
    167                         writejsdouble(fid,[modelname '.solidearth.settings.optim'],self.optim);
    168160                end % }}}
    169161                function self = extrude(self,md) % {{{

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

    r26087 r26165  
    1818md.solidearth.settings.cross_section_shape=1;    % for square-edged x-section
    1919md.solidearth.settings.computesealevelchange=0;  %do not compute sea level, only deformation
    20 md.solidearth.requested_outputs={'Sealevel','SealevelUGrd'};
     20md.solidearth.requested_outputs={'Sealevel','BedGRD'};
    2121
    2222%Loading history
    … …  
    5050field_names     ={'UGrd'};
    5151field_tolerances={1e-13};
    52 field_values={md.results.TransientSolution(2).SealevelUGrd};
     52field_values={md.results.TransientSolution(2).BedGRD};

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

    r26099 r26165  
    33%mesh earth:
    44md=model;
    5 md.mesh=gmshplanet('radius',6.371012*10^3,'resolution',700.); %700 km resolution mesh
     5md.mesh=gmshplanet('radius',6.371012*10^3,'resolution',500.); %700 km resolution mesh
    66
    77%Geometry for the bed, arbitrary thickness of 1000:
    … …  
    2424pos=find(late < -80);
    2525md.masstransport.spcthickness(md.mesh.elements(pos,:))= md.masstransport.spcthickness(md.mesh.elements(pos,:))-100;
     26posant=pos;
    2627%greenland
    2728pos=find(late>70 & late<80 & longe>-60 & longe<-30);
    2829md.masstransport.spcthickness(md.mesh.elements(pos,:))= md.masstransport.spcthickness(md.mesh.elements(pos,:))-100;
     30posgre=pos;
    2931
    3032%elastic loading from love numbers:
    … …  
    3436%mask:  {{{
    3537mask=gmtmask(md.mesh.lat,md.mesh.long);
     38oceanmask=-ones(md.mesh.numberofvertices,1);
     39pos=find(mask==0); oceanmask(pos)=1;
     40
    3641icemask=ones(md.mesh.numberofvertices,1);
    37 pos=find(mask==0);
    38 icemask(pos)=-1;
    39 pos=find(sum(mask(md.mesh.elements),2)<3);
    40 icemask(md.mesh.elements(pos,:))=-1;
     42icemask(md.mesh.elements(posant))=-1;
     43icemask(md.mesh.elements(posgre))=-1;
     44
    4145md.mask.ice_levelset=icemask;
    42 md.mask.ocean_levelset=-icemask;
    43 
    44 
    45 %make sure wherever there is an ice load, that the mask is set to ice:
    46 pos=find(md.masstransport.spcthickness(1:end-1));
    47 md.mask.ice_levelset(md.mesh.elements(pos,:))=-1;
     46md.mask.ocean_levelset=oceanmask;
    4847% }}}
    4948
    … …  
    6665md.materials=materials('hydro');
    6766
    68 
    6967%Miscellaneous
    7068md.miscellaneous.name='test2002';
    7169
    7270%Solution parameters
     71md.cluster.np=3;
    7372md.solidearth.settings.reltol=NaN;
    7473md.solidearth.settings.abstol=1e-3;
    … …  
    8382md.transient.isthermal=0;
    8483md.transient.ismasstransport=1;
    85 md.transient.isoceantransport=1;
    8684md.transient.isslc=1;
    8785md.solidearth.requested_outputs={'Sealevel'};
Note: See TracChangeset for help on using the changeset viewer.