Index: /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp (revision 21467) +++ /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp (revision 21468) @@ -235,6 +235,6 @@ basalelement ->JacobianDeterminant(&Jdet,xyz_list,gauss); - epl_transmitivity = EplTransmitivity(basalelement,gauss,epl_thick_input); - epl_storing = EplStoring(basalelement,gauss,epl_thick_input); + epl_transmitivity = EplTransmitivity(basalelement,gauss,epl_thick_input,epl_head_input,base_input); + epl_storing = EplStoring(basalelement,gauss,epl_thick_input,epl_head_input,base_input); /*Diffusivity*/ @@ -348,17 +348,14 @@ basalelement ->NodalFunctions(basis,gauss); - epl_storing = EplStoring(basalelement,gauss,epl_thick_input); + epl_storing = EplStoring(basalelement,gauss,epl_thick_input,epl_head_input,base_input); /*Loading term*/ water_input->GetInputValue(&water_load,gauss); scalar = Jdet*gauss->weight*(water_load); if(dt!=0.) scalar = scalar*dt; - for(int i=0;ivalues[i]+=scalar*basis[i]; - } + for(int i=0;ivalues[i]+=scalar*basis[i]; /*Transient and transfer terms*/ if(dt!=0.){ - old_wh_input ->GetInputValue(&water_head,gauss); - + old_wh_input->GetInputValue(&water_head,gauss); /*Dealing with the epl part of the transfer term*/ transfer=GetHydrologyPVectorTransfer(basalelement,gauss,sed_head_input); @@ -425,6 +422,5 @@ /*Use the dof list to index into the solution vector: */ - /*If the EPL is not active we revert to the initialisation vallue*/ - /* For now we keep OldValue but it is probably not the best*/ + /*If the EPL is not active we revert to the bedrock elevation*/ if(active_element){ for(int i=0;iGetInputListOnVertices(&eplHeads[0],EplHeadOldEnum); + basalelement->GetInputListOnVertices(&eplHeads[0],BaseEnum); for(int i=0;i(eplHeads[i])) _error_("NaN found in solution vector"); @@ -519,6 +515,4 @@ Element* element=(Element*)femmodel->elements->GetObjectByOffset(j); - /* element->parameters->FindParam(&iseplthickcomp,HydrologydcEplThickCompEnum); */ - /* if(iseplthickcomp==0) return; */ switch(domaintype){ @@ -576,5 +570,5 @@ /*Compute first the effective pressure in the EPL*/ EPL_N=gravity*((rho_ice*ice_thickness[i])-(rho_water*(eplhead[i]-bed[i]))); - if(EPL_N<0.0)EPL_N=0.0; + // if(EPL_N<0.0)EPL_N=0.0; /*Get then the square of the gradient of EPL heads*/ EPLgrad2 = (epl_slopeX[i]*epl_slopeX[i])+(epl_slopeY[i]*epl_slopeY[i]); @@ -688,14 +682,14 @@ else if(old_active[i]>0.){ vec_mask->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL); - /*If epl thickness gets under colapse thickness, close the layer*/ + /* If epl thickness gets under colapse thickness, close the layer */ if(epl_thickness[i]SetValue(basalelement->nodes[i]->Sid(),0.,INS_VAL); recurence->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL); } - /* /\*If epl head gets under base elevation, close the layer*\/ */ - /* else if(eplhead[i]SetValue(basalelement->nodes[i]->Sid(),0.,INS_VAL); */ - /* recurence->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL); */ - /* } */ + /*If epl head gets under base elevation, close the layer*/ + else if(eplhead[i]SetValue(basalelement->nodes[i]->Sid(),0.,INS_VAL); + recurence->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL); + } } /*Increase of the efficient system is needed if the epl head reach the maximum value (sediment max value for now)*/ @@ -724,7 +718,8 @@ } /*}}}*/ -IssmDouble HydrologyDCEfficientAnalysis::EplStoring(Element* element,Gauss* gauss, Input* epl_thick_input){/*{{{*/ +IssmDouble HydrologyDCEfficientAnalysis::EplStoring(Element* element,Gauss* gauss, Input* epl_thick_input, Input* epl_head_input, Input* base_input){/*{{{*/ IssmDouble epl_storing; - IssmDouble epl_thick; + IssmDouble water_sheet,storing; + IssmDouble epl_thickness,prestep_head,base_elev; IssmDouble rho_freshwater = element->GetMaterialParameter(MaterialsRhoFreshwaterEnum); IssmDouble g = element->GetMaterialParameter(ConstantsGEnum); @@ -733,18 +728,40 @@ IssmDouble water_compressibility = element->GetMaterialParameter(HydrologydcWaterCompressibilityEnum); - epl_thick_input->GetInputValue(&epl_thick,gauss); + epl_thick_input->GetInputValue(&epl_thickness,gauss); + epl_head_input->GetInputValue(&prestep_head,gauss); + base_input->GetInputValue(&base_elev,gauss); + water_sheet=max(0.0,(prestep_head-base_elev)); + + storing=rho_freshwater*g*epl_porosity*epl_thickness*(water_compressibility+(epl_compressibility/epl_porosity)); + + //porosity for unconfined region + if (water_sheet<=0.9*epl_thickness){ + epl_storing=epl_porosity; + } + //continuity ramp + else if((water_sheet0.9*epl_thickness)){ + epl_storing=(epl_thickness-water_sheet)*(epl_porosity-storing)/(0.1*epl_thickness)+storing; + } + //storing coefficient for confined + else{ + epl_storing=storing; + } + //return storing; + return epl_storing; +}/*}}}*/ + +IssmDouble HydrologyDCEfficientAnalysis::EplTransmitivity(Element* element,Gauss* gauss, Input* epl_thick_input, Input* epl_head_input, Input* base_input){/*{{{*/ + IssmDouble epl_transmitivity; + IssmDouble water_sheet; + IssmDouble epl_thickness,base_elev,prestep_head; + IssmDouble epl_conductivity = element->GetMaterialParameter(HydrologydcEplConductivityEnum); + epl_thick_input->GetInputValue(&epl_thickness,gauss); + epl_head_input->GetInputValue(&prestep_head,gauss); + base_input->GetInputValue(&base_elev,gauss); + + water_sheet=max(0.0,(prestep_head-base_elev)); - epl_storing=rho_freshwater*g*epl_porosity*epl_thick*(water_compressibility+(epl_compressibility/epl_porosity)); - return epl_storing; -}/*}}}*/ - -IssmDouble HydrologyDCEfficientAnalysis::EplTransmitivity(Element* element,Gauss* gauss, Input* epl_thick_input){/*{{{*/ - IssmDouble epl_transmitivity; - IssmDouble epl_thick; - IssmDouble epl_conductivity = element->GetMaterialParameter(HydrologydcEplConductivityEnum); - epl_thick_input->GetInputValue(&epl_thick,gauss); - - epl_transmitivity=epl_conductivity*epl_thick; - + //epl_transmitivity=epl_conductivity*epl_thickness; + epl_transmitivity=epl_conductivity*min(water_sheet,epl_thickness); return epl_transmitivity; }/*}}}*/ Index: /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.h (revision 21467) +++ /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.h (revision 21468) @@ -38,6 +38,6 @@ IssmDouble GetHydrologyKMatrixTransfer(Element* element); IssmDouble GetHydrologyPVectorTransfer(Element* element, Gauss* gauss, Input* sed_head_input); - IssmDouble EplStoring(Element* element,Gauss* gauss, Input* epl_thick_input); - IssmDouble EplTransmitivity(Element* element,Gauss* gauss, Input* epl_thick_input); + IssmDouble EplStoring(Element* element,Gauss* gauss, Input* epl_thick_input, Input* epl_head_input, Input* base_input); + IssmDouble EplTransmitivity(Element* element,Gauss* gauss, Input* epl_thick_input, Input* epl_head_input, Input* base_input); void HydrologyEPLGetMask(Vector* vec_mask, Vector* recurence, int* eplzigzag_counter, Element* element); void HydrologyEPLGetActive(Vector* active_vec, Element* element); Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp (revision 21467) +++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp (revision 21468) @@ -170,5 +170,4 @@ /*{{{*//*Retrieve the EPL head slopes and compute EPL Thickness*/ if(VerboseSolution()) _printf0_("computing EPL Head slope...\n"); - //inefanalysis->ElementizeEplMask(femmodel); femmodel->SetCurrentConfiguration(L2ProjectionEPLAnalysisEnum); femmodel->UpdateConstraintsL2ProjectionEPLx(&L2Count); @@ -182,6 +181,6 @@ effanalysis->ComputeEPLThickness(femmodel); //updating mask after the computation of the epl thickness (Allow to close too thin EPL) - femmodel->HydrologyEPLupdateDomainx(&ThickCount); - inefanalysis->ElementizeEplMask(femmodel); + /* femmodel->HydrologyEPLupdateDomainx(&ThickCount); */ + /* inefanalysis->ElementizeEplMask(femmodel); */ /*}}}*/