Index: /issm/trunk-jpl/src/c/analyses/HydrologyShaktiAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyShaktiAnalysis.cpp (revision 26489) +++ /issm/trunk-jpl/src/c/analyses/HydrologyShaktiAnalysis.cpp (revision 26490) @@ -181,4 +181,5 @@ IssmDouble Jdet; IssmDouble* xyz_list = NULL; + IssmDouble gap,bed,thickness,head,g,rho_ice,rho_water,A,B,n,head_old;/***/ /*Fetch number of nodes and dof for this finite element*/ @@ -201,4 +202,16 @@ element->FindParam(&dt,TimesteppingTimeStepEnum); + /***Get all inputs and parameters*/ + element->FindParam(&rho_water,MaterialsRhoFreshwaterEnum); + element->FindParam(&rho_ice,MaterialsRhoIceEnum); + element->FindParam(&g,ConstantsGEnum); + Input* B_input = element->GetInput(MaterialsRheologyBEnum); _assert_(B_input); + Input* n_input = element->GetInput(MaterialsRheologyNEnum); _assert_(n_input); + Input* gap_input = element->GetInput(HydrologyGapHeightEnum); _assert_(gap_input); + Input* thickness_input = element->GetInput(ThicknessEnum); _assert_(thickness_input); + Input* head_input = element->GetInput(HydrologyHeadEnum); _assert_(head_input); + Input* base_input = element->GetInput(BaseEnum); _assert_(base_input); + Input* head_old_input = element->GetInput(HydrologyHeadOldEnum); _assert_(head_old_input); + /* Start looping on the number of gaussian points: */ Gauss* gauss=element->NewGauss(1); @@ -209,8 +222,31 @@ element->NodalFunctions(basis,gauss); + /***/ + base_input->GetInputValue(&bed,gauss); + thickness_input->GetInputValue(&thickness,gauss); + gap_input->GetInputValue(&gap,gauss); + head_input->GetInputValue(&head,gauss); + head_old_input->GetInputValue(&head_old,gauss); + + /***Get ice A parameter*/ + B_input->GetInputValue(&B,gauss); + n_input->GetInputValue(&n,gauss); + A=pow(B,-n); + + /***Get water and ice pressures*/ + IssmDouble pressure_ice = rho_ice*g*thickness; _assert_(pressure_ice>0.); + IssmDouble pressure_water = rho_water*g*(head-bed); + if(pressure_water>pressure_ice) pressure_water = pressure_ice; + + /***Get water pressure from previous time step to use in lagged creep term*/ + IssmDouble pressure_water_old = rho_water*g*(head_old-bed); + if(pressure_water_old>pressure_ice) pressure_water_old = pressure_ice; + + for(int i=0;ivalues[i*numnodes+j] += conductivity*gauss->weight*Jdet*(dbasis[0*numnodes+i]*dbasis[0*numnodes+j] + dbasis[1*numnodes+i]*dbasis[1*numnodes+j]) - + gauss->weight*Jdet*storage/dt*basis[i]*basis[j]; + + gauss->weight*Jdet*storage/dt*basis[i]*basis[j] + +gauss->weight*Jdet*A*(n)*(pow(fabs(pressure_ice-pressure_water),(n-1))*rho_water*g)*gap*basis[i]*basis[j]; } } @@ -232,8 +268,9 @@ IssmDouble Jdet,meltrate,G,dh[2],B,A,n; IssmDouble gap,bed,thickness,head,ieb,head_old; - IssmDouble lr,br,vx,vy,beta; + IssmDouble lr,br,vx,vy,beta,lc; IssmDouble alpha2,frictionheat; IssmDouble PMPheat,dpressure_water[2],dbed[2]; IssmDouble* xyz_list = NULL; + IssmDouble dgap[2]; /***/ /*Fetch number of nodes and dof for this finite element*/ @@ -261,4 +298,6 @@ Input* br_input = element->GetInput(HydrologyBumpHeightEnum); _assert_(br_input); Input* headold_input = element->GetInput(HydrologyHeadOldEnum); _assert_(headold_input); +// Input* dgapxx_input = element->GetInput(HydrologyGapXXEnum); /***/ +// Input* dgapyy_input = element->GetInput(HydrologyGapYYEnum); /***/ /*Get conductivity from inputs*/ @@ -289,5 +328,5 @@ lr_input->GetInputValue(&lr,gauss); br_input->GetInputValue(&br,gauss); - headold_input->GetInputValue(&head_old,gauss); + headold_input->GetInputValue(&head_old,gauss); /*Get ice A parameter*/ @@ -302,4 +341,12 @@ beta = 0.; + /* TEST Compute lc term */ +// if(gapGetAlpha2(&alpha2,gauss); @@ -324,12 +371,35 @@ _assert_(meltrate>0.); - for(int i=0;ivalues[i]+=Jdet*gauss->weight* - ( - meltrate*(1/rho_water-1/rho_ice) - +A*pow(fabs(pressure_ice - pressure_water),n-1)*(pressure_ice - pressure_water)*gap - -beta*sqrt(vx*vx+vy*vy) - +ieb - +storage*head_old/dt - )*basis[i]; + +// for(int i=0;ivalues[i]+=Jdet*gauss->weight* +// ( +// meltrate*(1/rho_water-1/rho_ice) +// +A*pow(fabs(pressure_ice - pressure_water),n-1)*(pressure_ice - pressure_water)*gap +// -beta*sqrt(vx*vx+vy*vy) +// +ieb +// +storage*head_old/dt +// )*basis[i]; + + /* With weighted creep term*/ + for(int i=0;ivalues[i]+=Jdet*gauss->weight* + ( + meltrate*(1/rho_water-1/rho_ice) + +A*pow(fabs(pressure_ice - pressure_water),n-1)*(pressure_ice + rho_water*g*bed)*gap + +(n-1)*A*pow(fabs(pressure_ice - pressure_water),n-1)*(rho_water*g*head)*gap + -beta*sqrt(vx*vx+vy*vy) + +ieb + +storage*head_old/dt + )*basis[i]; + + +// /*Test with linear creep rate*/ + // for(int i=0;ivalues[i]+=Jdet*gauss->weight* + // ( + // meltrate*(1/rho_water-1/rho_ice) + // +(pressure_ice - pressure_water)*gap/pow(10,13) + // -beta*sqrt(vx*vx+vy*vy) + // +ieb + // +storage*head_old/dt + // )*basis[i]; } /*Clean up and return*/ @@ -382,13 +452,13 @@ values[i]=solution[doflist[i]]; - /*make sure that p_water hrho_ice*thickness[i]/rho_water+bed[i]){ - values[i] = rho_ice*thickness[i]/rho_water+bed[i]; - } + /*overburden cap: make sure that p_water hrho_ice*thickness[i]/rho_water+bed[i]){ +// values[i] = rho_ice*thickness[i]/rho_water+bed[i]; +// } /*Make sure that negative pressure is not allowed*/ - // if(values[i]GetAlpha2(&alpha2,gauss); @@ -555,4 +633,20 @@ +beta*sqrt(vx*vx+vy*vy) )); + + /* TEST with gap height "diffusivity" for melting walls - not yet completed*/ +// newgap += gauss->weight*Jdet*(gap+dt*( +// meltrate/rho_ice +// -A*pow(fabs(pressure_ice-pressure_water),n-1)*(pressure_ice-pressure_water)*lc +// +beta*sqrt(vx*vx+vy*vy) +// +(meltrate*gap/rho_ice*(dgapxx+dgapyy)) +// )); + + + /*TEST with linear creep rate*/ +// newgap += gauss->weight*Jdet*(gap+dt*( + // meltrate/rho_ice + // -(pressure_ice-pressure_water)*gap/pow(10,13) + // +beta*sqrt(vx*vx+vy*vy) + // )); totalweights +=gauss->weight*Jdet;