Index: /issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp (revision 28049) +++ /issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp (revision 28050) @@ -151,4 +151,6 @@ iomodel->FetchDataToInput(inputs,elements,"md.initialization.watercolumn",HydrologySheetThicknessEnum); iomodel->FetchDataToInput(inputs,elements,"md.initialization.hydraulic_potential",HydraulicPotentialEnum); + iomodel->FetchDataToInput(inputs,elements,"md.hydrology.rheology_B_base",HydrologyRheologyBBaseEnum); + if(iomodel->domaintype==Domain2DhorizontalEnum){ iomodel->FetchDataToInput(inputs,elements,"md.initialization.vx",VxEnum); @@ -249,12 +251,12 @@ IssmDouble g = element->FindParam(ConstantsGEnum); IssmDouble e_v = element->FindParam(HydrologyEnglacialVoidRatioEnum); - Input* hr_input = element->GetInput(HydrologyBumpHeightEnum);_assert_(hr_input); + Input* hr_input = element->GetInput(HydrologyBumpHeightEnum); _assert_(hr_input); Input* k_input = element->GetInput(HydrologySheetConductivityEnum);_assert_(k_input); - Input* phi_input = element->GetInput(HydraulicPotentialEnum); _assert_(phi_input); - Input* h_input = element->GetInput(HydrologySheetThicknessEnum); _assert_(h_input); - Input* H_input = element->GetInput(ThicknessEnum); _assert_(H_input); - Input* b_input = element->GetInput(BedEnum); _assert_(b_input); - Input* B_input = element->GetInput(MaterialsRheologyBEnum); _assert_(B_input); - Input* n_input = element->GetInput(MaterialsRheologyNEnum); _assert_(n_input); + Input* phi_input = element->GetInput(HydraulicPotentialEnum); _assert_(phi_input); + Input* h_input = element->GetInput(HydrologySheetThicknessEnum); _assert_(h_input); + Input* H_input = element->GetInput(ThicknessEnum); _assert_(H_input); + Input* b_input = element->GetInput(BedEnum); _assert_(b_input); + Input* B_input = element->GetInput(HydrologyRheologyBBaseEnum); _assert_(B_input); + Input* n_input = element->GetInput(MaterialsRheologyNEnum); _assert_(n_input); /* Start looping on the number of gaussian points: */ @@ -275,7 +277,4 @@ b_input->GetInputValue(&b,gauss); H_input->GetInputValue(&H,gauss); - - /*Hard code B*/ - B = Cuffey(273.15-2); /*Get norm of gradient of hydraulic potential and make sure it is >0*/ @@ -364,15 +363,15 @@ IssmDouble g = element->FindParam(ConstantsGEnum); IssmDouble e_v = element->FindParam(HydrologyEnglacialVoidRatioEnum); - Input* hr_input = element->GetInput(HydrologyBumpHeightEnum);_assert_(hr_input); - Input* h_input = element->GetInput(HydrologySheetThicknessEnum);_assert_(h_input); - Input* H_input = element->GetInput(ThicknessEnum); _assert_(H_input); - Input* b_input = element->GetInput(BedEnum); _assert_(b_input); - Input* G_input = element->GetInput(BasalforcingsGeothermalfluxEnum);_assert_(G_input); + Input* hr_input = element->GetInput(HydrologyBumpHeightEnum); _assert_(hr_input); + Input* h_input = element->GetInput(HydrologySheetThicknessEnum); _assert_(h_input); + Input* H_input = element->GetInput(ThicknessEnum); _assert_(H_input); + Input* b_input = element->GetInput(BedEnum); _assert_(b_input); + Input* G_input = element->GetInput(BasalforcingsGeothermalfluxEnum); _assert_(G_input); Input* melt_input = element->GetInput(BasalforcingsGroundediceMeltingRateEnum);_assert_(melt_input); Input* RO_input = NULL; - Input* B_input = element->GetInput(MaterialsRheologyBEnum); _assert_(B_input); - Input* n_input = element->GetInput(MaterialsRheologyNEnum); _assert_(n_input); - Input* phiold_input = element->GetInput(HydraulicPotentialOldEnum); _assert_(phiold_input); - Input* phi_input = element->GetInput(HydraulicPotentialEnum); _assert_(phi_input); + Input* B_input = element->GetInput(HydrologyRheologyBBaseEnum); _assert_(B_input); + Input* n_input = element->GetInput(MaterialsRheologyNEnum); _assert_(n_input); + Input* phiold_input = element->GetInput(HydraulicPotentialOldEnum); _assert_(phiold_input); + Input* phi_input = element->GetInput(HydraulicPotentialEnum); _assert_(phi_input); /*Build friction element, needed later: */ @@ -397,7 +396,4 @@ melt_input->GetInputValue(&melt,gauss); - /*Hard code B*/ - B = Cuffey(273.15-2); - /*Get basal velocity*/ friction->GetBasalSlidingSpeeds(&vx, &vy ,gauss); @@ -457,9 +453,9 @@ element->InputUpdateFromSolutionOneDof(solution,HydraulicPotentialEnum); - /*Compute Hydrology Vx and Vy for time stepping purposes (These inputs do not affect GlaDS)*/ + /*Compute Hydrology Vx and Vy for time stepping purposes, and Sheet Discharge as an optional output (These inputs do not affect GlaDS)*/ /*Intermediaries*/ IssmDouble dphi[3],h,k,phi; - IssmDouble h_r; + IssmDouble h_r; IssmDouble oceanLS,iceLS; IssmDouble* xyz_list = NULL; @@ -468,7 +464,8 @@ int numvertices = element->GetNumberOfVertices(); - /*Initialize new sheet thickness*/ + /*Initialize water sheet velocity and discharge*/ IssmDouble* vx = xNew(numvertices); IssmDouble* vy = xNew(numvertices); + IssmDouble* d = xNew(numvertices); /*Set to 0 if inactive element*/ @@ -476,8 +473,11 @@ for(int iv=0;ivAddInput(HydrologyWaterVxEnum,vx,P1DGEnum); element->AddInput(HydrologyWaterVyEnum,vy,P1DGEnum); + element->AddInput(HydrologySheetDischargeEnum,d,P1DGEnum); xDelete(vx); xDelete(vy); + xDelete(d); return; } @@ -499,5 +499,5 @@ Input *iceLS_input = element->GetInput(MaskIceLevelsetEnum); _assert_(iceLS_input); - /* Start looping on the number of gaussian points: */ + /* Start looping on the number of gaussian points: */ Gauss* gauss=element->NewGauss(); for(int iv=0;ivGetInputValue(&phi,gauss); h_input->GetInputValue(&h,gauss); - hr_input->GetInputValue(&h_r,gauss); + hr_input->GetInputValue(&h_r,gauss); k_input->GetInputValue(&k,gauss); oceanLS_input->GetInputValue(&oceanLS,gauss); iceLS_input->GetInputValue(&iceLS,gauss); - /*Set sheet thickness to zero if floating or no ice*/ + /*Set to zero if floating or no ice*/ if(oceanLS<0. || iceLS>0.){ vx[iv] = 0.; vy[iv] = 0.; + d[iv] = 0.; } else{ @@ -526,17 +527,18 @@ /*If omega is zero, use standard model, otherwise transition model*/ IssmDouble nu = mu_water/rho_water; - IssmDouble coeff; - if(istransition==1 && omega>=AEPS){ - IssmDouble hratio = fabs(h/h_r); - IssmDouble coarg = 1. + 4.*pow(hratio,3-2*alpha)*omega*k*pow(h,3)*normgradphi/nu; - coeff = nu/2./omega*pow(hratio,2*alpha-3) * (-1 + pow(coarg, 0.5))/normgradphi/max(AEPS,h); // divide by h to get speed instead of discharge - } - else { - coeff = k*pow(h,alpha)*pow(normgradphi,beta-2.)/max(AEPS,h); // divide by h to get speed instead of discharge - } - - - vx[iv] = -coeff*dphi[0]; - vy[iv] = -coeff*dphi[1]; + IssmDouble coeff; + if(istransition==1 && omega>=AEPS){ + IssmDouble hratio = fabs(h/h_r); + IssmDouble coarg = 1. + 4.*pow(hratio,3-2*alpha)*omega*k*pow(h,3)*normgradphi/nu; + coeff = nu/2./omega*pow(hratio,2*alpha-3) * (-1 + pow(coarg, 0.5))/normgradphi; // coeff gives discharge; divide by h to get speed instead of discharge + } + else { + coeff = k*pow(h,alpha)*pow(normgradphi,beta-2.); // coeff gives discharge; divide by h to get speed instead of discharge + } + + vx[iv] = -coeff/max(AEPS,h)*dphi[0]; + vy[iv] = -coeff/max(AEPS,h)*dphi[1]; + + d[iv] = coeff*normgradphi; } } @@ -544,4 +546,5 @@ element->AddInput(HydrologyWaterVxEnum,vx,P1DGEnum); element->AddInput(HydrologyWaterVyEnum,vy,P1DGEnum); + element->AddInput(HydrologySheetDischargeEnum,d,P1DGEnum); /*Clean up and return*/ @@ -549,4 +552,5 @@ xDelete(vx); xDelete(vy); + xDelete(d); delete gauss; }/*}}}*/ @@ -646,16 +650,15 @@ IssmDouble rho_water = element->FindParam(MaterialsRhoFreshwaterEnum); IssmDouble g = element->FindParam(ConstantsGEnum); - Input* hr_input = element->GetInput(HydrologyBumpHeightEnum);_assert_(hr_input); - Input* vx_input = element->GetInput(VxBaseEnum);_assert_(vx_input); - Input* vy_input = element->GetInput(VyBaseEnum);_assert_(vy_input); - Input* H_input = element->GetInput(ThicknessEnum); _assert_(H_input); - Input* b_input = element->GetInput(BedEnum); _assert_(b_input); - Input* hold_input = element->GetInput(HydrologySheetThicknessOldEnum);_assert_(hold_input); - Input* B_input = element->GetInput(MaterialsRheologyBEnum); _assert_(B_input); - Input* n_input = element->GetInput(MaterialsRheologyNEnum); _assert_(n_input); + Input* hr_input = element->GetInput(HydrologyBumpHeightEnum); _assert_(hr_input); + Input* vx_input = element->GetInput(VxBaseEnum); _assert_(vx_input); + Input* vy_input = element->GetInput(VyBaseEnum); _assert_(vy_input); + Input* H_input = element->GetInput(ThicknessEnum); _assert_(H_input); + Input* b_input = element->GetInput(BedEnum); _assert_(b_input); + Input* hold_input = element->GetInput(HydrologySheetThicknessOldEnum);_assert_(hold_input); + Input* B_input = element->GetInput(HydrologyRheologyBBaseEnum); _assert_(B_input); + Input* n_input = element->GetInput(MaterialsRheologyNEnum); _assert_(n_input); Input* phi_input = element->GetInput(HydraulicPotentialEnum); _assert_(phi_input); - Input* oceanLS_input = element->GetInput(MaskOceanLevelsetEnum); _assert_(oceanLS_input); - Input* iceLS_input = element->GetInput(MaskIceLevelsetEnum); _assert_(iceLS_input); - + Input* oceanLS_input = element->GetInput(MaskOceanLevelsetEnum); _assert_(oceanLS_input); + Input* iceLS_input = element->GetInput(MaskIceLevelsetEnum); _assert_(iceLS_input); /* Start looping on the number of gaussian points: */ @@ -677,7 +680,4 @@ iceLS_input->GetInputValue(&iceLS,gauss); - /*Hard code B*/ - B = Cuffey(273.15-2); - /*Set sheet thickness to zero if floating or no ice*/ if(oceanLS<0. || iceLS>0.){ Index: /issm/trunk-jpl/src/c/classes/Loads/Channel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Channel.cpp (revision 28049) +++ /issm/trunk-jpl/src/c/classes/Loads/Channel.cpp (revision 28050) @@ -354,5 +354,5 @@ /*Initialize Element matrix and return if necessary*/ Tria* tria=(Tria*)element; - if(!tria->IsIceInElement()) return NULL; + if(!tria->IsIceOnlyInElement()) return NULL; _assert_(tria->FiniteElement()==P1Enum); int index1=tria->GetVertexIndex(vertices[0]); @@ -396,9 +396,9 @@ Input* H_input = element->GetInput(ThicknessEnum); _assert_(H_input); Input* b_input = element->GetInput(BedEnum); _assert_(b_input); - Input* B_input = element->GetInput(MaterialsRheologyBEnum); _assert_(B_input); + Input* B_input = element->GetInput(HydrologyRheologyBBaseEnum); _assert_(B_input); Input* n_input = element->GetInput(MaterialsRheologyNEnum); _assert_(n_input); Input* ks_input = element->GetInput(HydrologySheetConductivityEnum); _assert_(ks_input); Input* kc_input = element->GetInput(HydrologyChannelConductivityEnum); _assert_(kc_input); - Input* hr_input = element->GetInput(HydrologyBumpHeightEnum);_assert_(hr_input); + Input* hr_input = element->GetInput(HydrologyBumpHeightEnum); _assert_(hr_input); Input* phi_input = element->GetInput(HydraulicPotentialEnum); _assert_(phi_input); @@ -432,7 +432,4 @@ b_input->GetInputValue(&b,gauss); H_input->GetInputValue(&H,gauss); - - /*Hard code B*/ - B = Cuffey(273.15-2); /*Get values for a few potentials*/ @@ -493,7 +490,15 @@ } - /*Closing rate term, see Gagliardini and Werder 2018 eq. A2 (v = v1*phi_i + v2(phi_{i+1}))*/ - A=pow(B,-n); - v1 = 2./pow(n,n)*A*S*(pow(fabs(phi_0 - phi),n-1.)*( - n)); + /*Closing rate term*/ + /*See Gagliardini and Werder 2018 eq. A2 (v = v1*phi_i + v2(phi_{i+1}))*/ + A = pow(B,-n); + if(phi_0-phi<0){ + v1 = 0.; + } + else{ + v1 = 2./pow(n,n)*A*S*(pow(fabs(phi_0-phi),n-1.)*( - n)); + + } + for(int i=0;iIsIceInElement()) return NULL; + if(!tria->IsIceOnlyInElement()) return NULL; _assert_(tria->FiniteElement()==P1Enum); int index1=tria->GetVertexIndex(vertices[0]); @@ -550,10 +555,10 @@ Input* H_input = element->GetInput(ThicknessEnum); _assert_(H_input); Input* b_input = element->GetInput(BedEnum); _assert_(b_input); - Input* B_input = element->GetInput(MaterialsRheologyBEnum); _assert_(B_input); + Input* B_input = element->GetInput(HydrologyRheologyBBaseEnum); _assert_(B_input); Input* n_input = element->GetInput(MaterialsRheologyNEnum); _assert_(n_input); Input* ks_input = element->GetInput(HydrologySheetConductivityEnum); _assert_(ks_input); Input* kc_input = element->GetInput(HydrologyChannelConductivityEnum); _assert_(kc_input); Input* phi_input = element->GetInput(HydraulicPotentialEnum); _assert_(phi_input); - Input* hr_input = element->GetInput(HydrologyBumpHeightEnum);_assert_(hr_input); + Input* hr_input = element->GetInput(HydrologyBumpHeightEnum); _assert_(hr_input); /*Get tangent vector*/ @@ -584,6 +589,4 @@ hr_input->GetInputValue(&h_r,gauss); - /*Hard code B*/ - B = Cuffey(273.15-2); /*Get values for a few potentials*/ @@ -625,6 +628,11 @@ /*Compute closing rate*/ /*See Gagliardini and Werder 2018 eq. A2 (v = v2(phi_i) + v1*phi_{i+1})*/ - A=pow(B,-n); - v2 = 2./pow(n,n)*A*this->S*(pow(fabs(phi_0 - phi),n-1.)*(phi_0 +(n-1.)*phi)); + A = pow(B,-n); + if(phi_0-phi<0){ + v2 = 0.; + } + else{ + v2 = 2./pow(n,n)*A*this->S*(pow(fabs(phi_0 - phi),n-1.)*(phi_0 +(n-1.)*phi)); + } for(int i=0;iboundary){ + if(this->boundary || !tria->IsIceOnlyInElement()){ this->S = 0.; return; @@ -700,10 +708,10 @@ Input* H_input = element->GetInput(ThicknessEnum); _assert_(H_input); Input* b_input = element->GetInput(BedEnum); _assert_(b_input); - Input* B_input = element->GetInput(MaterialsRheologyBEnum); _assert_(B_input); + Input* B_input = element->GetInput(HydrologyRheologyBBaseEnum); _assert_(B_input); Input* n_input = element->GetInput(MaterialsRheologyNEnum); _assert_(n_input); Input* ks_input = element->GetInput(HydrologySheetConductivityEnum); _assert_(ks_input); Input* kc_input = element->GetInput(HydrologyChannelConductivityEnum); _assert_(kc_input); Input* phi_input = element->GetInput(HydraulicPotentialEnum); _assert_(phi_input); - Input* hr_input = element->GetInput(HydrologyBumpHeightEnum);_assert_(hr_input); + Input* hr_input = element->GetInput(HydrologyBumpHeightEnum); _assert_(hr_input); /*Get tangent vector*/ @@ -728,7 +736,4 @@ - /*Hard code B*/ - B = Cuffey(273.15-2); - /*Get values for a few potentials*/ phi_0 = rho_water*g*b + rho_ice*g*H; @@ -754,5 +759,5 @@ /*Ice rate factor*/ - A=pow(B,-n); + A = pow(B,-n); IssmDouble C = C_W*c_t*rho_water; @@ -777,4 +782,10 @@ + C*Qprime*pow(Snew,alpha_c-1.)*dPw ) - 2./pow(n,n)*A*pow(fabs(N),n-1.)*N; + if(N<0){ + alpha = 1./(rho_ice*L)*( + fabs(Qprime*pow(Snew,alpha_c-1.)*dphids) + + C*Qprime*pow(Snew,alpha_c-1.)*dPw + ); + } IssmDouble beta = 1./(rho_ice*L)*( fabs(lc*qc*dphids) + C*fFactor*dPw ); @@ -786,8 +797,5 @@ /*Constrain the cross section to be between 0 and 500 m^2*/ if(this->S<0.) this->S = 0.; - if(this->S>100.) this->S = 100.; - - /*Do not allow channels to grow in areas with no sheet thickness*/ - if(H<200.) this->S = 0.; + if(this->S>500.) this->S = 500.; count++; Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 28049) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 28050) @@ -899,5 +899,7 @@ HydrologyNeumannfluxEnum, HydrologyReynoldsEnum, + HydrologyRheologyBBaseEnum, HydrologySheetConductivityEnum, + HydrologySheetDischargeEnum, HydrologySheetThicknessEnum, HydrologySheetThicknessOldEnum, Index: /issm/trunk-jpl/src/m/classes/hydrologyglads.m =================================================================== --- /issm/trunk-jpl/src/m/classes/hydrologyglads.m (revision 28049) +++ /issm/trunk-jpl/src/m/classes/hydrologyglads.m (revision 28050) @@ -14,4 +14,5 @@ sheet_alpha = NaN; sheet_beta = NaN; + rheology_B_base = NaN; %Channels @@ -83,5 +84,5 @@ md = checkfield(md,'fieldname','hydrology.sheet_alpha', 'numel', [1], '>', 0); md = checkfield(md,'fieldname','hydrology.sheet_beta', 'numel', [1], '>', 0); - + md = checkfield(md,'fieldname','hydrology.rheology_B_base','size',[md.mesh.numberofvertices 1],'>=',0,'NaN',1,'Inf',1); %Channels md = checkfield(md,'fieldname','hydrology.ischannels','numel',[1],'values',[0 1]); @@ -113,4 +114,5 @@ fielddisplay(self,'bump_height','typical bump height (h_r) [m]'); fielddisplay(self,'omega','transition parameter (omega) []'); + fielddisplay(self,'rheology_B_base','Ice rheology factor B at base of ice (B) [Pa s^(-1/3)]'); disp(sprintf(' CHANNELS')); fielddisplay(self,'ischannels','Do we allow for channels? 1: yes, 0: no'); @@ -143,4 +145,5 @@ WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','sheet_alpha','format','Double'); WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','sheet_beta','format','Double'); + WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','rheology_B_base','format','DoubleMat','mattype',1); %Channels Index: /issm/trunk-jpl/src/m/classes/hydrologyglads.py =================================================================== --- /issm/trunk-jpl/src/m/classes/hydrologyglads.py (revision 28049) +++ /issm/trunk-jpl/src/m/classes/hydrologyglads.py (revision 28050) @@ -22,4 +22,5 @@ self.sheet_alpha = np.nan; self.sheet_beta = np.nan; + self.rheology_B_base = np.nan; # Channels @@ -59,4 +60,5 @@ s += '{}\n'.format(fielddisplay(self, 'bump_height', 'typical bump height (h_r) [m]')) s += '{}\n'.format(fielddisplay(self, 'omega', 'transition parameter (omega) []')) #TH + s += '{}\n'.format(fielddisplay(self, 'rheology_B_base', 'ice rheology factor B at base of ice (B) [Pa s^(-1/3)]')) #SE s = '\t--CHANNELS\n' s += '{}\n'.format(fielddisplay(self, 'ischannels', 'Do we allow for channels? 1: yes, 0: no')) @@ -129,4 +131,5 @@ md = checkfield(md,'fieldname','hydrology.sheet_alpha', 'numel', [1], '>', 0); md = checkfield(md,'fieldname','hydrology.sheet_beta', 'numel', [1], '>', 0); + md = checkfield(md,'fieldname','hydrology.rheology_B_base', 'size', [md.mesh.numberofvertices], '>=', 0, 'np.nan', 1, 'Inf', 1) # Channels @@ -162,4 +165,5 @@ WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','sheet_alpha','format','Double'); WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','sheet_beta','format','Double'); + WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','rheology_B_base','format','DoubleMat', 'mattype', 1); # Channels Index: /issm/trunk-jpl/test/NightlyRun/test355.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test355.m (revision 28049) +++ /issm/trunk-jpl/test/NightlyRun/test355.m (revision 28050) @@ -59,4 +59,5 @@ md.hydrology.sheet_conductivity= 1e-3 * ones(md.mesh.numberofvertices,1); md.hydrology.channel_conductivity= 5.e-2 * ones(md.mesh.numberofvertices,1); +md.hydrology.rheology_B_base = cuffey(273.15 - 2)*ones(md.mesh.numberofvertices,1); % BCs for hydrology @@ -78,5 +79,5 @@ 2e-10,2e-08,4e-07,... 3e-10,2e-08,4e-07,... - 4e-10,1e-08,4e-07}; + 4e-10,2e-08,4e-07}; field_values={... md.results.TransientSolution(1).HydrologySheetThickness, ... Index: /issm/trunk-jpl/test/NightlyRun/test355.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test355.py (revision 28049) +++ /issm/trunk-jpl/test/NightlyRun/test355.py (revision 28050) @@ -68,4 +68,5 @@ md.hydrology.sheet_conductivity = 1.e-3 * np.ones((md.mesh.numberofvertices)) md.hydrology.channel_conductivity = 5.e-2 * np.ones((md.mesh.numberofvertices)) +md.hydrology.rheology_B_base = 8.3788e7 * np.ones((md.mesh.numberofvertices)) # BCs for hydrology