Index: /issm/trunk-jpl/src/c/analyses/SealevelchangeAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/SealevelchangeAnalysis.cpp (revision 26164) +++ /issm/trunk-jpl/src/c/analyses/SealevelchangeAnalysis.cpp (revision 26165) @@ -145,5 +145,4 @@ parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.glfraction",SolidearthSettingsGlfractionEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.cross_section_shape",SolidearthSettingsCrossSectionShapeEnum)); - parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.optim",SolidearthSettingsOptimEnum)); parameters->AddObject(new DoubleParam(CumBslcEnum,0.0)); parameters->AddObject(new DoubleParam(CumBslcIceEnum,0.0)); Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 26164) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 26165) @@ -383,21 +383,12 @@ virtual IssmDouble GetArea3D(void)=0; virtual IssmDouble GetAreaSpherical(void)=0; - virtual IssmDouble OceanAverage(IssmDouble* Sg, SealevelMasks* masks)=0; - virtual void SealevelchangeMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks)=0; - virtual void SealevelchangeMomentOfInertiaElement(IssmDouble* dI_list, IssmDouble* loads, IssmDouble* sealevelloads)=0; - virtual IssmDouble SealevelchangeBarystaticIce(IssmDouble* Sgi, SealevelMasks* masks,Vector* barystatic_contribution,IssmDouble* partition,IssmDouble oceanarea)=0; - virtual IssmDouble SealevelchangeBarystaticHydro(IssmDouble* Sgi, SealevelMasks* masks,Vector* barystatic_contribution, IssmDouble* partition, IssmDouble oceanarea)=0; - virtual void SealevelchangeBarystaticBottomPressure(IssmDouble* Sgi, SealevelMasks* masks)=0; - virtual void SealevelchangeGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, - IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze)=0; - virtual void SealevelchangeSal(IssmDouble* Sgo, IssmDouble* Sg_old,SealevelMasks* mask)=0; - virtual void DeformationFromSurfaceLoads(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg,SealevelMasks* masks)=0; - virtual void GiaDeflection(Vector* wg,Vector* dwgdt,Matlitho* litho, IssmDouble* x,IssmDouble* y)=0; - - virtual void SealevelchangeGeometryOptim(IssmDouble* lat,IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae)=0; - virtual void SealevelchangeBarystaticLoads(Vector* loads, BarystaticContributions* barycontrib, SealevelMasks* masks)=0; - virtual void SealevelchangeConvolution(Vector* sealevelloads, Vector* oceanareas, IssmDouble* allsealevelloads, IssmDouble* allloads,IssmDouble* rotationaxismotionvector)=0; - virtual void SealevelchangeDeformationConvolution(IssmDouble* sealevelloads, IssmDouble* loads, IssmDouble* rotationvector)=0; - virtual void SealevelchangeShift(Vector* loads, IssmDouble offset, SealevelMasks* masks)=0; + virtual void SealevelchangeMomentOfInertia(IssmDouble* dI_list, IssmDouble* loads, IssmDouble* sealevelloads)=0; + virtual void GiaDeflection(Vector* wg,Vector* dwgdt,Matlitho* litho, IssmDouble* x,IssmDouble* y)=0; + + virtual void SealevelchangeGeometry(IssmDouble* lat,IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae)=0; + virtual void SealevelchangeBarystaticLoads(Vector* loads, BarystaticContributions* barycontrib, SealevelMasks* masks)=0; + virtual void SealevelchangeConvolution(Vector* sealevelloads, Vector* oceanareas, IssmDouble* allsealevelloads, IssmDouble* allloads,IssmDouble* rotationaxismotionvector)=0; + virtual void SealevelchangeDeformationConvolution(IssmDouble* sealevelloads, IssmDouble* loads, IssmDouble* rotationvector)=0; + virtual void SealevelchangeShift(Vector* loads, IssmDouble offset, SealevelMasks* masks)=0; #endif Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 26164) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 26165) @@ -219,18 +219,9 @@ #endif #ifdef _HAVE_SEALEVELCHANGE_ - IssmDouble OceanAverage(IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet!");}; - void SetSealevelMasks(SealevelMasks* masks){_error_("not implemented yet!");}; - void SealevelchangeMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks){_error_("not implemented yet!");}; - void SealevelchangeMomentOfInertiaElement(IssmDouble* dI_list, IssmDouble* loads, IssmDouble* sealevelloads){_error_("not implemented yet!");}; - void SealevelchangeGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, - IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze){_error_("not implemented yet!");}; - IssmDouble SealevelchangeBarystaticIce(IssmDouble* Sgi, SealevelMasks* masks, Vector* barystatic_contribution, IssmDouble* partition, IssmDouble oceanarea){_error_("not implemented yet!");}; - IssmDouble SealevelchangeBarystaticHydro(IssmDouble* Sgi, SealevelMasks* masks, Vector* barystatic_contribution, IssmDouble* partition,IssmDouble oceanarea){_error_("not implemented yet!");}; - void SealevelchangeBarystaticBottomPressure(IssmDouble* Sgi,SealevelMasks* masks){_error_("not implemented yet!");}; - void SealevelchangeSal(IssmDouble* Sgo, IssmDouble* Sg_old, SealevelMasks* masks){_error_("not implemented yet!");}; - void DeformationFromSurfaceLoads(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet!");}; - void GiaDeflection(Vector* wg,Vector* dwgdt,Matlitho* litho, IssmDouble* x,IssmDouble* y){_error_("not implemented yet");}; - - 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");}; + void SetSealevelMasks(SealevelMasks* masks){_error_("not implemented yet!");}; + void SealevelchangeMomentOfInertia(IssmDouble* dI_list, IssmDouble* loads, IssmDouble* sealevelloads){_error_("not implemented yet!");}; + void GiaDeflection(Vector* wg,Vector* dwgdt,Matlitho* litho, IssmDouble* x,IssmDouble* y){_error_("not implemented yet");}; + + 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");}; void SealevelchangeBarystaticLoads(Vector* loads, BarystaticContributions* barycontrib, SealevelMasks* masks){_error_("not implemented yet");}; void SealevelchangeConvolution(Vector* sealevelloads, Vector* oceanareas, IssmDouble* allsealevelloads, IssmDouble* allloads,IssmDouble* rotationaxismotionvector){_error_("not implemented yet");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 26164) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 26165) @@ -174,18 +174,8 @@ #endif #ifdef _HAVE_SEALEVELCHANGE_ - void SealevelchangeMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks){_error_("not implemented yet!");}; - void SealevelchangeMomentOfInertiaElement(IssmDouble* dI_list, IssmDouble* loads, IssmDouble* sealevelloads){_error_("not implemented yet!");}; - void SetSealevelMasks(SealevelMasks* masks){_error_("not implemented yet!");}; - void SealevelchangeGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, - IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze){_error_("not implemented yet!");}; - IssmDouble SealevelchangeBarystaticIce(IssmDouble* Sgi, SealevelMasks* masks, Vector* barystatic_contribution, IssmDouble* partition, IssmDouble oceanarea){_error_("not implemented yet!");}; - IssmDouble SealevelchangeBarystaticHydro(IssmDouble* Sgi, SealevelMasks* masks, Vector* barystatic_contribution, IssmDouble* partition,IssmDouble oceanarea){_error_("not implemented yet!");}; - void SealevelchangeBarystaticBottomPressure(IssmDouble* Sgi,SealevelMasks* masks){_error_("not implemented yet!");}; - void SealevelchangeSal(IssmDouble* Sgo,IssmDouble* Sg_old, SealevelMasks* masks){_error_("not implemented yet!");}; - void DeformationFromSurfaceLoads(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet!");}; - IssmDouble OceanAverage(IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet!");}; - void GiaDeflection(Vector* wg,Vector* dwgdt,Matlitho* litho, IssmDouble* x,IssmDouble* y){_error_("not implemented yet");}; - - 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");}; + void SealevelchangeMomentOfInertia(IssmDouble* dI_list, IssmDouble* loads, IssmDouble* sealevelloads){_error_("not implemented yet!");}; + void SetSealevelMasks(SealevelMasks* masks){_error_("not implemented yet!");}; + void GiaDeflection(Vector* wg,Vector* dwgdt,Matlitho* litho, IssmDouble* x,IssmDouble* y){_error_("not implemented yet");}; + 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");}; void SealevelchangeBarystaticLoads(Vector* loads, BarystaticContributions* barycontrib, SealevelMasks* masks){_error_("not implemented yet");}; void SealevelchangeConvolution(Vector* sealevelloads, Vector* oceanareas, IssmDouble* allsealevelloads, IssmDouble* allloads,IssmDouble* rotationaxismotionvector){_error_("not implemented yet");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Tetra.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 26164) +++ /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 26165) @@ -181,17 +181,8 @@ #endif #ifdef _HAVE_SEALEVELCHANGE_ - void SetSealevelMasks(SealevelMasks* masks){_error_("not implemented yet!");}; - void SealevelchangeMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks){_error_("not implemented yet!");}; - void SealevelchangeMomentOfInertiaElement(IssmDouble* dI_list, IssmDouble* loads, IssmDouble* sealevelloads){_error_("not implemented yet!");}; - void SealevelchangeGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, - IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze){_error_("not implemented yet!");}; - IssmDouble SealevelchangeBarystaticIce(IssmDouble* Sgi, SealevelMasks* masks, Vector* barystatic_contribution, IssmDouble* partition, IssmDouble oceanarea){_error_("not implemented yet!");}; - IssmDouble SealevelchangeBarystaticHydro(IssmDouble* Sgi, SealevelMasks* masks, Vector* barystatic_contribution, IssmDouble* partition,IssmDouble oceanarea){_error_("not implemented yet!");}; - void SealevelchangeBarystaticBottomPressure(IssmDouble* Sgi,SealevelMasks* masks){_error_("not implemented yet!");}; - void SealevelchangeSal(IssmDouble* Sgo, IssmDouble* Sg_old, SealevelMasks* masks){_error_("not implemented yet!");}; - void DeformationFromSurfaceLoads(IssmDouble* Up ,IssmDouble* North, IssmDouble* East, IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet!");}; - IssmDouble OceanAverage(IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet!");}; - void GiaDeflection(Vector* wg,Vector* dwgdt, Matlitho* litho, IssmDouble* x,IssmDouble* y){_error_("not implemented yet");}; - 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");}; + void SetSealevelMasks(SealevelMasks* masks){_error_("not implemented yet!");}; + void SealevelchangeMomentOfInertia(IssmDouble* dI_list, IssmDouble* loads, IssmDouble* sealevelloads){_error_("not implemented yet!");}; + void GiaDeflection(Vector* wg,Vector* dwgdt, Matlitho* litho, IssmDouble* x,IssmDouble* y){_error_("not implemented yet");}; + 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");}; void SealevelchangeBarystaticLoads(Vector* loads, BarystaticContributions* barycontrib, SealevelMasks* masks){_error_("not implemented yet");}; void SealevelchangeConvolution(Vector* sealevelloads, Vector* oceanareas, IssmDouble* allsealevelloads, IssmDouble* allloads,IssmDouble* rotationaxismotionvector){_error_("not implemented yet");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 26164) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 26165) @@ -5556,175 +5556,4 @@ #ifdef _HAVE_SEALEVELCHANGE_ //old code -IssmDouble Tria::OceanAverage(IssmDouble* Sg, SealevelMasks* masks){ /*{{{*/ - - if(masks->isoceanin[this->lid]){ - - IssmDouble area; - - /*Get area of element:*/ - this->Element::GetInputValue(&area,AreaEnum); - - /*Average Sg over vertices:*/ - IssmDouble Sg_avg=0; for(int i=0;ivertices[i]->Sid()]/NUMVERTICES; - - /*return: */ - return area*Sg_avg; - } - else return 0; - -} -/*}}}*/ -void Tria::SealevelchangeMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks){/*{{{*/ - /*early return if we are not on an ice cap OR ocean:*/ - if(!masks->isiceonly[this->lid] && !masks->isoceanin[this->lid]){ - dI_list[0] = 0.0; // this is important!!! - dI_list[1] = 0.0; // this is important!!! - dI_list[2] = 0.0; // this is important!!! - return; - } - - /*Compute area of element:*/ - IssmDouble area,planetarea; - this->Element::GetInputValue(&area,AreaEnum); - - /*recover earth area: */ - this->parameters->FindParam(&planetarea,SolidearthPlanetAreaEnum); - - /*Compute lat,long,radius of elemental centroid: */ - bool spherical=true; - IssmDouble llr_list[NUMVERTICES][3]; - IssmDouble late,longe,re; - /* Where is the centroid of this element?:{{{*/ - ::GetVerticesCoordinates(&llr_list[0][0],this->vertices,NUMVERTICES,spherical); - - IssmDouble minlong=400; - IssmDouble maxlong=-20; - for (int i=0;imaxlong)maxlong=llr_list[i][1]; - if(llr_list[i][1]180){ - if (llr_list[0][1]==0)llr_list[0][1]=360; - if (llr_list[1][1]==0)llr_list[1][1]=360; - if (llr_list[2][1]==0)llr_list[2][1]=360; - } - - // correction at the north pole - if(llr_list[0][0]==0)llr_list[0][1]=(llr_list[1][1]+llr_list[2][1])/2.0; - if(llr_list[1][0]==0)llr_list[1][1]=(llr_list[0][1]+llr_list[2][1])/2.0; - if(llr_list[2][0]==0)llr_list[2][1]=(llr_list[0][1]+llr_list[1][1])/2.0; - - //correction at the south pole - if(llr_list[0][0]==180)llr_list[0][1]=(llr_list[1][1]+llr_list[2][1])/2.0; - if(llr_list[1][0]==180)llr_list[1][1]=(llr_list[0][1]+llr_list[2][1])/2.0; - if(llr_list[2][0]==180)llr_list[2][1]=(llr_list[0][1]+llr_list[1][1])/2.0; - - late=(llr_list[0][0]+llr_list[1][0]+llr_list[2][0])/3.0; - longe=(llr_list[0][1]+llr_list[1][1]+llr_list[2][1])/3.0; - - late=90.-late; - if(longe>180.)longe=(longe-180.)-180.; - - late=late/180.*M_PI; - longe=longe/180.*M_PI; - /*}}}*/ - re=(llr_list[0][2]+llr_list[1][2]+llr_list[2][2])/3.0; - - if(masks->isoceanin[this->lid]){ - IssmDouble rho_water, S; - - /*recover material parameters: */ - rho_water=FindParam(MaterialsRhoSeawaterEnum); - - /*From Sg_old, recover water sea level change:*/ - S=0; for(int i=0;ivertices[i]->Sid()]/NUMVERTICES; - - /* Perturbation terms for moment of inertia (moi_list): - * computed analytically (see Wu & Peltier, eqs 10 & 32) - * also consistent with my GMD formulation! - * ALL in geographic coordinates - * */ - dI_list[0] = -4*M_PI*(rho_water*S*area)*pow(re,4)*(sin(late)*cos(late)*cos(longe))/planetarea; - dI_list[1] = -4*M_PI*(rho_water*S*area)*pow(re,4)*(sin(late)*cos(late)*sin(longe))/planetarea; - dI_list[2] = +4*M_PI*(rho_water*S*area)*pow(re,4)*(1-pow(sin(late),2))/planetarea; - } - if(masks->isiceonly[this->lid]){ - - IssmDouble rho_ice,I; - bool computerigid= false; - - bool notfullygrounded=false; - - bool scaleoceanarea= false; - int glfraction=1; - IssmDouble phi=1.0; - /*{{{*/ - if (masks->notfullygrounded[this->lid]) notfullygrounded=true; //used later on. - - /*recover some parameters:*/ - rho_ice=FindParam(MaterialsRhoIceEnum); - this->parameters->FindParam(&computerigid,SolidearthSettingsRigidEnum); - this->parameters->FindParam(&scaleoceanarea,SolidearthSettingsOceanAreaScalingEnum); - this->parameters->FindParam(&glfraction,SolidearthSettingsGlfractionEnum); - - /*Get area of element: precomputed in the sealevelchange_geometry:*/ - this->Element::GetInputValue(&area,AreaEnum); - - /*Compute fraction of the element that is grounded: */ - if(notfullygrounded){ - IssmDouble xyz_list[NUMVERTICES][3]; - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - - 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 - if(glfraction==0)phi=1; - #ifdef _ISSM_DEBUG_ - this->AddInput(SealevelBarystaticMaskEnum,&phi,P0Enum); - #endif - } - else phi=1.0; - - /*Retrieve surface load for ice: */ - Input* deltathickness_input=this->GetInput(DeltaIceThicknessEnum); - if (!deltathickness_input)_error_("delta thickness input needed to compute sea level change!"); - - /*/Average ice thickness over grounded area of the element only: {{{*/ - if(!notfullygrounded)deltathickness_input->GetInputAverage(&I); - else{ - IssmDouble total_weight=0; - bool mainlyfloating = true; - int point1; - IssmDouble fraction1,fraction2; - - /*Recover portion of element that is grounded*/ - this->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating); - Gauss* gauss = this->NewGauss(point1,fraction1,fraction2,mainlyfloating,2); - - /* Start looping on the number of gaussian points and average over these gaussian points: */ - total_weight=0; - I=0; - while(gauss->next()){ - IssmDouble Ig=0; - deltathickness_input->GetInputValue(&Ig,gauss); - I+=Ig*gauss->weight; - total_weight+=gauss->weight; - } - if(total_weight) I=I/total_weight; - delete gauss; - } - /*}}}*/ - - /*}}}*/ - /*convert to kg/m^2*/ - I=I*phi; - - dI_list[0] = -4*M_PI*(rho_ice*I*area)*pow(re,4)*(sin(late)*cos(late)*cos(longe))/planetarea; - dI_list[1] = -4*M_PI*(rho_ice*I*area)*pow(re,4)*(sin(late)*cos(late)*sin(longe))/planetarea; - dI_list[2] = +4*M_PI*(rho_ice*I*area)*pow(re,4)*(1-pow(sin(late),2))/planetarea; - } - - return; -}/*}}}*/ void Tria::SetSealevelMasks(SealevelMasks* masks){ /*{{{*/ @@ -5740,636 +5569,4 @@ if ((gr_input->GetInputMin())<0) masks->notfullygrounded[this->lid]=true; else masks->notfullygrounded[this->lid]=false; -} -/*}}}*/ -void Tria::SealevelchangeGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze){ /*{{{*/ - - /*diverse:*/ - int gsize; - bool spherical=true; - IssmDouble llr_list[NUMVERTICES][3]; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble area,planetarea,planetradius; - IssmDouble I; //change in ice thickness or water level(Farrel and Clarke, Equ. 4) - IssmDouble rho_earth; - IssmDouble late,longe,re; - IssmDouble lati,longi,ri; - IssmDouble constant; - IssmDouble x_element,y_element,z_element,x,y,z,dx,dy,dz,N_azim,E_azim; - - #ifdef _HAVE_RESTRICT_ - IssmDouble* __restrict__ G=NULL; - IssmDouble* __restrict__ GU=NULL; - IssmDouble* __restrict__ GN=NULL; - IssmDouble* __restrict__ GE=NULL; - IssmDouble* __restrict__ G_elastic_precomputed=NULL; - IssmDouble* __restrict__ G_rigid_precomputed=NULL; - IssmDouble* __restrict__ U_elastic_precomputed=NULL; - IssmDouble* __restrict__ H_elastic_precomputed=NULL; - IssmDouble* __restrict__ indices=NULL; - #else - IssmDouble* G=NULL; - IssmDouble* GU=NULL; - IssmDouble* GN=NULL; - IssmDouble* GE=NULL; - IssmDouble* G_elastic_precomputed=NULL; - IssmDouble* G_rigid_precomputed=NULL; - IssmDouble* U_elastic_precomputed=NULL; - IssmDouble* H_elastic_precomputed=NULL; - IssmDouble* indices=NULL; - #endif - - /*elastic green function:*/ - int index; - int M,Mtest; - - /*Computational flags:*/ - bool computerigid = false; - bool computeelastic = false; - int horiz; - - /*recover parameters: */ - rho_earth=FindParam(MaterialsEarthDensityEnum); - this->parameters->FindParam(&computerigid,SolidearthSettingsRigidEnum); - this->parameters->FindParam(&computeelastic,SolidearthSettingsElasticEnum); - this->parameters->FindParam(&gsize,MeshNumberofverticesEnum); - this->parameters->FindParam(&planetarea,SolidearthPlanetAreaEnum); - this->parameters->FindParam(&planetradius,SolidearthPlanetRadiusEnum); - this->parameters->FindParam(&horiz,SolidearthSettingsHorizEnum); - - /*compute area and add to inputs:*/ - area=GetAreaSpherical(); - this->inputs->SetDoubleInput(AreaEnum,this->lid,area); - this->AddInput(SealevelAreaEnum,&area,P0Enum); - if(!computerigid)return; - - /*recover precomputed green function kernels:*/ - DoubleVecParam* parameter = static_cast(this->parameters->FindParamObject(SealevelchangeGRigidEnum)); _assert_(parameter); - parameter->GetParameterValueByPointer((IssmDouble**)&G_rigid_precomputed,&M); - _assert_(M>0); - - /*allocate indices:*/ - indices=xNew(gsize); - G=xNewZeroInit(gsize); - - if(computeelastic){ - parameter = static_cast(this->parameters->FindParamObject(SealevelchangeGElasticEnum)); _assert_(parameter); - parameter->GetParameterValueByPointer((IssmDouble**)&G_elastic_precomputed,&Mtest); _assert_(Mtest==M); - - parameter = static_cast(this->parameters->FindParamObject(SealevelchangeHElasticEnum)); _assert_(parameter); - parameter->GetParameterValueByPointer((IssmDouble**)&H_elastic_precomputed,&Mtest); _assert_(Mtest==M); - - parameter = static_cast(this->parameters->FindParamObject(SealevelchangeUElasticEnum)); _assert_(parameter); - parameter->GetParameterValueByPointer((IssmDouble**)&U_elastic_precomputed,&Mtest); _assert_(Mtest==M); - - GU=xNewZeroInit(gsize); - if(horiz){ - GN=xNewZeroInit(gsize); - GE=xNewZeroInit(gsize); - } - } - - /* Where is the centroid of this element:*/ - - /*retrieve coordinates: */ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - - /*figure out gravity center of our element (Cartesian): */ - x_element=(xyz_list[0][0]+xyz_list[1][0]+xyz_list[2][0])/3.0; - y_element=(xyz_list[0][1]+xyz_list[1][1]+xyz_list[2][1])/3.0; - z_element=(xyz_list[0][2]+xyz_list[1][2]+xyz_list[2][2])/3.0; - - /*compute gravity center in lat,long: */ - late= asin(z_element/planetradius); - longe = atan2(y_element,x_element); - - constant=3/rho_earth*area/planetarea; - - for(int i=0;iM_PI)delLambda=2*M_PI-delLambda; - alpha=2.*asin(sqrt(pow(sin(delPhi/2),2)+cos(lati)*cos(late)*pow(sin(delLambda/2),2))); - indices[i]=alpha/M_PI*reCast(M-1); - index=reCast(indices[i]); - _assert_(index>=0); _assert_(indexinputs->SetArrayInput(SealevelchangeGEnum,this->lid,G,gsize); - if(computeelastic){ - this->inputs->SetArrayInput(SealevelchangeGUEnum,this->lid,GU,gsize); - if(horiz){ - this->inputs->SetArrayInput(SealevelchangeGNEnum,this->lid,GN,gsize); - this->inputs->SetArrayInput(SealevelchangeGEEnum,this->lid,GE,gsize); - } - } - this->inputs->SetDoubleInput(AreaEnum,this->lid,area); - this->AddInput(SealevelAreaEnum,&area,P0Enum); - - /*Free allocations:*/ - #ifdef _HAVE_RESTRICT_ - delete indices; - delete G; - if(computeelastic){ - delete GU; - if(horiz){ - delete GN; - delete GE; - } - } - #else - xDelete(indices); - xDelete(G); - if(computeelastic){ - xDelete(GU); - if(horiz){ - xDelete(GN); - xDelete(GE); - } - } - #endif - - return; -} -/*}}}*/ -IssmDouble Tria::SealevelchangeBarystaticIce(IssmDouble* Sgi, SealevelMasks* masks, Vector* barystatic_contribution, IssmDouble* partition, IssmDouble oceanarea){ /*{{{*/ - - /*diverse:*/ - int gsize; - IssmDouble area; - IssmDouble phi=1.0; //WARNING: do not touch this, default is entire elemnt contributes barystatic - IssmDouble I; //change in ice thickness or water level(Farrel and Clarke, Equ. 4) - bool notfullygrounded=false; - bool scaleoceanarea= false; - bool computerigid= false; - int glfraction=1; - - /*output: */ - IssmDouble bslcice=0; - - /*elastic green function:*/ - IssmDouble* G=NULL; - - /*ice properties: */ - IssmDouble rho_ice,rho_water; - - /*constants:*/ - IssmDouble constant=0; - - /*early return if we are not on an ice cap:*/ - if(!masks->isiceonly[this->lid]){ - #ifdef _ISSM_DEBUG_ - constant=0; this->AddInput(SealevelBarystaticMaskEnum,&constant,P0Enum); - #endif - bslcice=0; - return bslcice; - } - - /*early return if we are fully floating:*/ - if (masks->isfullyfloating[this->lid]){ - constant=0; - #ifdef _ISSM_DEBUG_ - this->AddInput(SealevelBarystaticMaskEnum,&constant,P0Enum); - #endif - bslcice=0; - return bslcice; - } - - /*If we are here, we are on ice that is fully grounded or half-way to floating:*/ - if (masks->notfullygrounded[this->lid]) notfullygrounded=true; //used later on. - - /*Inform mask: */ - constant=1; - #ifdef _ISSM_DEBUG_ - this->AddInput(SealevelBarystaticMaskEnum,&constant,P0Enum); - #endif - - /*recover some parameters:*/ - rho_ice=FindParam(MaterialsRhoIceEnum); - rho_water=FindParam(MaterialsRhoSeawaterEnum); - this->parameters->FindParam(&computerigid,SolidearthSettingsRigidEnum); - this->parameters->FindParam(&scaleoceanarea,SolidearthSettingsOceanAreaScalingEnum); - this->parameters->FindParam(&glfraction,SolidearthSettingsGlfractionEnum); - - /*retrieve precomputed G:*/ - if(computerigid)this->inputs->GetArrayPtr(SealevelchangeGEnum,this->lid,&G,&gsize); - - /*Get area of element: precomputed in the sealevelchange_geometry:*/ - this->Element::GetInputValue(&area,AreaEnum); - - /*Compute fraction of the element that is grounded: */ - if(notfullygrounded){ - IssmDouble xyz_list[NUMVERTICES][3]; - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - - 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 - if(glfraction==0)phi=1; - #ifdef _ISSM_DEBUG_ - this->AddInput(SealevelBarystaticMaskEnum,&phi,P0Enum); - #endif - } - else phi=1.0; - - /*Retrieve surface load for ice: */ - Input* deltathickness_input=this->GetInput(DeltaIceThicknessEnum); - if (!deltathickness_input)_error_("delta thickness input needed to compute sea level change!"); - - /*/Average ice thickness over grounded area of the element only: {{{*/ - if(!notfullygrounded)deltathickness_input->GetInputAverage(&I); - else{ - IssmDouble total_weight=0; - bool mainlyfloating = true; - int point1; - IssmDouble fraction1,fraction2; - - /*Recover portion of element that is grounded*/ - this->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating); - Gauss* gauss = this->NewGauss(point1,fraction1,fraction2,mainlyfloating,2); - - /* Start looping on the number of gaussian points and average over these gaussian points: */ - total_weight=0; - I=0; - while(gauss->next()){ - IssmDouble Ig=0; - deltathickness_input->GetInputValue(&Ig,gauss); - I+=Ig*gauss->weight; - total_weight+=gauss->weight; - } - I=I/total_weight; - delete gauss; - } - /*}}}*/ - - /*Compute barystatic contribution:*/ - _assert_(oceanarea>0.); - if(scaleoceanarea) oceanarea=3.619e+14; // use true ocean area, m^2 - bslcice = rho_ice*area*phi*I/(oceanarea*rho_water); - _assert_(!xIsNan(bslcice)); - - if(computerigid){ - /*convert from m to kg/m^2:*/ - I=I*rho_ice*phi; - - /*convolve:*/ - for(int i=0;i(partition[this->Sid()])+1; - barystatic_contribution->SetValue(id,bslcice,ADD_VAL); - } - /*return :*/ - return bslcice; -} -/*}}}*/ -IssmDouble Tria::SealevelchangeBarystaticHydro(IssmDouble* Sgi, SealevelMasks* masks, Vector* barystatic_contribution, IssmDouble* partition, IssmDouble oceanarea){ /*{{{*/ - - /*diverse:*/ - int gsize; - IssmDouble area; - IssmDouble phi=1.0; //WARNING: do not touch this, default is entire elemnt contributes barystatic - IssmDouble W; //change in water height thickness (Farrel and Clarke, Equ. 4) - bool notfullygrounded=false; - bool scaleoceanarea= false; - bool computeelastic= false; - - /*elastic green function:*/ - IssmDouble* G=NULL; - - /*ice properties: */ - IssmDouble rho_water; - IssmDouble rho_freshwater; - - /*output:*/ - IssmDouble bslchydro = 0; - - /*early return if we are on an ice cap. Nop, we may have hydro - * and ice on the same masscon:*/ - /*if(masks->isiceonly[this->lid]){ - bslchydro=0; - return bslchydro; - }*/ - - /*early return if we are fully floating:*/ - if (masks->isfullyfloating[this->lid]){ - bslchydro=0; - return bslchydro; - } - - /*If we are here, we are on earth, not on ice: */ - - /*recover parameters: */ - rho_water=FindParam(MaterialsRhoSeawaterEnum); - rho_freshwater=FindParam(MaterialsRhoFreshwaterEnum); - this->parameters->FindParam(&computeelastic,SolidearthSettingsElasticEnum); - this->parameters->FindParam(&scaleoceanarea,SolidearthSettingsOceanAreaScalingEnum); - - /*retrieve precomputed G:*/ - if(computeelastic)this->inputs->GetArrayPtr(SealevelchangeGEnum,this->lid,&G,&gsize); - - /*Get area of element: precomputed in the sealevelchange_geometry:*/ - this->Element::GetInputValue(&area,AreaEnum); - - /*Retrieve water height at vertices: */ - Input* deltathickness_input=this->GetInput(DeltaTwsEnum); - if (!deltathickness_input)_error_("DeltaTwsEnum input needed to compute sea level change!"); - deltathickness_input->GetInputAverage(&W); - - /*Compute barystatic component:*/ - _assert_(oceanarea>0.); - if(scaleoceanarea) oceanarea=3.619e+14; // use true ocean area, m^2 - bslchydro = rho_freshwater*area*phi*W/(oceanarea*rho_water); - _assert_(!xIsNan(bslchydro)); - - if(computeelastic){ - /*convert from m to kg/m^2:*/ - W=W*rho_freshwater*phi; - - /*convolve:*/ - for(int i=0;i(partition[this->Sid()])+1; - barystatic_contribution->SetValue(id,bslchydro,ADD_VAL); - } - /*output:*/ - return bslchydro; -} -/*}}}*/ -void Tria::SealevelchangeBarystaticBottomPressure(IssmDouble* Sgi,SealevelMasks* masks){ /*{{{*/ - - /*diverse:*/ - int gsize; - IssmDouble area; - IssmDouble BP; //change in bottom pressure (Farrel and Clarke, Equ. 4) - IssmDouble constant; - bool computeelastic= false; - - /*elastic green function:*/ - IssmDouble* G=NULL; - - /*water properties: */ - IssmDouble rho_water; - - /*exit now?:*/ - this->parameters->FindParam(&computeelastic,SolidearthSettingsElasticEnum); - if(!computeelastic)return; - - /*we are here to compute fingerprints originating fromn bottom pressure changes:*/ - if(!masks->isoceanin[this->lid])return; - - /*Inform mask: */ - #ifdef _ISSM_DEBUG_ - constant=1; this->AddInput(SealevelBarystaticMaskEnum,&constant,P0Enum); - #endif - - /*recover material parameters: */ - rho_water=FindParam(MaterialsRhoSeawaterEnum); - - /*retrieve precomputed G:*/ - this->inputs->GetArrayPtr(SealevelchangeGEnum,this->lid,&G,&gsize); - - /*Get area of element: precomputed in the sealevelchange_geometry:*/ - this->Element::GetInputValue(&area,AreaEnum); - - /*Retrieve bottom pressure change and average over the element: */ - Input* bottompressure_change_input=this->GetInput(DeltaBottomPressureEnum); - if (!bottompressure_change_input)_error_("DeltaBottomPressureEnum pressure input needed to compute sea level change fingerprint!"); - bottompressure_change_input->GetInputAverage(&BP); - - /*convert from m to kg/m^2:*/ - BP=BP*rho_water; - - /*convolve:*/ - for(int i=0;iisoceanin[this->lid]){ - constant=0; - #ifdef _ISSM_DEBUG_ - this->AddInput(SealevelBarystaticOceanMaskEnum,&constant,P0Enum); - #endif - return; - } - constant=1; - #ifdef _ISSM_DEBUG_ - this->AddInput(SealevelBarystaticOceanMaskEnum,&constant,P0Enum); - #endif - - /*how many dofs are we working with here? */ - this->parameters->FindParam(&gsize,MeshNumberofverticesEnum); - - /*retrieve precomputed G:*/ - this->inputs->GetArrayPtr(SealevelchangeGEnum,this->lid,&G,&dummy); _assert_(dummy==gsize); - - /*From Sg_old, recover water sea level change:*/ - S=0; for(int i=0;ivertices[i]->Sid()]/NUMVERTICES; - - /*convert to kg/m^2: */ - S=S*rho_water; - - for(int i=0;iisiceonly[this->lid] && !masks->isoceanin[this->lid]) return; - - /*recover parameters:*/ - this->parameters->FindParam(&computeelastic,SolidearthSettingsElasticEnum); - this->parameters->FindParam(&horiz,SolidearthSettingsHorizEnum); - this->parameters->FindParam(&bp_compute_fingerprints,SolidearthSettingsComputeBpGrdEnum); - - /*early return if elastic not requested:*/ - if(!computeelastic) return; - - /*recover material parameters: */ - rho_ice=FindParam(MaterialsRhoIceEnum); - rho_water=FindParam(MaterialsRhoSeawaterEnum); - - /*recover elastic Green's functions for displacement:*/ - this->inputs->GetArrayPtr(SealevelchangeGUEnum,this->lid,&GU,&gsize); - if(horiz){ - this->inputs->GetArrayPtr(SealevelchangeGEEnum,this->lid,&GE,&gsize); - this->inputs->GetArrayPtr(SealevelchangeGNEnum,this->lid,&GN,&gsize); - } - - - if(masks->isoceanin[this->lid]){ - /*From Sg, recover water sea level change:*/ - S=0; for(int i=0;ivertices[i]->Sid()]/NUMVERTICES; - - /*convert to kg/m^2:*/ - S=rho_water*S; - - /*If bottom pressures are available, retrieve them to load the bedrock:*/ - if(bp_compute_fingerprints){ - Input* bottompressure_change_input=this->GetInput(DeltaBottomPressureEnum); - if (!bottompressure_change_input)_error_("bottom pressure input needed to compute sea level change fingerprint!"); - bottompressure_change_input->GetInputAverage(&BP); - - /*convert from m to kg/m^2:*/ - BP=BP*rho_water; - - S+=BP; - } - - for(int i=0;iisiceonly[this->lid]){ - - if (masks->isfullyfloating[this->lid]) return; - - - if (masks->notfullygrounded[this->lid]) notfullygrounded=true; //used later on. - - /*recover some parameters:*/ - rho_ice=FindParam(MaterialsRhoIceEnum); - rho_water=FindParam(MaterialsRhoSeawaterEnum); - this->parameters->FindParam(&computerigid,SolidearthSettingsRigidEnum); - this->parameters->FindParam(&scaleoceanarea,SolidearthSettingsOceanAreaScalingEnum); - this->parameters->FindParam(&glfraction,SolidearthSettingsGlfractionEnum); - - /*Get area of element: precomputed in the sealevelchange_geometry:*/ - this->Element::GetInputValue(&area,AreaEnum); - - /*Compute fraction of the element that is grounded: */ - if(notfullygrounded){ - IssmDouble xyz_list[NUMVERTICES][3]; - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - - 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 - if(glfraction==0)phi=1; - #ifdef _ISSM_DEBUG_ - this->AddInput(SealevelBarystaticMaskEnum,&phi,P0Enum); - #endif - } - else phi=1.0; - - /*Retrieve surface load for ice: */ - Input* deltathickness_input=this->GetInput(DeltaIceThicknessEnum); - if (!deltathickness_input)_error_("delta thickness input needed to compute sea level change!"); - - /*/Average ice thickness over grounded area of the element only: {{{*/ - if(!notfullygrounded)deltathickness_input->GetInputAverage(&I); - else{ - IssmDouble total_weight=0; - bool mainlyfloating = true; - int point1; - IssmDouble fraction1,fraction2; - - /*Recover portion of element that is grounded*/ - this->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating); - Gauss* gauss = this->NewGauss(point1,fraction1,fraction2,mainlyfloating,2); - - /* Start looping on the number of gaussian points and average over these gaussian points: */ - total_weight=0; - I=0; - while(gauss->next()){ - IssmDouble Ig=0; - deltathickness_input->GetInputValue(&Ig,gauss); - I+=Ig*gauss->weight; - total_weight+=gauss->weight; - } - I=I/total_weight; - delete gauss; - } - /*}}}*/ - - /*convert to kg/m^2*/ - I=I*rho_ice*phi; - - for(int i=0;i* barystatic_contribution, IssmDouble* partition, IssmDouble oceanarea); - IssmDouble SealevelchangeBarystaticHydro(IssmDouble* Sgi, SealevelMasks* masks, Vector* barystatic_contribution, IssmDouble* partition,IssmDouble oceanarea); - void SealevelchangeBarystaticBottomPressure(IssmDouble* Sgi,SealevelMasks* masks); - void SealevelchangeSal(IssmDouble* Sgo,IssmDouble* Sg_old,SealevelMasks* masks); - void DeformationFromSurfaceLoads(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg,SealevelMasks* masks); void GiaDeflection(Vector* wg,Vector* dwgdt,Matlitho* litho, IssmDouble* x,IssmDouble* y); - IssmDouble OceanAverage(IssmDouble* Sg, SealevelMasks* masks); void SetSealevelMasks(SealevelMasks* masks); - - void SealevelchangeGeometryOptim(IssmDouble* lat,IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae); + void SealevelchangeGeometry(IssmDouble* lat,IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae); void SealevelchangeBarystaticLoads(Vector* loads, BarystaticContributions* barycontrib, SealevelMasks* masks); void SealevelchangeConvolution(Vector* sealevelloads, Vector* oceanareas, IssmDouble* allsealevelloads, IssmDouble* allloads,IssmDouble* rotationaxismotionvector); Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 26164) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 26165) @@ -4748,408 +4748,4 @@ /*}}}*/ #endif -#ifdef _HAVE_SEALEVELCHANGE_ -void FemModel::SealevelchangeBarystatic(Vector* pRSLgi, IssmDouble* poceanarea, IssmDouble* pbslc,IssmDouble* pbslcice, IssmDouble* pbslchydro, IssmDouble** pbslcice_partition,IssmDouble** pbslchydro_partition,SealevelMasks* masks) { /*{{{*/ - - /*serialized vectors:*/ - IssmDouble bslcice = 0.; - IssmDouble bslcice_cpu = 0.; - IssmDouble bslchydro = 0.; - IssmDouble bslchydro_cpu = 0.; - IssmDouble area = 0.; - IssmDouble oceanarea = 0.; - IssmDouble oceanarea_cpu = 0.; - int bp_compute_fingerprints= 0; - bool isoceantransport=false; - - Vector* bslcice_partition=NULL; - IssmDouble* bslcice_partition_serial=NULL; - IssmDouble* partitionice=NULL; - int npartice,nel; - - Vector* bslchydro_partition=NULL; - IssmDouble* bslchydro_partition_serial=NULL; - IssmDouble* partitionhydro=NULL; - bool istws=0; - int nparthydro; - - int npartocean; - Vector* bslcocean_partition=NULL; - IssmDouble* bslcocean_partition_serial=NULL; - IssmDouble* partitionocean=NULL; - - /*Initialize temporary vector that will be used to sum barystatic components - * on all local elements, prior to assembly:*/ - int gsize = this->nodes->NumberOfDofs(GsetEnum); - IssmDouble* RSLgi=xNewZeroInit(gsize); - int* indices=xNew(gsize); - for(int i=0;ielements->objects){ - i +=1; - Element* element = xDynamicCast(object); - element->GetInputValue(&area,AreaEnum); - if (masks->isoceanin[i]) oceanarea_cpu += area; - } - ISSM_MPI_Reduce (&oceanarea_cpu,&oceanarea,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); - ISSM_MPI_Bcast(&oceanarea,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); - _assert_(oceanarea>0.); - - /*Initialize partition vectors to retrieve barystatic contributions: */ - this->parameters->FindParam(&npartice,SolidearthNpartIceEnum); - if(npartice){ - this->parameters->FindParam(&partitionice,&nel,NULL,SolidearthPartitionIceEnum); - bslcice_partition= new Vector(npartice); - } - - this->parameters->FindParam(&nparthydro,SolidearthNpartHydroEnum); - if(nparthydro){ - this->parameters->FindParam(&partitionhydro,&nel,NULL,SolidearthPartitionHydroEnum); - bslchydro_partition= new Vector(nparthydro); - } - - this->parameters->FindParam(&npartocean,SolidearthNpartOceanEnum); - if(npartocean){ - this->parameters->FindParam(&partitionocean,&nel,NULL,SolidearthPartitionOceanEnum); - bslchydro_partition= new Vector(npartocean); - } - /*For later: - npartbarystatic=npartice; - if(nparthydro>npartbarystatic)npartbarystatic=nparthydro; - if(npartocean>npartbarystatic)npartbarystatic=npartocean; - bslc_partition=new Matrix(IssmDouble>(npartbarystatic,3); - - bslc_cpu[0]=0; bslc_cpu[1]=0; bslc_cpu[2]=0; - for(Object* & object : this->elements->objects){ - Element* element = xDynamicCast(object); - element->SealevelchangeBarystaticLoads(&bslc_cpu[0], localloads,masks, bslcice_partition,partitionice,oceanarea); - } - MPI Bcast localloads -> loads - - for(Object* & object : this->elements->objects){ - Element* element = xDynamicCast(object); - element->SealevelchangeConvolution(loads); - } - */ - - - - - - /*Call the barystatic sea level change core for ice : */ - bslcice_cpu=0; - for(Object* & object : this->elements->objects){ - Element* element = xDynamicCast(object); - bslcice_cpu+=element->SealevelchangeBarystaticIce(RSLgi,masks, bslcice_partition,partitionice,oceanarea); - } - - /*Call the barystatic sea level change core for hydro: */ - bslchydro_cpu=0; //make sure to initialize this, so we have a total barystatic contribution computed at 0. - this->parameters->FindParam(&istws,TransientIshydrologyEnum); - if(istws){ - for(int i=0;iSize();i++){ - Element* element=xDynamicCast(elements->GetObjectByOffset(i)); - bslchydro_cpu+=element->SealevelchangeBarystaticHydro(RSLgi,masks, bslchydro_partition,partitionhydro,oceanarea); - } - } - - /*Call the barystatic sea level change core for bottom pressures: */ - this->parameters->FindParam(&bp_compute_fingerprints,SolidearthSettingsComputeBpGrdEnum); - this->parameters->FindParam(&isoceantransport,TransientIsoceantransportEnum); - if(bp_compute_fingerprints && isoceantransport){ - for(int i=0;iSize();i++){ - Element* element=xDynamicCast(elements->GetObjectByOffset(i)); - element->SealevelchangeBarystaticBottomPressure(RSLgi,masks); - } - } - - /*Plug values once and assemble: */ - pRSLgi->SetValues(gsize,indices,RSLgi,ADD_VAL); - pRSLgi->Assemble(); - - /*Sum all barystatic components from all cpus:*/ - ISSM_MPI_Reduce (&bslcice_cpu,&bslcice,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); - ISSM_MPI_Bcast(&bslcice,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); - _assert_(!xIsNan(bslcice)); - - ISSM_MPI_Reduce (&bslchydro_cpu,&bslchydro,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); - ISSM_MPI_Bcast(&bslchydro,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); - _assert_(!xIsNan(bslchydro)); - - /*Take care of partition vectors:*/ - if(bslcice_partition){ - bslcice_partition->Assemble(); - bslcice_partition_serial=bslcice_partition->ToMPISerial(); - } - if(bslchydro_partition){ - bslchydro_partition->Assemble(); - bslchydro_partition_serial=bslchydro_partition->ToMPISerial(); - } - - - /*Free ressources:*/ - xDelete(indices); - xDelete(RSLgi); - if(bslchydro_partition)delete bslchydro_partition; - if(bslcice_partition)delete bslcice_partition; - if(partitionhydro)xDelete(partitionhydro); - if(partitionice)xDelete(partitionice); - - /*Assign output pointers:*/ - *poceanarea = oceanarea; - *pbslcice = bslcice; - *pbslchydro = bslchydro; - *pbslc=bslchydro+bslcice; - *pbslcice_partition=bslcice_partition_serial; - *pbslchydro_partition=bslchydro_partition_serial; - -} -/*}}}*/ -void FemModel::SealevelchangeSal(Vector* pRSLgo, Vector* pRSLg_old, SealevelMasks* masks, bool verboseconvolution){/*{{{*/ - - /*serialized vectors:*/ - IssmDouble* RSLg_old=NULL; - - IssmDouble* RSLgo = NULL; - int* indices = NULL; - int gsize; - - bool computerigid = true; - bool computeelastic= true; - - /*recover computational flags: */ - this->parameters->FindParam(&computerigid,SolidearthSettingsRigidEnum); - - /*Initialize temporary vector that will be used to sum barystatic components on all local elements, prior - * to assembly:*/ - gsize = this->nodes->NumberOfDofs(GsetEnum); - RSLgo=xNewZeroInit(gsize); - indices=xNew(gsize); for (int i=0;iToMPISerial(); - - /*Call the sal sea level change core only if required: */ - if(computerigid){ - for(Object* & object : this->elements->objects){ - Element* element = xDynamicCast(object); - element->SealevelchangeSal(RSLgo,RSLg_old,masks); - } - } - pRSLgo->SetValues(gsize,indices,RSLgo,ADD_VAL); - pRSLgo->Assemble(); - - /*Free ressources:*/ - xDelete(indices); - xDelete(RSLgo); - xDelete(RSLg_old); - -} -/*}}}*/ -void FemModel::SealevelchangeRotationalFeedback(Vector* pRSLgo_rot, Vector* pRSLg_old, IssmDouble* pIxz, IssmDouble* pIyz, IssmDouble* pIzz, SealevelMasks* masks){/*{{{*/ - - /*serialized vectors:*/ - bool spherical=true; - IssmDouble* RSLg_old=NULL; - IssmDouble* tide_love_h = NULL; - IssmDouble* tide_love_k = NULL; - IssmDouble* load_love_k = NULL; - IssmDouble tide_love_k2secular; - IssmDouble moi_e, moi_p, omega, g; - IssmDouble m1, m2, m3; - IssmDouble lati, longi, radi, value; - IssmDouble *latitude = NULL; - IssmDouble *longitude = NULL; - IssmDouble *radius = NULL; - - /*Serialize vectors from previous iteration:*/ - RSLg_old=pRSLg_old->ToMPISerial(); - - IssmDouble moi_list[3]={0,0,0}; - IssmDouble moi_list_cpu[3]={0,0,0}; - for(Object* & object : this->elements->objects){ - Element* element = xDynamicCast(object); - element->SealevelchangeMomentOfInertia(&moi_list[0],RSLg_old,masks ); - moi_list_cpu[0] += moi_list[0]; - moi_list_cpu[1] += moi_list[1]; - moi_list_cpu[2] += moi_list[2]; - } - ISSM_MPI_Reduce (&moi_list_cpu[0],&moi_list[0],1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); - ISSM_MPI_Bcast(&moi_list[0],1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); - // - ISSM_MPI_Reduce (&moi_list_cpu[1],&moi_list[1],1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); - ISSM_MPI_Bcast(&moi_list[1],1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); - // - ISSM_MPI_Reduce (&moi_list_cpu[2],&moi_list[2],1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); - ISSM_MPI_Bcast(&moi_list[2],1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); - - /*pull out some useful parameters: */ - parameters->FindParam(&load_love_k,NULL,NULL,LoadLoveKEnum); - parameters->FindParam(&tide_love_h,NULL,NULL,TidalLoveHEnum); - parameters->FindParam(&tide_love_k,NULL,NULL,TidalLoveKEnum); - parameters->FindParam(&tide_love_k2secular,TidalLoveK2SecularEnum); - parameters->FindParam(&moi_e,RotationalEquatorialMoiEnum); - parameters->FindParam(&moi_p,RotationalPolarMoiEnum); - parameters->FindParam(&omega,RotationalAngularVelocityEnum); - - /*compute perturbation terms for angular velocity vector: */ - m1 = 1/(1-tide_love_k[2]/tide_love_k2secular) * (1+load_love_k[2])/(moi_p-moi_e) * moi_list[0]; - m2 = 1/(1-tide_love_k[2]/tide_love_k2secular) * (1+load_love_k[2])/(moi_p-moi_e) * moi_list[1]; - m3 = -(1+load_love_k[2])/moi_p * moi_list[2]; // term associated with fluid number (3-order-of-magnitude smaller) is negelected - - /*recover lat,long and radius vectors from vertices: */ - VertexCoordinatesx(&latitude,&longitude,&radius,this->vertices,spherical); - - /* Green's function (1+k_2-h_2/g): checked against Glenn Milne's thesis Chapter 3 (eqs: 3.3-4, 3.10-11) - * Perturbation terms for angular velocity vector (m1, m2, m3): checked against Mitrovica (2005 Appendix) & Jensen et al (2013 Appendix A3) - * Sea level rotational feedback: checked against GMD eqs 8-9 (only first order terms, i.e., degree 2 order 0 & 1 considered) - * all DONE in Geographic coordinates: theta \in [-90,90], lambda \in [-180 180] - */ - for(Object* & object : vertices->objects){ - Vertex* vertex=xDynamicCast(object); - int sid=vertex->Sid(); - - lati=latitude[sid]/180*PI; - longi=longitude[sid]/180*PI; - radi=radius[sid]; - - /*only first order terms are considered now: */ - value=((1.0+tide_love_k[2]-tide_love_h[2])/9.81)*pow(omega*radi,2.0)* - (-m3/6.0 + 0.5*m3*cos(2.0*lati) - 0.5*sin(2.*lati)*(m1*cos(longi)+m2*sin(longi))); - - pRSLgo_rot->SetValue(sid,value,INS_VAL); //INS_VAL ensures that you don't add several times - } - - /*Assemble mesh velocity*/ - pRSLgo_rot->Assemble(); - - /*Assign output pointers:*/ - if(pIxz)*pIxz=moi_list[0]; - if(pIyz)*pIyz=moi_list[1]; - if(pIzz)*pIzz=moi_list[2]; - xDelete(latitude); - xDelete(longitude); - xDelete(tide_love_h); - xDelete(tide_love_k); - xDelete(load_love_k); - - xDelete(radius); - - /*Free ressources:*/ - xDelete(RSLg_old); - -} -/*}}}*/ -void FemModel::SealevelchangeDeformation(Vector* pgeoid, Vector* pUp, Vector* pNorth, Vector* pEast, Vector* pRSLg, SealevelMasks* masks){/*{{{*/ - - /*serialized vectors:*/ - IssmDouble* RSLg=NULL; - - IssmDouble* Up = NULL; - IssmDouble* North = NULL; - IssmDouble* East = NULL; - int* indices = NULL; - int gsize; - int horiz; - - /*retrieve parameters:*/ - this->parameters->FindParam(&horiz,SolidearthSettingsHorizEnum); - - /*Serialize vectors from previous iteration:*/ - RSLg=pRSLg->ToMPISerial(); - - /*Initialize temporary vector that will be used to sum barystatic components on all local elements, prior - * to assembly:*/ - gsize = this->nodes->NumberOfDofs(GsetEnum); - Up=xNewZeroInit(gsize); - if(horiz){ - North=xNewZeroInit(gsize); - East=xNewZeroInit(gsize); - } - indices=xNew(gsize); for (int i=0;ielements->objects){ - Element* element = xDynamicCast(object); - element->DeformationFromSurfaceLoads(Up,North,East,RSLg,masks); - } - - pUp->SetValues(gsize,indices,Up,ADD_VAL); - pUp->Assemble(); - - - /*Add RSL to Up to find the geoid: */ - pUp->Copy(pgeoid); pgeoid->AXPY(pRSLg,1); - - if (horiz){ - pNorth->SetValues(gsize,indices,North,ADD_VAL); - pNorth->Assemble(); - pEast->SetValues(gsize,indices,East,ADD_VAL); - pEast->Assemble(); - } - - /*Free ressources:*/ - xDelete(Up); - if(horiz){ - xDelete(North); - xDelete(East); - } - xDelete(indices); - xDelete(RSLg); -} -/*}}}*/ -IssmDouble FemModel::SealevelchangeOceanAverage(Vector* RSLg,SealevelMasks* masks, IssmDouble oceanarea) { /*{{{*/ - - IssmDouble* RSLg_serial=NULL; - IssmDouble oceanvalue,oceanvalue_cpu; - - /*Serialize vectors from previous iteration:*/ - RSLg_serial=RSLg->ToMPISerial(); - - /*Initialize:*/ - oceanvalue_cpu=0; - - /*Go through elements, and add contribution from each element and divide by overall ocean area:*/ - for(Object* & object : this->elements->objects){ - Element* element = xDynamicCast(object); - oceanvalue_cpu += element->OceanAverage(RSLg_serial,masks); - } - - ISSM_MPI_Reduce (&oceanvalue_cpu,&oceanvalue,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); - ISSM_MPI_Bcast(&oceanvalue,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); - - /*Free ressources:*/ - xDelete(RSLg_serial); - - return oceanvalue/oceanarea; -} -/*}}}*/ -void FemModel::IvinsDeformation(Vector* wg,Vector* dwgdt, IssmDouble* x, IssmDouble* y){ /*{{{*/ - - /*Find the litho material to be used by all the elements:*/ - Matlitho* matlitho=NULL; - for (Object* & object: this->materials->objects){ - Material* material=xDynamicCast(object); - if(material->ObjectEnum()==MatlithoEnum){ - matlitho=xDynamicCast(material); - break; - } - } - _assert_(matlitho); - - - /*Go through elements, and add contribution from each element to the deflection vector wg:*/ - for(Object* & object : this->elements->objects){ - Element* element = xDynamicCast(object); - element->GiaDeflection(wg,dwgdt, matlitho, x,y); - } - - /*Assemble parallel vector:*/ - dwgdt->Assemble(); - wg->Assemble(); -} -/*}}}*/ -#endif void FemModel::HydrologyEPLupdateDomainx(IssmDouble* pEplcount){ /*{{{*/ Index: /issm/trunk-jpl/src/c/classes/FemModel.h =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.h (revision 26164) +++ /issm/trunk-jpl/src/c/classes/FemModel.h (revision 26165) @@ -163,12 +163,4 @@ void EsaGeodetic3D(Vector* pUp, Vector* pNorth, Vector* pEast, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz); #endif - #ifdef _HAVE_SEALEVELCHANGE_ - void SealevelchangeBarystatic(Vector* pRSLgi, IssmDouble* poceanarea, IssmDouble* pbslc,IssmDouble* pbslcice, IssmDouble* pbslchydro, IssmDouble** pbslcice_partition,IssmDouble** pbslchydro_partition, SealevelMasks* masks); - void SealevelchangeSal(Vector* pSgo, Vector* pSg_old, SealevelMasks* masks,bool verboseconvolution); - void SealevelchangeRotationalFeedback(Vector* pRSLgo_rot, Vector* pRSLg_old, IssmDouble* pIxz, IssmDouble* pIyz, IssmDouble* pIzz, SealevelMasks* masks); - void SealevelchangeDeformation(Vector* pgeoid,Vector* pUp, Vector* pNorth, Vector* pEast, Vector* pSg_old, SealevelMasks* masks); - IssmDouble SealevelchangeOceanAverage(Vector* Sg,SealevelMasks* masks, IssmDouble oceanarea); - void IvinsDeformation(Vector* wg,Vector* dwgdt, IssmDouble* x, IssmDouble* y); - #endif void HydrologyEPLupdateDomainx(IssmDouble* pEplcount); void HydrologyIDSupdateDomainx(IssmDouble* pIDScount); Index: /issm/trunk-jpl/src/c/cores/cores.h =================================================================== --- /issm/trunk-jpl/src/c/cores/cores.h (revision 26164) +++ /issm/trunk-jpl/src/c/cores/cores.h (revision 26165) @@ -63,10 +63,6 @@ #endif void grd_core(FemModel* femmodel); -void grd_core_optim(FemModel* femmodel); void solidearthexternal_core(FemModel* femmodel); void dynstr_core(FemModel* femmodel); -Vector* sealevelchange_core_barystatic(FemModel* femmodel,SealevelMasks* mask, IssmDouble* poceanarea); -Vector* sealevelchange_core_sal(FemModel* femmodel,SealevelMasks* masks, Vector* RSLg_barystatic,IssmDouble oceanarea); -void sealevelchange_core_deformation(Vector** pN_radial, Vector** pU_radial, Vector** pU_north,Vector** pU_east,FemModel* femmodel,Vector* RSLg, SealevelMasks* masks); void couplerinput_core(FemModel* femmodel); void coupleroutput_core(FemModel* femmodel); Index: /issm/trunk-jpl/src/c/cores/sealevelchange_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/sealevelchange_core.cpp (revision 26164) +++ /issm/trunk-jpl/src/c/cores/sealevelchange_core.cpp (revision 26165) @@ -21,4 +21,5 @@ void RotationAxisMotion(IssmDouble* m, FemModel* femmodel,IssmDouble* loads, IssmDouble* sealeveloads); void ConserveOceanMass(FemModel* femmodel,Vector* sealevelloads, IssmDouble offset, SealevelMasks* masks); +void ivins_deformation_core(FemModel* femmodel); /*}}}*/ @@ -31,9 +32,7 @@ /*Parameters, variables:*/ bool save_results; - int optim=0; /*Retrieve parameters:*/ femmodel->parameters->FindParam(&save_results,SaveResultsEnum); - femmodel->parameters->FindParam(&optim,SolidearthSettingsOptimEnum); /*Verbose: */ @@ -53,6 +52,5 @@ /*Run geodetic:*/ - if(!optim) grd_core(femmodel); - else grd_core_optim(femmodel); + grd_core(femmodel); /*Run steric core for sure:*/ @@ -208,35 +206,52 @@ }; /*}}}*/ -void grd_core(FemModel* femmodel){ /*{{{*/ - - /*Gravity rotation deformation core GRD: */ - - /*variables:*/ - Vector *RSLg = NULL; - Vector *RSLg_barystatic = NULL; - Vector *U_grd = NULL; - Vector *N_grd = NULL; - Vector *U_north_grd = NULL; - Vector *U_east_grd = NULL; - Vector *bedrock = NULL; - Vector *bedrockeast = NULL; - Vector *bedrocknorth = NULL; - Vector *geoid= NULL; +void grd_core(FemModel* femmodel) { /*{{{*/ + + /*variables:{{{*/ + int nel; + BarystaticContributions* barycontrib=NULL; SealevelMasks* masks=NULL; - - /*parameters:*/ + GenericParam* barycontribparam=NULL; + IssmDouble rotationaxismotionvector[3]={0}; + + Vector* loads=NULL; + IssmDouble* allloads=NULL; + Vector* sealevelloads=NULL; + Vector* oldsealevelloads=NULL; + IssmDouble* allsealevelloads=NULL; + Vector* oceanareas=NULL; + IssmDouble oceanarea; + IssmDouble oceanaverage; + bool scaleoceanarea=false; + IssmDouble rho_water; + + IssmDouble eps_rel; + IssmDouble eps_abs; + int max_nonlinear_iterations; + int iterations=0; + int step; + IssmDouble time; + bool converged=false; + int modelid,earthid; int horiz; - IssmDouble oceanarea; int count,frequency,iscoupling; int grd=0; + int grdmodel; + int computesealevel=0; + + /*}}}*/ /*Verbose: */ if(VerboseSolution()) _printf0_(" computing GRD patterns\n"); - - /*retrieve parameters:*/ + + /*retrieve parameters:{{{*/ femmodel->parameters->FindParam(&grd,SolidearthSettingsGRDEnum); femmodel->parameters->FindParam(&frequency,SolidearthSettingsRunFrequencyEnum); femmodel->parameters->FindParam(&count,SealevelchangeRunCountEnum); + femmodel->parameters->FindParam(&computesealevel,SolidearthSettingsComputesealevelchangeEnum); + femmodel->parameters->FindParam(&max_nonlinear_iterations,SolidearthSettingsMaxiterEnum); + femmodel->parameters->FindParam(&grdmodel,GrdModelEnum); + /*}}}*/ /*only run if grd was requested, if we are the earth, and we have reached @@ -251,63 +266,116 @@ } + /*branch directly to Ivins deformation core if requested:*/ + if(grdmodel==IvinsEnum){ + ivins_deformation_core(femmodel); + return; + } + + /*retrieve parameters: {{{*/ + femmodel->parameters->FindParam(&scaleoceanarea,SolidearthSettingsOceanAreaScalingEnum); + barycontribparam = xDynamicCast*>(femmodel->parameters->FindParamObject(BarystaticContributionsEnum)); + barycontrib=barycontribparam->GetParameterValue(); + femmodel->parameters->FindParam(&rho_water,MaterialsRhoSeawaterEnum); + femmodel->parameters->FindParam(&eps_rel,SolidearthSettingsReltolEnum); + femmodel->parameters->FindParam(&eps_abs,SolidearthSettingsAbstolEnum); + femmodel->parameters->FindParam(&horiz,SolidearthSettingsHorizEnum); + /*}}}*/ + + /*initialize matrices and vectors:*/ + femmodel->parameters->FindParam(&nel,MeshNumberofelementsEnum); + loads=new Vector(nel); + oceanareas=new Vector(nel); + sealevelloads=new Vector(nel); + sealevelloads->Set(0);sealevelloads->Assemble(); + /*call masks core: */ masks=sealevel_masks(femmodel); - - /*call barystatic core (generalized eustatic - Farrel and Clark, Eq 4, 1st, 3rd and 4rd terms on the RHS) */ - RSLg_barystatic=sealevelchange_core_barystatic(femmodel,masks,&oceanarea); - - /*call self attraction and loading module (ocean loading tems - 2nd and 5th terms on the RHS of Farrel and Clark) */ - RSLg=sealevelchange_core_sal(femmodel,masks,RSLg_barystatic,oceanarea); - - /*compute bedrock motion and derive geoid: */ - sealevelchange_core_deformation(&N_grd,&U_grd,&U_north_grd,&U_east_grd,femmodel,RSLg,masks); + if(VerboseSolution()) _printf0_(" starting GRD convolutions\n"); + + /*buildup loads: */ + for(Object* & object : femmodel->elements->objects){ + Element* element = xDynamicCast(object); + element->SealevelchangeBarystaticLoads(loads, barycontrib,masks); + } + loads->Assemble(); + + //broadcast loads + allloads=loads->ToMPISerial(); + + //compute rotation axis motion: + RotationAxisMotion(&rotationaxismotionvector[0],femmodel,allloads,NULL); + + /*skip computation of sea level if requested, which means sea level loads should be zeroed */ + if(!computesealevel){ + allsealevelloads=xNewZeroInit(nel); + goto deformation; + } + + if(VerboseSolution()) _printf0_(" converging GRD convolutions\n"); + for(;;){ + + oldsealevelloads=sealevelloads->Duplicate(); sealevelloads->Copy(oldsealevelloads); + + /*convolve load and sealevel loads on oceans:*/ + for(Object* & object : femmodel->elements->objects){ + Element* element = xDynamicCast(object); + element->SealevelchangeConvolution(sealevelloads, oceanareas, allsealevelloads, allloads,rotationaxismotionvector); + } + sealevelloads->Assemble(); + + /*compute ocean areas:*/ + if(!allsealevelloads){ //first time in the loop + oceanareas->Assemble(); oceanareas->Sum(&oceanarea); _assert_(oceanarea>0.); + if(scaleoceanarea) oceanarea=3.619e+14; // use true ocean area, m^2 + } + + //Conserve ocean mass: + oceanaverage=SealevelloadsOceanAverage(sealevelloads,oceanareas,oceanarea); + ConserveOceanMass(femmodel,sealevelloads,barycontrib->Total()/oceanarea - oceanaverage,masks); + + //broadcast sea level loads + allsealevelloads=sealevelloads->ToMPISerial(); + + //compute rotation axis motion: + RotationAxisMotion(&rotationaxismotionvector[0],femmodel,allloads,allsealevelloads); + + //convergence? + slcconvergence(&converged,sealevelloads,oldsealevelloads,eps_rel,eps_abs); + if (converged)break; + + //early return? + if(iterations>=max_nonlinear_iterations)break; + else iterations++; + } + + deformation: + + if(VerboseSolution()) _printf0_(" deformation GRD convolutions\n"); + + /*convolve loads and sea level loads to get the deformation:*/ + for(Object* & object : femmodel->elements->objects){ + Element* element = xDynamicCast(object); + element->SealevelchangeDeformationConvolution(allsealevelloads, allloads, rotationaxismotionvector); + } + + if(VerboseSolution()) _printf0_(" updating GRD fields\n"); /*Update bedrock motion and geoid:*/ - GetVectorFromInputsx(&geoid,femmodel,SealevelEnum,VertexSIdEnum); - GetVectorFromInputsx(&bedrock,femmodel,BedEnum,VertexSIdEnum); + if(computesealevel){ + femmodel->inputs->Shift(SealevelGRDEnum,barycontrib->Total()/rho_water/oceanarea- oceanaverage/rho_water); + + /*cumulate barystatic contributions and save to results: */ + barycontrib->Cumulate(femmodel->parameters); + barycontrib->Save(femmodel->results,femmodel->parameters,oceanarea); + } + + femmodel->inputs->AXPY(1,SealevelGRDEnum,SealevelEnum); + femmodel->inputs->AXPY(1,BedGRDEnum,BedEnum); if(horiz){ - GetVectorFromInputsx(&bedrockeast,femmodel,BedEastEnum,VertexSIdEnum); - GetVectorFromInputsx(&bedrocknorth,femmodel,BedNorthEnum,VertexSIdEnum); - } - - geoid->AXPY(N_grd,1); - bedrock->AXPY(U_grd,1); - if(horiz){ - bedrockeast->AXPY(U_east_grd,1); - bedrocknorth->AXPY(U_north_grd,1); - } - - /*get some of the updates into elements:*/ - InputUpdateFromVectorx(femmodel,U_grd,SealevelUGrdEnum,VertexSIdEnum); - InputUpdateFromVectorx(femmodel,N_grd,SealevelEnum,VertexSIdEnum); - InputUpdateFromVectorx(femmodel,N_grd,SealevelNGrdEnum,VertexSIdEnum); - InputUpdateFromVectorx(femmodel,bedrock,BedEnum,VertexSIdEnum); - if(RSLg)InputUpdateFromVectorx(femmodel,RSLg,SealevelRSLEnum,VertexSIdEnum); - if(RSLg_barystatic)InputUpdateFromVectorx(femmodel,RSLg_barystatic,SealevelRSLBarystaticEnum,VertexSIdEnum); - if (horiz){ - InputUpdateFromVectorx(femmodel,U_north_grd,SealevelUNorthEsaEnum,VertexSIdEnum); - InputUpdateFromVectorx(femmodel,U_east_grd,SealevelUEastEsaEnum,VertexSIdEnum); - } - - /*reset counter to 1:*/ - femmodel->parameters->SetParam(1,SealevelchangeRunCountEnum); //reset counter. - - /*free ressources:{{{*/ - delete RSLg; - delete RSLg_barystatic; - delete U_grd; - delete N_grd; - delete bedrock; - delete geoid; - if(horiz){ - delete U_north_grd; - delete U_east_grd; - delete bedrockeast; - delete bedrocknorth; - } - delete masks; - /*}}}*/ - -} + femmodel->inputs->AXPY(1,BedEastGRDEnum,BedEastEnum); + femmodel->inputs->AXPY(1,BedNorthGRDEnum, BedNorthEnum); + } + +} /*}}}*/ void dynstr_core(FemModel* femmodel){ /*{{{*/ @@ -397,167 +465,53 @@ } }; /*}}}*/ - -void grd_core_optim(FemModel* femmodel) { /*{{{*/ - - /*variables:{{{*/ - int nel; - BarystaticContributions* barycontrib=NULL; - SealevelMasks* masks=NULL; - GenericParam* barycontribparam=NULL; - IssmDouble rotationaxismotionvector[3]={0}; - - Vector* loads=NULL; - IssmDouble* allloads=NULL; - Vector* sealevelloads=NULL; - Vector* oldsealevelloads=NULL; - IssmDouble* allsealevelloads=NULL; - Vector* oceanareas=NULL; - IssmDouble oceanarea; - IssmDouble oceanaverage; - bool scaleoceanarea=false; - IssmDouble rho_water; - - IssmDouble eps_rel; - IssmDouble eps_abs; - int max_nonlinear_iterations; - int iterations=0; - int step; - IssmDouble time; - bool converged=false; - - int modelid,earthid; - int horiz; - int count,frequency,iscoupling; - int grd=0; - int computesealevel=0; - - /*}}}*/ - - /*Verbose: */ - if(VerboseSolution()) _printf0_(" computing GRD patterns\n"); - - /*retrieve parameters:{{{*/ - femmodel->parameters->FindParam(&grd,SolidearthSettingsGRDEnum); - femmodel->parameters->FindParam(&frequency,SolidearthSettingsRunFrequencyEnum); - femmodel->parameters->FindParam(&count,SealevelchangeRunCountEnum); - femmodel->parameters->FindParam(&computesealevel,SolidearthSettingsComputesealevelchangeEnum); - femmodel->parameters->FindParam(&max_nonlinear_iterations,SolidearthSettingsMaxiterEnum); - /*}}}*/ - - /*only run if grd was requested, if we are the earth, and we have reached - * the necessary number of time steps dictated by :*/ - if(!grd) return; - if(count!=frequency)return; - femmodel->parameters->FindParam(&iscoupling,IsSlcCouplingEnum); - if(iscoupling){ - femmodel->parameters->FindParam(&modelid,ModelIdEnum); - femmodel->parameters->FindParam(&earthid,EarthIdEnum); - if(modelid!=earthid)return; - } - /*retrieve parameters: {{{*/ - femmodel->parameters->FindParam(&scaleoceanarea,SolidearthSettingsOceanAreaScalingEnum); - barycontribparam = xDynamicCast*>(femmodel->parameters->FindParamObject(BarystaticContributionsEnum)); - barycontrib=barycontribparam->GetParameterValue(); - femmodel->parameters->FindParam(&rho_water,MaterialsRhoSeawaterEnum); - femmodel->parameters->FindParam(&eps_rel,SolidearthSettingsReltolEnum); - femmodel->parameters->FindParam(&eps_abs,SolidearthSettingsAbstolEnum); - femmodel->parameters->FindParam(&horiz,SolidearthSettingsHorizEnum); - /*}}}*/ - - /*initialize matrices and vectors:*/ - femmodel->parameters->FindParam(&nel,MeshNumberofelementsEnum); - loads=new Vector(nel); - oceanareas=new Vector(nel); - sealevelloads=new Vector(nel); - sealevelloads->Set(0);sealevelloads->Assemble(); - - /*call masks core: */ - masks=sealevel_masks(femmodel); - if(VerboseSolution()) _printf0_(" starting GRD convolutions\n"); - - /*buildup loads: */ +void ivins_deformation_core(FemModel* femmodel){ /*{{{*/ + + int gsize; + Vector *bedup = NULL; + Vector *beduprate= NULL; + IssmDouble *xx = NULL; + IssmDouble *yy = NULL; + + if(VerboseSolution()) _printf0_(" computing vertical deformation using Ivins model. \n"); + + /*find size of vectors:*/ + gsize = femmodel->nodes->NumberOfDofs(GsetEnum); + + /*Find the litho material to be used by all the elements:*/ + Matlitho* matlitho=NULL; + for (Object* & object: femmodel->materials->objects){ + Material* material=xDynamicCast(object); + if(material->ObjectEnum()==MatlithoEnum){ + matlitho=xDynamicCast(material); + break; + } + } + + /*initialize vectors:*/ + bedup = new Vector(gsize); + beduprate = new Vector(gsize); + + /*retrieve geometric information: */ + VertexCoordinatesx(&xx,&yy,NULL,femmodel->vertices); + + /*Go through elements, and add contribution from each element to the deflection vector wg:*/ for(Object* & object : femmodel->elements->objects){ Element* element = xDynamicCast(object); - element->SealevelchangeBarystaticLoads(loads, barycontrib,masks); - } - loads->Assemble(); - - //broadcast loads - allloads=loads->ToMPISerial(); - - //compute rotation axis motion: - RotationAxisMotion(&rotationaxismotionvector[0],femmodel,allloads,NULL); - - /*skip computation of sea level if requested, which means sea level loads should be zeroed */ - if(!computesealevel){ - allsealevelloads=xNewZeroInit(nel); - goto deformation; - } - - if(VerboseSolution()) _printf0_(" converging GRD convolutions\n"); - for(;;){ - - oldsealevelloads=sealevelloads->Duplicate(); sealevelloads->Copy(oldsealevelloads); - - /*convolve load and sealevel loads on oceans:*/ - for(Object* & object : femmodel->elements->objects){ - Element* element = xDynamicCast(object); - element->SealevelchangeConvolution(sealevelloads, oceanareas, allsealevelloads, allloads,rotationaxismotionvector); - } - sealevelloads->Assemble(); - - /*compute ocean areas:*/ - if(!allsealevelloads){ //first time in the loop - oceanareas->Assemble(); oceanareas->Sum(&oceanarea); _assert_(oceanarea>0.); - if(scaleoceanarea) oceanarea=3.619e+14; // use true ocean area, m^2 - } - - //Conserve ocean mass: - oceanaverage=SealevelloadsOceanAverage(sealevelloads,oceanareas,oceanarea); - ConserveOceanMass(femmodel,sealevelloads,barycontrib->Total()/oceanarea - oceanaverage,masks); - - //broadcast sea level loads - allsealevelloads=sealevelloads->ToMPISerial(); - - //compute rotation axis motion: - RotationAxisMotion(&rotationaxismotionvector[0],femmodel,allloads,allsealevelloads); - - //convergence? - slcconvergence(&converged,sealevelloads,oldsealevelloads,eps_rel,eps_abs); - if (converged)break; - - //early return? - if(iterations>=max_nonlinear_iterations)break; - else iterations++; - } - - deformation: - - if(VerboseSolution()) _printf0_(" deformation GRD convolutions\n"); - - /*convolve loads and sea level loads to get the deformation:*/ - for(Object* & object : femmodel->elements->objects){ - Element* element = xDynamicCast(object); - element->SealevelchangeDeformationConvolution(allsealevelloads, allloads, rotationaxismotionvector); - } - - if(VerboseSolution()) _printf0_(" updating GRD fields\n"); - - /*Update bedrock motion and geoid:*/ - if(computesealevel){ - femmodel->inputs->Shift(SealevelGRDEnum,barycontrib->Total()/rho_water/oceanarea- oceanaverage/rho_water); - - /*cumulate barystatic contributions and save to results: */ - barycontrib->Cumulate(femmodel->parameters); - barycontrib->Save(femmodel->results,femmodel->parameters,oceanarea); - } - - femmodel->inputs->AXPY(1,SealevelGRDEnum,SealevelEnum); + element->GiaDeflection(bedup,beduprate, matlitho, xx,yy); + } + + /*Assemble parallel vector:*/ + beduprate->Assemble(); + bedup->Assemble(); + + /*Save results:*/ + InputUpdateFromVectorx(femmodel,bedup,BedGRDEnum,VertexSIdEnum); femmodel->inputs->AXPY(1,BedGRDEnum,BedEnum); - if(horiz){ - femmodel->inputs->AXPY(1,BedEastGRDEnum,BedEastEnum); - femmodel->inputs->AXPY(1,BedNorthGRDEnum, BedNorthEnum); - } - + + /*Free ressources: */ + xDelete(xx); + xDelete(yy); + delete beduprate; + delete bedup; } /*}}}*/ @@ -612,5 +566,4 @@ femmodel->parameters->FindParam(&geometrydone,SealevelchangeGeometryDoneEnum); femmodel->parameters->FindParam(&grdmodel,GrdModelEnum); - femmodel->parameters->FindParam(&optim,SolidearthSettingsOptimEnum); /*early return?:*/ @@ -635,6 +588,5 @@ for(Object* & object : femmodel->elements->objects){ Element* element=xDynamicCast(object); - if(optim) element->SealevelchangeGeometryOptim(latitude,longitude,radius,xx,yy,zz,xxe,yye,zze,areae); - else element->SealevelchangeGeometry(latitude,longitude,radius,xx,yy,zz,xxe,yye,zze); + element->SealevelchangeGeometry(latitude,longitude,radius,xx,yy,zz,xxe,yye,zze,areae); } @@ -658,300 +610,4 @@ }/*}}}*/ - -//GRD: -Vector* sealevelchange_core_barystatic(FemModel* femmodel,SealevelMasks* masks, IssmDouble* poceanarea){ /*{{{*/ - - /*Barystatic core of the SLR solution (terms that are constant with respect to sea-level)*/ - - Vector *RSLgi = NULL; - IssmDouble RSLgi_oceanaverage = 0; - - /*parameters: */ - int gsize; - IssmDouble oceanarea; - int step; - IssmDouble time; - - /*barystatic contribution:*/ - IssmDouble bslc; - IssmDouble bslcice; - IssmDouble* bslcice_partition=NULL; - IssmDouble bslchydro; - IssmDouble* bslchydro_partition=NULL; - IssmDouble cumbslc; - IssmDouble cumbslcice; - IssmDouble cumbslchydro; - IssmDouble* cumbslcice_partition=NULL; - int npartice; - IssmDouble* cumbslchydro_partition=NULL; - int nparthydro; - int computesealevel=0; - - /*early return if we are not computing sea level, but rather deformation: */ - femmodel->parameters->FindParam(&computesealevel,SolidearthSettingsComputesealevelchangeEnum); - if (!computesealevel)return NULL; - - if(VerboseSolution()) _printf0_(" computing bslc components on ice\n"); - - /*Figure out size of g-set deflection vector and allocate solution vector: */ - gsize = femmodel->nodes->NumberOfDofs(GsetEnum); - - /*some parameters:*/ - femmodel->parameters->FindParam(&step,StepEnum); - femmodel->parameters->FindParam(&time,TimeEnum); - femmodel->parameters->FindParam(&cumbslc,CumBslcEnum); - femmodel->parameters->FindParam(&cumbslcice,CumBslcIceEnum); - femmodel->parameters->FindParam(&cumbslchydro,CumBslcHydroEnum); - femmodel->parameters->FindParam(&npartice,SolidearthNpartIceEnum); - femmodel->parameters->FindParam(&nparthydro,SolidearthNpartHydroEnum); - if(npartice) femmodel->parameters->FindParam(&cumbslcice_partition,&npartice,NULL,CumBslcIcePartitionEnum); - if(nparthydro) femmodel->parameters->FindParam(&cumbslchydro_partition,&nparthydro,NULL,CumBslcHydroPartitionEnum); - - /*Initialize:*/ - RSLgi = new Vector(gsize); - - /*call the bslc main module: */ - femmodel->SealevelchangeBarystatic(RSLgi,&oceanarea,&bslc, &bslcice, &bslchydro, &bslcice_partition, &bslchydro_partition,masks); //this computes - - /*we need to average RSLgi over the ocean: RHS term 4 in Eq.4 of Farrel and clarke. Only the elements can do that: */ - RSLgi_oceanaverage=femmodel->SealevelchangeOceanAverage(RSLgi,masks, oceanarea); - - /*RSLg is the sum of the pure bslc component (term 3) and the contribution from the perturbation to the graviation potential due to the - * presence of ice (terms 1 and 4 in Eq.4 of Farrel and Clarke):*/ - RSLgi->Shift(-bslc-RSLgi_oceanaverage); - - /*save bslc and cumulated bslc value for results:{{{ */ - femmodel->results->AddResult(new GenericExternalResult(femmodel->results->Size()+1,BslcEnum,-bslc,step,time)); - femmodel->results->AddResult(new GenericExternalResult(femmodel->results->Size()+1,BslcIceEnum,-bslcice,step,time)); - femmodel->results->AddResult(new GenericExternalResult(femmodel->results->Size()+1,BslcHydroEnum,-bslchydro,step,time)); - - //cumulative barystatic contribution: - cumbslc=cumbslc-bslc; - cumbslcice=cumbslcice-bslcice; - cumbslchydro=cumbslchydro-bslchydro; - - femmodel->parameters->SetParam(cumbslc,CumBslcEnum); - femmodel->parameters->SetParam(cumbslcice,CumBslcIceEnum); - femmodel->parameters->SetParam(cumbslchydro,CumBslcHydroEnum); - - femmodel->results->AddResult(new GenericExternalResult(femmodel->results->Size()+1,CumBslcEnum,cumbslc,step,time)); - femmodel->results->AddResult(new GenericExternalResult(femmodel->results->Size()+1,CumBslcIceEnum,cumbslcice,step,time)); - femmodel->results->AddResult(new GenericExternalResult(femmodel->results->Size()+1,CumBslcHydroEnum,cumbslchydro,step,time)); - - //cumulative barystatic contributions by partition: - if(npartice){ - for(int i=0;iparameters->SetParam(cumbslcice_partition,npartice,1,CumBslcIcePartitionEnum); - femmodel->results->AddResult(new GenericExternalResult(femmodel->results->Size()+1,CumBslcIcePartitionEnum,cumbslcice_partition,npartice,1,step,time)); - } - - if(nparthydro){ - for(int i=0;iparameters->SetParam(cumbslchydro_partition,nparthydro,1,CumBslcHydroPartitionEnum); - femmodel->results->AddResult(new GenericExternalResult(femmodel->results->Size()+1,CumBslcHydroPartitionEnum,cumbslchydro_partition,nparthydro,1,step,time)); - } - /*}}}*/ - - /*Assign output pointers and return: */ - *poceanarea=oceanarea; - return RSLgi; -}/*}}}*/ -Vector* sealevelchange_core_sal(FemModel* femmodel, SealevelMasks* masks, Vector* RSLg_barystatic,IssmDouble oceanarea){ /*{{{*/ - - /*this core computes the contributions from Eq.4 of Farrel and Clarke, rhs terms 2 and 5. - sal core of the SLR solution */ - - Vector *RSLg = NULL; - Vector *RSLg_old = NULL; - - Vector *RSLgo = NULL; //ocean convolution of the perturbation to gravity potential. - Vector *RSLgo_rot= NULL; // rotational feedback - IssmDouble RSLgo_oceanaverage = 0; //average of RSLgo over the ocean. - - /*parameters: */ - int count; - bool save_results; - int gsize; - bool converged=true; - bool rotation=true; - bool verboseconvolution=true; - int max_nonlinear_iterations; - IssmDouble eps_rel; - IssmDouble eps_abs; - IssmDouble Ixz, Iyz, Izz; - int computesealevel=0; - - /*early return if we are not computing sea level, but rather deformation: */ - femmodel->parameters->FindParam(&computesealevel,SolidearthSettingsComputesealevelchangeEnum); - if (!computesealevel)return NULL; - - if(VerboseSolution()) _printf0_(" converging on ocean components\n"); - - /*Recover some parameters: */ - femmodel->parameters->FindParam(&max_nonlinear_iterations,SolidearthSettingsMaxiterEnum); - femmodel->parameters->FindParam(&eps_rel,SolidearthSettingsReltolEnum); - femmodel->parameters->FindParam(&eps_abs,SolidearthSettingsAbstolEnum); - - /*computational flag: */ - femmodel->parameters->FindParam(&rotation,SolidearthSettingsRotationEnum); - - /*Figure out size of g-set deflection vector and allocate solution vector: */ - gsize = femmodel->nodes->NumberOfDofs(GsetEnum); - - /*Initialize:*/ - RSLg = new Vector(gsize); - RSLg->Assemble(); - RSLg_barystatic->Copy(RSLg); //first initialize RSLg with the barystatic component computed in sealevelchange_core_barystatic. - - RSLg_old = new Vector(gsize); - RSLg_old->Assemble(); - - count=1; - converged=false; - - /*Start loop: */ - for(;;){ - - //save pointer to old sea level - delete RSLg_old; RSLg_old=RSLg; - - /*Initialize solution vector: */ - RSLg = new Vector(gsize); RSLg->Assemble(); - RSLgo = new Vector(gsize); RSLgo->Assemble(); - - /*call the sal module: */ - femmodel->SealevelchangeSal(RSLgo, RSLg_old, masks, verboseconvolution); - - /*assemble solution vector: */ - RSLgo->Assemble(); - - if(rotation){ - - /*call rotational feedback module: */ - RSLgo_rot = new Vector(gsize); RSLgo_rot->Assemble(); - femmodel->SealevelchangeRotationalFeedback(RSLgo_rot,RSLg_old,&Ixz,&Iyz,&Izz, masks); - RSLgo_rot->Assemble(); - - /*save changes in inertia tensor as results: */ - femmodel->results->AddResult(new GenericExternalResult(femmodel->results->Size()+1,SealevelInertiaTensorXZEnum,Ixz)); - femmodel->results->AddResult(new GenericExternalResult(femmodel->results->Size()+1,SealevelInertiaTensorYZEnum,Iyz)); - femmodel->results->AddResult(new GenericExternalResult(femmodel->results->Size()+1,SealevelInertiaTensorZZEnum,Izz)); - - RSLgo->AXPY(RSLgo_rot,1); - } - - /*we need to average RSLgo over the ocean: RHS term 5 in Eq.4 of Farrel and clarke. Only the elements can do that: */ - RSLgo_oceanaverage=femmodel->SealevelchangeOceanAverage(RSLgo,masks, oceanarea); - - /*RSLg is the sum of the barystatic term, and the ocean terms: */ - RSLg_barystatic->Copy(RSLg); RSLg->AXPY(RSLgo,1); - RSLg->Shift(-RSLgo_oceanaverage); - - /*convergence criterion:*/ - slcconvergence(&converged,RSLg,RSLg_old,eps_rel,eps_abs); - - - /*free ressources: */ - delete RSLgo; - delete RSLgo_rot; - - /*Increase count: */ - count++; - if(converged==true){ - break; - } - if(count>=max_nonlinear_iterations){ - _printf0_(" maximum number of nonlinear iterations (" << max_nonlinear_iterations << ") exceeded\n"); - converged=true; - break; - } - - /*some minor verbosing adjustment:*/ - if(count>1)verboseconvolution=false; - - } - if(VerboseConvergence()) _printf0_("\n total number of iterations: " << count-1 << "\n"); - - - delete RSLg_old; - - return RSLg; -} /*}}}*/ -void sealevelchange_core_deformation(Vector** pN_grd, Vector** pU_grd, Vector** pU_north_grd,Vector** pU_east_grd,FemModel* femmodel,Vector* RSLg, SealevelMasks* masks){ /*{{{*/ - - Vector *U_grd = NULL; - Vector *dUdt_grd = NULL; - Vector *N_grd = NULL; - Vector *U_north_grd = NULL; - Vector *U_east_grd = NULL; - - /*parameters: */ - int gsize; - bool spherical=true; - - IssmDouble *latitude = NULL; - IssmDouble *longitude = NULL; - IssmDouble *radius = NULL; - IssmDouble *xx = NULL; - IssmDouble *yy = NULL; - IssmDouble *zz = NULL; - int horiz; - int grdmodel; - - if(VerboseSolution()) _printf0_(" computing vertical and horizontal geodetic signatures\n"); - - /*retrieve some parameters:*/ - femmodel->parameters->FindParam(&horiz,SolidearthSettingsHorizEnum); - femmodel->parameters->FindParam(&grdmodel,GrdModelEnum); - - /*find size of vectors:*/ - gsize = femmodel->nodes->NumberOfDofs(GsetEnum); - - /*intialize vectors:*/ - U_grd = new Vector(gsize); - if(grdmodel==IvinsEnum) dUdt_grd = new Vector(gsize); - N_grd = new Vector(gsize); - if (horiz){ - U_north_grd = new Vector(gsize); - U_east_grd = new Vector(gsize); - } - - /*retrieve geometric information: */ - VertexCoordinatesx(&latitude,&longitude,&radius,femmodel->vertices,spherical); - VertexCoordinatesx(&xx,&yy,&zz,femmodel->vertices); - - /*call the deformation module: */ - switch(grdmodel){ - case NoneEnum: - //do nothing: - break; - case IvinsEnum: - femmodel->IvinsDeformation(U_grd,dUdt_grd,xx,yy); - break; - case ElasticEnum: - femmodel->SealevelchangeDeformation(N_grd, U_grd,U_north_grd,U_east_grd,RSLg, masks); - break; - default: - _error_("Grd model " << EnumToStringx(grdmodel) << " not supported yet"); - } - - /*Assign output pointers:*/ - *pU_grd=U_grd; - *pN_grd=N_grd; - if(horiz){ - *pU_east_grd=U_east_grd; - *pU_north_grd=U_north_grd; - } - - /*Free ressources: */ - xDelete(longitude); - xDelete(latitude); - xDelete(xx); - xDelete(yy); - xDelete(zz); - xDelete(radius); - if(grdmodel==IvinsEnum)delete dUdt_grd; -} -/*}}}*/ /*Support routines:*/ @@ -1284,5 +940,5 @@ for(Object* & object : femmodel->elements->objects){ Element* element = xDynamicCast(object); - element->SealevelchangeMomentOfInertiaElement(&moi_list[0],loads,sealevelloads); + element->SealevelchangeMomentOfInertia(&moi_list[0],loads,sealevelloads); moi_list_cpu[0] += moi_list[0]; moi_list_cpu[1] += moi_list[1]; Index: /issm/trunk-jpl/src/c/shared/Enum/Enum.vim =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 26164) +++ /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 26165) @@ -350,5 +350,4 @@ syn keyword cConstant SolidearthSettingsAbstolEnum syn keyword cConstant SolidearthSettingsCrossSectionShapeEnum -syn keyword cConstant SolidearthSettingsOptimEnum syn keyword cConstant RotationalAngularVelocityEnum syn keyword cConstant SolidearthSettingsElasticEnum @@ -1452,5 +1451,4 @@ syn keyword cType Cfsurfacesquare syn keyword cType Channel -syn keyword cType classes syn keyword cType Constraint syn keyword cType Constraints @@ -1459,6 +1457,6 @@ syn keyword cType ControlInput syn keyword cType Covertree +syn keyword cType DataSetParam syn keyword cType DatasetInput -syn keyword cType DataSetParam syn keyword cType Definition syn keyword cType DependentObject @@ -1473,6 +1471,6 @@ syn keyword cType ElementInput syn keyword cType ElementMatrix +syn keyword cType ElementVector syn keyword cType Elements -syn keyword cType ElementVector syn keyword cType ExponentialVariogram syn keyword cType ExternalResult @@ -1481,10 +1479,9 @@ syn keyword cType Friction syn keyword cType Gauss -syn keyword cType GaussianVariogram -syn keyword cType gaussobjects syn keyword cType GaussPenta syn keyword cType GaussSeg syn keyword cType GaussTetra syn keyword cType GaussTria +syn keyword cType GaussianVariogram syn keyword cType GenericExternalResult syn keyword cType GenericOption @@ -1501,5 +1498,4 @@ syn keyword cType IssmDirectApplicInterface syn keyword cType IssmParallelDirectApplicInterface -syn keyword cType krigingobjects syn keyword cType Load syn keyword cType Loads @@ -1512,5 +1508,4 @@ syn keyword cType Matice syn keyword cType Matlitho -syn keyword cType matrixobjects syn keyword cType MatrixParam syn keyword cType Misfit @@ -1525,6 +1520,6 @@ syn keyword cType Observations syn keyword cType Option +syn keyword cType OptionUtilities syn keyword cType Options -syn keyword cType OptionUtilities syn keyword cType Param syn keyword cType Parameters @@ -1540,11 +1535,11 @@ syn keyword cType Regionaloutput syn keyword cType Results +syn keyword cType RiftStruct syn keyword cType Riftfront -syn keyword cType RiftStruct syn keyword cType SealevelMasks syn keyword cType Seg syn keyword cType SegInput +syn keyword cType SegRef syn keyword cType Segment -syn keyword cType SegRef syn keyword cType SpcDynamic syn keyword cType SpcStatic @@ -1565,4 +1560,8 @@ syn keyword cType Vertex syn keyword cType Vertices +syn keyword cType classes +syn keyword cType gaussobjects +syn keyword cType krigingobjects +syn keyword cType matrixobjects syn keyword cType AdjointBalancethickness2Analysis syn keyword cType AdjointBalancethicknessAnalysis Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 26164) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 26165) @@ -344,5 +344,4 @@ SolidearthSettingsAbstolEnum, SolidearthSettingsCrossSectionShapeEnum, - SolidearthSettingsOptimEnum, RotationalAngularVelocityEnum, SolidearthSettingsElasticEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 26164) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 26165) @@ -352,5 +352,4 @@ case SolidearthSettingsAbstolEnum : return "SolidearthSettingsAbstol"; case SolidearthSettingsCrossSectionShapeEnum : return "SolidearthSettingsCrossSectionShape"; - case SolidearthSettingsOptimEnum : return "SolidearthSettingsOptim"; case RotationalAngularVelocityEnum : return "RotationalAngularVelocity"; case SolidearthSettingsElasticEnum : return "SolidearthSettingsElastic"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 26164) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 26165) @@ -358,5 +358,4 @@ else if (strcmp(name,"SolidearthSettingsAbstol")==0) return SolidearthSettingsAbstolEnum; else if (strcmp(name,"SolidearthSettingsCrossSectionShape")==0) return SolidearthSettingsCrossSectionShapeEnum; - else if (strcmp(name,"SolidearthSettingsOptim")==0) return SolidearthSettingsOptimEnum; else if (strcmp(name,"RotationalAngularVelocity")==0) return RotationalAngularVelocityEnum; else if (strcmp(name,"SolidearthSettingsElastic")==0) return SolidearthSettingsElasticEnum; @@ -383,9 +382,9 @@ else if (strcmp(name,"SolidearthSettingsReltol")==0) return SolidearthSettingsReltolEnum; else if (strcmp(name,"SealevelchangeRequestedOutputs")==0) return SealevelchangeRequestedOutputsEnum; + else if (strcmp(name,"SolidearthSettingsRigid")==0) return SolidearthSettingsRigidEnum; else stage=4; } if(stage==4){ - if (strcmp(name,"SolidearthSettingsRigid")==0) return SolidearthSettingsRigidEnum; - else if (strcmp(name,"SolidearthSettingsRotation")==0) return SolidearthSettingsRotationEnum; + if (strcmp(name,"SolidearthSettingsRotation")==0) return SolidearthSettingsRotationEnum; else if (strcmp(name,"SealevelchangeRunCount")==0) return SealevelchangeRunCountEnum; else if (strcmp(name,"SealevelchangeTransitions")==0) return SealevelchangeTransitionsEnum; @@ -506,9 +505,9 @@ else if (strcmp(name,"TransientIsgroundingline")==0) return TransientIsgroundinglineEnum; else if (strcmp(name,"TransientIshydrology")==0) return TransientIshydrologyEnum; + else if (strcmp(name,"TransientIsmasstransport")==0) return TransientIsmasstransportEnum; else stage=5; } if(stage==5){ - if (strcmp(name,"TransientIsmasstransport")==0) return TransientIsmasstransportEnum; - else if (strcmp(name,"TransientIsoceantransport")==0) return TransientIsoceantransportEnum; + if (strcmp(name,"TransientIsoceantransport")==0) return TransientIsoceantransportEnum; else if (strcmp(name,"TransientIsmovingfront")==0) return TransientIsmovingfrontEnum; else if (strcmp(name,"TransientIsoceancoupling")==0) return TransientIsoceancouplingEnum; @@ -629,9 +628,9 @@ else if (strcmp(name,"EplHead")==0) return EplHeadEnum; else if (strcmp(name,"EplHeadOld")==0) return EplHeadOldEnum; + else if (strcmp(name,"EplHeadSlopeX")==0) return EplHeadSlopeXEnum; else stage=6; } if(stage==6){ - if (strcmp(name,"EplHeadSlopeX")==0) return EplHeadSlopeXEnum; - else if (strcmp(name,"EplHeadSlopeY")==0) return EplHeadSlopeYEnum; + if (strcmp(name,"EplHeadSlopeY")==0) return EplHeadSlopeYEnum; else if (strcmp(name,"EplHeadSubstep")==0) return EplHeadSubstepEnum; else if (strcmp(name,"EplHeadTransient")==0) return EplHeadTransientEnum; @@ -752,9 +751,9 @@ else if (strcmp(name,"RadarIcePeriod")==0) return RadarIcePeriodEnum; else if (strcmp(name,"RadarPowerMacGregor")==0) return RadarPowerMacGregorEnum; + else if (strcmp(name,"RadarPowerWolff")==0) return RadarPowerWolffEnum; else stage=7; } if(stage==7){ - if (strcmp(name,"RadarPowerWolff")==0) return RadarPowerWolffEnum; - else if (strcmp(name,"RheologyBAbsGradient")==0) return RheologyBAbsGradientEnum; + if (strcmp(name,"RheologyBAbsGradient")==0) return RheologyBAbsGradientEnum; else if (strcmp(name,"RheologyBInitialguess")==0) return RheologyBInitialguessEnum; else if (strcmp(name,"RheologyBInitialguessMisfit")==0) return RheologyBInitialguessMisfitEnum; @@ -875,9 +874,9 @@ else if (strcmp(name,"SmbReini")==0) return SmbReiniEnum; else if (strcmp(name,"SmbRunoff")==0) return SmbRunoffEnum; + else if (strcmp(name,"SmbRunoffSubstep")==0) return SmbRunoffSubstepEnum; else stage=8; } if(stage==8){ - if (strcmp(name,"SmbRunoffSubstep")==0) return SmbRunoffSubstepEnum; - else if (strcmp(name,"SmbRunoffTransient")==0) return SmbRunoffTransientEnum; + if (strcmp(name,"SmbRunoffTransient")==0) return SmbRunoffTransientEnum; else if (strcmp(name,"SmbS0gcm")==0) return SmbS0gcmEnum; else if (strcmp(name,"SmbS0p")==0) return SmbS0pEnum; @@ -998,9 +997,9 @@ else if (strcmp(name,"Outputdefinition22")==0) return Outputdefinition22Enum; else if (strcmp(name,"Outputdefinition23")==0) return Outputdefinition23Enum; + else if (strcmp(name,"Outputdefinition24")==0) return Outputdefinition24Enum; else stage=9; } if(stage==9){ - if (strcmp(name,"Outputdefinition24")==0) return Outputdefinition24Enum; - else if (strcmp(name,"Outputdefinition25")==0) return Outputdefinition25Enum; + if (strcmp(name,"Outputdefinition25")==0) return Outputdefinition25Enum; else if (strcmp(name,"Outputdefinition26")==0) return Outputdefinition26Enum; else if (strcmp(name,"Outputdefinition27")==0) return Outputdefinition27Enum; @@ -1121,9 +1120,9 @@ else if (strcmp(name,"CalvingHab")==0) return CalvingHabEnum; else if (strcmp(name,"CalvingLevermann")==0) return CalvingLevermannEnum; + else if (strcmp(name,"CalvingVonmises")==0) return CalvingVonmisesEnum; else stage=10; } if(stage==10){ - if (strcmp(name,"CalvingVonmises")==0) return CalvingVonmisesEnum; - else if (strcmp(name,"Cfdragcoeffabsgrad")==0) return CfdragcoeffabsgradEnum; + if (strcmp(name,"Cfdragcoeffabsgrad")==0) return CfdragcoeffabsgradEnum; else if (strcmp(name,"Cfsurfacelogvel")==0) return CfsurfacelogvelEnum; else if (strcmp(name,"Cfsurfacesquare")==0) return CfsurfacesquareEnum; @@ -1244,9 +1243,9 @@ else if (strcmp(name,"IntParam")==0) return IntParamEnum; else if (strcmp(name,"IntVecParam")==0) return IntVecParamEnum; + else if (strcmp(name,"Inputs")==0) return InputsEnum; else stage=11; } if(stage==11){ - if (strcmp(name,"Inputs")==0) return InputsEnum; - else if (strcmp(name,"Internal")==0) return InternalEnum; + if (strcmp(name,"Internal")==0) return InternalEnum; else if (strcmp(name,"Intersect")==0) return IntersectEnum; else if (strcmp(name,"InversionVzObs")==0) return InversionVzObsEnum; @@ -1367,9 +1366,9 @@ else if (strcmp(name,"SMBd18opdd")==0) return SMBd18opddEnum; else if (strcmp(name,"SMBforcing")==0) return SMBforcingEnum; + else if (strcmp(name,"SMBgcm")==0) return SMBgcmEnum; else stage=12; } if(stage==12){ - if (strcmp(name,"SMBgcm")==0) return SMBgcmEnum; - else if (strcmp(name,"SMBgemb")==0) return SMBgembEnum; + if (strcmp(name,"SMBgemb")==0) return SMBgembEnum; else if (strcmp(name,"SMBgradients")==0) return SMBgradientsEnum; else if (strcmp(name,"SMBgradientscomponents")==0) return SMBgradientscomponentsEnum; Index: /issm/trunk-jpl/src/m/classes/solidearthsettings.m =================================================================== --- /issm/trunk-jpl/src/m/classes/solidearthsettings.m (revision 26164) +++ /issm/trunk-jpl/src/m/classes/solidearthsettings.m (revision 26165) @@ -22,5 +22,4 @@ grdmodel = 0; %grd model (0 by default, 1 for elastic, 2 for Ivins) cross_section_shape = 0; %cross section only used when grd model is Ivins - optim = 0; %new optimized version of the GRD code. end methods @@ -66,7 +65,4 @@ self.cross_section_shape=1; %square as default (see iedge in GiaDeflectionCorex) - %optim? - self.optim=0; - %no grd model by default: self.grdmodel=0; @@ -87,5 +83,4 @@ md = checkfield(md,'fieldname','solidearth.settings.grdmodel','values',[1 2]); md = checkfield(md,'fieldname','solidearth.settings.cross_section_shape','numel',[1],'values',[1,2]); - md = checkfield(md,'fieldname','solidearth.settings.optim','values',[0,1]); %checks on computational flags @@ -131,5 +126,4 @@ fielddisplay(self,'grdmodel','type of deformation model, 1 for elastic, 2 for visco-elastic from Ivins'); fielddisplay(self,'cross_section_shape','1: square-edged (default). 2: elliptical. See iedge in GiaDeflectionCore'); - fielddisplay(self,'optim','use optimized version of the GRD code? (default 0)'); end % }}} function marshall(self,prefix,md,fid) % {{{ @@ -150,5 +144,4 @@ WriteData(fid,prefix,'object',self,'fieldname','grdmodel','name','md.solidearth.settings.grdmodel','format','Integer'); WriteData(fid,prefix,'object',self,'fieldname','cross_section_shape','name','md.solidearth.settings.cross_section_shape','format','Integer'); - WriteData(fid,prefix,'object',self,'fieldname','optim','name','md.solidearth.settings.optim','format','Integer'); end % }}} function savemodeljs(self,fid,modelname) % {{{ @@ -165,5 +158,4 @@ writejsdouble(fid,[modelname '.solidearth.settings.glfraction'],self.glfraction); writejsdouble(fid,[modelname '.solidearth.settings.cross_section_shape'],self.cross_section_shape); - writejsdouble(fid,[modelname '.solidearth.settings.optim'],self.optim); end % }}} function self = extrude(self,md) % {{{ Index: /issm/trunk-jpl/test/NightlyRun/test2001.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test2001.m (revision 26164) +++ /issm/trunk-jpl/test/NightlyRun/test2001.m (revision 26165) @@ -18,5 +18,5 @@ md.solidearth.settings.cross_section_shape=1; % for square-edged x-section md.solidearth.settings.computesealevelchange=0; %do not compute sea level, only deformation -md.solidearth.requested_outputs={'Sealevel','SealevelUGrd'}; +md.solidearth.requested_outputs={'Sealevel','BedGRD'}; %Loading history @@ -50,3 +50,3 @@ field_names ={'UGrd'}; field_tolerances={1e-13}; -field_values={md.results.TransientSolution(2).SealevelUGrd}; +field_values={md.results.TransientSolution(2).BedGRD}; Index: /issm/trunk-jpl/test/NightlyRun/test2002.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test2002.m (revision 26164) +++ /issm/trunk-jpl/test/NightlyRun/test2002.m (revision 26165) @@ -3,5 +3,5 @@ %mesh earth: md=model; -md.mesh=gmshplanet('radius',6.371012*10^3,'resolution',700.); %700 km resolution mesh +md.mesh=gmshplanet('radius',6.371012*10^3,'resolution',500.); %700 km resolution mesh %Geometry for the bed, arbitrary thickness of 1000: @@ -24,7 +24,9 @@ pos=find(late < -80); md.masstransport.spcthickness(md.mesh.elements(pos,:))= md.masstransport.spcthickness(md.mesh.elements(pos,:))-100; +posant=pos; %greenland pos=find(late>70 & late<80 & longe>-60 & longe<-30); md.masstransport.spcthickness(md.mesh.elements(pos,:))= md.masstransport.spcthickness(md.mesh.elements(pos,:))-100; +posgre=pos; %elastic loading from love numbers: @@ -34,16 +36,13 @@ %mask: {{{ mask=gmtmask(md.mesh.lat,md.mesh.long); +oceanmask=-ones(md.mesh.numberofvertices,1); +pos=find(mask==0); oceanmask(pos)=1; + icemask=ones(md.mesh.numberofvertices,1); -pos=find(mask==0); -icemask(pos)=-1; -pos=find(sum(mask(md.mesh.elements),2)<3); -icemask(md.mesh.elements(pos,:))=-1; +icemask(md.mesh.elements(posant))=-1; +icemask(md.mesh.elements(posgre))=-1; + md.mask.ice_levelset=icemask; -md.mask.ocean_levelset=-icemask; - - -%make sure wherever there is an ice load, that the mask is set to ice: -pos=find(md.masstransport.spcthickness(1:end-1)); -md.mask.ice_levelset(md.mesh.elements(pos,:))=-1; +md.mask.ocean_levelset=oceanmask; % }}} @@ -66,9 +65,9 @@ md.materials=materials('hydro'); - %Miscellaneous md.miscellaneous.name='test2002'; %Solution parameters +md.cluster.np=3; md.solidearth.settings.reltol=NaN; md.solidearth.settings.abstol=1e-3; @@ -83,5 +82,4 @@ md.transient.isthermal=0; md.transient.ismasstransport=1; -md.transient.isoceantransport=1; md.transient.isslc=1; md.solidearth.requested_outputs={'Sealevel'};