Index: /issm/trunk-jpl/src/c/analyses/HydrologyShaktiAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyShaktiAnalysis.cpp (revision 26695) +++ /issm/trunk-jpl/src/c/analyses/HydrologyShaktiAnalysis.cpp (revision 26696) @@ -202,5 +202,5 @@ IssmDouble Jdet; IssmDouble* xyz_list = NULL; - IssmDouble gap,bed,thickness,head,g,rho_ice,rho_water,A,B,n,head_old;/***/ + IssmDouble gap,bed,thickness,head,g,rho_ice,rho_water,A,B,n; /*Fetch number of nodes and dof for this finite element*/ @@ -223,5 +223,5 @@ element->FindParam(&dt,TimesteppingTimeStepEnum); - /***Get all inputs and parameters*/ + /*Get all inputs and parameters*/ element->FindParam(&rho_water,MaterialsRhoFreshwaterEnum); element->FindParam(&rho_ice,MaterialsRhoIceEnum); @@ -233,5 +233,4 @@ 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: */ @@ -248,19 +247,15 @@ gap_input->GetInputValue(&gap,gauss); head_input->GetInputValue(&head,gauss); - head_old_input->GetInputValue(&head_old,gauss); - - /***Get ice A parameter*/ + + /*Get ice A parameter*/ B_input->GetInputValue(&B,gauss); n_input->GetInputValue(&n,gauss); A=pow(B,-n); - /***Get water and ice pressures*/ + /*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; @@ -293,6 +288,6 @@ IssmDouble PMPheat,dpressure_water[2],dbed[2]; IssmDouble* xyz_list = NULL; - IssmDouble dgapxx; /***/ - IssmDouble dgapyy; /***/ +// IssmDouble dgapxx; /***/ +// IssmDouble dgapyy; /***/ /*Fetch number of nodes and dof for this finite element*/ @@ -320,6 +315,6 @@ Input* br_input = element->GetInput(HydrologyBumpHeightEnum); _assert_(br_input); Input* headold_input = element->GetInput(HydrologyHeadOldEnum); _assert_(headold_input); - Input* dgapxx_input = element->GetInput(HydrologyGapHeightXXEnum); /***/ - Input* dgapyy_input = element->GetInput(HydrologyGapHeightYYEnum); /***/ +// Input* dgapxx_input = element->GetInput(HydrologyGapHeightXXEnum); /***/ +// Input* dgapyy_input = element->GetInput(HydrologyGapHeightYYEnum); /***/ /*Get conductivity from inputs*/ @@ -351,6 +346,6 @@ br_input->GetInputValue(&br,gauss); headold_input->GetInputValue(&head_old,gauss); - dgapxx_input->GetInputValue(&dgapxx,gauss); /***/ - dgapyy_input->GetInputValue(&dgapyy,gauss); /***/ +// dgapxx_input->GetInputValue(&dgapxx,gauss); /***/ +// dgapyy_input->GetInputValue(&dgapyy,gauss); /***/ /*Get ice A parameter*/ @@ -365,12 +360,4 @@ beta = 0.; - /* TEST Compute lc term */ -// if(gapGetAlpha2(&alpha2,gauss); @@ -378,5 +365,5 @@ frictionheat=alpha2*(vx*vx+vy*vy); - /* *** TEST take out frictional heat for large gap height *** */ + /* Take out frictional heat for large gap height */ if(gap>br) frictionheat=0.; @@ -399,26 +386,4 @@ _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]; - - /* 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]; - - /* With weighted creep term */ for(int i=0;ivalues[i]+=Jdet*gauss->weight* ( @@ -443,15 +408,5 @@ // )*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*/ @@ -597,6 +552,6 @@ IssmDouble q = 0.; IssmDouble channelization = 0.; - IssmDouble dgapxx; /***/ - IssmDouble dgapyy; /***/ +// IssmDouble dgapxx; /***/ +// IssmDouble dgapyy; /***/ /*Retrieve all inputs and parameters*/ @@ -616,6 +571,6 @@ Input* lr_input = element->GetInput(HydrologyBumpSpacingEnum); _assert_(lr_input); Input* br_input = element->GetInput(HydrologyBumpHeightEnum); _assert_(br_input); - Input* dgapxx_input = element->GetInput(HydrologyGapHeightXXEnum); /***/ - Input* dgapyy_input = element->GetInput(HydrologyGapHeightYYEnum); /***/ + // Input* dgapxx_input = element->GetInput(HydrologyGapHeightXXEnum); /***/ + // Input* dgapyy_input = element->GetInput(HydrologyGapHeightYYEnum); /***/ /*Get conductivity from inputs*/ @@ -643,6 +598,6 @@ lr_input->GetInputValue(&lr,gauss); br_input->GetInputValue(&br,gauss); - dgapxx_input->GetInputValue(&dgapxx,gauss); /***/ - dgapyy_input->GetInputValue(&dgapyy,gauss); /***/ +// dgapxx_input->GetInputValue(&dgapxx,gauss); /***/ +// dgapyy_input->GetInputValue(&dgapyy,gauss); /***/ /*Get ice A parameter*/ @@ -657,12 +612,4 @@ beta = 0.; - /* TEST Compute lc term*/ -// if(gapGetAlpha2(&alpha2,gauss); @@ -670,5 +617,5 @@ frictionheat=alpha2*(vx*vx+vy*vy); - /* *** TEST take out frictional heat for large gap height *** */ + /* Take out frictional heat for large gap height */ if(gap>br) frictionheat=0.; @@ -702,10 +649,4 @@ - /*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;