Index: /issm/trunk-jpl/src/c/Makefile.am =================================================================== --- /issm/trunk-jpl/src/c/Makefile.am (revision 26221) +++ /issm/trunk-jpl/src/c/Makefile.am (revision 26222) @@ -249,4 +249,5 @@ ./cores/ResetBoundaryConditions.cpp \ ./cores/WrapperCorePointerFromSolutionEnum.cpp \ + ./cores/WrapperPreCorePointerFromSolutionEnum.cpp \ ./cores/CorePointerFromSolutionEnum.cpp \ ./cores/ad_core.cpp \ @@ -557,5 +558,6 @@ issm_sources += \ ./cores/sealevelchange_core.cpp \ - ./analyses/SealevelchangeAnalysis.cpp + ./analyses/SealevelchangeAnalysis.cpp\ + ./classes/SealevelGeometry.cpp #gia ivins physics (only if have fortran) Index: /issm/trunk-jpl/src/c/analyses/SealevelchangeAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/SealevelchangeAnalysis.cpp (revision 26221) +++ /issm/trunk-jpl/src/c/analyses/SealevelchangeAnalysis.cpp (revision 26222) @@ -143,5 +143,4 @@ parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.compute_bp_grd",SolidearthSettingsComputeBpGrdEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.planetradius",SolidearthPlanetRadiusEnum)); - parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.glfraction",SolidearthSettingsGlfractionEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.cross_section_shape",SolidearthSettingsCrossSectionShapeEnum)); parameters->AddObject(new DoubleParam(CumBslcEnum,0.0)); @@ -429,5 +428,5 @@ /*run sea level change core geometry only once, after the Model Processor is done:*/ - sealevelchange_geometry(femmodel); + sealevelchange_initialgeometry(femmodel); }/*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Dakota/IssmParallelDirectApplicInterface.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Dakota/IssmParallelDirectApplicInterface.cpp (revision 26221) +++ /issm/trunk-jpl/src/c/classes/Dakota/IssmParallelDirectApplicInterface.cpp (revision 26222) @@ -55,4 +55,5 @@ bool control_analysis = false; void (*solutioncore)(FemModel*) = NULL; + void (*solutionprecore)(FemModel*) = NULL; bool nodakotacore = true; @@ -94,6 +95,7 @@ responses=xNewZeroInit(numFns); - /*Hack:*/ - for(int i=0;inummodels;i++) if(femmodel_init->analysis_type_list[i]==SealevelchangeAnalysisEnum) sealevelchange_geometry(femmodel_init); + /*Launch cores that are not used during the uncertainty quantification: */ + WrapperPreCorePointerFromSolutionEnum(&solutionprecore,femmodel_init->parameters,solution_type); + if(solutionprecore)solutionprecore(femmodel_init); /*Make a copy of femmodel, so we start this new evaluation run for this specific sample with a brand Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 26221) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 26222) @@ -2034,5 +2034,5 @@ bool Element::IsIceOnlyInElement(){/*{{{*/ Input* input=this->GetInput(MaskIceLevelsetEnum); _assert_(input); - return (input->GetInputMax()<0.); + return (input->GetInputMax()<=0.); } /*}}}*/ @@ -2045,4 +2045,9 @@ Input* input=this->GetInput(MaskOceanLevelsetEnum); _assert_(input); return (input->GetInputMin()<0.); +} +/*}}}*/ +bool Element::IsOceanOnlyInElement(){/*{{{*/ + Input* input=this->GetInput(MaskOceanLevelsetEnum); _assert_(input); + return (input->GetInputMax()<=0.); } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 26221) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 26222) @@ -32,4 +32,5 @@ class IoModel; class SealevelMasks; +class SealevelGeometry; class Gauss; class ElementVector; @@ -146,4 +147,5 @@ bool IsIceOnlyInElement(); bool IsOceanInElement(); + bool IsOceanOnlyInElement(); bool IsLandInElement(); void Ismip6FloatingiceMeltingRate(); @@ -240,4 +242,5 @@ virtual void ElementSizes(IssmDouble* phx,IssmDouble* phy,IssmDouble* phz)=0; virtual void ElementCoordinates(Vector* vxe,Vector* vye,Vector* vze, Vector* vareae, bool spherical=false)=0; + virtual void ElementCoordinates(Vector* vlonge,Vector* vlate,Vector* vareae)=0; virtual int FiniteElement(void)=0; virtual IssmDouble FloatingArea(bool scaled)=0; @@ -248,6 +251,9 @@ virtual void GetGroundedPart(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* mainlyfloating)=0; virtual IssmDouble GetGroundedPortion(IssmDouble* xyz_list)=0; - virtual void GetFractionGeometry(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* pmainlynegative, IssmDouble* gl)=0; - virtual IssmDouble GetFractionArea(IssmDouble* xyz_list, IssmDouble* gl)=0; + virtual void GetFractionGeometry(IssmDouble* weights, IssmDouble* pphi, int* ppoint1,IssmDouble* pfraction1,IssmDouble* pfraction2, bool* ptrapezeisnegative, IssmDouble* gl)=0; + virtual void GetNodalWeightsAndAreaAndCentroidsFromLeveset(IssmDouble* loadweights, IssmDouble* ploadarea, IssmDouble* platbar, IssmDouble* plongbar, IssmDouble late, IssmDouble longe, IssmDouble area, int levelsetenum)=0; + virtual void GetNodalWeightsAndAreaAndCentroidsFromLeveset(IssmDouble* loadweights, IssmDouble* ploadarea, IssmDouble* platbar, IssmDouble* plongbar, IssmDouble late, IssmDouble longe, IssmDouble area, int levelset1enum, int levelset2enum)=0; + virtual void GetBarycenterFromLevelset(IssmDouble* platbar, IssmDouble* plongbar,IssmDouble phi,IssmDouble fraction1,IssmDouble fraction2,IssmDouble late, IssmDouble longe, int point1,int istrapeze1, IssmDouble planetradius)=0; + virtual IssmDouble GetIcefrontArea(){_error_("not implemented");}; virtual void GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum)=0; @@ -379,16 +385,19 @@ #endif #ifdef _HAVE_SEALEVELCHANGE_ - virtual void LevelsetAverage(IssmDouble* paverage, IssmDouble* parea, IssmDouble* field_on_localvertices, int levelsetenum){_error_("not implemented");}; - virtual void SetSealevelMasks(SealevelMasks* masks)=0; virtual IssmDouble GetArea3D(void)=0; virtual IssmDouble GetAreaSpherical(void)=0; - virtual void SealevelchangeMomentOfInertia(IssmDouble* dI_list, IssmDouble* loads, IssmDouble* sealevelloads)=0; + virtual IssmDouble GetTriangleAreaSpherical(IssmDouble xyz_list[3][3])=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; + virtual void SealevelchangeGeometryFractionKernel(SealevelGeometry* slgeom)=0; + virtual void SealevelchangeMomentOfInertiaCentroid(IssmDouble* dI_list, IssmDouble* loads, IssmDouble* sealevelloads, SealevelGeometry* slgeom)=0; + virtual void SealevelchangeShift(Vector* loads, Vector* subelementloads, IssmDouble offset, SealevelGeometry* slgeom)=0; + virtual void SealevelchangeGeometryInitial(IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae)=0; + virtual void SealevelchangeGeometryCentroidLoads(SealevelGeometry* slgeom, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae)=0; + virtual void SealevelchangeGeometrySubElementLoads(SealevelGeometry* slgeom, IssmDouble* areae)=0; + virtual void SealevelchangeBarystaticLoads(Vector* loads, Vector* subelementiceloads, Vector* subelementhydroloads, Vector* subelementbploads, BarystaticContributions* barycontrib, SealevelGeometry* slgeom)=0; + virtual void SealevelchangeConvolution(Vector* sealevelloads, Vector* subelementsealevelloads, Vector* oceanareas, Vector* subelementoceanareas, IssmDouble* allsealevelloads, IssmDouble* allloads,IssmDouble* allsubelementiceloads, IssmDouble* allsubelementhydroloads, IssmDouble* allsubelementbploads, IssmDouble* allsubelementsealevelloads, IssmDouble* rotationvector,SealevelGeometry* slgeom)=0; + virtual void SealevelchangeDeformationConvolution(IssmDouble* sealevelloads, IssmDouble* subelementsealevelloads, IssmDouble* loads, IssmDouble* subelementiceloads, IssmDouble* subelementhydroloads, IssmDouble* subelementbploads, IssmDouble* rotationvector,SealevelGeometry* slgeom)=0; + virtual void SealevelchangeMomentOfInertiaSubElement(IssmDouble* dI_list, IssmDouble* subelementiceloads, IssmDouble* subelementhydroloads, IssmDouble* subelementbploads, IssmDouble* subelementsealevelloads, SealevelGeometry* slgeom)=0; #endif Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 26221) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 26222) @@ -71,4 +71,5 @@ void ElementResponse(IssmDouble* presponse,int response_enum); void ElementCoordinates(Vector* vxe,Vector* vye,Vector* vze,Vector* vareae,bool spherical=false){_error_("not implemented yet");}; + void ElementCoordinates(Vector* vlonge,Vector* vlate,Vector* vareae){_error_("not implemented yet");}; void ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz); int FiniteElement(void); @@ -78,4 +79,5 @@ IssmDouble GetAreaSpherical(void){_error_("not implemented yet!");}; void GetAreaCoordinates(IssmDouble *area_coordinates,IssmDouble* xyz_zero,IssmDouble* xyz_list,int numpoints); + IssmDouble GetTriangleAreaSpherical(IssmDouble xyz_list[3][3]){_error_("not implemented yet");}; Element* GetBasalElement(void); Penta* GetBasalPenta(void); @@ -83,6 +85,8 @@ void GetGroundedPart(int* point1,IssmDouble* fraction1, IssmDouble* fraction2,bool* mainlyfloating); IssmDouble GetGroundedPortion(IssmDouble* xyz_list); - void GetFractionGeometry(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* pmainlynegative, IssmDouble* gl){_error_("not implemented yet");}; - IssmDouble GetFractionArea(IssmDouble* xyz_list, IssmDouble* gl){_error_("not implemented yet");}; + void GetFractionGeometry(IssmDouble* weights, IssmDouble* pphi, int* ppoint1,IssmDouble* pfraction1,IssmDouble* pfraction2, bool* ptrapezeisnegative, IssmDouble* gl){_error_("not implemented yet");}; + void GetNodalWeightsAndAreaAndCentroidsFromLeveset(IssmDouble* loadweights, IssmDouble* ploadarea, IssmDouble* platbar, IssmDouble* plongbar, IssmDouble late, IssmDouble longe, IssmDouble area, int levelsetenum){_error_("not implemented yet");}; + void GetNodalWeightsAndAreaAndCentroidsFromLeveset(IssmDouble* loadweights, IssmDouble* ploadarea, IssmDouble* platbar, IssmDouble* plongbar, IssmDouble late, IssmDouble longe, IssmDouble area, int levelset1enum, int levelset2enum){_error_("not implemented yet");}; + void GetBarycenterFromLevelset(IssmDouble* platbar, IssmDouble* plongbar,IssmDouble phi,IssmDouble fraction1,IssmDouble fraction2,IssmDouble late, IssmDouble longe, int point1,int istrapeze1, IssmDouble planetradius){_error_("not implemented yet");}; IssmDouble GetIcefrontArea(); void GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum); @@ -218,13 +222,15 @@ #endif #ifdef _HAVE_SEALEVELCHANGE_ - 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");}; - void SealevelchangeDeformationConvolution(IssmDouble* sealevelloads, IssmDouble* loads, IssmDouble* rotationvector){_error_("not implemented yet");}; - void SealevelchangeShift(Vector* loads, IssmDouble offset, SealevelMasks* masks){_error_("not implemented yet");}; + void SealevelchangeGeometryFractionKernel(SealevelGeometry* slgeom){_error_("not implemented yet");}; + void SealevelchangeMomentOfInertiaCentroid(IssmDouble* dI_list, IssmDouble* loads, IssmDouble* sealevelloads, SealevelGeometry* slgeom){_error_("not implemented yet");}; + void SealevelchangeShift(Vector* loads, Vector* subelementloads, IssmDouble offset, SealevelGeometry* slgeom){_error_("not implemented yet");}; + void SealevelchangeGeometryInitial(IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae){_error_("not implemented yet");}; + void SealevelchangeGeometryCentroidLoads(SealevelGeometry* slgeom, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae){_error_("not implemented yet");}; + void SealevelchangeGeometrySubElementLoads(SealevelGeometry* slgeom, IssmDouble* areae){_error_("not implemented yet");}; + void SealevelchangeBarystaticLoads(Vector* loads, Vector* subelementiceloads, Vector* subelementhydroloads, Vector* subelementbploads, BarystaticContributions* barycontrib, SealevelGeometry* slgeom){_error_("not implemented yet");}; + void SealevelchangeConvolution(Vector* sealevelloads, Vector* subelementsealevelloads, Vector* oceanareas, Vector* subelementoceanareas, IssmDouble* allsealevelloads, IssmDouble* allloads,IssmDouble* allsubelementiceloads, IssmDouble* allsubelementhydroloads, IssmDouble* allsubelementbploads, IssmDouble* allsubelementsealevelloads, IssmDouble* rotationvector,SealevelGeometry* slgeom){_error_("not implemented yet");}; + void SealevelchangeDeformationConvolution(IssmDouble* sealevelloads, IssmDouble* subelementsealevelloads, IssmDouble* loads, IssmDouble* subelementiceloads, IssmDouble* subelementhydroloads, IssmDouble* subelementbploads, IssmDouble* rotationvector,SealevelGeometry* slgeom){_error_("not implemented yet");}; + void SealevelchangeMomentOfInertiaSubElement(IssmDouble* dI_list, IssmDouble* subelementiceloads, IssmDouble* subelementhydroloads, IssmDouble* subelementbploads, IssmDouble* subelementsealevelloads, SealevelGeometry* slgeom){_error_("not implemented yet");}; #endif Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 26221) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 26222) @@ -55,4 +55,5 @@ void ElementResponse(IssmDouble* presponse,int response_enum){_error_("not implemented yet");}; void ElementCoordinates(Vector* vxe,Vector* vye,Vector* vze,Vector* vareae,bool spherical=false){_error_("not implemented yet");}; + void ElementCoordinates(Vector* vlonge,Vector* vlate,Vector* vareae){_error_("not implemented yet");}; void ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz){_error_("not implemented yet");}; int FiniteElement(void); @@ -63,7 +64,9 @@ void GetGroundedPart(int* point1,IssmDouble* fraction1, IssmDouble* fraction2,bool* mainlyfloating){_error_("not implemented yet");}; IssmDouble GetGroundedPortion(IssmDouble* xyz_list); - void GetFractionGeometry(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* pmainlynegative, IssmDouble* gl){_error_("not implemented yet");}; - IssmDouble GetFractionArea(IssmDouble* xyz_list, IssmDouble* gl){_error_("not implemented yet");}; - void GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum); + void GetFractionGeometry(IssmDouble* weights, IssmDouble* pphi, int* ppoint1,IssmDouble* pfraction1,IssmDouble* pfraction2, bool* ptrapezeisnegative, IssmDouble* gl){_error_("not implemented yet");}; + void GetNodalWeightsAndAreaAndCentroidsFromLeveset(IssmDouble* loadweights, IssmDouble* ploadarea, IssmDouble* platbar, IssmDouble* plongbar, IssmDouble late, IssmDouble longe, IssmDouble area, int levelsetenum){_error_("not implemented yet");}; + void GetNodalWeightsAndAreaAndCentroidsFromLeveset(IssmDouble* loadweights, IssmDouble* ploadarea, IssmDouble* platbar, IssmDouble* plongbar, IssmDouble late, IssmDouble longe, IssmDouble area, int levelset1enum, int levelset2enum){_error_("not implemented yet");}; + void GetBarycenterFromLevelset(IssmDouble* platbar, IssmDouble* plongbar,IssmDouble phi,IssmDouble fraction1,IssmDouble fraction2,IssmDouble late, IssmDouble longe, int point1,int istrapeze1, IssmDouble planetradius){_error_("not implemented yet");}; + void GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum); Input* GetInput(int enumtype); Input* GetInput(int enumtype,IssmDouble time); @@ -167,4 +170,5 @@ IssmDouble GetArea3D(void){_error_("not implemented yet!");}; IssmDouble GetAreaSpherical(void){_error_("not implemented yet!");}; + IssmDouble GetTriangleAreaSpherical(IssmDouble xyz_list[3][3]){_error_("not implemented yet");}; #ifdef _HAVE_ESA_ @@ -173,13 +177,15 @@ #endif #ifdef _HAVE_SEALEVELCHANGE_ - 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");}; - void SealevelchangeDeformationConvolution(IssmDouble* sealevelloads, IssmDouble* loads, IssmDouble* rotationvector){_error_("not implemented yet");}; - void SealevelchangeShift(Vector* loads, IssmDouble offset, SealevelMasks* masks){_error_("not implemented yet");}; - + void SealevelchangeGeometryFractionKernel(SealevelGeometry* slgeom){_error_("not implemented yet");}; + void SealevelchangeMomentOfInertiaCentroid(IssmDouble* dI_list, IssmDouble* loads, IssmDouble* sealevelloads, SealevelGeometry* slgeom){_error_("not implemented yet");}; + void SealevelchangeShift(Vector* loads, Vector* subelementloads, IssmDouble offset, SealevelGeometry* slgeom){_error_("not implemented yet");}; + void SealevelchangeGeometryInitial(IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae){_error_("not implemented yet");}; + void SealevelchangeGeometryCentroidLoads(SealevelGeometry* slgeom, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae){_error_("not implemented yet");}; + void SealevelchangeGeometrySubElementLoads(SealevelGeometry* slgeom, IssmDouble* areae){_error_("not implemented yet");}; + void SealevelchangeBarystaticLoads(Vector* loads, Vector* subelementiceloads, Vector* subelementhydroloads, Vector* subelementbploads, BarystaticContributions* barycontrib, SealevelGeometry* slgeom){_error_("not implemented yet");}; + void SealevelchangeConvolution(Vector* sealevelloads, Vector* subelementsealevelloads, Vector* oceanareas, Vector* subelementoceanareas, IssmDouble* allsealevelloads, IssmDouble* allloads,IssmDouble* allsubelementiceloads, IssmDouble* allsubelementhydroloads, IssmDouble* allsubelementbploads, IssmDouble* allsubelementsealevelloads, IssmDouble* rotationvector,SealevelGeometry* slgeom){_error_("not implemented yet");}; + void SealevelchangeDeformationConvolution(IssmDouble* sealevelloads, IssmDouble* subelementsealevelloads, IssmDouble* loads, IssmDouble* subelementiceloads, IssmDouble* subelementhydroloads, IssmDouble* subelementbploads, IssmDouble* rotationvector,SealevelGeometry* slgeom){_error_("not implemented yet");}; + void SealevelchangeMomentOfInertiaSubElement(IssmDouble* dI_list, IssmDouble* subelementiceloads, IssmDouble* subelementhydroloads, IssmDouble* subelementbploads, IssmDouble* subelementsealevelloads, SealevelGeometry* slgeom){_error_("not implemented yet");}; #endif Index: /issm/trunk-jpl/src/c/classes/Elements/Tetra.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 26221) +++ /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 26222) @@ -54,4 +54,5 @@ void ElementResponse(IssmDouble* presponse,int response_enum){_error_("not implemented yet");}; void ElementCoordinates(Vector* vxe,Vector* vye,Vector* vze,Vector* vareae,bool spherical=false){_error_("not implemented yet");}; + void ElementCoordinates(Vector* vlonge,Vector* vlate,Vector* vareae){_error_("not implemented yet");}; void ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz); void FaceOnBaseIndices(int* pindex1,int* pindex2,int* pindex3); @@ -63,10 +64,13 @@ IssmDouble GetArea3D(void){_error_("not implemented yet!");}; IssmDouble GetAreaSpherical(void){_error_("not implemented yet!");}; + IssmDouble GetTriangleAreaSpherical(IssmDouble xyz_list[3][3]){_error_("not implemented yet");}; Element* GetBasalElement(void){_error_("not implemented yet");}; int GetElementType(void); void GetGroundedPart(int* point1,IssmDouble* fraction1, IssmDouble* fraction2,bool* mainlyfloating){_error_("not implemented yet");}; IssmDouble GetGroundedPortion(IssmDouble* xyz_list){_error_("not implemented yet");}; - void GetFractionGeometry(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* pmainlynegative, IssmDouble* gl){_error_("not implemented yet");}; - IssmDouble GetFractionArea(IssmDouble* xyz_list, IssmDouble* gl){_error_("not implemented yet");}; + void GetFractionGeometry(IssmDouble* weights, IssmDouble* pphi, int* ppoint1,IssmDouble* pfraction1,IssmDouble* pfraction2, bool* ptrapezeisnegative, IssmDouble* gl){_error_("not implemented yet");}; + void GetNodalWeightsAndAreaAndCentroidsFromLeveset(IssmDouble* loadweights, IssmDouble* ploadarea, IssmDouble* platbar, IssmDouble* plongbar, IssmDouble late, IssmDouble longe, IssmDouble area, int levelsetenum){_error_("not implemented yet");}; + void GetNodalWeightsAndAreaAndCentroidsFromLeveset(IssmDouble* loadweights, IssmDouble* ploadarea, IssmDouble* platbar, IssmDouble* plongbar, IssmDouble late, IssmDouble longe, IssmDouble area, int levelset1enum, int levelset2enum){_error_("not implemented yet");}; + void GetBarycenterFromLevelset(IssmDouble* platbar, IssmDouble* plongbar,IssmDouble phi,IssmDouble fraction1,IssmDouble fraction2,IssmDouble late, IssmDouble longe, int point1,int istrapeze1, IssmDouble planetradius){_error_("not implemented yet");}; void GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum){_error_("not implemented yet");}; Input* GetInput(int enumtype); @@ -180,13 +184,15 @@ #endif #ifdef _HAVE_SEALEVELCHANGE_ - 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");}; - void SealevelchangeDeformationConvolution(IssmDouble* sealevelloads, IssmDouble* loads, IssmDouble* rotationvector){_error_("not implemented yet");}; - void SealevelchangeShift(Vector* loads, IssmDouble offset, SealevelMasks* masks){_error_("not implemented yet");}; - + void SealevelchangeGeometryFractionKernel(SealevelGeometry* slgeom){_error_("not implemented yet");}; + void SealevelchangeMomentOfInertiaCentroid(IssmDouble* dI_list, IssmDouble* loads, IssmDouble* sealevelloads, SealevelGeometry* slgeom){_error_("not implemented yet");}; + void SealevelchangeShift(Vector* loads, Vector* subelementloads, IssmDouble offset, SealevelGeometry* slgeom){_error_("not implemented yet");}; + void SealevelchangeGeometryInitial(IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae){_error_("not implemented yet");}; + void SealevelchangeGeometryCentroidLoads(SealevelGeometry* slgeom, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae){_error_("not implemented yet");}; + void SealevelchangeGeometrySubElementLoads(SealevelGeometry* slgeom, IssmDouble* areae){_error_("not implemented yet");}; + void SealevelchangeBarystaticLoads(Vector* loads, Vector* subelementiceloads, Vector* subelementhydroloads, Vector* subelementbploads, BarystaticContributions* barycontrib, SealevelGeometry* slgeom){_error_("not implemented yet");}; + void SealevelchangeConvolution(Vector* sealevelloads, Vector* subelementsealevelloads, Vector* oceanareas, Vector* subelementoceanareas, IssmDouble* allsealevelloads, IssmDouble* allloads,IssmDouble* allsubelementiceloads, IssmDouble* allsubelementhydroloads, IssmDouble* allsubelementbploads, IssmDouble* allsubelementsealevelloads, IssmDouble* rotationvector,SealevelGeometry* slgeom){_error_("not implemented yet");}; + void SealevelchangeDeformationConvolution(IssmDouble* sealevelloads, IssmDouble* subelementsealevelloads, IssmDouble* loads, IssmDouble* subelementiceloads, IssmDouble* subelementhydroloads, IssmDouble* subelementbploads, IssmDouble* rotationvector,SealevelGeometry* slgeom){_error_("not implemented yet");}; + void SealevelchangeMomentOfInertiaSubElement(IssmDouble* dI_list, IssmDouble* subelementiceloads, IssmDouble* subelementhydroloads, IssmDouble* subelementbploads, IssmDouble* subelementsealevelloads, SealevelGeometry* slgeom){_error_("not implemented yet");}; #endif Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 26221) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 26222) @@ -1312,4 +1312,29 @@ } /*}}}*/ +void Tria::ElementCoordinates(Vector* vlonge,Vector* vlate,Vector* vareae){ /*{{{*/ + + IssmDouble planetradius; + + /*Look for x,y,z coordinates:*/ + IssmDouble xyz_list[NUMVERTICES][3]; + ::GetVerticesCoordinates(&xyz_list[0][0],this->vertices,NUMVERTICES); + IssmDouble area,xe,ye,ze,late,longe; + + /*Find centroid:*/ + xe=(xyz_list[0][0]+xyz_list[1][0]+xyz_list[2][0])/3.0; + ye=(xyz_list[0][1]+xyz_list[1][1]+xyz_list[2][1])/3.0; + ze=(xyz_list[0][2]+xyz_list[1][2]+xyz_list[2][2])/3.0; + + late= asin(ze/sqrt(pow(xe,2.0)+pow(ye,2.0)+pow(ze,2.0))); + longe= atan2(ye,xe); + area=this->GetAreaSpherical(); + + vlonge->SetValue(this->sid,longe,INS_VAL); + vlate->SetValue(this->sid,late,INS_VAL); + vareae->SetValue(this->sid,area,INS_VAL); + + return; +} +/*}}}*/ void Tria::ElementResponse(IssmDouble* presponse,int response_enum){/*{{{*/ @@ -1723,78 +1748,20 @@ } /*}}}*/ -IssmDouble Tria::GetFractionArea(IssmDouble* xyz_list, IssmDouble* gl){/*{{{*/ - - /*Computeportion of the element that is on a levelset: */ - bool mainlyfloating = true; - int domaintype,index1,index2; - const IssmPDouble epsilon = 1.e-15; - IssmDouble phi,s1,s2; - - /*Recover parameters and values*/ - parameters->FindParam(&domaintype,DomainTypeEnum); - - /*Be sure that values are not zero*/ - if(gl[0]==0.) gl[0]=gl[0]+epsilon; - if(gl[1]==0.) gl[1]=gl[1]+epsilon; - if(gl[2]==0.) gl[2]=gl[2]+epsilon; - - if(domaintype==Domain2DverticalEnum){ - this->EdgeOnBaseIndices(&index1,&index2); - if(gl[index1]>0 && gl[index2]>0) phi=1; // All grounded - else if(gl[index1]<0 && gl[index2]<0) phi=0; // All floating - else if(gl[index1]<0 && gl[index2]>0){ //index2 grounded - phi=1./(1.-gl[index1]/gl[index2]); - } - else if(gl[index2]<0 && gl[index1]>0){ //index1 grounded - phi=1./(1.-gl[index2]/gl[index1]); - } - - } - else if(domaintype==Domain2DhorizontalEnum || domaintype==Domain3DEnum || domaintype==Domain3DsurfaceEnum){ - /*Check that not all nodes are grounded or floating*/ - if(gl[0]>0 && gl[1]>0 && gl[2]>0){ // All grounded - phi=1; - } - else if(gl[0]<0 && gl[1]<0 && gl[2]<0){ //All floating - phi=0; - } - else{ - /*Figure out if two nodes are floating or grounded*/ - if(gl[0]*gl[1]*gl[2]>0) mainlyfloating=false; - - if(gl[0]*gl[1]>0){ //Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2 - s1=gl[2]/(gl[2]-gl[1]); - s2=gl[2]/(gl[2]-gl[0]); - } - else if(gl[1]*gl[2]>0){ //Nodes 1 and 2 are similar, so points must be found on segment 0-1 and 0-2 - s1=gl[0]/(gl[0]-gl[1]); - s2=gl[0]/(gl[0]-gl[2]); - } - else if(gl[0]*gl[2]>0){ //Nodes 0 and 2 are similar, so points must be found on segment 1-0 and 1-2 - s1=gl[1]/(gl[1]-gl[0]); - s2=gl[1]/(gl[1]-gl[2]); - } - else _error_("case not possible"); - if(mainlyfloating){ - phi = (1-s1*s2); - } - else{ - phi = s1*s2; - } - } - } - else _error_("mesh type "<=0.); - return 1-phi; -} -/*}}}*/ -void Tria::GetFractionGeometry(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* pmainlynegative, IssmDouble* gl){/*{{{*/ +void Tria::GetFractionGeometry(IssmDouble* weights, IssmDouble* pphi, int* ppoint1,IssmDouble* pfraction1,IssmDouble* pfraction2, bool* ptrapezeisnegative, IssmDouble* gl){/*{{{*/ /*Computeportion of the element that is grounded*/ - bool negative=false; + bool trapezeisnegative=true; //default value int point; const IssmPDouble epsilon= 1.e-15; - IssmDouble f1,f2; + IssmDouble f1,f2,phi; + + /*Weights: */ + Gauss* gauss=NULL; + IssmDouble loadweights_g[NUMVERTICES]; + IssmDouble total_weight=0; + + _assert_(!xIsNan(gl[0])); + _assert_(!xIsNan(gl[1])); + _assert_(!xIsNan(gl[2])); /*Be sure that values are not zero*/ @@ -1815,5 +1782,5 @@ } else{ - if(gl[0]*gl[1]*gl[2]<0) negative=true; + if(gl[0]*gl[1]*gl[2]<0) trapezeisnegative=false; //no matter what configuration, there has to be two positive vertices, which means the trapeze is positive. if(gl[0]*gl[1]>0){ //Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2 @@ -1834,10 +1801,467 @@ else _error_("case not possible"); } - *point1=point; - *fraction1=f1; - *fraction2=f2; - *pmainlynegative=negative; -} -/*}}}*/ + if(trapezeisnegative) phi=1-f1*f2; + else phi=f1*f2; + + + /*Compute weights:*/ + gauss = this->NewGauss(point,f1,f2,trapezeisnegative,2); + + total_weight=0; + for(int i=0;inext()){ + TriaRef::GetNodalFunctions(&loadweights_g[0], gauss,P1Enum); + for(int i=0;iweight; + total_weight+=gauss->weight; + } + //normalize to phi. + if(total_weight)for(int i=0;i barycenter_trapeze = (centroid_el - barycenter_triangle * (1-phi_trapeze) )/phi_trapeze + for (int i=0;i<3;i++) barycenter[i] =(centroid[i] -barycenter[i]*(1.0-phi))/phi; + + } + + //recompute planetradius from the barycenter onwards: + planetradius=sqrt( pow(barycenter[0],2.0)+ pow(barycenter[1],2.0)+ pow(barycenter[2],2.0)); + + *platbar=asin(barycenter[2]/planetradius)*180.0/M_PI; + *plongbar=atan2(barycenter[1],barycenter[0])*180.0/M_PI; + +} /*}}}*/ +void Tria::GetNodalWeightsAndAreaAndCentroidsFromLeveset(IssmDouble* loadweights, IssmDouble* ploadarea, IssmDouble* platbar, IssmDouble* plongbar, IssmDouble late, IssmDouble longe, IssmDouble area, int levelsetenum){ /*{{{*/ + + IssmDouble phi; + IssmDouble fraction1,fraction2; + bool istrapeze1; + bool flip1=false; + int point1; + IssmDouble levelset[NUMVERTICES]; + IssmDouble planetradius; + + this->parameters->FindParam(&planetradius,SolidearthPlanetRadiusEnum); + + //figure out if we are flipping the levelsets: + if(levelsetenum<0){ + levelsetenum=-levelsetenum; + flip1=true; + } + //figure out area where we have loads + Element::GetInputListOnVertices(&levelset[0],levelsetenum); + if(flip1)for(int i=0;iGetFractionGeometry(loadweights,&phi,&point1,&fraction1,&fraction2,&istrapeze1,levelset); + + this->GetBarycenterFromLevelset(platbar,plongbar, phi, fraction1, fraction2, late, longe, point1,istrapeze1,planetradius); + + /*assign output pointers:*/ + *ploadarea=phi*area; + +} /*}}}*/ +void Tria::GetNodalWeightsAndAreaAndCentroidsFromLeveset(IssmDouble* loadweights, IssmDouble* ploadarea, IssmDouble* platbar, IssmDouble* plongbar, IssmDouble late, IssmDouble longe, IssmDouble area, int levelset1enum, int levelset2enum){ /*{{{*/ + + + bool istrapeze1, istrapeze2; + IssmDouble phi1,phi2, f11,f12,f21,f22; + int point1, point2, i0,i1,i2,j0,j1,j2; + IssmDouble weights1[3],weights2[3]; + IssmDouble levelset1[3]; + IssmDouble levelset2[3]; + + bool flip1=false; + bool flip2=false; + + IssmDouble xyz0[3][3]; + IssmDouble xyz1[3][3]={0}; + IssmDouble xyz2[3][3]={0}; + IssmDouble xyz3[3][3]={0}; + IssmDouble xyz[8][3]={0}; + IssmDouble f1o; + IssmDouble w[8][NUMVERTICES]={0}; + IssmDouble areasub=0; + IssmDouble area1=0; + IssmDouble area2=0; + IssmDouble area3=0; + + int tria0[3]={0,1,2}; + int tria1[3]={-1}; + int tria2[3]={-1}; + int tria3[3]={-1}; + + IssmDouble w1[3][3]={0}; + IssmDouble w2[3][3]={0}; + IssmDouble w3[3][3]={0}; + + IssmDouble barycenter[3]={0}; + IssmDouble planetradius; + + //figure out if we are flipping the levelsets: + if(levelset1enum<0){ + levelset1enum=-levelset1enum; + flip1=true; + } + if(levelset2enum<0){ + levelset2enum=-levelset2enum; + flip2=true; + } + + //recover levelsets: + Element::GetInputListOnVertices(&levelset1[0],levelset1enum); + if(flip1)for(int i=0;i=0 on all vertices + if (levelset1[0]>=0 && levelset1[1]>=0 && levelset1[2]>=0) { + for (int i=0;i=0 && levelset2[1]>=0 && levelset2[2]>=0) { + for (int i=0;iparameters->FindParam(&planetradius,SolidearthPlanetRadiusEnum); + + + //If just one levelset is all negative, just take the partitioning of the other, no interaction between them + if (levelset1[0]<=0 && levelset1[1]<=0 && levelset1[2]<=0) { + this->GetFractionGeometry(loadweights,&phi2,&point2,&f21,&f22,&istrapeze2,levelset2); + this->GetBarycenterFromLevelset(platbar,plongbar, phi2, f21, f22, late, longe, point2,istrapeze2,planetradius); + *ploadarea=area*phi2; + return; + } + if (levelset2[0]<=0 && levelset2[1]<=0 && levelset2[2]<=0) { + this->GetFractionGeometry(loadweights,&phi1,&point1,&f11,&f12,&istrapeze1,levelset1); + this->GetBarycenterFromLevelset(platbar,plongbar, phi1, f11, f12, late, longe, point1,istrapeze1,planetradius); + *ploadarea=area*phi1; + return; + } + + + this->GetFractionGeometry(&weights1[0],&phi1,&point1,&f11,&f12,&istrapeze1,levelset1); + this->GetFractionGeometry(&weights2[0],&phi2,&point2,&f21,&f22,&istrapeze2,levelset2); + + + //Early return if levelsets are not independent + if (istrapeze1==istrapeze2 && point1==point2 && phi1==phi2){ + //the two levelsets are redundant: levelset1 = positivescalar * levelset2 + this->GetBarycenterFromLevelset(platbar,plongbar, phi1, f11, f12, late, longe, point1,istrapeze1,planetradius); + *ploadarea=area*phi1; + for (int i=0;i point_f12} where the levelsets intersect (only used when f1o>=0 and f1o<=1) + if(point2==i0) f1o= f22*(f11-f21)/(f11*f22-f12*f21); + else if(point2==i1) f1o=f21*(1.0-f22-f11)/(f21*(f12-f11)-f12*f22); + else f1o= (f22*(1.0-f21-f11)+f21*f11)/(f22*(f12-f11) +f21*f11); + + + //interpolant weights of each point. Any field F[0,1,2] provided at the original vertices [0,1,2] will be equal on point k to sum_i (F[i] * w[k][i]) + w[0][0]=1; + w[1][1]=1; + w[2][2]=1; + w[3][i0]=1.0-f11; w[3][i1]=f11; + w[4][i0]=1.0-f12; w[4][i2]=f12; + w[5][j0]=1.0-f21; w[5][j1]=f21; + w[6][j0]=1.0-f22; w[6][j2]=f22; + for (int j=0;j<3;j++) w[7][j]=w[3][j]*(1.0-f1o) + w[4][j]*f1o; //we interpolate the intersection point between point_f11 and point_f12 at fraction f1o + + for (int k=0;k<8;k++){ + for (int i=0;i1.0-f21){ /*{{{*/ + if (!istrapeze1 && !istrapeze2){ + tria1[0]=5; tria1[1]= 7; tria1[2]= 4; + } + else if (!istrapeze1 && istrapeze2){ + tria1[0]=i0; tria1[1]= 3; tria1[2]= 7; + tria2[0]=i0; tria2[1]= 7; tria2[2]= 5; + } + else if (istrapeze1 && !istrapeze2){ + tria1[0]=7; tria1[1]= 6; tria1[2]= 4; + tria2[0]=4; tria2[1]= 6; tria2[2]= i2; + } + else { //istrapeze1 && istrapeze2 + tria1[0]=3; tria1[1]= i1; tria1[2]= 7; + tria2[0]=7; tria2[1]= i1; tria2[2]= 6; + } /*}}}*/ + } + else if (f12<=1.0-f21){ /*{{{*/ + if (!istrapeze1 && !istrapeze2){ + } + else if (!istrapeze1 && istrapeze2){ + tria1[0]=i0; tria1[1]= 3; tria1[2]= 4; + } + else if (istrapeze1 && !istrapeze2){ + tria1[0]=5; tria1[1]= 6; tria1[2]= i2; + } + else { //istrapeze1 && istrapeze2 + tria1[0]=3; tria1[1]= i1; tria1[2]= 4; + tria2[0]=4; tria2[1]= i1; tria2[2]= 5; + tria3[0]=5; tria3[1]= i1; tria3[2]= 6; + } /*}}}*/ + } + }/*}}}*/ + else if(point2==i1){ /*{{{*/ + if (f11>1.0-f22){ /*{{{*/ + if (!istrapeze1 && !istrapeze2){ + tria1[0]=6; tria1[1]= 3; tria1[2]= 7; + } + else if (!istrapeze1 && istrapeze2){ + tria1[0]=i0; tria1[1]= 6; tria1[2]= 7; + tria2[0]=i0; tria2[1]= 7; tria2[2]= 4; + } + else if (istrapeze1 && !istrapeze2){ + tria1[0]=3; tria1[1]= i1; tria1[2]= 7; + tria2[0]=7; tria2[1]= i1; tria2[2]= 5; + } + else { //istrapeze1 && istrapeze2 + tria1[0]=7; tria1[1]= 5; tria1[2]= 4; + tria2[0]=4; tria2[1]= 5; tria2[2]= i2; + } /*}}}*/ + } + else if (f11<=1.0-f22){ /*{{{*/ + if (!istrapeze1 && !istrapeze2){ + } + else if (!istrapeze1 && istrapeze2){ + tria1[0]=i0; tria1[1]= 3; tria1[2]= 4; + } + else if (istrapeze1 && !istrapeze2){ + tria1[0]=6; tria1[1]= i1; tria1[2]= 5; + } + else { //istrapeze1 && istrapeze2 + tria1[0]=3; tria1[1]= 6; tria1[2]= 4; + tria2[0]=4; tria2[1]= 6; tria2[2]= 5; + tria3[0]=4; tria3[1]= 5; tria3[2]= i2; + } /*}}}*/ + } + + }/*}}}*/ + else{ /*{{{*/ + if (f11<=f21 && f12>=f22){ /*{{{*/ + if (!istrapeze1 && !istrapeze2){ + tria1[0]=i0; tria1[1]= 3; tria1[2]= 7; + tria2[0]=i0; tria2[1]= 7; tria2[2]= 6; + } + else if (!istrapeze1 && istrapeze2){ + tria1[0]=6; tria1[1]= 7; tria1[2]= 4; + } + else if (istrapeze1 && !istrapeze2){ + tria1[0]=3; tria1[1]= 5; tria1[2]= 7; + } + else { //istrapeze1 && istrapeze2 + tria1[0]=7; tria1[1]= 5; tria1[2]= i1; + tria2[0]=7; tria2[1]= i1; tria2[2]= 4; + tria3[0]=4; tria3[1]= i1; tria3[2]= i2; + } /*}}}*/ + } + else if (f11>=f21 && f12<=f22){ /*{{{*/ + if (!istrapeze1 && !istrapeze2){ + tria1[0]=i0; tria1[1]= 5; tria1[2]= 7; + tria2[0]=i0; tria2[1]= 4; tria2[2]= 7; + }else if (!istrapeze1 && istrapeze2){ + tria1[0]=5; tria1[1]= 3; tria1[2]= 7; + }else if (istrapeze1 && !istrapeze2){ + tria1[0]=4; tria1[1]= 7; tria1[2]= 6; + }else { //istrapeze1 && istrapeze2 + tria1[0]=3; tria1[1]= i1; tria1[2]= 7; + tria2[0]=7; tria2[1]= i1; tria2[2]= 6; + tria3[0]=6; tria3[1]= i1; tria3[2]= i2; + } /*}}}*/ + } + else if (f11<=f21 && f12<=f22){ /*{{{*/ + if (!istrapeze1 && !istrapeze2){ + tria1[0]=i0; tria1[1]= 3; tria1[2]= 4; + }else if (!istrapeze1 && istrapeze2){ + }else if (istrapeze1 && !istrapeze2){ + tria1[0]=3; tria1[1]= 5; tria1[2]= 4; + tria2[0]=4; tria2[1]= 5; tria2[2]= 6; + }else { //istrapeze1 && istrapeze2 + tria1[0]=6; tria1[1]= 5; tria1[2]= i1; + tria2[0]=6; tria2[1]= i1; tria2[2]= i2; + } /*}}}*/ + } + else if (f11>=f21 && f12>=f22){ /*{{{*/ + if (!istrapeze1 && !istrapeze2){ + tria1[0]=i0; tria1[1]= 5; tria1[2]= 6; + }else if (!istrapeze1 && istrapeze2){ + tria1[0]=5; tria1[1]= 3; tria1[2]= 6; + tria2[0]=6; tria2[1]= 3; tria2[2]= 4; + }else if (istrapeze1 && !istrapeze2){ + }else { //istrapeze1 && istrapeze2 + tria1[0]=4; tria1[1]= 3; tria1[2]= i1; + tria2[0]=i1; tria2[1]= i2; tria2[2]= 4; + } /*}}}*/ + } + } /*}}}*/ + + if(tria1[0]>-1){ + for (int i=0;i-1){ + for (int i=0;i-1){ + for (int i=0;i0){ + for (int j=0;j<3;j++){ + for (int i=0;iisiceonly[this->lid]=this->IsIceOnlyInElement(); - masks->isoceanin[this->lid]=this->IsOceanInElement(); - - /*are we fully floating:*/ - Input* gr_input=this->GetInput(MaskOceanLevelsetEnum); _assert_(gr_input); - if (gr_input->GetInputMax()<=0)masks->isfullyfloating[this->lid]=true; - else masks->isfullyfloating[this->lid]=false; - - /*are we not fully grounded: */ - if ((gr_input->GetInputMin())<0) masks->notfullygrounded[this->lid]=true; - else masks->notfullygrounded[this->lid]=false; -} -/*}}}*/ void Tria::GiaDeflection(Vector* wg,Vector* dwgdt, Matlitho* litho, IssmDouble* x, IssmDouble* y){/*{{{*/ @@ -5676,7 +6084,7 @@ } /*}}}*/ -void Tria::SealevelchangeGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae){ /*{{{*/ - - /*diverse:*/ +void Tria::SealevelchangeGeometryInitial(IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae){ /*{{{*/ + + /*Declarations:{{{*/ int nel; IssmDouble area,planetarea,planetradius; @@ -5684,8 +6092,11 @@ IssmDouble I; //change in ice thickness or water level(Farrel and Clarke, Equ. 4) IssmDouble rho_earth; + IssmDouble NewtonG; + IssmDouble g; IssmDouble lati,longi; + IssmDouble latitude[NUMVERTICES]; + IssmDouble longitude[NUMVERTICES]; IssmDouble x,y,z,dx,dy,dz,N_azim,E_azim; - int sidlist[NUMVERTICES]; - int sid; + IssmDouble xyz_list[NUMVERTICES][3]; #ifdef _HAVE_RESTRICT_ @@ -5698,6 +6109,4 @@ IssmDouble* __restrict__ U_elastic_precomputed=NULL; IssmDouble* __restrict__ H_elastic_precomputed=NULL; - IssmDouble* __restrict__ lates=NULL; - IssmDouble* __restrict__ longes=NULL; #else IssmDouble* G=NULL; @@ -5709,6 +6118,4 @@ IssmDouble* U_elastic_precomputed=NULL; IssmDouble* H_elastic_precomputed=NULL; - IssmDouble* lates=NULL; - IssmDouble* longes=NULL; #endif @@ -5728,12 +6135,11 @@ IssmDouble* load_love_k = NULL; IssmDouble tide_love_k2secular; - IssmDouble moi_e, moi_p, omega, g; + IssmDouble moi_e, moi_p, omega; IssmDouble Grotm1[3]; IssmDouble Grotm2[3]; IssmDouble Grotm3[3]; IssmDouble pre; - - - /*recover parameters: */ + /*}}}*/ + /*Recover parameters:{{{ */ rho_earth=FindParam(MaterialsEarthDensityEnum); this->parameters->FindParam(&computerigid,SolidearthSettingsRigidEnum); @@ -5743,4 +6149,5 @@ this->parameters->FindParam(&planetarea,SolidearthPlanetAreaEnum); this->parameters->FindParam(&planetradius,SolidearthPlanetRadiusEnum); + this->parameters->FindParam(&NewtonG,ConstantsNewtonGravityEnum); this->parameters->FindParam(&horiz,SolidearthSettingsHorizEnum); @@ -5754,16 +6161,12 @@ parameters->FindParam(&omega,RotationalAngularVelocityEnum); } + /*}}}*/ /*early return:*/ if(!computerigid)return; - /*recover precomputed green function kernels:*/ + /*Recover precomputed green function kernels:{{{*/ DoubleVecParam* parameter = static_cast(this->parameters->FindParamObject(SealevelchangeGRigidEnum)); _assert_(parameter); parameter->GetParameterValueByPointer((IssmDouble**)&G_rigid_precomputed,&M); - - /*allocate arrays:*/ - G=xNewZeroInit(3*nel); - longes=xNewZeroInit(nel); - lates=xNewZeroInit(nel); if(computeelastic){ @@ -5777,4 +6180,21 @@ parameter->GetParameterValueByPointer((IssmDouble**)&U_elastic_precomputed,&M); + + } + /*}}}*/ + /*Compute lat long of all vertices in the element:{{{*/ + ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); + for(int i=0;i(3*nel); + if(computeelastic){ GU=xNewZeroInit(3*nel); if(horiz){ @@ -5783,16 +6203,4 @@ } } - - /*compute centroids of all elements:*/ - for(int i=0;iGetVerticesSidList(&sidlist[0]); for(int e=0;eM_PI)delLambda=2*M_PI-delLambda; alpha=2.*asin(sqrt(pow(sin(delPhi/2),2)+cos(lati)*cos(late)*pow(sin(delLambda/2),2))); index=reCast( alpha/M_PI*reCast(M-1) ); + _assert_(index>0 && indexAddInput(SealevelGrotm1Enum,&Grotm1[0],P1Enum); - this->AddInput(SealevelGrotm2Enum,&Grotm2[0],P1Enum); - this->AddInput(SealevelGrotm3Enum,&Grotm3[0],P1Enum); } @@ -5869,6 +6263,26 @@ } } - - /*Free allocations:*/ + /*}}}*/ + /*Compute rotation kernel:{{{*/ + if(computerotation){ + + /*What is the gravity at planet's surface: */ + g=4.0/3.0*M_PI*rho_earth*NewtonG*planetradius; + + for(int i=0;i<3;i++){ + lati=latitude[i]; + longi=longitude[i]; + + pre=((1.0+tide_love_k[2]-tide_love_h[2])/g)*pow(omega*planetradius,2.0); + Grotm1[i]= - pre* 0.5*sin(2.*lati)*cos(longi); + Grotm2[i]= - pre* 0.5*sin(2.*lati)*sin(longi); + Grotm3[i]= - pre* (1.0/6.0 - 0.5*cos(2.0*lati)); + } + this->AddInput(SealevelGrotm1Enum,&Grotm1[0],P1Enum); + this->AddInput(SealevelGrotm2Enum,&Grotm2[0],P1Enum); + this->AddInput(SealevelGrotm3Enum,&Grotm3[0],P1Enum); + } + /*}}}*/ + /*Free allocations:{{{*/ #ifdef _HAVE_RESTRICT_ delete G; @@ -5880,6 +6294,4 @@ } } - delete lates; - delete longes; #else xDelete(G); @@ -5891,33 +6303,513 @@ } } - delete lates; - delete longes; #endif - + /*}}}*/ + + return; -} -/*}}}*/ -void Tria::SealevelchangeBarystaticLoads(Vector* loads, BarystaticContributions* barycontrib, SealevelMasks* masks){ /*{{{*/ - - //Compute sea level barystatic loads (ice, water column or bottom pressure, see Farrel and Clarke, Equ. 4) - - /*diverse:*/ + +} +/*}}}*/ +void Tria::SealevelchangeGeometryCentroidLoads(SealevelGeometry* slgeom, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae){ /*{{{*/ + + /* Classic buildup of load weights, centroids and areas *for elements which are fully inside a mask. + * At the same time, we'll tag the elements that are fractionally only inside a mask*/ + + IssmDouble loadweights[3]={0}; IssmDouble area; - IssmDouble phi_ice=1.0; //WARNING: do not touch this, default is entire elemnt contributes barystatic - IssmDouble phi_water=1.0; //WARNING: do not touch this, default is entire elemnt contributes barystatic - IssmDouble phi_bp=1.0; //WARNING: do not touch this, default is entire elemnt contributes barystatic - IssmDouble I=0; //Do not change this! - IssmDouble W=0; //Do not change this! - IssmDouble BP=0; //Do not change this! - - bool notfullygrounded=false; - bool computerigid= false; - int glfraction=1; - int npartice,nparthydro,npartocean; + IssmDouble loadweightsocean[3]; //to keep memory of these loads, no need to recompute for bottom pressure. + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble planetradius; + IssmDouble late,longe; /*flags:*/ bool isocean=false; - bool ishydro=false; + bool isoceanonly=false; bool isice=false; + bool isiceonly=false; + bool computeice=false; + bool computebp=false; + bool computehydro=false; + + /*constants:*/ + IssmDouble constant=0; + + /*recover parameters:*/ + this->parameters->FindParam(&computeice,TransientIsmasstransportEnum); + this->parameters->FindParam(&computebp,TransientIsoceantransportEnum); + this->parameters->FindParam(&computehydro,TransientIshydrologyEnum); + this->parameters->FindParam(&planetradius,SolidearthPlanetRadiusEnum); + + /*get vertex information:*/ + ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); + + /*answer mask questions:*/ + isiceonly=this->IsIceOnlyInElement(); + isice=this->IsIceInElement(); + isoceanonly=this->IsOceanOnlyInElement(); + isocean=this->IsOceanInElement(); + slgeom->isoceanin[this->lid]=isocean; //keep track for later. + area=areae[this->sid]; + + /*set barycentre for all elements, to be updated for fractional loads in the next routine: */ + //late= asin(zze[this->sid]/planetradius)*180.0/M_PI; + late= asin(zze[this->sid]/sqrt( pow(xxe[this->sid],2.0)+ pow(yye[this->sid],2.0)+ pow(zze[this->sid],2.0)))*180.0/M_PI; + longe= atan2(yye[this->sid],xxe[this->sid])*180.0/M_PI; + slgeom->longe[this->lid]=longe; + slgeom->late[this->lid]=late; + + /*compute fractional areas and load weights for ocean:*/ + if(isoceanonly){ + slgeom->LoadArea[SLGEOM_OCEAN][this->lid]=area; + for(int i=0;iLoadWeigths[SLGEOM_OCEAN][i][this->lid]=1.0/3.0; + + #ifdef _ISSM_DEBUG_ /*{{{*/ + /*Inform mask: */ + constant=1.0; + for(int i=0;iAddInput(SealevelBarystaticOceanMaskEnum,&constant,P0Enum); + this->AddInput(SealevelBarystaticOceanWeightsEnum,loadweightsocean,P1DGEnum); + this->AddInput(SealevelBarystaticOceanAreaEnum,&area,P0Enum); + #endif /*}}}*/ + } + else if(!isocean){ + slgeom->LoadArea[SLGEOM_OCEAN][this->lid]=0; + for(int i=0;iLoadWeigths[SLGEOM_OCEAN][i][this->lid]=0.0; + #ifdef _ISSM_DEBUG_ /*{{{*/ + /*Inform mask: */ + constant=0.0; + for(int i=0;iAddInput(SealevelBarystaticOceanMaskEnum,&constant,P0Enum); + this->AddInput(SealevelBarystaticOceanWeightsEnum,loadweightsocean,P1DGEnum); + this->AddInput(SealevelBarystaticOceanAreaEnum,&constant,P0Enum); + #endif /*}}}*/ + } + else{ + slgeom->issubelement[SLGEOM_OCEAN][this->lid]=true; + slgeom->nsubel[SLGEOM_OCEAN]++; + } + + /*early return if we are not on an ice sheet , and we are not requesting + *hydrology or bottom pressure loads :*/ + if(!computebp && !computehydro){ + if(!isice || isoceanonly) { + #ifdef _ISSM_DEBUG_ + constant=0; + this->AddInput(SealevelBarystaticIceMaskEnum,&constant,P0Enum); + this->AddInput(SealevelBarystaticIceAreaEnum,&constant,P0Enum); + this->AddInput(SealevelBarystaticIceWeightsEnum,loadweights,P1DGEnum); + this->AddInput(SealevelBarystaticHydroMaskEnum,&constant,P0Enum); + this->AddInput(SealevelBarystaticHydroWeightsEnum,loadweights,P1DGEnum); + this->AddInput(SealevelBarystaticHydroAreaEnum,&constant,P0Enum); + this->AddInput(SealevelBarystaticBpMaskEnum,&constant,P0Enum); + this->AddInput(SealevelBarystaticBpWeightsEnum,loadweights,P1DGEnum); + this->AddInput(SealevelBarystaticBpAreaEnum,&constant,P0Enum); + #endif + for(int i=0;iLoadWeigths[SLGEOM_ICE][i][this->lid]=0; + slgeom->LoadWeigths[SLGEOM_WATER][i][this->lid]=0; + } + slgeom->LoadArea[SLGEOM_ICE][this->lid]=0; + slgeom->LoadArea[SLGEOM_WATER][this->lid]=0; + return; + } + } + + /*early return if we are fully floating and we are not doing bottom pressure loads:*/ + if(!computebp){ + if (isoceanonly){ + #ifdef _ISSM_DEBUG_ + constant=0; + this->AddInput(SealevelBarystaticIceMaskEnum,&constant,P0Enum); + this->AddInput(SealevelBarystaticIceWeightsEnum,loadweights,P1DGEnum); + this->AddInput(SealevelBarystaticIceAreaEnum,&constant,P0Enum); + this->AddInput(SealevelBarystaticHydroMaskEnum,&constant,P0Enum); + this->AddInput(SealevelBarystaticHydroWeightsEnum,loadweights,P1DGEnum); + this->AddInput(SealevelBarystaticHydroAreaEnum,&constant,P0Enum); + this->AddInput(SealevelBarystaticBpMaskEnum,&constant,P0Enum); + this->AddInput(SealevelBarystaticBpWeightsEnum,loadweights,P1DGEnum); + this->AddInput(SealevelBarystaticBpAreaEnum,&constant,P0Enum); + #endif + for(int i=0;iLoadWeigths[SLGEOM_ICE][i][this->lid]=0; + slgeom->LoadWeigths[SLGEOM_WATER][i][this->lid]=0; + } + slgeom->LoadArea[SLGEOM_ICE][this->lid]=0; + slgeom->LoadArea[SLGEOM_WATER][this->lid]=0; + return; + } + } + + /*early return if we are not on the ocean and we are not doing ice mass transport of + * hydrology:*/ + if(!computeice && !computehydro){ + if(!isocean){ + #ifdef _ISSM_DEBUG_ + constant=0; + this->AddInput(SealevelBarystaticIceMaskEnum,&constant,P0Enum); + this->AddInput(SealevelBarystaticIceWeightsEnum,loadweights,P1DGEnum); + this->AddInput(SealevelBarystaticIceAreaEnum,&constant,P0Enum); + this->AddInput(SealevelBarystaticHydroMaskEnum,&constant,P0Enum); + this->AddInput(SealevelBarystaticHydroWeightsEnum,loadweights,P1DGEnum); + this->AddInput(SealevelBarystaticHydroAreaEnum,&constant,P0Enum); + this->AddInput(SealevelBarystaticBpMaskEnum,&constant,P0Enum); + this->AddInput(SealevelBarystaticBpWeightsEnum,loadweights,P1DGEnum); + this->AddInput(SealevelBarystaticBpAreaEnum,&constant,P0Enum); + #endif + for(int i=0;iLoadWeigths[SLGEOM_ICE][i][this->lid]=0; + slgeom->LoadWeigths[SLGEOM_WATER][i][this->lid]=0; + } + slgeom->LoadArea[SLGEOM_ICE][this->lid]=0; + slgeom->LoadArea[SLGEOM_WATER][this->lid]=0; + return; + } + } + + /*Deal with ice loads if we are on grounded ice:*/ + if(isice && !isoceanonly && computeice){ + if(isiceonly){ + slgeom->LoadArea[SLGEOM_ICE][this->lid]=area; + for(int i=0;iLoadWeigths[SLGEOM_ICE][i][this->lid]=1.0/3.0; + + #ifdef _ISSM_DEBUG_ /*{{{*/ + /*Inform mask: */ + constant=1.0; + for(int i=0;iAddInput(SealevelBarystaticIceMaskEnum,&constant,P0Enum); + this->AddInput(SealevelBarystaticIceWeightsEnum,loadweights,P1DGEnum); + this->AddInput(SealevelBarystaticIceAreaEnum,&area,P0Enum); + #endif /*}}}*/ + } + else{ + slgeom->issubelement[SLGEOM_ICE][this->lid]=true; + slgeom->nsubel[SLGEOM_ICE]++; + } + } + + /*Deal with water loads if we are on ground:*/ + if(!isoceanonly && computehydro){ + + if(!isocean){ + slgeom->LoadArea[SLGEOM_WATER][this->lid]=area; + for(int i=0;iLoadWeigths[SLGEOM_WATER][i][this->lid]=1.0/3.0; + + #ifdef _ISSM_DEBUG_ /*{{{*/ + /*Inform mask: */ + constant=1.0; + for(int i=0;iAddInput(SealevelBarystaticHydroMaskEnum,&constant,P0Enum); + this->AddInput(SealevelBarystaticHydroWeightsEnum,loadweights,P1DGEnum); + this->AddInput(SealevelBarystaticHydroAreaEnum,&area,P0Enum); + #endif /*}}}*/ + } + else{ + slgeom->issubelement[SLGEOM_WATER][this->lid]=true; + slgeom->nsubel[SLGEOM_WATER]++; + } + } + +} +/*}}}*/ +void Tria::SealevelchangeGeometrySubElementLoads(SealevelGeometry* slgeom, IssmDouble* areae){ /*{{{*/ + + /* Classic buildup of load weights, centroids and areas *for elements which are fully inside a mask. + * At the same time, we'll tag the elements that are fractionally only inside a mask*/ + + IssmDouble loadweights[3]={0}; + IssmDouble area,loadarea; + IssmDouble loadareaocean; + IssmDouble loadweightsocean[3]; //to keep memory of these loads, no need to recompute for bottom pressure. + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble latbar=0; + IssmDouble longbar=0; + IssmDouble constant; + IssmDouble nanconstant=NAN; + + /*get vertex and area information:*/ + ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); + area=areae[this->sid]; + + #ifdef _ISSM_DEBUG_ + this->AddInput(SealevelBarystaticIceLatbarEnum,&nanconstant,P0Enum); + this->AddInput(SealevelBarystaticIceLongbarEnum,&nanconstant,P0Enum); + this->AddInput(SealevelBarystaticHydroLatbarEnum,&nanconstant,P0Enum); + this->AddInput(SealevelBarystaticHydroLongbarEnum,&nanconstant,P0Enum); + this->AddInput(SealevelBarystaticOceanLatbarEnum,&nanconstant,P0Enum); + this->AddInput(SealevelBarystaticOceanLongbarEnum,&nanconstant,P0Enum); + #endif + + if(slgeom->issubelement[SLGEOM_OCEAN][this->lid]){ + int intj=slgeom->subelementmapping[SLGEOM_OCEAN][this->lid]; + + this->GetNodalWeightsAndAreaAndCentroidsFromLeveset(&loadweightsocean[0],&loadareaocean,&latbar, &longbar, slgeom->late[this->lid], slgeom->longe[this->lid], area, MaskOceanLevelsetEnum); + slgeom->LoadArea[SLGEOM_OCEAN][this->lid]=loadareaocean; + slgeom->vareae_subel[SLGEOM_OCEAN]->SetValue(intj,loadareaocean,INS_VAL); + slgeom->vlatbarycentre[SLGEOM_OCEAN]->SetValue(intj,latbar,INS_VAL); + slgeom->vlongbarycentre[SLGEOM_OCEAN]->SetValue(intj,longbar,INS_VAL); + + for(int i=0;iLoadWeigths[SLGEOM_OCEAN][i][this->lid]=loadweightsocean[i]; + #ifdef _ISSM_DEBUG_ /*{{{*/ + /*Inform mask: */ + constant=loadareaocean/area; + this->AddInput(SealevelBarystaticOceanMaskEnum,&constant,P0Enum); + this->AddInput(SealevelBarystaticOceanWeightsEnum,loadweightsocean,P1DGEnum); + this->AddInput(SealevelBarystaticOceanAreaEnum,&loadareaocean,P0Enum); + + this->AddInput(SealevelBarystaticOceanLatbarEnum,&latbar,P0Enum); + this->AddInput(SealevelBarystaticOceanLongbarEnum,&longbar,P0Enum); + #endif /*}}}*/ + } + if(slgeom->issubelement[SLGEOM_ICE][this->lid]){ + int intj=slgeom->subelementmapping[SLGEOM_ICE][this->lid]; + + this->GetNodalWeightsAndAreaAndCentroidsFromLeveset(&loadweights[0],&loadarea,&latbar, &longbar, slgeom->late[this->lid], slgeom->longe[this->lid], area, -MaskOceanLevelsetEnum,MaskIceLevelsetEnum); + + slgeom->LoadArea[SLGEOM_ICE][this->lid]=loadarea; + slgeom->vareae_subel[SLGEOM_ICE]->SetValue(intj,loadarea,INS_VAL); + slgeom->vlatbarycentre[SLGEOM_ICE]->SetValue(intj,latbar,INS_VAL); + slgeom->vlongbarycentre[SLGEOM_ICE]->SetValue(intj,longbar,INS_VAL); + + for(int i=0;iLoadWeigths[SLGEOM_ICE][i][this->lid]=loadweights[i]; + + #ifdef _ISSM_DEBUG_ + /*Inform mask: */ + constant=loadarea/area; + this->AddInput(SealevelBarystaticIceMaskEnum,&constant,P0Enum); + this->AddInput(SealevelBarystaticIceWeightsEnum,loadweights,P1DGEnum); + this->AddInput(SealevelBarystaticIceAreaEnum,&loadarea,P0Enum); + + this->AddInput(SealevelBarystaticIceLatbarEnum,&latbar,P0Enum); + this->AddInput(SealevelBarystaticIceLongbarEnum,&longbar,P0Enum); + + #endif + } + if(slgeom->issubelement[SLGEOM_WATER][this->lid]){ + int intj=slgeom->subelementmapping[SLGEOM_WATER][this->lid]; + + this->GetNodalWeightsAndAreaAndCentroidsFromLeveset(&loadweights[0],&loadarea,&latbar, &longbar, slgeom->late[this->lid], slgeom->longe[this->lid], area, -MaskOceanLevelsetEnum); + + slgeom->LoadArea[SLGEOM_WATER][this->lid]=loadarea; + slgeom->vareae_subel[SLGEOM_WATER]->SetValue(intj,loadarea,INS_VAL); + slgeom->vlatbarycentre[SLGEOM_WATER]->SetValue(intj,latbar,INS_VAL); + slgeom->vlongbarycentre[SLGEOM_WATER]->SetValue(intj,longbar,INS_VAL); + + for(int i=0;iLoadWeigths[SLGEOM_WATER][i][this->lid]=loadweights[i]; + + #ifdef _ISSM_DEBUG_ + /*Inform mask: */ + constant=loadarea/area; + this->AddInput(SealevelBarystaticHydroMaskEnum,&constant,P0Enum); + this->AddInput(SealevelBarystaticHydroWeightsEnum,loadweights,P1DGEnum); + this->AddInput(SealevelBarystaticHydroAreaEnum,&loadarea,P0Enum); + + this->AddInput(SealevelBarystaticHydroLatbarEnum,&latbar,P0Enum); + this->AddInput(SealevelBarystaticHydroLongbarEnum,&longbar,P0Enum); + + #endif + } + +} +/*}}}*/ +void Tria::SealevelchangeGeometryFractionKernel(SealevelGeometry* slgeom){ /*{{{*/ + + /*Declarations:{{{*/ + int nel; + IssmDouble planetarea,planetradius; + IssmDouble constant,ratioe; + IssmDouble rho_earth; + IssmDouble lati,longi; + IssmDouble latitude[NUMVERTICES]; + IssmDouble longitude[NUMVERTICES]; + IssmDouble x,y,z,dx,dy,dz,N_azim,E_azim; + IssmDouble xyz_list[NUMVERTICES][3]; + + #ifdef _HAVE_RESTRICT_ + 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__ Gsubel=NULL; + IssmDouble** __restrict__ GUsubel=NULL; + IssmDouble** __restrict__ GNsubel=NULL; + IssmDouble** __restrict__ GEsubel=NULL; + + #else + IssmDouble* G_elastic_precomputed=NULL; + IssmDouble* G_rigid_precomputed=NULL; + IssmDouble* U_elastic_precomputed=NULL; + IssmDouble* H_elastic_precomputed=NULL; + IssmDouble** Gsubel=NULL; + IssmDouble** GUsubel=NULL; + IssmDouble** GNsubel=NULL; + IssmDouble** GEsubel=NULL; + #endif + + /*elastic green function:*/ + int index; + int M; + + /*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(&nel,MeshNumberofelementsEnum); + this->parameters->FindParam(&planetarea,SolidearthPlanetAreaEnum); + this->parameters->FindParam(&planetradius,SolidearthPlanetRadiusEnum); + this->parameters->FindParam(&horiz,SolidearthSettingsHorizEnum); + /*}}}*/ + + /*early return:*/ + 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); + + if(computeelastic){ + parameter = static_cast(this->parameters->FindParamObject(SealevelchangeGElasticEnum)); _assert_(parameter); + parameter->GetParameterValueByPointer((IssmDouble**)&G_elastic_precomputed,&M); + + parameter = static_cast(this->parameters->FindParamObject(SealevelchangeHElasticEnum)); _assert_(parameter); + parameter->GetParameterValueByPointer((IssmDouble**)&H_elastic_precomputed,&M); + + parameter = static_cast(this->parameters->FindParamObject(SealevelchangeUElasticEnum)); _assert_(parameter); + parameter->GetParameterValueByPointer((IssmDouble**)&U_elastic_precomputed,&M); + } + /*}}}*/ + /*Compute lat long of all vertices in the element:{{{*/ + ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); + for(int i=0;i(SLGEOM_NUMLOADS); + if(computeelastic){ + GUsubel=xNew(SLGEOM_NUMLOADS); + if(horiz){ + GNsubel=xNew(SLGEOM_NUMLOADS); + GEsubel=xNew(SLGEOM_NUMLOADS); + } + } + + //Allocate: + for(int l=0;lnbar[l]; + Gsubel[l]=xNewZeroInit(3*nbar); + if(computeelastic){ + GUsubel[l]=xNewZeroInit(3*nbar); + if(horiz){ + GNsubel[l]=xNewZeroInit(3*nbar); + GEsubel[l]=xNewZeroInit(3*nbar); + } + } + + for(int e=0;earea_subel[l][e]; + for (int i=0;i<3;i++){ + IssmDouble alpha; + IssmDouble delPhi,delLambda; + /*recover info for this element and vertex:*/ + IssmDouble late= slgeom->latbarycentre[l][e]; + IssmDouble longe= slgeom->longbarycentre[l][e]; + late=late/180*M_PI; + longe=longe/180*M_PI; + + lati=latitude[i]; + longi=longitude[i]; + + /*Compute alpha angle between centroid and current vertex and index into precomputed tables: */ + delPhi=fabs(lati-late); delLambda=fabs(longi-longe); if (delLambda>M_PI)delLambda=2*M_PI-delLambda; + alpha=2.*asin(sqrt(pow(sin(delPhi/2.0),2.0)+cos(lati)*cos(late)*pow(sin(delLambda/2.0),2.0))); + index=reCast( alpha/M_PI*reCast(M-1) ); + + /*Rigid earth gravitational perturbation: */ + Gsubel[l][i*nbar+e]+=G_rigid_precomputed[index]; + + if(computeelastic){ + Gsubel[l][i*nbar+e]+=G_elastic_precomputed[index]; + } + Gsubel[l][i*nbar+e]*=ratioe; + + /*Elastic components:*/ + if(computeelastic){ + GUsubel[l][i*nbar+e] = ratioe * U_elastic_precomputed[index]; + if(horiz){ + /*Compute azimuths, both north and east components: */ + x = xyz_list[i][0]; y = xyz_list[i][1]; z = xyz_list[i][2]; + if(lati==90){ + x=1e-12; y=1e-12; + } + if(lati==-90){ + x=1e-12; y=1e-12; + } + IssmDouble xbar=planetradius*cos(late)*cos(longe); + IssmDouble ybar=planetradius*cos(late)*sin(longe); + IssmDouble zbar=planetradius*sin(late); + + dx = xbar-x; dy = ybar-y; dz = zbar-z; + N_azim = (-z*x*dx-z*y*dy+(pow(x,2)+pow(y,2))*dz) /pow((pow(x,2)+pow(y,2))*(pow(x,2)+pow(y,2)+pow(z,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5); + E_azim = (-y*dx+x*dy) /pow((pow(x,2)+pow(y,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5); + + GNsubel[l][i*nbar+e] = ratioe*H_elastic_precomputed[index]*N_azim; + GEsubel[l][i*nbar+e] = ratioe*H_elastic_precomputed[index]*E_azim; + } + } + } + } + } + + /*Save all these arrayins for each element:*/ + for (int l=0;linputs->SetArrayInput(slgeom->GEnum(l),this->lid,Gsubel[l],slgeom->nbar[l]*3); + if(computeelastic){ + this->inputs->SetArrayInput(slgeom->GUEnum(l),this->lid,GUsubel[l],slgeom->nbar[l]*3); + if(horiz){ + this->inputs->SetArrayInput(slgeom->GNEnum(l),this->lid,GNsubel[l],slgeom->nbar[l]*3); + this->inputs->SetArrayInput(slgeom->GEEnum(l),this->lid,GEsubel[l],slgeom->nbar[l]*3); + } + } + } + /*}}}*/ + /*Free memory:{{{*/ + for (int l=0;l(Gsubel[l]); + if(computeelastic){ + xDelete(GUsubel[l]); + if(horiz){ + xDelete(GNsubel[l]); + xDelete(GEsubel[l]); + } + } + } + xDelete(Gsubel); + if(computeelastic){ + xDelete(GUsubel); + if(horiz){ + xDelete(GNsubel); + xDelete(GEsubel); + } + } + /*}}}*/ + return; + +} +/*}}}*/ +void Tria::SealevelchangeBarystaticLoads(Vector* loads, Vector* subelementiceloads, Vector* subelementhydroloads, Vector* subelementbploads, BarystaticContributions* barycontrib, SealevelGeometry* slgeom){ /*{{{*/ + + /*Inputs:*/ + IssmDouble I[NUMVERTICES]; + IssmDouble W[NUMVERTICES]; + IssmDouble BP[NUMVERTICES]; /*output: */ @@ -5925,53 +6817,10 @@ IssmDouble bslchydro=0; IssmDouble bslcbp=0; + IssmDouble BPavg=0; + IssmDouble Iavg=0; + IssmDouble Wavg=0; /*ice properties: */ IssmDouble rho_ice,rho_water,rho_freshwater; - - /*constants:*/ - IssmDouble constant=0; - - /*recover parameters:*/ - this->parameters->FindParam(&isice,TransientIsmasstransportEnum); - this->parameters->FindParam(&isocean,TransientIsoceantransportEnum); - this->parameters->FindParam(&ishydro,TransientIshydrologyEnum); - - /*early return if we are not on an ice cap, and we are not requesting - *hydrology or bottom pressure loads :*/ - if(!isocean && !ishydro){ - if(!masks->isiceonly[this->lid]){ - #ifdef _ISSM_DEBUG_ - constant=0; this->AddInput(SealevelBarystaticMaskEnum,&constant,P0Enum); - #endif - return; - } - } - - /*early return if we are fully floating and we are not doing bottom pressure loads:*/ - if(!isocean){ - if (masks->isfullyfloating[this->lid]){ - constant=0; - #ifdef _ISSM_DEBUG_ - this->AddInput(SealevelBarystaticMaskEnum,&constant,P0Enum); - #endif - return; - } - } - - /*early return if we are not on the ocean and we are not doing ice mass transport of - * hydrology:*/ - if(!isice && !ishydro){ - if(!masks->isoceanin[this->lid]){ - constant=0; - #ifdef _ISSM_DEBUG_ - this->AddInput(SealevelBarystaticMaskEnum,&constant,P0Enum); - #endif - return; - } - - } - /*if we are an ice shelf, are we fully grounded or not? (used later):*/ - if (masks->notfullygrounded[this->lid]) notfullygrounded=true; //used later on. - /*recover some parameters:*/ @@ -5979,117 +6828,55 @@ this->parameters->FindParam(&rho_water,MaterialsRhoSeawaterEnum); this->parameters->FindParam(&rho_freshwater,MaterialsRhoFreshwaterEnum); - this->parameters->FindParam(&computerigid,SolidearthSettingsRigidEnum); - this->parameters->FindParam(&glfraction,SolidearthSettingsGlfractionEnum); - this->parameters->FindParam(&npartice,SolidearthNpartIceEnum); - this->parameters->FindParam(&nparthydro,SolidearthNpartHydroEnum); - - /*Get area of element: precomputed in the sealevelchange_geometry:*/ - this->Element::GetInputValue(&area,AreaEnum); - - /*Deal with ice loads if we are on grounded ice:*/ - if(masks->isiceonly[this->lid] && !masks->isfullyfloating[this->lid] && isice){ - - /*Compute fraction of the element that is grounded: {{{*/ - if(notfullygrounded){ - IssmDouble xyz_list[NUMVERTICES][3]; - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - - phi_ice=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_ice=1; -#ifdef _ISSM_DEBUG_ - this->AddInput(SealevelBarystaticMaskEnum,&phi_ice,P0Enum); -#endif - } - else phi_ice=1.0; - /*}}}*/ - - /*Inform mask: */ - constant+=100; //1 for ice. - #ifdef _ISSM_DEBUG_ - this->AddInput(SealevelBarystaticMaskEnum,&constant,P0Enum); - #endif - - /*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; - } - /*}}}*/ - - /*Compute barystatic contribution in kg:*/ - bslcice = -rho_ice*phi_ice*area*I; - _assert_(!xIsNan(bslcice)); - - /*Transfer thickness change into kg/m^2:*/ - I=I*rho_ice*phi_ice; - } - - /*Deal with water loads if we are on ground:*/ - if(!masks->isfullyfloating[this->lid] && ishydro){ - - constant+=10; //1 for water. - #ifdef _ISSM_DEBUG_ - this->AddInput(SealevelBarystaticMaskEnum,&constant,P0Enum); - #endif - - /*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 in kg:*/ - bslchydro = -rho_freshwater*area*phi_water*W; - _assert_(!xIsNan(bslchydro)); - - /*convert from m to kg/m^2:*/ - W=W*rho_freshwater*phi_water; - } - - /*Deal with ocean bottom pressures:*/ - if(masks->isoceanin[this->lid] && isocean){ - - constant+=1; //1 for bottom pressure. - #ifdef _ISSM_DEBUG_ - this->AddInput(SealevelBarystaticMaskEnum,&constant,P0Enum); - #endif - - /*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); - - /*Compute barystatic component in kg:*/ - bslcbp = -rho_water*area*phi_bp*BP; - - /*convert from m to kg/m^2:*/ - BP=BP*rho_water*phi_bp; - } - - /*Plug all loads into total load vector:*/ - loads->SetValue(this->sid,I+W+BP,INS_VAL); + + /*Retrieve inputs:*/ + Element::GetInputListOnVertices(&I[0],DeltaIceThicknessEnum); + Element::GetInputListOnVertices(&W[0],DeltaTwsEnum); + Element::GetInputListOnVertices(&BP[0],DeltaBottomPressureEnum); + + for(int i=0;iLoadWeigths[SLGEOM_ICE][i][this->lid]; + Wavg+=W[i]*slgeom->LoadWeigths[SLGEOM_WATER][i][this->lid]; + BPavg+=BP[i]*slgeom->LoadWeigths[SLGEOM_OCEAN][i][this->lid]; + } + + /*convert from m to kg/m^2:*/ + Iavg*=rho_ice; + Wavg*=rho_freshwater; + BPavg*=rho_water; + + #ifdef _ISSM_DEBUG_ + this->AddInput(SealevelBarystaticIceLoadEnum,&Iavg,P0Enum); + this->AddInput(SealevelBarystaticHydroLoadEnum,&Wavg,P0Enum); + this->AddInput(SealevelBarystaticBpLoadEnum,&BPavg,P0Enum); + #endif + + /*Compute barystatic component in kg:*/ + bslcice = -slgeom->LoadArea[SLGEOM_ICE][this->lid]*Iavg; + bslchydro = -slgeom->LoadArea[SLGEOM_WATER][this->lid]*Wavg; + bslcbp = -slgeom->LoadArea[SLGEOM_OCEAN][this->lid]*BPavg; + _assert_(!xIsNan(bslcice)); + _assert_(!xIsNan(bslchydro)); + _assert_(!xIsNan(bslcbp)); + + /*Plug values into subelement load vector:*/ + if(slgeom->issubelement[SLGEOM_ICE][this->lid]){ + int intj=slgeom->subelementmapping[SLGEOM_ICE][this->lid]; + subelementiceloads->SetValue(intj,Iavg,INS_VAL); + Iavg=0; //avoid double counting centroid loads and subelement loads! + } + if(slgeom->issubelement[SLGEOM_WATER][this->lid]){ + int intj=slgeom->subelementmapping[SLGEOM_WATER][this->lid]; + subelementhydroloads->SetValue(intj,Wavg,INS_VAL); + Wavg=0; + } + if(slgeom->issubelement[SLGEOM_OCEAN][this->lid]){ + int intj=slgeom->subelementmapping[SLGEOM_OCEAN][this->lid]; + subelementbploads->SetValue(intj,BPavg,INS_VAL); + BPavg=0; + } + /*Plug remaining values into centroi load vector:*/ + loads->SetValue(this->sid,Iavg+Wavg+BPavg,INS_VAL); + /*Keep track of barystatic contributions:*/ barycontrib->Set(this->Sid(),bslcice,bslchydro,bslcbp); @@ -6097,10 +6884,14 @@ } /*}}}*/ -void Tria::SealevelchangeConvolution(Vector* sealevelloads, Vector* oceanareas, IssmDouble* allsealevelloads, IssmDouble* allloads,IssmDouble* rotationvector){ /*{{{*/ +void Tria::SealevelchangeConvolution(Vector* sealevelloads, Vector* subelementsealevelloads, Vector* oceanareas, Vector* subelementoceanareas, IssmDouble* allsealevelloads, IssmDouble* allloads,IssmDouble* allsubelementiceloads, IssmDouble* allsubelementhydroloads, IssmDouble* allsubelementbploads, IssmDouble* allsubelementsealevelloads, IssmDouble* rotationvector,SealevelGeometry* slgeom){ /*{{{*/ /*sal green function:*/ IssmDouble* G=NULL; + IssmDouble* GsubelIce=NULL; + IssmDouble* GsubelHydro=NULL; + IssmDouble* GsubelOcean=NULL; IssmDouble SealevelGRD[NUMVERTICES]={0}; - IssmDouble oceanaverage,oceanarea=0; + IssmDouble oceanaverage=0; + IssmDouble oceanarea=0; IssmDouble rho_water; @@ -6108,5 +6899,5 @@ bool rotation= false; int size; - int nel; + int nel,nbar; IssmDouble Grotm1[3]; IssmDouble Grotm2[3]; @@ -6128,4 +6919,7 @@ if(sal){ this->inputs->GetArrayPtr(SealevelchangeGEnum,this->lid,&G,&size); + this->inputs->GetArrayPtr(SealevelchangeGsubelIceEnum,this->lid,&GsubelIce,&size); + this->inputs->GetArrayPtr(SealevelchangeGsubelHydroEnum,this->lid,&GsubelHydro,&size); + this->inputs->GetArrayPtr(SealevelchangeGsubelOceanEnum,this->lid,&GsubelOcean,&size); if(allsealevelloads){ @@ -6134,4 +6928,16 @@ SealevelGRD[i]+=G[i*nel+e]*(allsealevelloads[e]+allloads[e]); } + nbar=slgeom->nbar[SLGEOM_ICE]; + for (int e=0;enbar[SLGEOM_WATER]; + for (int e=0;enbar[SLGEOM_OCEAN]; + for (int e=0;enbar[SLGEOM_ICE]; + for (int e=0;enbar[SLGEOM_WATER]; + for (int e=0;enbar[SLGEOM_OCEAN]; + for (int e=0;eLoadWeigths[SLGEOM_OCEAN][i][this->lid]; + } + #ifdef _ISSM_DEBUG_ + this->AddInput(SealevelBarystaticOceanLoadEnum,&oceanaverage,P0Enum); + #endif + oceanarea=slgeom->LoadArea[SLGEOM_OCEAN][this->lid]; + /*add ocean average in the global sealevelloads vector:*/ - sealevelloads->SetValue(this->sid,oceanaverage*rho_water,INS_VAL); - if(!allsealevelloads)oceanareas->SetValue(this->sid,oceanarea,INS_VAL); + if(slgeom->issubelement[SLGEOM_OCEAN][this->lid]){ + int intj=slgeom->subelementmapping[SLGEOM_OCEAN][this->lid]; + subelementsealevelloads->SetValue(intj,oceanaverage*rho_water,INS_VAL); + } + else sealevelloads->SetValue(this->sid,oceanaverage*rho_water,INS_VAL); + + if(!allsealevelloads){ + oceanareas->SetValue(this->sid,oceanarea,INS_VAL); + if(slgeom->issubelement[SLGEOM_OCEAN][this->lid]){ + int intj=slgeom->subelementmapping[SLGEOM_OCEAN][this->lid]; + subelementoceanareas->SetValue(intj,oceanarea,INS_VAL); + } + } + return; } /*}}}*/ -void Tria::SealevelchangeDeformationConvolution(IssmDouble* sealevelloads, IssmDouble* loads, IssmDouble* rotationvector){ /*{{{*/ +void Tria::SealevelchangeDeformationConvolution(IssmDouble* sealevelloads, IssmDouble* subelementsealevelloads, IssmDouble* loads, IssmDouble* subelementiceloads, IssmDouble* subelementhydroloads, IssmDouble* subelementbploads, IssmDouble* rotationvector,SealevelGeometry* slgeom){ /*{{{*/ IssmDouble Sealevel[3]={0,0,0}; @@ -6161,5 +6997,5 @@ IssmDouble SealevelN[3]={0,0,0}; IssmDouble SealevelE[3]={0,0,0}; - int nel; + int nel,nbar; bool sal = false; IssmDouble* G=NULL; @@ -6167,4 +7003,17 @@ IssmDouble* GE=NULL; IssmDouble* GN=NULL; + IssmDouble* GsubelIce=NULL; + IssmDouble* GsubelHydro=NULL; + IssmDouble* GsubelOcean=NULL; + IssmDouble* GUsubelIce=NULL; + IssmDouble* GUsubelHydro=NULL; + IssmDouble* GUsubelOcean=NULL; + IssmDouble* GNsubelIce=NULL; + IssmDouble* GNsubelHydro=NULL; + IssmDouble* GNsubelOcean=NULL; + IssmDouble* GEsubelIce=NULL; + IssmDouble* GEsubelHydro=NULL; + IssmDouble* GEsubelOcean=NULL; + int horiz; int size; @@ -6191,8 +7040,22 @@ if(sal){ this->inputs->GetArrayPtr(SealevelchangeGEnum,this->lid,&G,&size); - if(elastic) this->inputs->GetArrayPtr(SealevelchangeGUEnum,this->lid,&GU,&size); - if(horiz && elastic){ - this->inputs->GetArrayPtr(SealevelchangeGEEnum,this->lid,&GE,&size); - this->inputs->GetArrayPtr(SealevelchangeGNEnum,this->lid,&GN,&size); + this->inputs->GetArrayPtr(SealevelchangeGsubelIceEnum,this->lid,&GsubelIce,&size); + this->inputs->GetArrayPtr(SealevelchangeGsubelHydroEnum,this->lid,&GsubelHydro,&size); + this->inputs->GetArrayPtr(SealevelchangeGsubelOceanEnum,this->lid,&GsubelOcean,&size); + if(elastic){ + this->inputs->GetArrayPtr(SealevelchangeGUEnum,this->lid,&GU,&size); + this->inputs->GetArrayPtr(SealevelchangeGUsubelIceEnum,this->lid,&GUsubelIce,&size); + this->inputs->GetArrayPtr(SealevelchangeGUsubelHydroEnum,this->lid,&GUsubelHydro,&size); + this->inputs->GetArrayPtr(SealevelchangeGUsubelOceanEnum,this->lid,&GUsubelOcean,&size); + if(horiz){ + this->inputs->GetArrayPtr(SealevelchangeGEEnum,this->lid,&GE,&size); + this->inputs->GetArrayPtr(SealevelchangeGEsubelIceEnum,this->lid,&GEsubelIce,&size); + this->inputs->GetArrayPtr(SealevelchangeGEsubelOceanEnum,this->lid,&GEsubelOcean,&size); + this->inputs->GetArrayPtr(SealevelchangeGEsubelHydroEnum,this->lid,&GEsubelHydro,&size); + this->inputs->GetArrayPtr(SealevelchangeGNEnum,this->lid,&GN,&size); + this->inputs->GetArrayPtr(SealevelchangeGNsubelIceEnum,this->lid,&GNsubelIce,&size); + this->inputs->GetArrayPtr(SealevelchangeGNsubelHydroEnum,this->lid,&GNsubelHydro,&size); + this->inputs->GetArrayPtr(SealevelchangeGNsubelOceanEnum,this->lid,&GNsubelOcean,&size); + } } @@ -6201,8 +7064,33 @@ SealevelRSL[i]+=G[i*nel+e]*(sealevelloads[e]+loads[e]); } + nbar=slgeom->nbar[SLGEOM_ICE]; + for (int e=0;enbar[SLGEOM_WATER]; + for (int e=0;enbar[SLGEOM_OCEAN]; + for (int e=0;enbar[SLGEOM_ICE]; + for (int e=0;enbar[SLGEOM_WATER]; + for (int e=0;enbar[SLGEOM_OCEAN]; + for (int e=0;enbar[SLGEOM_ICE]; + for (int e=0;enbar[SLGEOM_WATER]; + for (int e=0;enbar[SLGEOM_OCEAN]; + for (int e=0;eGetInput(levelsetenum); - this->Element::GetInputListOnVertices(&levelsetvalues[0],levelsetenum); - - /*Early return if no vertices are inside the desired area:*/ - if(levelset->GetInputMin()>=0){ - *paverage=0; - *parea=0; - return; - } - - /*Get area of element:*/ - this->Element::GetInputValue(&area,AreaEnum); - - /*Are we fully in the desired area, in which case average over the entire area?:*/ - if(levelset->GetInputMax()<=0){ - for(int i=0;iGetFractionArea(&xyz_list[0][0],levelsetvalues); - area=phi*area; - - /*Average over the fraction area only, using the right gaussian points: */ - this->GetFractionGeometry(&point1,&fraction1,&fraction2,&fractiongeometryflag,levelsetvalues); - Gauss* gauss = this->NewGauss(point1,fraction1,fraction2,fractiongeometryflag,2); - - total_weight=0; - average=0; - while(gauss->next()){ - IssmDouble field_gauss=0; - TriaRef::GetInputValue(&field_gauss, field_on_localvertices, gauss,P1Enum); - average+=field_gauss*gauss->weight; - total_weight+=gauss->weight; - } - if(total_weight) average=average/total_weight; - - /*free ressources:*/ - delete gauss; - - *paverage=average; - *parea=area; - return; - -} -/*}}}*/ -void Tria::SealevelchangeMomentOfInertia(IssmDouble* dI_list, IssmDouble* loads, IssmDouble* sealevelloads){/*{{{*/ +void Tria::SealevelchangeMomentOfInertiaCentroid(IssmDouble* dI_list, IssmDouble* loads, IssmDouble* sealevelloads, SealevelGeometry* slgeom){/*{{{*/ IssmDouble S=0; /*Compute area of element:*/ - IssmDouble area,planetarea; + IssmDouble area,planetarea,re; + IssmDouble late,longe; this->Element::GetInputValue(&area,AreaEnum); - /*recover earth area: */ + /*recover parameters: */ 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; + this->parameters->FindParam(&re,SolidearthPlanetRadiusEnum); + late=slgeom->late[this->lid]/180*M_PI; + longe=slgeom->longe[this->lid]/180*M_PI; @@ -6350,5 +7150,5 @@ if(loads) S+=loads[this->Sid()]; if(sealevelloads) S+=sealevelloads[this->Sid()]; - + /* Perturbation terms for moment of inertia (moi_list): * computed analytically (see Wu & Peltier, eqs 10 & 32) @@ -6361,9 +7161,68 @@ return; }/*}}}*/ -void Tria::SealevelchangeShift(Vector* loads, IssmDouble offset, SealevelMasks* masks){ /*{{{*/ +void Tria::SealevelchangeMomentOfInertiaSubElement(IssmDouble* dI_list, IssmDouble* subelementiceloads, IssmDouble* subelementhydroloads, IssmDouble* subelementbploads, IssmDouble* subelementsealevelloads, SealevelGeometry* slgeom){/*{{{*/ + + IssmDouble SA=0; + IssmDouble* loads=NULL; + IssmDouble* sealevelloads=NULL; + IssmDouble late,longe,re; + int intj; + IssmDouble area; + IssmDouble planetarea; + + /*recover parameters: */ + this->parameters->FindParam(&planetarea,SolidearthPlanetAreaEnum); + this->parameters->FindParam(&re,SolidearthPlanetRadiusEnum); - if(masks->isoceanin[this->lid]){ - loads->SetValue(this->sid,offset,ADD_VAL); - } + /*Initalize:*/ + for(int i=0;i<3;i++)dI_list[i]=0; + + /*Go through our loads:*/ + for(int i=0;iissubelement[i][this->lid]){ + switch(i){ + case SLGEOM_ICE: + loads=subelementiceloads; + break; + case SLGEOM_WATER: + loads=subelementhydroloads; + break; + case SLGEOM_OCEAN: + loads=subelementbploads; + sealevelloads=subelementsealevelloads; + break; + } + intj=slgeom->subelementmapping[i][this->lid]; + late=slgeom->latbarycentre[i][intj]/180*M_PI; + longe=slgeom->longbarycentre[i][intj]/180*M_PI; + area=slgeom->area_subel[i][intj]; + + /*recover total load: */ + if(loads) SA+=loads[intj]*area; + if(sealevelloads) SA+=sealevelloads[intj]*area; + } + } + + /* 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*(SA)*pow(re,4)*(sin(late)*cos(late)*cos(longe))/planetarea; + dI_list[1] += -4*M_PI*(SA)*pow(re,4)*(sin(late)*cos(late)*sin(longe))/planetarea; + dI_list[2] += +4*M_PI*(SA)*pow(re,4)*(1-pow(sin(late),2))/planetarea; + + return; +}/*}}}*/ +void Tria::SealevelchangeShift(Vector* loads, Vector* subelementloads, IssmDouble offset, SealevelGeometry* slgeom){ /*{{{*/ + + if(slgeom->isoceanin[this->lid]){ + if(slgeom->issubelement[SLGEOM_OCEAN][this->lid]){ + int intj=slgeom->subelementmapping[SLGEOM_OCEAN][this->lid]; + subelementloads->SetValue(intj,offset,ADD_VAL); + } + else loads->SetValue(this->sid,offset,ADD_VAL); + } + } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 26221) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 26222) @@ -71,4 +71,5 @@ void CreateDistanceInputFromSegmentlist(IssmDouble* distances,int distanceenum); void ElementCoordinates(Vector* vxe,Vector* vye,Vector* vze, Vector* vareae, bool spherical=false); + void ElementCoordinates(Vector* vlonge,Vector* vlate,Vector* vareae); int EdgeOnBaseIndex(); void EdgeOnBaseIndices(int* pindex1,int* pindex); @@ -84,6 +85,9 @@ void GetGroundedPart(int* point1,IssmDouble* fraction1, IssmDouble* fraction2,bool* mainlyfloating); IssmDouble GetGroundedPortion(IssmDouble* xyz_list); - void GetFractionGeometry(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* pmainlynegative, IssmDouble* gl); - IssmDouble GetFractionArea(IssmDouble* xyz_list, IssmDouble* gl); + void GetFractionGeometry(IssmDouble* weights, IssmDouble* pphi, int* ppoint1,IssmDouble* pfraction1,IssmDouble* pfraction2, bool* ptrapezeisnegative, IssmDouble* gl); + void GetNodalWeightsAndAreaAndCentroidsFromLeveset(IssmDouble* loadweights, IssmDouble* ploadarea, IssmDouble* platbar, IssmDouble* plongbar, IssmDouble late, IssmDouble longe, IssmDouble area, int levelsetenum); + void GetNodalWeightsAndAreaAndCentroidsFromLeveset(IssmDouble* loadweights, IssmDouble* ploadarea, IssmDouble* platbar, IssmDouble* plongbar, IssmDouble late, IssmDouble longe, IssmDouble area, int levelset1enum, int levelset2enum); + void GetBarycenterFromLevelset(IssmDouble* platbar, IssmDouble* plongbar,IssmDouble phi,IssmDouble fraction1,IssmDouble fraction2,IssmDouble late, IssmDouble longe, int point1,int istrapeze1, IssmDouble planetradius); + IssmDouble GetTriangleAreaSpherical(IssmDouble xyz_list[3][3]); IssmDouble GetIcefrontArea(); void GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum); @@ -164,14 +168,15 @@ #endif #ifdef _HAVE_SEALEVELCHANGE_ - void LevelsetAverage(IssmDouble* paverage, IssmDouble* parea, IssmDouble* field_on_localvertices, int levelsetenum); - void SealevelchangeMomentOfInertia(IssmDouble* dI_list, IssmDouble* loads, IssmDouble* sealevelloads); - void SealevelchangeGeometry(IssmDouble* lat, IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze); void GiaDeflection(Vector* wg,Vector* dwgdt,Matlitho* litho, IssmDouble* x,IssmDouble* y); - void SetSealevelMasks(SealevelMasks* masks); - 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); - void SealevelchangeDeformationConvolution(IssmDouble* sealevelloads, IssmDouble* loads, IssmDouble* rotationvector); - void SealevelchangeShift(Vector* loads, IssmDouble offset, SealevelMasks* masks); + void SealevelchangeGeometryInitial(IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae); + void SealevelchangeGeometryFractionKernel(SealevelGeometry* slgeom); + void SealevelchangeGeometrySubElementLoads(SealevelGeometry* slgeom, IssmDouble* areae); + void SealevelchangeGeometryCentroidLoads(SealevelGeometry* slgeom, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae); + void SealevelchangeBarystaticLoads(Vector* loads, Vector* subelementiceloads, Vector* subelementhydroloads, Vector* subelementbploads, BarystaticContributions* barycontrib, SealevelGeometry* slgeom); + void SealevelchangeConvolution(Vector* sealevelloads, Vector* subelementsealevelloads, Vector* oceanareas, Vector* subelementoceanareas, IssmDouble* allsealevelloads, IssmDouble* allloads,IssmDouble* allsubelementiceloads, IssmDouble* allsubelementhydroloads, IssmDouble* allsubelementbploads, IssmDouble* allsubelementsealevelloads, IssmDouble* rotationvector,SealevelGeometry* slgeom); + void SealevelchangeDeformationConvolution(IssmDouble* sealevelloads, IssmDouble* subelementsealevelloads, IssmDouble* loads, IssmDouble* subelementiceloads, IssmDouble* subelementhydroloads, IssmDouble* subelementbploads, IssmDouble* rotationvector,SealevelGeometry* slgeom); + void SealevelchangeShift(Vector* loads, Vector* subelementloads, IssmDouble offset, SealevelGeometry* slgeom); + void SealevelchangeMomentOfInertiaCentroid(IssmDouble* dI_list, IssmDouble* loads, IssmDouble* sealevelloads, SealevelGeometry* slgeom); + void SealevelchangeMomentOfInertiaSubElement(IssmDouble* dI_list, IssmDouble* subelementiceloads, IssmDouble* subelementhydroloads, IssmDouble* subelementbploads, IssmDouble* subelementsealevelloads, SealevelGeometry* slgeom); #endif /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/SealevelGeometry.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/SealevelGeometry.cpp (revision 26222) +++ /issm/trunk-jpl/src/c/classes/SealevelGeometry.cpp (revision 26222) @@ -0,0 +1,136 @@ +/* + * \file SealevelGeometry.cpp + * \brief: Implementation of SealevelGeometry class + */ + +/*Headers: {{{*/ +#ifdef HAVE_CONFIG_H + #include +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +#include "./SealevelGeometry.h" + +using namespace std; +/*}}}*/ + +/*Object constructors and destructor*/ +SealevelGeometry::SealevelGeometry(int localnelin){ /*{{{*/ + localnel=localnelin; + for(int i=0;i(localnel); + } + vlatbarycentre[i]=NULL; //we don't know yet + vlongbarycentre[i]=NULL; + vareae_subel[i]=NULL; + latbarycentre[i]=NULL; //we don't know yet + longbarycentre[i]=NULL; + area_subel[i]=NULL; + + LoadArea[i]=xNewZeroInit(localnel); + issubelement[i]=xNewZeroInit(localnel); + subelementmapping[i]=NULL; + nsubel[i]=0; + nbar[i]=0; + } + late=xNew(localnel); + longe=xNew(localnel); + isoceanin=xNew(localnel); + +}; /*}}}*/ +SealevelGeometry::~SealevelGeometry(){ /*{{{*/ + for(int i=0;i(LoadWeigths[i][j]); + } + xDelete(LoadArea[i]); + xDelete(issubelement[i]); + xDelete(subelementmapping[i]); + delete vlatbarycentre[i]; + delete vlongbarycentre[i]; + delete vareae_subel[i]; + xDelete(latbarycentre[i]); + xDelete(longbarycentre[i]); + xDelete(area_subel[i]); + } + xDelete(late); + xDelete(longe); + xDelete(isoceanin); +}; /*}}}*/ + +void SealevelGeometry::InitializeMappingsAndBarycentres(void){ /*{{{*/ + + int dummy; + bool fromlocalsize=true; + int lower_row; + + + for (int i=0;i(localnel); + GetOwnershipBoundariesFromRange(&lower_row,&dummy,nsubel[i],IssmComm::GetComm()); + + int count=0; + for (int j=0;j(nsubel[i],fromlocalsize); + vlongbarycentre[i]=new Vector(nsubel[i],fromlocalsize); + vareae_subel[i]=new Vector(nsubel[i],fromlocalsize); + vlatbarycentre[i]->GetSize(&nbar[i]); + } + +} /*}}}*/ +void SealevelGeometry::Assemble(void){ /*{{{*/ + + /*Initialize barycentre vectors, now that we know their size: */ + for (int i=0;iAssemble(); + vlongbarycentre[i]->Assemble(); + vareae_subel[i]->Assemble(); + + latbarycentre[i]=vlatbarycentre[i]->ToMPISerial(); + longbarycentre[i]=vlongbarycentre[i]->ToMPISerial(); + area_subel[i]=vareae_subel[i]->ToMPISerial(); + } + + /*Also, we'll need the barycentre associated areas:*/ + + +} /*}}}*/ +int SealevelGeometry::GEnum(int l){ /*{{{*/ + + if(l==SLGEOM_OCEAN) return SealevelchangeGsubelOceanEnum; + if(l==SLGEOM_ICE) return SealevelchangeGsubelIceEnum; + if(l==SLGEOM_WATER) return SealevelchangeGsubelHydroEnum; + +} /*}}}*/ +int SealevelGeometry::GUEnum(int l){ /*{{{*/ + + if(l==SLGEOM_OCEAN) return SealevelchangeGUsubelOceanEnum; + if(l==SLGEOM_ICE) return SealevelchangeGUsubelIceEnum; + if(l==SLGEOM_WATER) return SealevelchangeGUsubelHydroEnum; + +} /*}}}*/ +int SealevelGeometry::GNEnum(int l){ /*{{{*/ + + if(l==SLGEOM_OCEAN) return SealevelchangeGNsubelOceanEnum; + if(l==SLGEOM_ICE) return SealevelchangeGNsubelIceEnum; + if(l==SLGEOM_WATER) return SealevelchangeGNsubelHydroEnum; + +} /*}}}*/ +int SealevelGeometry::GEEnum(int l){ /*{{{*/ + + if(l==SLGEOM_OCEAN) return SealevelchangeGEsubelOceanEnum; + if(l==SLGEOM_ICE) return SealevelchangeGEsubelIceEnum; + if(l==SLGEOM_WATER) return SealevelchangeGEsubelHydroEnum; + +} /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/SealevelGeometry.h =================================================================== --- /issm/trunk-jpl/src/c/classes/SealevelGeometry.h (revision 26222) +++ /issm/trunk-jpl/src/c/classes/SealevelGeometry.h (revision 26222) @@ -0,0 +1,47 @@ +/*!\file SealevelGeometry.h + * \brief: header file for SealevelMask geometry + */ + +#ifndef _SEALEVELGEOMETRY_H_ +#define _SEALEVELGEOMETRY_H_ + +/*Headers:*/ +#define SLGEOM_NUMLOADS 3 +#define SLGEOM_OCEAN 0 +#define SLGEOM_ICE 1 +#define SLGEOM_WATER 2 +#define MAXVERTICES 3 + +#include "../toolkits/toolkits.h" + +class SealevelGeometry{ + + public: + + int localnel; + IssmDouble* LoadWeigths[SLGEOM_NUMLOADS][MAXVERTICES]; + IssmDouble* LoadArea[SLGEOM_NUMLOADS]; + Vector* vlatbarycentre[SLGEOM_NUMLOADS]; + Vector* vlongbarycentre[SLGEOM_NUMLOADS]; + Vector* vareae_subel[SLGEOM_NUMLOADS]; + IssmDouble* latbarycentre[SLGEOM_NUMLOADS]; + IssmDouble* longbarycentre[SLGEOM_NUMLOADS]; + IssmDouble* area_subel[SLGEOM_NUMLOADS]; + IssmDouble* late; + IssmDouble* longe; + bool* isoceanin; + bool* issubelement[SLGEOM_NUMLOADS]; + int* subelementmapping[SLGEOM_NUMLOADS]; + int nsubel[SLGEOM_NUMLOADS]; + int nbar[SLGEOM_NUMLOADS]; + + SealevelGeometry(int localnel); + ~SealevelGeometry(); + void InitializeMappingsAndBarycentres(void); + void Assemble(void); + int GEnum(int l); + int GUEnum(int l); + int GNEnum(int l); + int GEEnum(int l); +}; +#endif /* _SEALEVELGEOMETRY_H_ */ Index: sm/trunk-jpl/src/c/classes/SealevelMasks.h =================================================================== --- /issm/trunk-jpl/src/c/classes/SealevelMasks.h (revision 26221) +++ (revision ) @@ -1,34 +1,0 @@ -/*!\file SealevelMasks.h - * \brief: header file for SealevelMask object - */ - -#ifndef _SEALEVELMASKS_H_ -#define _SEALEVELMASKS_H_ - -/*Headers:*/ - -class SealevelMasks { - - public: - - bool* isiceonly; - bool* isfullyfloating; - bool* notfullygrounded; - bool* isoceanin; - - /*SealevelMasks constructors, destructors :*/ - SealevelMasks(int localnel){ - /*allocate fields:*/ - this->isiceonly=xNew(localnel); - this->isfullyfloating=xNew(localnel); - this->notfullygrounded=xNew(localnel); - this->isoceanin=xNew(localnel); - }; - ~SealevelMasks(){ - xDelete(this->isiceonly); - xDelete(this->isfullyfloating); - xDelete(this->notfullygrounded); - xDelete(this->isoceanin); - }; -}; -#endif /* _SEALEVELMASKS_H_ */ Index: /issm/trunk-jpl/src/c/classes/classes.h =================================================================== --- /issm/trunk-jpl/src/c/classes/classes.h (revision 26221) +++ /issm/trunk-jpl/src/c/classes/classes.h (revision 26222) @@ -18,5 +18,5 @@ #include "./Massfluxatgate.h" #include "./Misfit.h" -#include "./SealevelMasks.h" +#include "./SealevelGeometry.h" #include "./BarystaticContributions.h" #include "./Nodalvalue.h" Index: /issm/trunk-jpl/src/c/cores/WrapperPreCorePointerFromSolutionEnum.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/WrapperPreCorePointerFromSolutionEnum.cpp (revision 26222) +++ /issm/trunk-jpl/src/c/cores/WrapperPreCorePointerFromSolutionEnum.cpp (revision 26222) @@ -0,0 +1,34 @@ +/*!\file: WrapperPreCorePointerFromSolutionEnum.cpp + * \brief: return solution core that is carried out once only for Dakota runs. + */ + +#ifdef HAVE_CONFIG_H + #include +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif +#include "./cores.h" +#include "../toolkits/toolkits.h" +#include "../classes/classes.h" +#include "../shared/shared.h" +#include "../modules/modules.h" +#include "../solutionsequences/solutionsequences.h" + +void WrapperPreCorePointerFromSolutionEnum(void (**psolutioncore)(FemModel*),Parameters* parameters,int solutiontype){ + + /*output: */ + void (*solutioncore)(FemModel*)=NULL; + + switch(solutiontype){ + + case TransientSolutionEnum: + solutioncore=&transient_precore; + break; + default: + break; + } + + /*Assign output pointer:*/ + *psolutioncore=solutioncore; + +} Index: /issm/trunk-jpl/src/c/cores/cores.h =================================================================== --- /issm/trunk-jpl/src/c/cores/cores.h (revision 26221) +++ /issm/trunk-jpl/src/c/cores/cores.h (revision 26222) @@ -9,5 +9,5 @@ class FemModel; class Parameters; -class SealevelMasks; +class SealevelGeometry; template class Matrix; template class Vector; @@ -45,4 +45,5 @@ void steadystate_core(FemModel* femmodel); void transient_core(FemModel* femmodel); +void transient_precore(FemModel* femmodel); void dakota_core(FemModel* femmodel); void ad_core(FemModel* femmodel); @@ -60,7 +61,8 @@ #ifdef _HAVE_SEALEVELCHANGE_ void sealevelchange_core(FemModel* femmodel); -void sealevelchange_geometry(FemModel* femmodel); +void sealevelchange_initialgeometry(FemModel* femmodel); +SealevelGeometry* sealevelchange_geometry(FemModel* femmodel); #endif -void grd_core(FemModel* femmodel); +void grd_core(FemModel* femmodel,SealevelGeometry* slgeom); void solidearthexternal_core(FemModel* femmodel); void dynstr_core(FemModel* femmodel); @@ -82,4 +84,5 @@ void CorePointerFromSolutionEnum(void (**psolutioncore)(FemModel*),Parameters* parameters,int solutiontype); void WrapperCorePointerFromSolutionEnum(void (**psolutioncore)(FemModel*),Parameters* parameters,int solutiontype,bool nodakotacore=false); +void WrapperPreCorePointerFromSolutionEnum(void (**psolutioncore)(FemModel*),Parameters* parameters,int solutiontype); void AdjointCorePointerFromSolutionEnum(void (**padjointcore)(FemModel*),int solutiontype); Index: /issm/trunk-jpl/src/c/cores/sealevelchange_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/sealevelchange_core.cpp (revision 26221) +++ /issm/trunk-jpl/src/c/cores/sealevelchange_core.cpp (revision 26222) @@ -2,4 +2,11 @@ * \brief: core of the sea-level change solution */ + +#ifdef HAVE_CONFIG_H + #include +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + #include "./cores.h" @@ -17,9 +24,9 @@ void TransferSealevel(FemModel* femmodel,int forcingenum); void slcconvergence(bool* pconverged, Vector* RSLg,Vector* RSLg_old,IssmDouble eps_rel,IssmDouble eps_abs); -IssmDouble SealevelloadsOceanAverage(Vector* sealevelloads,Vector* oceanareas, IssmDouble oceanarea); -SealevelMasks* sealevel_masks(FemModel* femmodel); -void RotationAxisMotion(IssmDouble* m, FemModel* femmodel,IssmDouble* loads, IssmDouble* sealeveloads); -void ConserveOceanMass(FemModel* femmodel,Vector* sealevelloads, IssmDouble offset, SealevelMasks* masks); +IssmDouble SealevelloadsOceanAverage(Vector* sealevelloads,Vector* subelementsealevelloads, Vector* oceanareas, Vector* subelementoceanareas, IssmDouble oceanarea); +void RotationAxisMotion(IssmDouble* m, FemModel* femmodel,IssmDouble* loads, IssmDouble* subelementiceloads, IssmDouble* subelementhydroloads, IssmDouble* subelementbploads, IssmDouble* sealevelloads, IssmDouble* subelementsealevelloads, SealevelGeometry* slgeom); +void ConserveOceanMass(FemModel* femmodel,Vector* sealevelloads, Vector* subelementsealevelloads, IssmDouble offset, SealevelGeometry* slgeom); void ivins_deformation_core(FemModel* femmodel); +IssmDouble* CombineLoads(IssmDouble* load,IssmDouble* subload,FemModel* femmodel, SealevelGeometry* slgeom,int loadtype,int nel); /*}}}*/ @@ -27,4 +34,6 @@ void sealevelchange_core(FemModel* femmodel){ /*{{{*/ + SealevelGeometry* slgeom=NULL; + /*Start profiler*/ femmodel->profiler->Start(SLRCORE); @@ -43,5 +52,5 @@ /*run geometry core: */ - sealevelchange_geometry(femmodel); + slgeom=sealevelchange_geometry(femmodel); /*any external forcings?:*/ @@ -52,5 +61,5 @@ /*Run geodetic:*/ - grd_core(femmodel); + grd_core(femmodel,slgeom); /*Run steric core for sure:*/ @@ -206,10 +215,9 @@ }; /*}}}*/ -void grd_core(FemModel* femmodel) { /*{{{*/ +void grd_core(FemModel* femmodel, SealevelGeometry* slgeom) { /*{{{*/ /*variables:{{{*/ int nel; BarystaticContributions* barycontrib=NULL; - SealevelMasks* masks=NULL; GenericParam* barycontribparam=NULL; IssmDouble rotationaxismotionvector[3]={0}; @@ -217,9 +225,18 @@ Vector* loads=NULL; IssmDouble* allloads=NULL; + Vector* subelementiceloads=NULL; + IssmDouble* allsubelementiceloads=NULL; + Vector* subelementhydroloads=NULL; + IssmDouble* allsubelementhydroloads=NULL; + Vector* subelementbploads=NULL; + IssmDouble* allsubelementbploads=NULL; Vector* sealevelloads=NULL; + IssmDouble* allsealevelloads=NULL; + Vector* subelementsealevelloads=NULL; + IssmDouble* allsubelementsealevelloads=NULL; Vector* oldsealevelloads=NULL; - IssmDouble* allsealevelloads=NULL; Vector* oceanareas=NULL; IssmDouble oceanarea; + Vector* subelementoceanareas=NULL; IssmDouble oceanaverage; bool scaleoceanarea=false; @@ -265,5 +282,4 @@ if(modelid!=earthid)return; } - /*branch directly to Ivins deformation core if requested:*/ if(grdmodel==IvinsEnum){ @@ -282,31 +298,44 @@ /*}}}*/ - /*initialize matrices and vectors:*/ + /*initialize loads and sea level loads:*/ femmodel->parameters->FindParam(&nel,MeshNumberofelementsEnum); + loads=new Vector(nel); - oceanareas=new Vector(nel); + subelementiceloads=new Vector(slgeom->nbar[SLGEOM_ICE]); + subelementhydroloads=new Vector(slgeom->nbar[SLGEOM_WATER]); + subelementbploads=new Vector(slgeom->nbar[SLGEOM_OCEAN]); + sealevelloads=new Vector(nel); sealevelloads->Set(0);sealevelloads->Assemble(); - - /*call masks core: */ - masks=sealevel_masks(femmodel); + subelementsealevelloads=new Vector(slgeom->nbar[SLGEOM_OCEAN]); + subelementoceanareas=new Vector(slgeom->nbar[SLGEOM_OCEAN]); + oceanareas=new Vector(nel); + if(VerboseSolution()) _printf0_(" starting GRD convolutions\n"); - + /*buildup loads: */ for(Object* & object : femmodel->elements->objects){ Element* element = xDynamicCast(object); - element->SealevelchangeBarystaticLoads(loads, barycontrib,masks); + element->SealevelchangeBarystaticLoads(loads, subelementiceloads, subelementhydroloads, subelementbploads, barycontrib,slgeom); } loads->Assemble(); + subelementiceloads->Assemble(); + subelementhydroloads->Assemble(); + subelementbploads->Assemble(); + //broadcast loads allloads=loads->ToMPISerial(); + allsubelementiceloads=subelementiceloads->ToMPISerial(); + allsubelementhydroloads=subelementhydroloads->ToMPISerial(); + allsubelementbploads=subelementbploads->ToMPISerial(); //compute rotation axis motion: - RotationAxisMotion(&rotationaxismotionvector[0],femmodel,allloads,NULL); + RotationAxisMotion(&rotationaxismotionvector[0],femmodel,allloads,allsubelementiceloads,allsubelementhydroloads,allsubelementbploads, NULL, NULL,slgeom); /*skip computation of sea level if requested, which means sea level loads should be zeroed */ if(!computesealevel){ allsealevelloads=xNewZeroInit(nel); + allsubelementsealevelloads=xNewZeroInit(slgeom->nbar[SLGEOM_OCEAN]); goto deformation; } @@ -314,5 +343,5 @@ if(VerboseSolution()) _printf0_(" converging GRD convolutions\n"); for(;;){ - + oldsealevelloads=sealevelloads->Duplicate(); sealevelloads->Copy(oldsealevelloads); @@ -320,23 +349,29 @@ for(Object* & object : femmodel->elements->objects){ Element* element = xDynamicCast(object); - element->SealevelchangeConvolution(sealevelloads, oceanareas, allsealevelloads, allloads,rotationaxismotionvector); - } + element->SealevelchangeConvolution(sealevelloads, subelementsealevelloads, oceanareas, subelementoceanareas, allsealevelloads, allloads,allsubelementiceloads, allsubelementhydroloads, allsubelementbploads, allsubelementsealevelloads, rotationaxismotionvector,slgeom); + } + sealevelloads->Assemble(); + subelementsealevelloads->Assemble(); /*compute ocean areas:*/ if(!allsealevelloads){ //first time in the loop - oceanareas->Assemble(); oceanareas->Sum(&oceanarea); _assert_(oceanarea>0.); + oceanareas->Assemble(); + subelementoceanareas->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,subelementsealevelloads, oceanareas,subelementoceanareas, oceanarea); - //Conserve ocean mass: - oceanaverage=SealevelloadsOceanAverage(sealevelloads,oceanareas,oceanarea); - ConserveOceanMass(femmodel,sealevelloads,barycontrib->Total()/oceanarea - oceanaverage,masks); + ConserveOceanMass(femmodel,sealevelloads,subelementsealevelloads,barycontrib->Total()/oceanarea - oceanaverage,slgeom); //broadcast sea level loads allsealevelloads=sealevelloads->ToMPISerial(); + allsubelementsealevelloads=subelementsealevelloads->ToMPISerial(); //compute rotation axis motion: - RotationAxisMotion(&rotationaxismotionvector[0],femmodel,allloads,allsealevelloads); + RotationAxisMotion(&rotationaxismotionvector[0],femmodel,allloads,allsubelementiceloads,allsubelementhydroloads,allsubelementbploads, allsealevelloads,allsubelementsealevelloads,slgeom); //convergence? @@ -356,14 +391,14 @@ for(Object* & object : femmodel->elements->objects){ Element* element = xDynamicCast(object); - element->SealevelchangeDeformationConvolution(allsealevelloads, allloads, rotationaxismotionvector); - } - + element->SealevelchangeDeformationConvolution(allsealevelloads, allsubelementsealevelloads,allloads, allsubelementiceloads, allsubelementhydroloads, allsubelementbploads,rotationaxismotionvector,slgeom); + } + 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: */ + femmodel->inputs->Shift(SealevelGRDEnum,barycontrib->Total()/rho_water/oceanarea- oceanaverage/rho_water); //given that we converged, no need to recompute ocean average + + //cumulate barystatic contributions and save to results: barycontrib->Cumulate(femmodel->parameters); barycontrib->Save(femmodel->results,femmodel->parameters,oceanarea); @@ -518,69 +553,27 @@ //Geometry: -SealevelMasks* sealevel_masks(FemModel* femmodel) { /*{{{*/ - - int grdmodel=0; - - /*early return?:*/ - femmodel->parameters->FindParam(&grdmodel,GrdModelEnum); - if(grdmodel==IvinsEnum) return NULL; - - if(VerboseSolution()) _printf0_(" computing sea level masks\n"); - - /*initialize SealevelMasks structure: */ - SealevelMasks* masks=new SealevelMasks(femmodel->elements->Size()); - - /*go through elements and fill the masks: */ - for(Object* & object : femmodel->elements->objects){ - Element* element=xDynamicCast(object); - element->SetSealevelMasks(masks); - } - - return masks; -}/*}}}*/ -void sealevelchange_geometry(FemModel* femmodel) { /*{{{*/ - - /*Geometry core where we compute indices into tables pre computed in the SealevelRiseAnalysis: */ +void sealevelchange_initialgeometry(FemModel* femmodel) { /*{{{*/ + + /*Geometry core where we compute geometrical kernels and weights:*/ /*parameters: */ - bool spherical=true; - IssmDouble *latitude = NULL; - IssmDouble *longitude = NULL; - IssmDouble *radius = NULL; - IssmDouble *xx = NULL; - IssmDouble *yy = NULL; - IssmDouble *zz = NULL; IssmDouble *xxe = NULL; IssmDouble *yye = NULL; IssmDouble *zze = NULL; IssmDouble* areae = NULL; - - int horiz; - bool geometrydone = false; - int optim; + int nel; int grdmodel=0; - /*retrieve parameters:*/ - femmodel->parameters->FindParam(&horiz,SolidearthSettingsHorizEnum); - femmodel->parameters->FindParam(&geometrydone,SealevelchangeGeometryDoneEnum); femmodel->parameters->FindParam(&grdmodel,GrdModelEnum); + nel=femmodel->elements->NumberOfElements(); /*early return?:*/ if(grdmodel==IvinsEnum) return; - if(geometrydone){ - if(VerboseSolution()) _printf0_(" geometrical offsets have already been computed, skipping \n"); - return; //don't need to run this again. - } - /*Verbose: */ - if(VerboseSolution()) _printf0_(" computing geometrical offsets into precomputed Green tables \n"); - - /*first, recover lat,long and radius vectors from vertices: */ - VertexCoordinatesx(&latitude,&longitude,&radius,femmodel->vertices,spherical); - if(horiz) VertexCoordinatesx(&xx,&yy,&zz,femmodel->vertices); - - /*first, recover x,y,z and areas from elements: */ + if(VerboseSolution()) _printf0_(" computing initial sea level geometrical kernels and weights.\n"); + + /*recover x,y,z and areas from elements: */ ElementCoordinatesx(&xxe,&yye,&zze,&areae,femmodel->elements); @@ -588,24 +581,88 @@ for(Object* & object : femmodel->elements->objects){ Element* element=xDynamicCast(object); - element->SealevelchangeGeometry(latitude,longitude,radius,xx,yy,zz,xxe,yye,zze,areae); - } - - /*Free ressources:*/ - if(horiz){ - xDelete(xx); - xDelete(yy); - xDelete(zz); - } - xDelete(xxe); - xDelete(yye); - xDelete(zze); - xDelete(areae); - xDelete(latitude); - xDelete(longitude); - xDelete(radius); - - /*Record the fact that we ran this module already: */ - femmodel->parameters->SetParam(true,SealevelchangeGeometryDoneEnum); - + element->SealevelchangeGeometryInitial(xxe,yye,zze,areae); + } + + femmodel->parameters->AddObject(new DoubleVecParam(XxeEnum,xxe,nel)); + femmodel->parameters->AddObject(new DoubleVecParam(YyeEnum,yye,nel)); + femmodel->parameters->AddObject(new DoubleVecParam(ZzeEnum,zze,nel)); + femmodel->parameters->AddObject(new DoubleVecParam(AreaeEnum,areae,nel)); + + + #ifdef _ISSM_DEBUG_ + femmodel->results->AddResult(new GenericExternalResult(femmodel->results->Size()+1,XxeEnum,xxe,nel,1,1,1)); + femmodel->results->AddResult(new GenericExternalResult(femmodel->results->Size()+1,YyeEnum,yye,nel,1,1,1)); + femmodel->results->AddResult(new GenericExternalResult(femmodel->results->Size()+1,ZzeEnum,zze,nel,1,1,1)); + femmodel->results->AddResult(new GenericExternalResult(femmodel->results->Size()+1,AreaeEnum,areae,nel,1,1,1)); + #endif + + return; + + +}/*}}}*/ +SealevelGeometry* sealevelchange_geometry(FemModel* femmodel) { /*{{{*/ + + /*Geometry core where we compute updates to the Green function kernels and weights, dependent + * on the evolution of levelsets: */ + + /*parameters: */ + IssmDouble *xxe = NULL; + IssmDouble *yye = NULL; + IssmDouble *zze = NULL; + IssmDouble* areae = NULL; + + int nel; + int grdmodel=0; + SealevelGeometry* slgeom=NULL; + + /*early return?:*/ + femmodel->parameters->FindParam(&grdmodel,GrdModelEnum); + if(grdmodel==IvinsEnum) return NULL; + + /*retrieve parameters:*/ + femmodel->parameters->FindParam(&xxe,&nel,XxeEnum); + femmodel->parameters->FindParam(&yye,&nel,YyeEnum); + femmodel->parameters->FindParam(&zze,&nel,ZzeEnum); + femmodel->parameters->FindParam(&areae,&nel,AreaeEnum); + + /*initialize SealevelMasks structure: */ + slgeom=new SealevelGeometry(femmodel->elements->Size()); + + /*Verbose: */ + if(VerboseSolution()) _printf0_(" computing sea level geometrical kernel and weight updates.\n"); + + + /*Run sealevel geometry routine for elements with full loading:*/ + for(Object* & object : femmodel->elements->objects){ + Element* element=xDynamicCast(object); + element->SealevelchangeGeometryCentroidLoads(slgeom,xxe,yye,zze,areae); + } + + /*Initialize fractional loading mapping: */ + slgeom->InitializeMappingsAndBarycentres(); + + + /*Run sealevel geometry routine for elements with fractional loading:*/ + for(Object* & object : femmodel->elements->objects){ + Element* element=xDynamicCast(object); + element->SealevelchangeGeometrySubElementLoads(slgeom,areae); + } + + /*Assemble barycentres of fraction loading elements:*/ + slgeom->Assemble(); + + /*Create fractional green function kernels: */ + for(Object* & object : femmodel->elements->objects){ + Element* element=xDynamicCast(object); + element->SealevelchangeGeometryFractionKernel(slgeom); + } + + + femmodel->parameters->AddObject(new DoubleVecParam(XxeEnum,xxe,nel)); + femmodel->parameters->AddObject(new DoubleVecParam(YyeEnum,yye,nel)); + femmodel->parameters->AddObject(new DoubleVecParam(ZzeEnum,zze,nel)); + femmodel->parameters->AddObject(new DoubleVecParam(AreaeEnum,areae,nel)); + + return slgeom; }/*}}}*/ @@ -903,18 +960,27 @@ } /*}}}*/ -IssmDouble SealevelloadsOceanAverage(Vector* sealevelloads,Vector* oceanareas, IssmDouble oceanarea){ /*{{{*/ +IssmDouble SealevelloadsOceanAverage(Vector* sealevelloads,Vector* subelementsealevelloads, Vector* oceanareas, Vector* subelementoceanareas, IssmDouble oceanarea){ /*{{{*/ IssmDouble sealevelloadsaverage; + IssmDouble subelementsealevelloadsaverage; Vector* sealevelloadsvolume=sealevelloads->Duplicate(); + Vector* subelementsealevelloadsvolume=subelementsealevelloads->Duplicate(); + sealevelloadsvolume->PointwiseMult(sealevelloads,oceanareas); + subelementsealevelloadsvolume->PointwiseMult(subelementsealevelloads,subelementoceanareas); + sealevelloadsvolume->Sum(&sealevelloadsaverage); + subelementsealevelloadsvolume->Sum(&subelementsealevelloadsaverage); delete sealevelloadsvolume; - - return sealevelloadsaverage/oceanarea; + delete subelementsealevelloadsvolume; + + //return (sealevelloadsaverage+subelementsealevelloadsaverage)/oceanarea; + return (sealevelloadsaverage)/oceanarea; } /*}}}*/ -void RotationAxisMotion(IssmDouble* m, FemModel* femmodel,IssmDouble* loads, IssmDouble* sealevelloads){ /*{{{*/ +void RotationAxisMotion(IssmDouble* m, FemModel* femmodel,IssmDouble* loads, IssmDouble* subelementiceloads, IssmDouble* subelementhydroloads, IssmDouble* subelementbploads, IssmDouble* sealevelloads, IssmDouble* subelementsealevelloads, SealevelGeometry* slgeom){ /*{{{*/ IssmDouble moi_list[3]={0,0,0}; + IssmDouble moi_list_sub[3]={0,0,0}; IssmDouble moi_list_cpu[3]={0,0,0}; IssmDouble* tide_love_h = NULL; @@ -940,9 +1006,14 @@ for(Object* & object : femmodel->elements->objects){ Element* element = xDynamicCast(object); - element->SealevelchangeMomentOfInertia(&moi_list[0],loads,sealevelloads); - moi_list_cpu[0] += moi_list[0]; - moi_list_cpu[1] += moi_list[1]; - moi_list_cpu[2] += moi_list[2]; - } + + element->SealevelchangeMomentOfInertiaCentroid(&moi_list[0],loads,sealevelloads,slgeom); + + element->SealevelchangeMomentOfInertiaSubElement(&moi_list_sub[0],subelementiceloads, subelementhydroloads, subelementbploads, subelementsealevelloads,slgeom); + + moi_list_cpu[0] += moi_list[0]+moi_list_sub[0]; + moi_list_cpu[1] += moi_list[1]+moi_list_sub[1]; + moi_list_cpu[2] += moi_list[2]+moi_list_sub[2]; + } + ISSM_MPI_Reduce (&moi_list_cpu[0],&moi_list[0],1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); @@ -965,12 +1036,42 @@ m[2]=m3; } /*}}}*/ -void ConserveOceanMass(FemModel* femmodel,Vector* sealevelloads, IssmDouble offset, SealevelMasks* masks){ /*{{{*/ +void ConserveOceanMass(FemModel* femmodel,Vector* sealevelloads, Vector* subelementsealevelloads, IssmDouble offset, SealevelGeometry* slgeom){ /*{{{*/ /*Shift sealevel loads by ocean average, only on ocean! :*/ for(Object* & object : femmodel->elements->objects){ Element* element = xDynamicCast(object); - element->SealevelchangeShift(sealevelloads,offset,masks); + element->SealevelchangeShift(sealevelloads,subelementsealevelloads,offset,slgeom); } sealevelloads->Assemble(); + subelementsealevelloads->Assemble(); } /*}}}*/ +IssmDouble* CombineLoads(IssmDouble* load,IssmDouble* subload,FemModel* femmodel, SealevelGeometry* slgeom,int loadtype,int nel){ /*{{{*/ + + int* indices=xNew(nel); + for(int i=0;i* vloadcopy=new Vector(nel); + IssmDouble* loadcopy=xNew(nel); + + vloadcopy->SetValues(nel,indices,load,INS_VAL); + vloadcopy->Assemble(); + + + if(subload){ + for (int i=0;ielements->Size();i++){ + if (slgeom->issubelement[loadtype][i]){ + int se= slgeom->subelementmapping[loadtype][i]; + IssmDouble subloadi=subload[se]; + Element* element=dynamic_cast(femmodel->elements->GetObjectByOffset(i)); + vloadcopy->SetValue(element->Sid(),subloadi,ADD_VAL); + } + } + } + vloadcopy->Assemble(); + loadcopy=vloadcopy->ToMPISerial(); + + return loadcopy; + +} /*}}}*/ + Index: /issm/trunk-jpl/src/c/cores/transient_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/transient_core.cpp (revision 26221) +++ /issm/trunk-jpl/src/c/cores/transient_core.cpp (revision 26222) @@ -24,8 +24,8 @@ /*parameters: */ IssmDouble finaltime,dt,yts; - bool isoceancoupling,iscontrol,isautodiff; + bool iscontrol,isautodiff; int timestepping; int output_frequency,checkpoint_frequency; - int amr_frequency,amr_restart; + int amr_frequency; char **requested_outputs = NULL; @@ -45,19 +45,10 @@ femmodel->parameters->FindParam(&output_frequency,SettingsOutputFrequencyEnum); femmodel->parameters->FindParam(×tepping,TimesteppingTypeEnum); - femmodel->parameters->FindParam(&isoceancoupling,TransientIsoceancouplingEnum); femmodel->parameters->FindParam(&amr_frequency,TransientAmrFrequencyEnum); femmodel->parameters->FindParam(&iscontrol,InversionIscontrolEnum); femmodel->parameters->FindParam(&isautodiff,AutodiffIsautodiffEnum); - #if defined(_HAVE_BAMG_) && !defined(_HAVE_AD_) - if(amr_frequency){ - femmodel->parameters->FindParam(&amr_restart,AmrRestartEnum); - if(amr_restart) femmodel->ReMesh(); - } - #endif - - #if defined(_HAVE_OCEAN_ ) - if(isoceancoupling) OceanExchangeDatax(femmodel,true); - #endif + /*call modules that are not dependent on time stepping:*/ + transient_precore(femmodel); while(time < finaltime - (yts*DBL_EPSILON)){ //make sure we run up to finaltime. @@ -253,4 +244,26 @@ } }/*}}}*/ +void transient_precore(FemModel* femmodel){/*{{{*/ + + bool isoceancoupling,isslc; + int amr_frequency,amr_restart; + + femmodel->parameters->FindParam(&isoceancoupling,TransientIsoceancouplingEnum); + femmodel->parameters->FindParam(&amr_frequency,TransientAmrFrequencyEnum); + femmodel->parameters->FindParam(&isslc,TransientIsslcEnum); + + #if defined(_HAVE_BAMG_) && !defined(_HAVE_AD_) + if(amr_frequency){ + femmodel->parameters->FindParam(&amr_restart,AmrRestartEnum); + if(amr_restart) femmodel->ReMesh(); + } + #endif + + #if defined(_HAVE_OCEAN_ ) + if(isoceancoupling) OceanExchangeDatax(femmodel,true); + #endif + + if(isslc) sealevelchange_initialgeometry(femmodel); +}/*}}}*/ #ifdef _HAVE_CODIPACK_ Index: /issm/trunk-jpl/src/c/modules/ElementCoordinatesx/ElementCoordinatesx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ElementCoordinatesx/ElementCoordinatesx.cpp (revision 26221) +++ /issm/trunk-jpl/src/c/modules/ElementCoordinatesx/ElementCoordinatesx.cpp (revision 26222) @@ -9,5 +9,5 @@ #include "../../toolkits/toolkits.h" -void ElementCoordinatesx( IssmDouble** pxe, IssmDouble** pye, IssmDouble** pze, IssmDouble** pareae, Elements* elements,bool spherical) { +void ElementCoordinatesx( IssmDouble** pxe, IssmDouble** pye, IssmDouble** pze, IssmDouble** pareae, Elements* elements,bool spherical) { /*{{{*/ /*figure out how many vertices we have: */ @@ -53,3 +53,42 @@ else xDelete(areae); -} +} /*}}}*/ +void ElementCoordinatesx( IssmDouble** plonge, IssmDouble** plate, IssmDouble** pareae, Elements* elements) { /*{{{*/ + + /*figure out how many vertices we have: */ + int numberofelements=elements->NumberOfElements(); + + Vector* vlonge=new Vector(numberofelements); + Vector* vlate=new Vector(numberofelements); + Vector* vareae=new Vector(numberofelements); + + /*march through our elements: */ + for(Object* & object : elements->objects){ + Element* element=(Element*)object; + element->ElementCoordinates(vlonge,vlate,vareae); + } + + /*Assemble*/ + vlonge->Assemble(); + vlate->Assemble(); + vareae->Assemble(); + + /*serialize: */ + IssmDouble* longe=vlonge->ToMPISerial(); + IssmDouble* late=vlate->ToMPISerial(); + IssmDouble* areae=vareae->ToMPISerial(); + + /*Free ressources: */ + delete vlonge; + delete vlate; + delete vareae; + + /*output: */ + if(plonge) *plonge=longe; + else xDelete(longe); + if(plate) *plate=late; + else xDelete(late); + if(pareae) *pareae=areae; + else xDelete(areae); + +} /*}}}*/ Index: /issm/trunk-jpl/src/c/modules/ElementCoordinatesx/ElementCoordinatesx.h =================================================================== --- /issm/trunk-jpl/src/c/modules/ElementCoordinatesx/ElementCoordinatesx.h (revision 26221) +++ /issm/trunk-jpl/src/c/modules/ElementCoordinatesx/ElementCoordinatesx.h (revision 26222) @@ -9,4 +9,5 @@ /* local prototypes: */ void ElementCoordinatesx( IssmDouble** pxe, IssmDouble** pye, IssmDouble** pze,IssmDouble** pareae, Elements* elements,bool spherical=false); +void ElementCoordinatesx( IssmDouble** plonge, IssmDouble** plate, IssmDouble** pareae, Elements* elements); #endif /* _ELEMENT_COORDINATESX_H */ Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 26221) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 26222) @@ -331,4 +331,5 @@ parameters->AddObject(iomodel->CopyConstantObject("md.materials.thermal_exchange_velocity",MaterialsThermalExchangeVelocityEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.constants.g",ConstantsGEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.constants.gravitational_constant",ConstantsNewtonGravityEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.materials.rheology_law",MaterialsRheologyLawEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.materials.earth_density",MaterialsEarthDensityEnum)); @@ -373,4 +374,5 @@ } parameters->AddObject(iomodel->CopyConstantObject("md.materials.earth_density",MaterialsEarthDensityEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.constants.gravitational_constant",ConstantsNewtonGravityEnum)); /*Free rssources:*/ xDelete(nature); Index: /issm/trunk-jpl/src/c/shared/Enum/Enum.vim =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 26221) +++ /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 26222) @@ -106,4 +106,5 @@ syn keyword cConstant ConfigurationTypeEnum syn keyword cConstant ConstantsGEnum +syn keyword cConstant ConstantsNewtonGravityEnum syn keyword cConstant ConstantsReferencetemperatureEnum syn keyword cConstant ConstantsYtsEnum @@ -365,5 +366,4 @@ syn keyword cConstant SolidearthSettingsComputesealevelchangeEnum syn keyword cConstant SolidearthSettingsGRDEnum -syn keyword cConstant SolidearthSettingsGlfractionEnum syn keyword cConstant SolidearthSettingsRunFrequencyEnum syn keyword cConstant SealevelchangeHElasticEnum @@ -506,4 +506,8 @@ syn keyword cConstant TransientRequestedOutputsEnum syn keyword cConstant VelocityEnum +syn keyword cConstant XxeEnum +syn keyword cConstant YyeEnum +syn keyword cConstant ZzeEnum +syn keyword cConstant AreaeEnum syn keyword cConstant WorldCommEnum syn keyword cConstant ParametersENDEnum @@ -745,5 +749,26 @@ syn keyword cConstant SealevelGRDEnum syn keyword cConstant SealevelBarystaticMaskEnum +syn keyword cConstant SealevelBarystaticIceMaskEnum +syn keyword cConstant SealevelBarystaticIceWeightsEnum +syn keyword cConstant SealevelBarystaticIceAreaEnum +syn keyword cConstant SealevelBarystaticIceLatbarEnum +syn keyword cConstant SealevelBarystaticIceLongbarEnum +syn keyword cConstant SealevelBarystaticIceLoadEnum +syn keyword cConstant SealevelBarystaticHydroMaskEnum +syn keyword cConstant SealevelBarystaticHydroWeightsEnum +syn keyword cConstant SealevelBarystaticHydroAreaEnum +syn keyword cConstant SealevelBarystaticHydroLatbarEnum +syn keyword cConstant SealevelBarystaticHydroLongbarEnum +syn keyword cConstant SealevelBarystaticHydroLoadEnum +syn keyword cConstant SealevelBarystaticBpMaskEnum +syn keyword cConstant SealevelBarystaticBpWeightsEnum +syn keyword cConstant SealevelBarystaticBpAreaEnum +syn keyword cConstant SealevelBarystaticBpLoadEnum syn keyword cConstant SealevelBarystaticOceanMaskEnum +syn keyword cConstant SealevelBarystaticOceanWeightsEnum +syn keyword cConstant SealevelBarystaticOceanAreaEnum +syn keyword cConstant SealevelBarystaticOceanLatbarEnum +syn keyword cConstant SealevelBarystaticOceanLongbarEnum +syn keyword cConstant SealevelBarystaticOceanLoadEnum syn keyword cConstant SealevelNEsaEnum syn keyword cConstant SealevelNEsaRateEnum @@ -769,4 +794,16 @@ syn keyword cConstant SealevelchangeGEEnum syn keyword cConstant SealevelchangeGNEnum +syn keyword cConstant SealevelchangeGsubelOceanEnum +syn keyword cConstant SealevelchangeGUsubelOceanEnum +syn keyword cConstant SealevelchangeGEsubelOceanEnum +syn keyword cConstant SealevelchangeGNsubelOceanEnum +syn keyword cConstant SealevelchangeGsubelIceEnum +syn keyword cConstant SealevelchangeGUsubelIceEnum +syn keyword cConstant SealevelchangeGEsubelIceEnum +syn keyword cConstant SealevelchangeGNsubelIceEnum +syn keyword cConstant SealevelchangeGsubelHydroEnum +syn keyword cConstant SealevelchangeGUsubelHydroEnum +syn keyword cConstant SealevelchangeGEsubelHydroEnum +syn keyword cConstant SealevelchangeGNsubelHydroEnum syn keyword cConstant SedimentHeadEnum syn keyword cConstant SedimentHeadOldEnum @@ -1540,5 +1577,5 @@ syn keyword cType RiftStruct syn keyword cType Riftfront -syn keyword cType SealevelMasks +syn keyword cType SealevelGeometry syn keyword cType Seg syn keyword cType SegInput Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 26221) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 26222) @@ -100,4 +100,5 @@ ConfigurationTypeEnum, ConstantsGEnum, + ConstantsNewtonGravityEnum, ConstantsReferencetemperatureEnum, ConstantsYtsEnum, @@ -359,5 +360,4 @@ SolidearthSettingsComputesealevelchangeEnum, SolidearthSettingsGRDEnum, - SolidearthSettingsGlfractionEnum, SolidearthSettingsRunFrequencyEnum, SealevelchangeHElasticEnum, @@ -500,4 +500,8 @@ TransientRequestedOutputsEnum, VelocityEnum, + XxeEnum, + YyeEnum, + ZzeEnum, + AreaeEnum, WorldCommEnum, /*}}}*/ @@ -741,5 +745,26 @@ SealevelGRDEnum, SealevelBarystaticMaskEnum, + SealevelBarystaticIceMaskEnum, + SealevelBarystaticIceWeightsEnum, + SealevelBarystaticIceAreaEnum, + SealevelBarystaticIceLatbarEnum, + SealevelBarystaticIceLongbarEnum, + SealevelBarystaticIceLoadEnum, + SealevelBarystaticHydroMaskEnum, + SealevelBarystaticHydroWeightsEnum, + SealevelBarystaticHydroAreaEnum, + SealevelBarystaticHydroLatbarEnum, + SealevelBarystaticHydroLongbarEnum, + SealevelBarystaticHydroLoadEnum, + SealevelBarystaticBpMaskEnum, + SealevelBarystaticBpWeightsEnum, + SealevelBarystaticBpAreaEnum, + SealevelBarystaticBpLoadEnum, SealevelBarystaticOceanMaskEnum, + SealevelBarystaticOceanWeightsEnum, + SealevelBarystaticOceanAreaEnum, + SealevelBarystaticOceanLatbarEnum, + SealevelBarystaticOceanLongbarEnum, + SealevelBarystaticOceanLoadEnum, SealevelNEsaEnum, SealevelNEsaRateEnum, @@ -765,4 +790,16 @@ SealevelchangeGEEnum, SealevelchangeGNEnum, + SealevelchangeGsubelOceanEnum, + SealevelchangeGUsubelOceanEnum, + SealevelchangeGEsubelOceanEnum, + SealevelchangeGNsubelOceanEnum, + SealevelchangeGsubelIceEnum, + SealevelchangeGUsubelIceEnum, + SealevelchangeGEsubelIceEnum, + SealevelchangeGNsubelIceEnum, + SealevelchangeGsubelHydroEnum, + SealevelchangeGUsubelHydroEnum, + SealevelchangeGEsubelHydroEnum, + SealevelchangeGNsubelHydroEnum, SedimentHeadEnum, SedimentHeadOldEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 26221) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 26222) @@ -108,4 +108,5 @@ case ConfigurationTypeEnum : return "ConfigurationType"; case ConstantsGEnum : return "ConstantsG"; + case ConstantsNewtonGravityEnum : return "ConstantsNewtonGravity"; case ConstantsReferencetemperatureEnum : return "ConstantsReferencetemperature"; case ConstantsYtsEnum : return "ConstantsYts"; @@ -367,5 +368,4 @@ case SolidearthSettingsComputesealevelchangeEnum : return "SolidearthSettingsComputesealevelchange"; case SolidearthSettingsGRDEnum : return "SolidearthSettingsGRD"; - case SolidearthSettingsGlfractionEnum : return "SolidearthSettingsGlfraction"; case SolidearthSettingsRunFrequencyEnum : return "SolidearthSettingsRunFrequency"; case SealevelchangeHElasticEnum : return "SealevelchangeHElastic"; @@ -508,4 +508,8 @@ case TransientRequestedOutputsEnum : return "TransientRequestedOutputs"; case VelocityEnum : return "Velocity"; + case XxeEnum : return "Xxe"; + case YyeEnum : return "Yye"; + case ZzeEnum : return "Zze"; + case AreaeEnum : return "Areae"; case WorldCommEnum : return "WorldComm"; case ParametersENDEnum : return "ParametersEND"; @@ -747,5 +751,26 @@ case SealevelGRDEnum : return "SealevelGRD"; case SealevelBarystaticMaskEnum : return "SealevelBarystaticMask"; + case SealevelBarystaticIceMaskEnum : return "SealevelBarystaticIceMask"; + case SealevelBarystaticIceWeightsEnum : return "SealevelBarystaticIceWeights"; + case SealevelBarystaticIceAreaEnum : return "SealevelBarystaticIceArea"; + case SealevelBarystaticIceLatbarEnum : return "SealevelBarystaticIceLatbar"; + case SealevelBarystaticIceLongbarEnum : return "SealevelBarystaticIceLongbar"; + case SealevelBarystaticIceLoadEnum : return "SealevelBarystaticIceLoad"; + case SealevelBarystaticHydroMaskEnum : return "SealevelBarystaticHydroMask"; + case SealevelBarystaticHydroWeightsEnum : return "SealevelBarystaticHydroWeights"; + case SealevelBarystaticHydroAreaEnum : return "SealevelBarystaticHydroArea"; + case SealevelBarystaticHydroLatbarEnum : return "SealevelBarystaticHydroLatbar"; + case SealevelBarystaticHydroLongbarEnum : return "SealevelBarystaticHydroLongbar"; + case SealevelBarystaticHydroLoadEnum : return "SealevelBarystaticHydroLoad"; + case SealevelBarystaticBpMaskEnum : return "SealevelBarystaticBpMask"; + case SealevelBarystaticBpWeightsEnum : return "SealevelBarystaticBpWeights"; + case SealevelBarystaticBpAreaEnum : return "SealevelBarystaticBpArea"; + case SealevelBarystaticBpLoadEnum : return "SealevelBarystaticBpLoad"; case SealevelBarystaticOceanMaskEnum : return "SealevelBarystaticOceanMask"; + case SealevelBarystaticOceanWeightsEnum : return "SealevelBarystaticOceanWeights"; + case SealevelBarystaticOceanAreaEnum : return "SealevelBarystaticOceanArea"; + case SealevelBarystaticOceanLatbarEnum : return "SealevelBarystaticOceanLatbar"; + case SealevelBarystaticOceanLongbarEnum : return "SealevelBarystaticOceanLongbar"; + case SealevelBarystaticOceanLoadEnum : return "SealevelBarystaticOceanLoad"; case SealevelNEsaEnum : return "SealevelNEsa"; case SealevelNEsaRateEnum : return "SealevelNEsaRate"; @@ -771,4 +796,16 @@ case SealevelchangeGEEnum : return "SealevelchangeGE"; case SealevelchangeGNEnum : return "SealevelchangeGN"; + case SealevelchangeGsubelOceanEnum : return "SealevelchangeGsubelOcean"; + case SealevelchangeGUsubelOceanEnum : return "SealevelchangeGUsubelOcean"; + case SealevelchangeGEsubelOceanEnum : return "SealevelchangeGEsubelOcean"; + case SealevelchangeGNsubelOceanEnum : return "SealevelchangeGNsubelOcean"; + case SealevelchangeGsubelIceEnum : return "SealevelchangeGsubelIce"; + case SealevelchangeGUsubelIceEnum : return "SealevelchangeGUsubelIce"; + case SealevelchangeGEsubelIceEnum : return "SealevelchangeGEsubelIce"; + case SealevelchangeGNsubelIceEnum : return "SealevelchangeGNsubelIce"; + case SealevelchangeGsubelHydroEnum : return "SealevelchangeGsubelHydro"; + case SealevelchangeGUsubelHydroEnum : return "SealevelchangeGUsubelHydro"; + case SealevelchangeGEsubelHydroEnum : return "SealevelchangeGEsubelHydro"; + case SealevelchangeGNsubelHydroEnum : return "SealevelchangeGNsubelHydro"; case SedimentHeadEnum : return "SedimentHead"; case SedimentHeadOldEnum : return "SedimentHeadOld"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 26221) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 26222) @@ -108,4 +108,5 @@ else if (strcmp(name,"ConfigurationType")==0) return ConfigurationTypeEnum; else if (strcmp(name,"ConstantsG")==0) return ConstantsGEnum; + else if (strcmp(name,"ConstantsNewtonGravity")==0) return ConstantsNewtonGravityEnum; else if (strcmp(name,"ConstantsReferencetemperature")==0) return ConstantsReferencetemperatureEnum; else if (strcmp(name,"ConstantsYts")==0) return ConstantsYtsEnum; @@ -136,9 +137,9 @@ else if (strcmp(name,"DamageStabilization")==0) return DamageStabilizationEnum; else if (strcmp(name,"DamageStressThreshold")==0) return DamageStressThresholdEnum; - else if (strcmp(name,"DamageStressUBound")==0) return DamageStressUBoundEnum; else stage=2; } if(stage==2){ - if (strcmp(name,"DebugProfiling")==0) return DebugProfilingEnum; + if (strcmp(name,"DamageStressUBound")==0) return DamageStressUBoundEnum; + else if (strcmp(name,"DebugProfiling")==0) return DebugProfilingEnum; else if (strcmp(name,"DomainDimension")==0) return DomainDimensionEnum; else if (strcmp(name,"DomainType")==0) return DomainTypeEnum; @@ -259,9 +260,9 @@ else if (strcmp(name,"LevelsetStabilization")==0) return LevelsetStabilizationEnum; else if (strcmp(name,"LockFileName")==0) return LockFileNameEnum; - else if (strcmp(name,"LoveAllowLayerDeletion")==0) return LoveAllowLayerDeletionEnum; else stage=3; } if(stage==3){ - if (strcmp(name,"LoveForcingType")==0) return LoveForcingTypeEnum; + if (strcmp(name,"LoveAllowLayerDeletion")==0) return LoveAllowLayerDeletionEnum; + else if (strcmp(name,"LoveForcingType")==0) return LoveForcingTypeEnum; else if (strcmp(name,"LoveFrequencies")==0) return LoveFrequenciesEnum; else if (strcmp(name,"LoveG0")==0) return LoveG0Enum; @@ -373,5 +374,4 @@ else if (strcmp(name,"SolidearthSettingsComputesealevelchange")==0) return SolidearthSettingsComputesealevelchangeEnum; else if (strcmp(name,"SolidearthSettingsGRD")==0) return SolidearthSettingsGRDEnum; - else if (strcmp(name,"SolidearthSettingsGlfraction")==0) return SolidearthSettingsGlfractionEnum; else if (strcmp(name,"SolidearthSettingsRunFrequency")==0) return SolidearthSettingsRunFrequencyEnum; else if (strcmp(name,"SealevelchangeHElastic")==0) return SealevelchangeHElasticEnum; @@ -520,4 +520,8 @@ else if (strcmp(name,"TransientRequestedOutputs")==0) return TransientRequestedOutputsEnum; else if (strcmp(name,"Velocity")==0) return VelocityEnum; + else if (strcmp(name,"Xxe")==0) return XxeEnum; + else if (strcmp(name,"Yye")==0) return YyeEnum; + else if (strcmp(name,"Zze")==0) return ZzeEnum; + else if (strcmp(name,"Areae")==0) return AreaeEnum; else if (strcmp(name,"WorldComm")==0) return WorldCommEnum; else if (strcmp(name,"ParametersEND")==0) return ParametersENDEnum; @@ -625,12 +629,12 @@ else if (strcmp(name,"EffectivePressureTransient")==0) return EffectivePressureTransientEnum; else if (strcmp(name,"Enthalpy")==0) return EnthalpyEnum; - else if (strcmp(name,"EnthalpyPicard")==0) return EnthalpyPicardEnum; + else stage=6; + } + if(stage==6){ + if (strcmp(name,"EnthalpyPicard")==0) return EnthalpyPicardEnum; 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,"EplHeadSlopeY")==0) return EplHeadSlopeYEnum; + else if (strcmp(name,"EplHeadSlopeY")==0) return EplHeadSlopeYEnum; else if (strcmp(name,"EplHeadSubstep")==0) return EplHeadSubstepEnum; else if (strcmp(name,"EplHeadTransient")==0) return EplHeadTransientEnum; @@ -748,12 +752,12 @@ else if (strcmp(name,"Radar")==0) return RadarEnum; else if (strcmp(name,"RadarAttenuationMacGregor")==0) return RadarAttenuationMacGregorEnum; - else if (strcmp(name,"RadarAttenuationWolff")==0) return RadarAttenuationWolffEnum; + else stage=7; + } + if(stage==7){ + if (strcmp(name,"RadarAttenuationWolff")==0) return RadarAttenuationWolffEnum; 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,"RheologyBAbsGradient")==0) return RheologyBAbsGradientEnum; + else if (strcmp(name,"RheologyBAbsGradient")==0) return RheologyBAbsGradientEnum; else if (strcmp(name,"RheologyBInitialguess")==0) return RheologyBInitialguessEnum; else if (strcmp(name,"RheologyBInitialguessMisfit")==0) return RheologyBInitialguessMisfitEnum; @@ -765,5 +769,26 @@ else if (strcmp(name,"SealevelGRD")==0) return SealevelGRDEnum; else if (strcmp(name,"SealevelBarystaticMask")==0) return SealevelBarystaticMaskEnum; + else if (strcmp(name,"SealevelBarystaticIceMask")==0) return SealevelBarystaticIceMaskEnum; + else if (strcmp(name,"SealevelBarystaticIceWeights")==0) return SealevelBarystaticIceWeightsEnum; + else if (strcmp(name,"SealevelBarystaticIceArea")==0) return SealevelBarystaticIceAreaEnum; + else if (strcmp(name,"SealevelBarystaticIceLatbar")==0) return SealevelBarystaticIceLatbarEnum; + else if (strcmp(name,"SealevelBarystaticIceLongbar")==0) return SealevelBarystaticIceLongbarEnum; + else if (strcmp(name,"SealevelBarystaticIceLoad")==0) return SealevelBarystaticIceLoadEnum; + else if (strcmp(name,"SealevelBarystaticHydroMask")==0) return SealevelBarystaticHydroMaskEnum; + else if (strcmp(name,"SealevelBarystaticHydroWeights")==0) return SealevelBarystaticHydroWeightsEnum; + else if (strcmp(name,"SealevelBarystaticHydroArea")==0) return SealevelBarystaticHydroAreaEnum; + else if (strcmp(name,"SealevelBarystaticHydroLatbar")==0) return SealevelBarystaticHydroLatbarEnum; + else if (strcmp(name,"SealevelBarystaticHydroLongbar")==0) return SealevelBarystaticHydroLongbarEnum; + else if (strcmp(name,"SealevelBarystaticHydroLoad")==0) return SealevelBarystaticHydroLoadEnum; + else if (strcmp(name,"SealevelBarystaticBpMask")==0) return SealevelBarystaticBpMaskEnum; + else if (strcmp(name,"SealevelBarystaticBpWeights")==0) return SealevelBarystaticBpWeightsEnum; + else if (strcmp(name,"SealevelBarystaticBpArea")==0) return SealevelBarystaticBpAreaEnum; + else if (strcmp(name,"SealevelBarystaticBpLoad")==0) return SealevelBarystaticBpLoadEnum; else if (strcmp(name,"SealevelBarystaticOceanMask")==0) return SealevelBarystaticOceanMaskEnum; + else if (strcmp(name,"SealevelBarystaticOceanWeights")==0) return SealevelBarystaticOceanWeightsEnum; + else if (strcmp(name,"SealevelBarystaticOceanArea")==0) return SealevelBarystaticOceanAreaEnum; + else if (strcmp(name,"SealevelBarystaticOceanLatbar")==0) return SealevelBarystaticOceanLatbarEnum; + else if (strcmp(name,"SealevelBarystaticOceanLongbar")==0) return SealevelBarystaticOceanLongbarEnum; + else if (strcmp(name,"SealevelBarystaticOceanLoad")==0) return SealevelBarystaticOceanLoadEnum; else if (strcmp(name,"SealevelNEsa")==0) return SealevelNEsaEnum; else if (strcmp(name,"SealevelNEsaRate")==0) return SealevelNEsaRateEnum; @@ -789,4 +814,16 @@ else if (strcmp(name,"SealevelchangeGE")==0) return SealevelchangeGEEnum; else if (strcmp(name,"SealevelchangeGN")==0) return SealevelchangeGNEnum; + else if (strcmp(name,"SealevelchangeGsubelOcean")==0) return SealevelchangeGsubelOceanEnum; + else if (strcmp(name,"SealevelchangeGUsubelOcean")==0) return SealevelchangeGUsubelOceanEnum; + else if (strcmp(name,"SealevelchangeGEsubelOcean")==0) return SealevelchangeGEsubelOceanEnum; + else if (strcmp(name,"SealevelchangeGNsubelOcean")==0) return SealevelchangeGNsubelOceanEnum; + else if (strcmp(name,"SealevelchangeGsubelIce")==0) return SealevelchangeGsubelIceEnum; + else if (strcmp(name,"SealevelchangeGUsubelIce")==0) return SealevelchangeGUsubelIceEnum; + else if (strcmp(name,"SealevelchangeGEsubelIce")==0) return SealevelchangeGEsubelIceEnum; + else if (strcmp(name,"SealevelchangeGNsubelIce")==0) return SealevelchangeGNsubelIceEnum; + else if (strcmp(name,"SealevelchangeGsubelHydro")==0) return SealevelchangeGsubelHydroEnum; + else if (strcmp(name,"SealevelchangeGUsubelHydro")==0) return SealevelchangeGUsubelHydroEnum; + else if (strcmp(name,"SealevelchangeGEsubelHydro")==0) return SealevelchangeGEsubelHydroEnum; + else if (strcmp(name,"SealevelchangeGNsubelHydro")==0) return SealevelchangeGNsubelHydroEnum; else if (strcmp(name,"SedimentHead")==0) return SedimentHeadEnum; else if (strcmp(name,"SedimentHeadOld")==0) return SedimentHeadOldEnum; @@ -838,5 +875,8 @@ else if (strcmp(name,"SmbEC")==0) return SmbECEnum; else if (strcmp(name,"SmbECDt")==0) return SmbECDtEnum; - else if (strcmp(name,"SmbECini")==0) return SmbECiniEnum; + else stage=8; + } + if(stage==8){ + if (strcmp(name,"SmbECini")==0) return SmbECiniEnum; else if (strcmp(name,"SmbEla")==0) return SmbElaEnum; else if (strcmp(name,"SmbEvaporation")==0) return SmbEvaporationEnum; @@ -875,8 +915,5 @@ else if (strcmp(name,"SmbRunoffSubstep")==0) return SmbRunoffSubstepEnum; else if (strcmp(name,"SmbRunoffTransient")==0) return SmbRunoffTransientEnum; - else stage=8; - } - if(stage==8){ - if (strcmp(name,"SmbS0gcm")==0) return SmbS0gcmEnum; + else if (strcmp(name,"SmbS0gcm")==0) return SmbS0gcmEnum; else if (strcmp(name,"SmbS0p")==0) return SmbS0pEnum; else if (strcmp(name,"SmbS0t")==0) return SmbS0tEnum; @@ -961,5 +998,8 @@ else if (strcmp(name,"VxMesh")==0) return VxMeshEnum; else if (strcmp(name,"VxObs")==0) return VxObsEnum; - else if (strcmp(name,"VxShear")==0) return VxShearEnum; + else stage=9; + } + if(stage==9){ + if (strcmp(name,"VxShear")==0) return VxShearEnum; else if (strcmp(name,"VxSurface")==0) return VxSurfaceEnum; else if (strcmp(name,"VyAverage")==0) return VyAverageEnum; @@ -998,8 +1038,5 @@ else if (strcmp(name,"Outputdefinition20")==0) return Outputdefinition20Enum; else if (strcmp(name,"Outputdefinition21")==0) return Outputdefinition21Enum; - else stage=9; - } - if(stage==9){ - if (strcmp(name,"Outputdefinition22")==0) return Outputdefinition22Enum; + else if (strcmp(name,"Outputdefinition22")==0) return Outputdefinition22Enum; else if (strcmp(name,"Outputdefinition23")==0) return Outputdefinition23Enum; else if (strcmp(name,"Outputdefinition24")==0) return Outputdefinition24Enum; @@ -1084,5 +1121,8 @@ else if (strcmp(name,"Outputdefinition96")==0) return Outputdefinition96Enum; else if (strcmp(name,"Outputdefinition97")==0) return Outputdefinition97Enum; - else if (strcmp(name,"Outputdefinition98")==0) return Outputdefinition98Enum; + else stage=10; + } + if(stage==10){ + if (strcmp(name,"Outputdefinition98")==0) return Outputdefinition98Enum; else if (strcmp(name,"Outputdefinition99")==0) return Outputdefinition99Enum; else if (strcmp(name,"Outputdefinition9")==0) return Outputdefinition9Enum; @@ -1121,8 +1161,5 @@ else if (strcmp(name,"BuddJacka")==0) return BuddJackaEnum; else if (strcmp(name,"CalvingDev2")==0) return CalvingDev2Enum; - else stage=10; - } - if(stage==10){ - if (strcmp(name,"CalvingHab")==0) return CalvingHabEnum; + else if (strcmp(name,"CalvingHab")==0) return CalvingHabEnum; else if (strcmp(name,"CalvingLevermann")==0) return CalvingLevermannEnum; else if (strcmp(name,"CalvingVonmises")==0) return CalvingVonmisesEnum; @@ -1207,5 +1244,8 @@ else if (strcmp(name,"Gradient1")==0) return Gradient1Enum; else if (strcmp(name,"Gradient2")==0) return Gradient2Enum; - else if (strcmp(name,"Gradient3")==0) return Gradient3Enum; + else stage=11; + } + if(stage==11){ + if (strcmp(name,"Gradient3")==0) return Gradient3Enum; else if (strcmp(name,"Gradient4")==0) return Gradient4Enum; else if (strcmp(name,"GroundedArea")==0) return GroundedAreaEnum; @@ -1244,8 +1284,5 @@ else if (strcmp(name,"IntMatExternalResult")==0) return IntMatExternalResultEnum; else if (strcmp(name,"IntMatParam")==0) return IntMatParamEnum; - else stage=11; - } - if(stage==11){ - if (strcmp(name,"IntParam")==0) return IntParamEnum; + else if (strcmp(name,"IntParam")==0) return IntParamEnum; else if (strcmp(name,"IntVecParam")==0) return IntVecParamEnum; else if (strcmp(name,"Inputs")==0) return InputsEnum; @@ -1330,5 +1367,8 @@ else if (strcmp(name,"OceantransportSolution")==0) return OceantransportSolutionEnum; else if (strcmp(name,"OldGradient")==0) return OldGradientEnum; - else if (strcmp(name,"OneLayerP4z")==0) return OneLayerP4zEnum; + else stage=12; + } + if(stage==12){ + if (strcmp(name,"OneLayerP4z")==0) return OneLayerP4zEnum; else if (strcmp(name,"Open")==0) return OpenEnum; else if (strcmp(name,"Option")==0) return OptionEnum; @@ -1367,8 +1407,5 @@ else if (strcmp(name,"SIAApproximation")==0) return SIAApproximationEnum; else if (strcmp(name,"SMBcomponents")==0) return SMBcomponentsEnum; - else stage=12; - } - if(stage==12){ - if (strcmp(name,"SMBd18opdd")==0) return SMBd18opddEnum; + else if (strcmp(name,"SMBd18opdd")==0) return SMBd18opddEnum; else if (strcmp(name,"SMBforcing")==0) return SMBforcingEnum; else if (strcmp(name,"SMBgcm")==0) return SMBgcmEnum; @@ -1453,5 +1490,8 @@ else if (strcmp(name,"VertexSId")==0) return VertexSIdEnum; else if (strcmp(name,"Vertices")==0) return VerticesEnum; - else if (strcmp(name,"ViscousHeating")==0) return ViscousHeatingEnum; + else stage=13; + } + if(stage==13){ + if (strcmp(name,"ViscousHeating")==0) return ViscousHeatingEnum; else if (strcmp(name,"Water")==0) return WaterEnum; else if (strcmp(name,"XTaylorHood")==0) return XTaylorHoodEnum; @@ -1476,5 +1516,5 @@ else if (strcmp(name,"TriangleInterp")==0) return TriangleInterpEnum; else if (strcmp(name,"MaximumNumberOfDefinitions")==0) return MaximumNumberOfDefinitionsEnum; - else stage=13; + else stage=14; } /*If we reach this point, the string provided has not been found*/ Index: /issm/trunk-jpl/src/m/classes/constants.m =================================================================== --- /issm/trunk-jpl/src/m/classes/constants.m (revision 26221) +++ /issm/trunk-jpl/src/m/classes/constants.m (revision 26222) @@ -10,4 +10,5 @@ yts = 0.; referencetemperature = 0.; + gravitational_constant = 0.; end methods @@ -33,4 +34,7 @@ %the reference temperature for enthalpy model (cf Aschwanden) self.referencetemperature=223.15; + + %gravity: + self.gravitational_constant = 6.67259e-11; end % }}} @@ -41,4 +45,5 @@ md = checkfield(md,'fieldname','constants.yts','>',0,'size',[1 1]); md = checkfield(md,'fieldname','constants.referencetemperature','size',[1 1]); + md = checkfield(md,'fieldname','constants.gravitational_constant','size',[1 1]); end % }}} @@ -50,4 +55,5 @@ fielddisplay(self,'yts','number of seconds in a year [s/yr]'); fielddisplay(self,'referencetemperature','reference temperature used in the enthalpy model [K]'); + fielddisplay(self,'gravitational_constant','Newtonian constant of gravitation [m^3/kg/s^2]'); end % }}} @@ -56,4 +62,5 @@ WriteData(fid,prefix,'object',self,'fieldname','yts','format','Double'); WriteData(fid,prefix,'object',self,'fieldname','referencetemperature','format','Double'); + WriteData(fid,prefix,'object',self,'fieldname','gravitational_constant','format','Double'); end % }}} function savemodeljs(self,fid,modelname) % {{{ @@ -63,4 +70,5 @@ writejsdouble(fid,[modelname '.constants.yts'],self.yts); writejsdouble(fid,[modelname '.constants.referencetemperature'],self.referencetemperature); + writejsdouble(fid,[modelname '.constants.gravitational_constant'],self.gravitational_constant); end % }}} Index: /issm/trunk-jpl/src/m/classes/solidearthsettings.m =================================================================== --- /issm/trunk-jpl/src/m/classes/solidearthsettings.m (revision 26221) +++ /issm/trunk-jpl/src/m/classes/solidearthsettings.m (revision 26222) @@ -19,5 +19,4 @@ degacc = 0; %degree increment for resolution of Green tables horiz = 0; %compute horizontal deformation - glfraction = 1; %barystatic contribution full or fractional (default fractional) 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 @@ -56,7 +55,4 @@ self.runfrequency=1; - %fractional contribution: - self.glfraction=1; - %horizontal displacement? (not by default) self.horiz=0; @@ -80,5 +76,4 @@ md = checkfield(md,'fieldname','solidearth.settings.degacc','size',[1 1],'>=',1e-10); md = checkfield(md,'fieldname','solidearth.settings.horiz','NaN',1,'Inf',1,'values',[0 1]); - md = checkfield(md,'fieldname','solidearth.settings.glfraction','values',[0 1]); md = checkfield(md,'fieldname','solidearth.settings.grdmodel','values',[1 2]); md = checkfield(md,'fieldname','solidearth.settings.cross_section_shape','numel',[1],'values',[1,2]); @@ -123,5 +118,4 @@ fielddisplay(self,'rotation','earth rotational potential perturbation'); fielddisplay(self,'degacc','accuracy (default .01 deg) for numerical discretization of the Green''s functions'); - fielddisplay(self,'glfraction','contribute fractionally (default, 1) to barystatic sea level'); 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'); @@ -141,5 +135,4 @@ WriteData(fid,prefix,'object',self,'fieldname','isgrd','name','md.solidearth.settings.isgrd','format','Integer'); WriteData(fid,prefix,'object',self,'fieldname','compute_bp_grd','name','md.solidearth.settings.compute_bp_grd','format','Integer'); - WriteData(fid,prefix,'object',self,'fieldname','glfraction','name','md.solidearth.settings.glfraction','format','Integer'); 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'); @@ -156,5 +149,4 @@ writejsdouble(fid,[modelname '.solidearth.settings.run_frequency'],self.run_frequency); writejsdouble(fid,[modelname '.solidearth.settings.degacc'],self.degacc); - writejsdouble(fid,[modelname '.solidearth.settings.glfraction'],self.glfraction); writejsdouble(fid,[modelname '.solidearth.settings.cross_section_shape'],self.cross_section_shape); end % }}} Index: /issm/trunk-jpl/test/NightlyRun/test2002.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test2002.m (revision 26221) +++ /issm/trunk-jpl/test/NightlyRun/test2002.m (revision 26222) @@ -20,11 +20,16 @@ md.smb.mass_balance=zeros(md.mesh.numberofvertices,1); %antarctica -late=sum(md.mesh.lat(md.mesh.elements),2)/3; -longe=sum(md.mesh.long(md.mesh.elements),2)/3; +xe=md.mesh.x(md.mesh.elements)*[1;1;1]/3; +ye=md.mesh.y(md.mesh.elements)*[1;1;1]/3; +ze=md.mesh.z(md.mesh.elements)*[1;1;1]/3; +re=sqrt(xe.^2+ye.^2+ze.^2); + +late=asind(ze./re); +longe=atan2d(ye,xe); 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); +pos=find(late>60 & late<90 & longe>-75 & longe<-15); md.masstransport.spcthickness(md.mesh.elements(pos,:))= md.masstransport.spcthickness(md.mesh.elements(pos,:))-100; posgre=pos; @@ -51,5 +56,4 @@ md.timestepping.time_step=1; md.timestepping.final_time=1; - md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1); @@ -93,4 +97,5 @@ md.solidearth.settings.elastic=0; md.solidearth.settings.rotation=0; + md=solve(md,'Transient'); Seustatic=md.results.TransientSolution.Sealevel; Index: /issm/trunk-jpl/test/NightlyRun/test2011.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test2011.m (revision 26222) +++ /issm/trunk-jpl/test/NightlyRun/test2011.m (revision 26222) @@ -0,0 +1,115 @@ +%Test Name: EarthSlc + +%mesh earth: +md=model; +md.mesh=gmshplanet('radius',6.371012*10^3,'resolution',500.); %700 km resolution mesh + +%Geometry for the bed, arbitrary thickness of 1000: +md.geometry.bed=-ones(md.mesh.numberofvertices,1); +md.geometry.base=md.geometry.bed; +md.geometry.thickness=1000*ones(md.mesh.numberofvertices,1); +md.geometry.surface=md.geometry.bed+md.geometry.thickness; + + +%parameterize solidearth solution: +%solidearth loading: {{{ +md.masstransport.spcthickness=[md.geometry.thickness;0]; +md.dsl.global_average_thermosteric_sea_level=[0;0]; +md.dsl.sea_surface_height_above_geoid=zeros(md.mesh.numberofvertices+1,1); +md.dsl.sea_water_pressure_at_sea_floor=zeros(md.mesh.numberofvertices+1,1); +md.smb.mass_balance=zeros(md.mesh.numberofvertices,1); +%antarctica +xe=md.mesh.x(md.mesh.elements)*[1;1;1]/3; +ye=md.mesh.y(md.mesh.elements)*[1;1;1]/3; +ze=md.mesh.z(md.mesh.elements)*[1;1;1]/3; +re=sqrt(xe.^2+ye.^2+ze.^2); + +late=asind(ze./re); +longe=atan2d(ye,xe); +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>60 & late<90 & longe>-75 & longe<-15); +md.masstransport.spcthickness(md.mesh.elements(pos,:))= md.masstransport.spcthickness(md.mesh.elements(pos,:))-100; +posgre=pos; + +%elastic loading from love numbers: +md.solidearth.lovenumbers=lovenumbers('maxdeg',100); + +%}}} +%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); +icemask(md.mesh.elements(posant))=-1; +icemask(md.mesh.elements(posgre))=-1; + +md.mask.ice_levelset=icemask; +md.mask.ocean_levelset=oceanmask; +% }}} + +%time stepping: +md.timestepping.start_time=0; +md.timestepping.time_step=1; +md.timestepping.final_time=1; + +md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1); +md.basalforcings.floatingice_melting_rate=zeros(md.mesh.numberofvertices,1); +md.initialization.vx=zeros(md.mesh.numberofvertices,1); +md.initialization.vy=zeros(md.mesh.numberofvertices,1); +md.initialization.sealevel=zeros(md.mesh.numberofvertices,1); +md.initialization.bottompressure=zeros(md.mesh.numberofvertices,1); +md.initialization.dsl=zeros(md.mesh.numberofvertices,1); +md.initialization.str=0; + +%Materials: +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; +md.solidearth.settings.computesealevelchange=1; +md.solidearth.settings.isgrd=1; +md.solidearth.settings.ocean_area_scaling=0; +md.solidearth.settings.grdmodel=1; + +%Physics: +md.transient.issmb=0; +md.transient.isstressbalance=0; +md.transient.isthermal=0; +md.transient.ismasstransport=1; +md.transient.isslc=1; +md.solidearth.requested_outputs={'Sealevel','DeltaIceThickness','SealevelBarystaticIceMask','SealevelBarystaticHydroMask','SealevelBarystaticBpMask','Bed','SealevelBarystaticIceWeights','SealevelBarystaticOceanWeights','SealevelBarystaticIceArea','SealevelBarystaticOceanArea','SealevelBarystaticOceanMask','SealevelGRD','SealevelBarystaticOceanLongbar'}; +md.settings.results_on_nodes={'SealevelBarystaticIceWeights','SealevelBarystaticOceanWeights'}; + +% max number of iteration reverted back to 10 (i.e., the original default value) +md.solidearth.settings.maxiter=10; + +%eustatic + rigid + elastic + rotation run: +md.solidearth.settings.rigid=1; +md.solidearth.settings.elastic=1; +md.solidearth.settings.rotation=1; +md=solve(md,'tr'); + +%recover barystatic: +bslc=md.results.TransientSolution(1).BslcIce; + +%alternative way of computing barystatic: +area=md.results.TransientSolution(1).SealevelBarystaticIceArea; +weights=md.results.TransientSolution(1).SealevelBarystaticIceWeights; +dH=md.results.TransientSolution(1).DeltaIceThickness(md.mesh.elements); +dHavg=sum(dH.*weights); +oceanarea=sum(md.results.TransientSolution(1).SealevelBarystaticOceanArea); +bslc2=-sum(dHavg.*area)*md.materials.rho_ice/md.materials.rho_water/oceanarea; + +%Fields and tolerances to track changes +field_names={'BarystaticIce','BarystaticIce2','BarystaticIceDiff'}; +field_tolerances={1e-13,1e-13,1e-13}; +field_values={bslc,bslc2,bslc2-bslc2};