Index: /issm/trunk-jpl/src/c/classes/BarystaticContributions.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/BarystaticContributions.cpp (revision 26148) +++ /issm/trunk-jpl/src/c/classes/BarystaticContributions.cpp (revision 26149) @@ -148,8 +148,8 @@ ice->Sum(&sumice); hydro->Sum(&sumhydro); ocean->Sum(&sumocean); - results->AddResult(new GenericExternalResult(results->Size()+1,BslcEnum,-this->Total()/oceanarea/rho_water,step,time)); - results->AddResult(new GenericExternalResult(results->Size()+1,BslcIceEnum,-sumice/oceanarea/rho_water,step,time)); - results->AddResult(new GenericExternalResult(results->Size()+1,BslcHydroEnum,-sumice/oceanarea/rho_water,step,time)); - results->AddResult(new GenericExternalResult(results->Size()+1,BslcOceanEnum,-sumocean/oceanarea/rho_water,step,time)); + results->AddResult(new GenericExternalResult(results->Size()+1,BslcEnum,this->Total()/oceanarea/rho_water,step,time)); + results->AddResult(new GenericExternalResult(results->Size()+1,BslcIceEnum,sumice/oceanarea/rho_water,step,time)); + results->AddResult(new GenericExternalResult(results->Size()+1,BslcHydroEnum,sumice/oceanarea/rho_water,step,time)); + results->AddResult(new GenericExternalResult(results->Size()+1,BslcOceanEnum,sumocean/oceanarea/rho_water,step,time)); cumice->Sum(&sumice); cumhydro->Sum(&sumhydro); cumocean->Sum(&sumocean); Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 26148) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 26149) @@ -239,5 +239,5 @@ virtual void ElementResponse(IssmDouble* presponse,int response_enum)=0; virtual void ElementSizes(IssmDouble* phx,IssmDouble* phy,IssmDouble* phz)=0; - virtual void ElementCoordinates(Vector* vxe,Vector* vye,Vector* vze,bool spherical=false)=0; + virtual void ElementCoordinates(Vector* vxe,Vector* vye,Vector* vze, Vector* vareae, bool spherical=false)=0; virtual int FiniteElement(void)=0; virtual IssmDouble FloatingArea(bool scaled)=0; @@ -392,5 +392,5 @@ virtual void GiaDeflection(Vector* wg,Vector* dwgdt,Matlitho* litho, IssmDouble* x,IssmDouble* y)=0; - virtual void SealevelchangeGeometryOptim(IssmDouble* lat,IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze)=0; + virtual void SealevelchangeGeometryOptim(IssmDouble* lat,IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae)=0; virtual void SealevelchangeBarystaticLoads(Vector* loads, BarystaticContributions* barycontrib, SealevelMasks* masks)=0; virtual void SealevelchangeInitialConvolution(Vector* sealevelloads,Vector* oceanareas,IssmDouble* loads, IssmDouble* rotationvector, SealevelMasks* masks)=0; Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 26148) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 26149) @@ -70,5 +70,5 @@ void CreateInputTimeAverage(int transientinput_enum,int averagedinput_enum,IssmDouble start_time,IssmDouble end_time,int averaging_method); void ElementResponse(IssmDouble* presponse,int response_enum); - void ElementCoordinates(Vector* vxe,Vector* vye,Vector* vze,bool spherical=false){_error_("not implemented yet");}; + void ElementCoordinates(Vector* vxe,Vector* vye,Vector* vze,Vector* vareae,bool spherical=false){_error_("not implemented yet");}; void ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz); int FiniteElement(void); @@ -229,5 +229,5 @@ void GiaDeflection(Vector* wg,Vector* dwgdt,Matlitho* litho, IssmDouble* x,IssmDouble* y){_error_("not implemented yet");}; - void SealevelchangeGeometryOptim(IssmDouble* lat,IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze){_error_("not implemented yet");}; + void SealevelchangeGeometryOptim(IssmDouble* lat,IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae){_error_("not implemented yet");}; void SealevelchangeBarystaticLoads(Vector* loads, BarystaticContributions* barycontrib, SealevelMasks* masks){_error_("not implemented yet");}; void SealevelchangeInitialConvolution(Vector* sealevelloads,Vector* oceanareas,IssmDouble* loads, IssmDouble* rotationvector, SealevelMasks* masks){_error_("not implemented yet");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 26148) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 26149) @@ -54,5 +54,5 @@ void ControlToVectors(Vector* vector_control, Vector* vector_gradient,int control_enum,int control_interp){_error_("not implemented yet");}; void ElementResponse(IssmDouble* presponse,int response_enum){_error_("not implemented yet");}; - void ElementCoordinates(Vector* vxe,Vector* vye,Vector* vze,bool spherical=false){_error_("not implemented yet");}; + void ElementCoordinates(Vector* vxe,Vector* vye,Vector* vze,Vector* vareae,bool spherical=false){_error_("not implemented yet");}; void ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz){_error_("not implemented yet");}; int FiniteElement(void); @@ -184,5 +184,5 @@ void GiaDeflection(Vector* wg,Vector* dwgdt,Matlitho* litho, IssmDouble* x,IssmDouble* y){_error_("not implemented yet");}; - void SealevelchangeGeometryOptim(IssmDouble* lat,IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze){_error_("not implemented yet");}; + void SealevelchangeGeometryOptim(IssmDouble* lat,IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae){_error_("not implemented yet");}; void SealevelchangeBarystaticLoads(Vector* loads, BarystaticContributions* barycontrib, SealevelMasks* masks){_error_("not implemented yet");}; void SealevelchangeInitialConvolution(Vector* sealevelloads,Vector* oceanareas,IssmDouble* loads, IssmDouble* rotationvector, SealevelMasks* masks){_error_("not implemented yet");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Tetra.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 26148) +++ /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 26149) @@ -53,5 +53,5 @@ IssmDouble DragCoefficientAbsGradient(void){_error_("not implemented yet");}; void ElementResponse(IssmDouble* presponse,int response_enum){_error_("not implemented yet");}; - void ElementCoordinates(Vector* vxe,Vector* vye,Vector* vze,bool spherical=false){_error_("not implemented yet");}; + void ElementCoordinates(Vector* vxe,Vector* vye,Vector* vze,Vector* vareae,bool spherical=false){_error_("not implemented yet");}; void ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz); void FaceOnBaseIndices(int* pindex1,int* pindex2,int* pindex3); @@ -190,5 +190,5 @@ IssmDouble OceanAverage(IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet!");}; void GiaDeflection(Vector* wg,Vector* dwgdt, Matlitho* litho, IssmDouble* x,IssmDouble* y){_error_("not implemented yet");}; - void SealevelchangeGeometryOptim(IssmDouble* lat,IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze){_error_("not implemented yet");}; + void SealevelchangeGeometryOptim(IssmDouble* lat,IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae){_error_("not implemented yet");}; void SealevelchangeBarystaticLoads(Vector* loads, BarystaticContributions* barycontrib, SealevelMasks* masks){_error_("not implemented yet");}; void SealevelchangeInitialConvolution(Vector* sealevelloads,Vector* oceanareas,IssmDouble* loads, IssmDouble* rotationvector, SealevelMasks* masks){_error_("not implemented yet");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 26148) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 26149) @@ -1286,5 +1286,5 @@ } /*}}}*/ -void Tria::ElementCoordinates(Vector* vxe,Vector* vye,Vector* vze, bool spherical){/*{{{*/ +void Tria::ElementCoordinates(Vector* vxe,Vector* vye,Vector* vze, Vector* vareae, bool spherical){/*{{{*/ /*Look for x,y,z coordinates:*/ @@ -1298,7 +1298,13 @@ if(!spherical){ + IssmDouble area; vxe->SetValue(this->sid,xe,INS_VAL); vye->SetValue(this->sid,ye,INS_VAL); vze->SetValue(this->sid,ze,INS_VAL); + area=this->GetAreaSpherical(); + vareae->SetValue(this->sid,area,INS_VAL); + + /*in addition, put in in the inputs:*/ + this->inputs->SetDoubleInput(AreaEnum,this->lid,area); } else _error_("spherical coordinates not supported yet!"); @@ -6239,16 +6245,13 @@ //new code -void Tria::SealevelchangeGeometryOptim(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze){ /*{{{*/ +void Tria::SealevelchangeGeometryOptim(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae){ /*{{{*/ /*diverse:*/ int nel; - bool spherical=true; - IssmDouble llr_list[NUMVERTICES][3]; - IssmDouble xyz_list[NUMVERTICES][3]; IssmDouble area,planetarea,planetradius; + IssmDouble constant,ratioe; IssmDouble I; //change in ice thickness or water level(Farrel and Clarke, Equ. 4) IssmDouble rho_earth; IssmDouble lati,longi; - IssmDouble constant; IssmDouble x,y,z,dx,dy,dz,N_azim,E_azim; int sidlist[NUMVERTICES]; @@ -6321,8 +6324,5 @@ } - /*compute area and add to inputs:*/ - area=GetAreaSpherical(); - this->inputs->SetDoubleInput(AreaEnum,this->lid,area); - this->AddInput(SealevelAreaEnum,&area,P0Enum); + /*early return:*/ if(!computerigid)return; @@ -6353,7 +6353,4 @@ } - /*retrieve coordinates: */ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - /*compute centroids of all elements:*/ for(int i=0;iinputs->SetDoubleInput(AreaEnum,this->lid,area); - this->AddInput(SealevelAreaEnum,&area,P0Enum); /*Free allocations:*/ @@ -6472,10 +6469,16 @@ 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:*/ 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 I,W,BP; //change in ice thickness, water column or bottom pressure (Farrel and Clarke, Equ. 4) - bool notfullygrounded=false; + 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; @@ -6529,4 +6532,8 @@ if(!isice && !ishydro){ if(!masks->isoceanin[this->lid]){ + constant=0; + #ifdef _ISSM_DEBUG_ + this->AddInput(SealevelBarystaticMaskEnum,&constant,P0Enum); + #endif return; } @@ -6603,5 +6610,5 @@ /*Compute barystatic contribution in kg:*/ - bslcice = rho_ice*area*I; + bslcice = -rho_ice*phi_ice*area*I; _assert_(!xIsNan(bslcice)); @@ -6612,4 +6619,9 @@ /*Deal with water loads if we are on ground:*/ if(!masks->isfullyfloating[this->lid]){ + + constant+=10; //1 for water. + #ifdef _ISSM_DEBUG_ + this->AddInput(SealevelBarystaticMaskEnum,&constant,P0Enum); + #endif /*Retrieve water height at vertices: */ @@ -6619,5 +6631,5 @@ /*Compute barystatic component in kg:*/ - bslchydro = rho_freshwater*area*phi_water*W; + bslchydro = -rho_freshwater*area*phi_water*W; _assert_(!xIsNan(bslchydro)); @@ -6628,4 +6640,9 @@ /*Deal with ocean bottom pressures:*/ if(masks->isoceanin[this->lid]){ + + constant+=1; //1 for bottom pressure. + #ifdef _ISSM_DEBUG_ + this->AddInput(SealevelBarystaticMaskEnum,&constant,P0Enum); + #endif /*Retrieve bottom pressure change and average over the element: */ @@ -6635,8 +6652,8 @@ /*Compute barystatic component in kg:*/ - bslcbp = rho_water*area*BP; + bslcbp = -rho_water*area*phi_bp*BP; /*convert from m to kg/m^2:*/ - BP=BP*rho_water; + BP=BP*rho_water*phi_bp; } Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 26148) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 26149) @@ -70,5 +70,5 @@ void ControlToVectors(Vector* vector_control, Vector* vector_gradient,int control_enum,int control_interp); void CreateDistanceInputFromSegmentlist(IssmDouble* distances,int distanceenum); - void ElementCoordinates(Vector* vxe,Vector* vye,Vector* vze,bool spherical=false); + void ElementCoordinates(Vector* vxe,Vector* vye,Vector* vze, Vector* vareae, bool spherical=false); int EdgeOnBaseIndex(); void EdgeOnBaseIndices(int* pindex1,int* pindex); @@ -174,5 +174,5 @@ void SetSealevelMasks(SealevelMasks* masks); - void SealevelchangeGeometryOptim(IssmDouble* lat,IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze); + void SealevelchangeGeometryOptim(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 SealevelchangeInitialConvolution(Vector* sealevelloads,Vector* oceanareas,IssmDouble* loads, IssmDouble* rotationvector, SealevelMasks* masks); Index: /issm/trunk-jpl/src/c/cores/sealevelchange_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/sealevelchange_core.cpp (revision 26148) +++ /issm/trunk-jpl/src/c/cores/sealevelchange_core.cpp (revision 26149) @@ -608,4 +608,5 @@ IssmDouble *yye = NULL; IssmDouble *zze = NULL; + IssmDouble* areae = NULL; int horiz; @@ -634,12 +635,13 @@ /*first, recover lat,long and radius vectors from vertices: */ VertexCoordinatesx(&latitude,&longitude,&radius,femmodel->vertices,spherical); - ElementCoordinatesx(&xxe,&yye,&zze,femmodel->elements); if(horiz) VertexCoordinatesx(&xx,&yy,&zz,femmodel->vertices); - + + /*first, recover x,y,z and areas from elements: */ + ElementCoordinatesx(&xxe,&yye,&zze,&areae,femmodel->elements); /*Run sealevel geometry routine in elements:*/ for(Object* & object : femmodel->elements->objects){ Element* element=xDynamicCast(object); - if(optim) element->SealevelchangeGeometryOptim(latitude,longitude,radius,xx,yy,zz,xxe,yye,zze); + if(optim) element->SealevelchangeGeometryOptim(latitude,longitude,radius,xx,yy,zz,xxe,yye,zze,areae); else element->SealevelchangeGeometry(latitude,longitude,radius,xx,yy,zz,xxe,yye,zze); } @@ -654,4 +656,5 @@ xDelete(yye); xDelete(zze); + xDelete(areae); xDelete(latitude); xDelete(longitude); @@ -1251,10 +1254,13 @@ } /*}}}*/ IssmDouble SealevelloadsOceanAverage(Vector* sealevelloads,Vector* oceanareas, IssmDouble oceanarea){ /*{{{*/ + IssmDouble sealevelloadsaverage; + Vector* sealevelloadsvolume=sealevelloads->Duplicate(); sealevelloadsvolume->PointwiseMult(sealevelloads,oceanareas); sealevelloadsvolume->Sum(&sealevelloadsaverage); delete sealevelloadsvolume; - return sealevelloadsaverage/oceanarea; + + return sealevelloadsaverage/oceanarea;- } /*}}}*/ void RotationAxisMotion(IssmDouble* m, FemModel* femmodel,IssmDouble* loads, IssmDouble* sealevelloads){ /*{{{*/ Index: /issm/trunk-jpl/src/c/modules/ElementCoordinatesx/ElementCoordinatesx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ElementCoordinatesx/ElementCoordinatesx.cpp (revision 26148) +++ /issm/trunk-jpl/src/c/modules/ElementCoordinatesx/ElementCoordinatesx.cpp (revision 26149) @@ -9,5 +9,5 @@ #include "../../toolkits/toolkits.h" -void ElementCoordinatesx( IssmDouble** pxe, IssmDouble** pye, IssmDouble** pze,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: */ @@ -17,9 +17,10 @@ Vector* vye=new Vector(numberofelements); Vector* vze=new Vector(numberofelements); + Vector* vareae=new Vector(numberofelements); /*march through our elements: */ for(Object* & object : elements->objects){ Element* element=(Element*)object; - element->ElementCoordinates(vxe,vye,vze,spherical); + element->ElementCoordinates(vxe,vye,vze,vareae,spherical); } @@ -28,4 +29,5 @@ vye->Assemble(); vze->Assemble(); + vareae->Assemble(); /*serialize: */ @@ -33,4 +35,5 @@ IssmDouble* ye=vye->ToMPISerial(); IssmDouble* ze=vze->ToMPISerial(); + IssmDouble* areae=vareae->ToMPISerial(); /*Free ressources: */ @@ -38,4 +41,5 @@ delete vye; delete vze; + delete vareae; /*output: */ @@ -46,3 +50,6 @@ if(pze) *pze=ze; else xDelete(ze); + 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 26148) +++ /issm/trunk-jpl/src/c/modules/ElementCoordinatesx/ElementCoordinatesx.h (revision 26149) @@ -8,5 +8,5 @@ /* local prototypes: */ -void ElementCoordinatesx( IssmDouble** pxe, IssmDouble** pye, IssmDouble** pze,Elements* elements,bool spherical=false); +void ElementCoordinatesx( IssmDouble** pxe, IssmDouble** pye, IssmDouble** pze,IssmDouble** pareae, Elements* elements,bool spherical=false); #endif /* _ELEMENT_COORDINATESX_H */