Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 23990) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 23991) @@ -320,4 +320,6 @@ virtual int TensorInterpolation()=0; virtual IssmDouble TimeAdapt()=0; + virtual IssmDouble TotalCalvingFlux(bool scaled){_error_("not implemented");}; + virtual IssmDouble TotalCalvingFluxLevelset(bool scaled){_error_("not implemented");}; virtual IssmDouble TotalFloatingBmb(bool scaled)=0; virtual IssmDouble TotalGroundedBmb(bool scaled)=0; Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 23990) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 23991) @@ -4060,4 +4060,189 @@ } /*}}}*/ +IssmDouble Tria::TotalCalvingFlux(bool scaled){/*{{{*/ + + /*Make sure there is an ice front here*/ + if(!IsIceInElement() || !IsIcefront()) return 0; + + /*Scaled not implemented yet...*/ + _assert_(!scaled); + + /*Get domain type*/ + int domaintype; + parameters->FindParam(&domaintype,DomainTypeEnum); + + /*Get ice front coordinates*/ + IssmDouble *xyz_list = NULL; + IssmDouble* xyz_front = NULL; + this->GetVerticesCoordinates(&xyz_list); + this->GetIcefrontCoordinates(&xyz_front,xyz_list,MaskIceLevelsetEnum); + + /*Get normal vector*/ + IssmDouble normal[3]; + this->NormalSection(&normal[0],xyz_front); + //normal[0] = -normal[0]; + //normal[1] = -normal[1]; + + /*Get inputs*/ + IssmDouble flux = 0.; + IssmDouble calvingratex,calvingratey,thickness,Jdet; + IssmDouble rho_ice=FindParam(MaterialsRhoIceEnum); + Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input); + Input* calvingratex_input=NULL; + Input* calvingratey_input=NULL; + if(domaintype==Domain2DhorizontalEnum){ + calvingratex_input=inputs->GetInput(CalvingratexEnum); _assert_(calvingratex_input); + calvingratey_input=inputs->GetInput(CalvingrateyEnum); _assert_(calvingratey_input); + } + else{ + calvingratex_input=inputs->GetInput(CalvingratexAverageEnum); _assert_(calvingratex_input); + calvingratey_input=inputs->GetInput(CalvingrateyAverageEnum); _assert_(calvingratey_input); + } + + /*Start looping on Gaussian points*/ + Gauss* gauss=this->NewGauss(xyz_list,xyz_front,3); + for(int ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + thickness_input->GetInputValue(&thickness,gauss); + calvingratex_input->GetInputValue(&calvingratex,gauss); + calvingratey_input->GetInputValue(&calvingratey,gauss); + this->JacobianDeterminantSurface(&Jdet,xyz_front,gauss); + + flux += rho_ice*Jdet*gauss->weight*thickness*(calvingratex*normal[0] + calvingratey*normal[1]); + } + + return flux; +} +/*}}}*/ +IssmDouble Tria::TotalCalvingFluxLevelset(bool scaled){/*{{{*/ + + /*Make sure there is an ice front here*/ + if(!IsIceInElement() || !IsZeroLevelset(MaskIceLevelsetEnum)) return 0; + + /*Scaled not implemented yet...*/ + _assert_(!scaled); + + int domaintype,index1,index2; + const IssmPDouble epsilon = 1.e-15; + IssmDouble s1,s2; + IssmDouble gl[NUMVERTICES]; + IssmDouble xyz_front[2][3]; + + + IssmDouble *xyz_list = NULL; + this->GetVerticesCoordinates(&xyz_list); + + /*Recover parameters and values*/ + parameters->FindParam(&domaintype,DomainTypeEnum); + GetInputListOnVertices(&gl[0],MaskIceLevelsetEnum); + + /*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){ + _error_("not implemented"); + } + else if(domaintype==Domain2DhorizontalEnum || domaintype==Domain3DEnum || domaintype==Domain3DsurfaceEnum){ + int pt1 = 0; + int pt2 = 1; + if(gl[0]*gl[1]>0){ //Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2 + + /*Portion of the segments*/ + s1=gl[2]/(gl[2]-gl[1]); + s2=gl[2]/(gl[2]-gl[0]); + if(gl[2]<0.){ + pt1 = 1; pt2 = 0; + } + xyz_front[pt2][0]=xyz_list[3*2+0]+s1*(xyz_list[3*1+0]-xyz_list[3*2+0]); + xyz_front[pt2][1]=xyz_list[3*2+1]+s1*(xyz_list[3*1+1]-xyz_list[3*2+1]); + xyz_front[pt2][2]=xyz_list[3*2+2]+s1*(xyz_list[3*1+2]-xyz_list[3*2+2]); + xyz_front[pt1][0]=xyz_list[3*2+0]+s2*(xyz_list[3*0+0]-xyz_list[3*2+0]); + xyz_front[pt1][1]=xyz_list[3*2+1]+s2*(xyz_list[3*0+1]-xyz_list[3*2+1]); + xyz_front[pt1][2]=xyz_list[3*2+2]+s2*(xyz_list[3*0+2]-xyz_list[3*2+2]); + } + 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 + + /*Portion of the segments*/ + s1=gl[0]/(gl[0]-gl[1]); + s2=gl[0]/(gl[0]-gl[2]); + if(gl[0]<0.){ + pt1 = 1; pt2 = 0; + } + + xyz_front[pt1][0]=xyz_list[3*0+0]+s1*(xyz_list[3*1+0]-xyz_list[3*0+0]); + xyz_front[pt1][1]=xyz_list[3*0+1]+s1*(xyz_list[3*1+1]-xyz_list[3*0+1]); + xyz_front[pt1][2]=xyz_list[3*0+2]+s1*(xyz_list[3*1+2]-xyz_list[3*0+2]); + xyz_front[pt2][0]=xyz_list[3*0+0]+s2*(xyz_list[3*2+0]-xyz_list[3*0+0]); + xyz_front[pt2][1]=xyz_list[3*0+1]+s2*(xyz_list[3*2+1]-xyz_list[3*0+1]); + xyz_front[pt2][2]=xyz_list[3*0+2]+s2*(xyz_list[3*2+2]-xyz_list[3*0+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 + + /*Portion of the segments*/ + s1=gl[1]/(gl[1]-gl[0]); + s2=gl[1]/(gl[1]-gl[2]); + if(gl[1]<0.){ + pt1 = 1; pt2 = 0; + } + + xyz_front[pt2][0]=xyz_list[3*1+0]+s1*(xyz_list[3*0+0]-xyz_list[3*1+0]); + xyz_front[pt2][1]=xyz_list[3*1+1]+s1*(xyz_list[3*0+1]-xyz_list[3*1+1]); + xyz_front[pt2][2]=xyz_list[3*1+2]+s1*(xyz_list[3*0+2]-xyz_list[3*1+2]); + xyz_front[pt1][0]=xyz_list[3*1+0]+s2*(xyz_list[3*2+0]-xyz_list[3*1+0]); + xyz_front[pt1][1]=xyz_list[3*1+1]+s2*(xyz_list[3*2+1]-xyz_list[3*1+1]); + xyz_front[pt1][2]=xyz_list[3*1+2]+s2*(xyz_list[3*2+2]-xyz_list[3*1+2]); + } + else{ + _error_("case not possible"); + } + + } + else _error_("mesh type "<=0 && s1<=1.); + _assert_(s2>=0 && s2<=1.); + + /*Get normal vector*/ + IssmDouble normal[3]; + this->NormalSection(&normal[0],&xyz_front[0][0]); + normal[0] = -normal[0]; + normal[1] = -normal[1]; + + /*Get inputs*/ + IssmDouble flux = 0.; + IssmDouble calvingratex,calvingratey,thickness,Jdet; + IssmDouble rho_ice=FindParam(MaterialsRhoIceEnum); + Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input); + Input* calvingratex_input=NULL; + Input* calvingratey_input=NULL; + if(domaintype==Domain2DhorizontalEnum){ + calvingratex_input=inputs->GetInput(CalvingratexEnum); _assert_(calvingratex_input); + calvingratey_input=inputs->GetInput(CalvingrateyEnum); _assert_(calvingratey_input); + } + else{ + calvingratex_input=inputs->GetInput(CalvingratexAverageEnum); _assert_(calvingratex_input); + calvingratey_input=inputs->GetInput(CalvingrateyAverageEnum); _assert_(calvingratey_input); + } + + /*Start looping on Gaussian points*/ + Gauss* gauss=this->NewGauss(xyz_list,&xyz_front[0][0],3); + for(int ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + thickness_input->GetInputValue(&thickness,gauss); + calvingratex_input->GetInputValue(&calvingratex,gauss); + calvingratey_input->GetInputValue(&calvingratey,gauss); + this->JacobianDeterminantSurface(&Jdet,&xyz_front[0][0],gauss); + + flux += rho_ice*Jdet*gauss->weight*thickness*(calvingratex*normal[0] + calvingratey*normal[1]); + } + + return flux; +} +/*}}}*/ IssmDouble Tria::TotalFloatingBmb(bool scaled){/*{{{*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 23990) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 23991) @@ -137,4 +137,6 @@ int TensorInterpolation(); IssmDouble TimeAdapt(); + IssmDouble TotalCalvingFlux(bool scaled); + IssmDouble TotalCalvingFluxLevelset(bool scaled); IssmDouble TotalFloatingBmb(bool scaled); IssmDouble TotalGroundedBmb(bool scaled); Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 23990) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 23991) @@ -2202,4 +2202,6 @@ case MaxAbsVzEnum: this->MaxAbsVzx(&double_result); break; case MassFluxEnum: this->MassFluxx(&double_result); break; + case TotalCalvingFluxEnum: this->TotalCalvingFluxx(&double_result,false); break; + case TotalCalvingFluxLevelsetEnum: this->TotalCalvingFluxLevelsetx(&double_result,false); break; case TotalFloatingBmbEnum: this->TotalFloatingBmbx(&double_result,false); break; case TotalFloatingBmbScaledEnum: this->TotalFloatingBmbx(&double_result,true); break; @@ -2446,4 +2448,6 @@ case DragCoefficientAbsGradientEnum: DragCoefficientAbsGradientx(responses, elements,nodes, vertices, loads, materials, parameters); break; case BalancethicknessMisfitEnum: BalancethicknessMisfitx(responses); break; + case TotalCalvingFluxEnum: this->TotalCalvingFluxx(responses, false); break; + case TotalCalvingFluxLevelsetEnum: this->TotalCalvingFluxLevelsetx(responses, false); break; case TotalFloatingBmbEnum: this->TotalFloatingBmbx(responses, false); break; case TotalFloatingBmbScaledEnum: this->TotalFloatingBmbx(responses, true); break; @@ -2710,4 +2714,36 @@ } /*}}}*/ +void FemModel::TotalCalvingFluxx(IssmDouble* pM, bool scaled){/*{{{*/ + + IssmDouble local_calving_flux = 0; + IssmDouble total_calving_flux; + + for(int i=0;ielements->Size();i++){ + Element* element=xDynamicCast(this->elements->GetObjectByOffset(i)); + local_calving_flux+=element->TotalCalvingFlux(scaled); + } + ISSM_MPI_Reduce(&local_calving_flux,&total_calving_flux,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); + ISSM_MPI_Bcast(&total_calving_flux,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); + + /*Assign output pointers: */ + *pM=total_calving_flux; + +}/*}}}*/ +void FemModel::TotalCalvingFluxLevelsetx(IssmDouble* pM, bool scaled){/*{{{*/ + + IssmDouble local_calving_flux = 0; + IssmDouble total_calving_flux; + + for(int i=0;ielements->Size();i++){ + Element* element=xDynamicCast(this->elements->GetObjectByOffset(i)); + local_calving_flux+=element->TotalCalvingFluxLevelset(scaled); + } + ISSM_MPI_Reduce(&local_calving_flux,&total_calving_flux,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); + ISSM_MPI_Bcast(&total_calving_flux,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); + + /*Assign output pointers: */ + *pM=total_calving_flux; + +}/*}}}*/ void FemModel::TotalFloatingBmbx(IssmDouble* pFbmb, bool scaled){/*{{{*/ Index: /issm/trunk-jpl/src/c/classes/FemModel.h =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.h (revision 23990) +++ /issm/trunk-jpl/src/c/classes/FemModel.h (revision 23991) @@ -127,4 +127,6 @@ void StressIntensityFactorx(); void RignotMeltParameterizationx(); + void TotalCalvingFluxx(IssmDouble* pFbmb, bool scaled); + void TotalCalvingFluxLevelsetx(IssmDouble* pGbmb, bool scaled); void TotalFloatingBmbx(IssmDouble* pFbmb, bool scaled); void TotalGroundedBmbx(IssmDouble* pGbmb, bool scaled); Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 23990) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 23991) @@ -145,5 +145,4 @@ FrictionVoidRatioEnum, FrontalForcingsBasinIcefrontAreaEnum, - FrontalForcingsBasinIdEnum, FrontalForcingsNumberofBasinsEnum, FrontalForcingsParamEnum, @@ -532,4 +531,5 @@ FrictionWaterLayerEnum, FrictionfEnum, + FrontalForcingsBasinIdEnum, FrontalForcingsSubglacialDischargeEnum, FrontalForcingsThermalForcingEnum, @@ -1208,4 +1208,6 @@ ThermalSolutionEnum, ThicknessErrorEstimatorEnum, + TotalCalvingFluxEnum, + TotalCalvingFluxLevelsetEnum, TotalFloatingBmbEnum, TotalFloatingBmbScaledEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 23990) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 23991) @@ -153,5 +153,4 @@ case FrictionVoidRatioEnum : return "FrictionVoidRatio"; case FrontalForcingsBasinIcefrontAreaEnum : return "FrontalForcingsBasinIcefrontArea"; - case FrontalForcingsBasinIdEnum : return "FrontalForcingsBasinId"; case FrontalForcingsNumberofBasinsEnum : return "FrontalForcingsNumberofBasins"; case FrontalForcingsParamEnum : return "FrontalForcingsParam"; @@ -538,4 +537,5 @@ case FrictionWaterLayerEnum : return "FrictionWaterLayer"; case FrictionfEnum : return "Frictionf"; + case FrontalForcingsBasinIdEnum : return "FrontalForcingsBasinId"; case FrontalForcingsSubglacialDischargeEnum : return "FrontalForcingsSubglacialDischarge"; case FrontalForcingsThermalForcingEnum : return "FrontalForcingsThermalForcing"; @@ -1212,4 +1212,6 @@ case ThermalSolutionEnum : return "ThermalSolution"; case ThicknessErrorEstimatorEnum : return "ThicknessErrorEstimator"; + case TotalCalvingFluxEnum : return "TotalCalvingFlux"; + case TotalCalvingFluxLevelsetEnum : return "TotalCalvingFluxLevelset"; case TotalFloatingBmbEnum : return "TotalFloatingBmb"; case TotalFloatingBmbScaledEnum : return "TotalFloatingBmbScaled"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 23990) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 23991) @@ -156,5 +156,4 @@ else if (strcmp(name,"FrictionVoidRatio")==0) return FrictionVoidRatioEnum; else if (strcmp(name,"FrontalForcingsBasinIcefrontArea")==0) return FrontalForcingsBasinIcefrontAreaEnum; - else if (strcmp(name,"FrontalForcingsBasinId")==0) return FrontalForcingsBasinIdEnum; else if (strcmp(name,"FrontalForcingsNumberofBasins")==0) return FrontalForcingsNumberofBasinsEnum; else if (strcmp(name,"FrontalForcingsParam")==0) return FrontalForcingsParamEnum; @@ -260,9 +259,9 @@ else if (strcmp(name,"MaterialsMeltingpoint")==0) return MaterialsMeltingpointEnum; else if (strcmp(name,"MaterialsMixedLayerCapacity")==0) return MaterialsMixedLayerCapacityEnum; + else if (strcmp(name,"MaterialsMuWater")==0) return MaterialsMuWaterEnum; else stage=3; } if(stage==3){ - if (strcmp(name,"MaterialsMuWater")==0) return MaterialsMuWaterEnum; - else if (strcmp(name,"MaterialsRheologyLaw")==0) return MaterialsRheologyLawEnum; + if (strcmp(name,"MaterialsRheologyLaw")==0) return MaterialsRheologyLawEnum; else if (strcmp(name,"MaterialsRhoFreshwater")==0) return MaterialsRhoFreshwaterEnum; else if (strcmp(name,"MaterialsRhoIce")==0) return MaterialsRhoIceEnum; @@ -383,9 +382,9 @@ else if (strcmp(name,"SmbT0wet")==0) return SmbT0wetEnum; else if (strcmp(name,"SmbTdiff")==0) return SmbTdiffEnum; + else if (strcmp(name,"SmbThermoDeltaTScaling")==0) return SmbThermoDeltaTScalingEnum; else stage=4; } if(stage==4){ - if (strcmp(name,"SmbThermoDeltaTScaling")==0) return SmbThermoDeltaTScalingEnum; - else if (strcmp(name,"SmoothThicknessMultiplier")==0) return SmoothThicknessMultiplierEnum; + if (strcmp(name,"SmoothThicknessMultiplier")==0) return SmoothThicknessMultiplierEnum; else if (strcmp(name,"SolutionType")==0) return SolutionTypeEnum; else if (strcmp(name,"SteadystateMaxiter")==0) return SteadystateMaxiterEnum; @@ -506,9 +505,9 @@ else if (strcmp(name,"Domain2Dhorizontal")==0) return Domain2DhorizontalEnum; else if (strcmp(name,"Domain2Dvertical")==0) return Domain2DverticalEnum; + else if (strcmp(name,"Domain3D")==0) return Domain3DEnum; else stage=5; } if(stage==5){ - if (strcmp(name,"Domain3D")==0) return Domain3DEnum; - else if (strcmp(name,"DragCoefficientAbsGradient")==0) return DragCoefficientAbsGradientEnum; + if (strcmp(name,"DragCoefficientAbsGradient")==0) return DragCoefficientAbsGradientEnum; else if (strcmp(name,"DrivingStressX")==0) return DrivingStressXEnum; else if (strcmp(name,"DrivingStressY")==0) return DrivingStressYEnum; @@ -550,4 +549,5 @@ else if (strcmp(name,"FrictionWaterLayer")==0) return FrictionWaterLayerEnum; else if (strcmp(name,"Frictionf")==0) return FrictionfEnum; + else if (strcmp(name,"FrontalForcingsBasinId")==0) return FrontalForcingsBasinIdEnum; else if (strcmp(name,"FrontalForcingsSubglacialDischarge")==0) return FrontalForcingsSubglacialDischargeEnum; else if (strcmp(name,"FrontalForcingsThermalForcing")==0) return FrontalForcingsThermalForcingEnum; @@ -1239,13 +1239,15 @@ else if (strcmp(name,"ThermalSolution")==0) return ThermalSolutionEnum; else if (strcmp(name,"ThicknessErrorEstimator")==0) return ThicknessErrorEstimatorEnum; + else if (strcmp(name,"TotalCalvingFlux")==0) return TotalCalvingFluxEnum; + else if (strcmp(name,"TotalCalvingFluxLevelset")==0) return TotalCalvingFluxLevelsetEnum; else if (strcmp(name,"TotalFloatingBmb")==0) return TotalFloatingBmbEnum; else if (strcmp(name,"TotalFloatingBmbScaled")==0) return TotalFloatingBmbScaledEnum; else if (strcmp(name,"TotalGroundedBmb")==0) return TotalGroundedBmbEnum; - else if (strcmp(name,"TotalGroundedBmbScaled")==0) return TotalGroundedBmbScaledEnum; - else if (strcmp(name,"TotalSmb")==0) return TotalSmbEnum; else stage=11; } if(stage==11){ - if (strcmp(name,"TotalSmbScaled")==0) return TotalSmbScaledEnum; + if (strcmp(name,"TotalGroundedBmbScaled")==0) return TotalGroundedBmbScaledEnum; + else if (strcmp(name,"TotalSmb")==0) return TotalSmbEnum; + else if (strcmp(name,"TotalSmbScaled")==0) return TotalSmbScaledEnum; else if (strcmp(name,"TransientArrayParam")==0) return TransientArrayParamEnum; else if (strcmp(name,"TransientInput")==0) return TransientInputEnum;