Index: /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp (revision 21435) +++ /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp (revision 21436) @@ -207,6 +207,7 @@ /* Intermediaries */ IssmDouble D_scalar,Jdet,dt; - IssmDouble epl_thickness; IssmDouble transfer; + IssmDouble epl_transmitivity; + IssmDouble epl_storing; IssmDouble *xyz_list = NULL; @@ -225,21 +226,18 @@ Input* epl_thick_input = basalelement->GetInput(HydrologydcEplThicknessEnum); _assert_(epl_thick_input); - Input* sed_head_input = basalelement->GetInput(SedimentHeadEnum); - Input* epl_head_input = basalelement->GetInput(EplHeadEnum); - Input* thick_input = basalelement->GetInput(ThicknessEnum); - Input* base_input = basalelement->GetInput(BaseEnum); - - IssmDouble epl_specificstoring = EplSpecificStoring(basalelement); - IssmDouble epl_conductivity = basalelement->GetMaterialParameter(HydrologydcEplConductivityEnum); + Input* epl_head_input = basalelement->GetInput(EplHeadEnum); _assert_(epl_head_input); + Input* base_input = basalelement->GetInput(BaseEnum); _assert_(base_input); /* Start looping on the number of gaussian points: */ - Gauss* gauss=basalelement->NewGauss(2); + Gauss* gauss = basalelement->NewGauss(2); for(int ig=gauss->begin();igend();ig++){ gauss ->GaussPoint(ig); basalelement ->JacobianDeterminant(&Jdet,xyz_list,gauss); - epl_thick_input ->GetInputValue(&epl_thickness,gauss); + + epl_transmitivity = EplTransmitivity(basalelement,gauss,epl_head_input,epl_thick_input,base_input); + epl_storing = EplStoring(basalelement,gauss,epl_head_input,epl_thick_input,base_input); /*Diffusivity*/ - D_scalar=epl_conductivity*epl_thickness*gauss->weight*Jdet; + D_scalar=epl_transmitivity*gauss->weight*Jdet; if(dt!=0.) D_scalar=D_scalar*dt; D[0][0]=D_scalar; @@ -254,5 +252,5 @@ if(dt!=0.){ basalelement->NodalFunctions(&basis[0],gauss); - D_scalar=epl_specificstoring*epl_thickness*gauss->weight*Jdet; + D_scalar=epl_storing*gauss->weight*Jdet; TripleMultiply(basis,numnodes,1,0, &D_scalar,1,1,0, @@ -262,5 +260,5 @@ /*Transfer EPL part*/ - transfer=GetHydrologyKMatrixTransfer(basalelement,gauss,sed_head_input,epl_head_input,thick_input,base_input); + transfer=GetHydrologyKMatrixTransfer(basalelement); D_scalar=dt*transfer*gauss->weight*Jdet; TripleMultiply(basis,numnodes,1,0, @@ -314,5 +312,5 @@ IssmDouble dt,scalar,water_head; IssmDouble water_load,transfer; - IssmDouble epl_thickness; + IssmDouble epl_storing; IssmDouble Jdet; IssmDouble residual,connectivity; @@ -334,16 +332,15 @@ Input* epl_thick_input = basalelement->GetInput(HydrologydcEplThicknessEnum); _assert_(epl_thick_input); - Input* sed_head_input = basalelement->GetInput(SedimentHeadEnum); - Input* epl_head_input = basalelement->GetInput(EplHeadEnum); - Input* thick_input = basalelement->GetInput(ThicknessEnum); - Input* base_input = basalelement->GetInput(BaseEnum); + Input* sed_head_input = basalelement->GetInput(SedimentHeadEnum); _assert_(sed_head_input); + Input* epl_head_input = basalelement->GetInput(EplHeadEnum); _assert_(epl_head_input); Input* water_input = basalelement->GetInput(BasalforcingsGroundediceMeltingRateEnum); _assert_(water_input); - Input* residual_input = basalelement->GetInput(SedimentHeadResidualEnum); _assert_(residual_input); - if(dt!= 0.){old_wh_input = basalelement->GetInput(EplHeadOldEnum); _assert_(old_wh_input);} - - IssmDouble epl_specificstoring = EplSpecificStoring(basalelement); - + Input* residual_input = basalelement->GetInput(SedimentHeadResidualEnum); _assert_(residual_input); + Input* base_input = basalelement->GetInput(BaseEnum); _assert_(base_input); + + if(dt!= 0.){ + old_wh_input = basalelement->GetInput(EplHeadOldEnum); _assert_(old_wh_input); + } /* Start looping on the number of gaussian points: */ - Gauss* gauss=basalelement->NewGauss(2); + Gauss* gauss = basalelement->NewGauss(2); for(int ig=gauss->begin();igend();ig++){ gauss->GaussPoint(ig); @@ -351,4 +348,5 @@ basalelement ->NodalFunctions(basis,gauss); + epl_storing = EplStoring(basalelement,gauss,epl_head_input,epl_thick_input,base_input); /*Loading term*/ water_input->GetInputValue(&water_load,gauss); @@ -362,9 +360,8 @@ if(dt!=0.){ old_wh_input ->GetInputValue(&water_head,gauss); - epl_thick_input ->GetInputValue(&epl_thickness,gauss); /*Dealing with the epl part of the transfer term*/ - transfer=GetHydrologyPVectorTransfer(basalelement,gauss,sed_head_input,epl_head_input,thick_input,base_input); - scalar = Jdet*gauss->weight*((water_head*epl_specificstoring*epl_thickness)+(dt*transfer)); + transfer=GetHydrologyPVectorTransfer(basalelement,gauss,sed_head_input); + scalar = Jdet*gauss->weight*((water_head*epl_storing)+(dt*transfer)); for(int i=0;ivalues[i]+=scalar*basis[i]; } @@ -458,24 +455,5 @@ /*Intermediaries*/ -IssmDouble HydrologyDCEfficientAnalysis::EplSpecificStoring(Element* element){/*{{{*/ - IssmDouble rho_freshwater = element->GetMaterialParameter(MaterialsRhoFreshwaterEnum); - IssmDouble g = element->GetMaterialParameter(ConstantsGEnum); - IssmDouble epl_porosity = element->GetMaterialParameter(HydrologydcEplPorosityEnum); - IssmDouble epl_compressibility = element->GetMaterialParameter(HydrologydcEplCompressibilityEnum); - IssmDouble water_compressibility = element->GetMaterialParameter(HydrologydcWaterCompressibilityEnum); - return rho_freshwater*g*epl_porosity*(water_compressibility+(epl_compressibility/epl_porosity)); -}/*}}}*/ - -IssmDouble HydrologyDCEfficientAnalysis::SedimentStoring(Element* element){/*{{{*/ - IssmDouble rho_freshwater = element->GetMaterialParameter(MaterialsRhoFreshwaterEnum); - IssmDouble g = element->GetMaterialParameter(ConstantsGEnum); - IssmDouble sediment_porosity = element->GetMaterialParameter(HydrologydcSedimentPorosityEnum); - IssmDouble sediment_thickness = element->GetMaterialParameter(HydrologydcSedimentThicknessEnum); - IssmDouble sediment_compressibility = element->GetMaterialParameter(HydrologydcSedimentCompressibilityEnum); - IssmDouble water_compressibility = element->GetMaterialParameter(HydrologydcWaterCompressibilityEnum); - return rho_freshwater*g*sediment_porosity*sediment_thickness*(water_compressibility+(sediment_compressibility/sediment_porosity)); -}/*}}}*/ - -IssmDouble HydrologyDCEfficientAnalysis::GetHydrologyKMatrixTransfer(Element* element, Gauss* gauss, Input* sed_head_input, Input* epl_head_input, Input* thick_input, Input* base_input){/*{{{*/ +IssmDouble HydrologyDCEfficientAnalysis::GetHydrologyKMatrixTransfer(Element* element){/*{{{*/ int transfermethod; @@ -483,5 +461,4 @@ element->FindParam(&transfermethod,HydrologydcTransferFlagEnum); - /*Switch between the different transfer methods cases*/ switch(transfermethod){ @@ -500,5 +477,5 @@ }/*}}}*/ -IssmDouble HydrologyDCEfficientAnalysis::GetHydrologyPVectorTransfer(Element* element, Gauss* gauss, Input* sed_head_input, Input* epl_head_input, Input* thick_input, Input* base_input){/*{{{*/ +IssmDouble HydrologyDCEfficientAnalysis::GetHydrologyPVectorTransfer(Element* element, Gauss* gauss, Input* sed_head_input){/*{{{*/ int transfermethod; @@ -507,5 +484,4 @@ element->FindParam(&transfermethod,HydrologydcTransferFlagEnum); - /*Switch between the different transfer methods cases*/ switch(transfermethod){ @@ -537,10 +513,12 @@ femmodel->parameters->FindParam(&domaintype,DomainTypeEnum); + femmodel->parameters->FindParam(&iseplthickcomp,HydrologydcEplThickCompEnum); + if(iseplthickcomp==0) return; for(int j=0;jelements->Size();j++){ Element* element=(Element*)femmodel->elements->GetObjectByOffset(j); - element->parameters->FindParam(&iseplthickcomp,HydrologydcEplThickCompEnum); - if(iseplthickcomp==0) return; + /* element->parameters->FindParam(&iseplthickcomp,HydrologydcEplThickCompEnum); */ + /* if(iseplthickcomp==0) return; */ switch(domaintype){ @@ -738,6 +716,52 @@ } /*}}}*/ - -void HydrologyDCEfficientAnalysis::HydrologyEPLGetActive(Vector* active_vec, Element* element){ +IssmDouble HydrologyDCEfficientAnalysis::EplStoring(Element* element,Gauss* gauss, Input* epl_head_input, Input* epl_thick_input, Input* base_input){/*{{{*/ + IssmDouble epl_storing; + IssmDouble storing,water_sheet; + IssmDouble base_elev,prestep_head,epl_thick; + IssmDouble rho_freshwater = element->GetMaterialParameter(MaterialsRhoFreshwaterEnum); + IssmDouble g = element->GetMaterialParameter(ConstantsGEnum); + IssmDouble epl_porosity = element->GetMaterialParameter(HydrologydcEplPorosityEnum); + IssmDouble epl_compressibility = element->GetMaterialParameter(HydrologydcEplCompressibilityEnum); + IssmDouble water_compressibility = element->GetMaterialParameter(HydrologydcWaterCompressibilityEnum); + + base_input->GetInputValue(&base_elev,gauss); + epl_head_input->GetInputValue(&prestep_head,gauss); + epl_thick_input->GetInputValue(&epl_thick,gauss); + + water_sheet=max(0.0,(prestep_head-base_elev)); + + storing=rho_freshwater*g*epl_porosity*epl_thick*(water_compressibility+(epl_compressibility/epl_porosity)); + //porosity for unconfined region + if (water_sheet<=0.9*epl_thick){ + epl_storing=epl_porosity; + } + //continuity ramp + else if((water_sheet0.9*epl_thick)){ + epl_storing=(epl_thick-water_sheet)*(epl_porosity-storing)/(0.1*epl_thick)+storing; + } + //storing coefficient for confined + else{ + epl_storing=storing; + } + return epl_storing; +}/*}}}*/ + +IssmDouble HydrologyDCEfficientAnalysis::EplTransmitivity(Element* element,Gauss* gauss, Input* epl_head_input, Input* epl_thick_input, Input* base_input){/*{{{*/ + IssmDouble epl_transmitivity; + IssmDouble base_elev,prestep_head,epl_thick; + IssmDouble water_sheet; + IssmDouble epl_conductivity = element->GetMaterialParameter(HydrologydcEplConductivityEnum); + base_input->GetInputValue(&base_elev,gauss); + epl_head_input->GetInputValue(&prestep_head,gauss); + epl_thick_input->GetInputValue(&epl_thick,gauss); + + water_sheet=max(0.0,(prestep_head-base_elev)); + epl_transmitivity=epl_conductivity*min(water_sheet,epl_thick); + + return epl_transmitivity; +}/*}}}*/ + +void HydrologyDCEfficientAnalysis::HydrologyEPLGetActive(Vector* active_vec, Element* element){/*{{{*/ /*Constants*/ int domaintype; @@ -787,4 +811,5 @@ xDelete(active); } +/*}}}*/ void HydrologyDCEfficientAnalysis::GetHydrologyDCInefficientHmax(IssmDouble* ph_max,Element* element, Node* innode){/*{{{*/ Index: /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.h (revision 21435) +++ /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.h (revision 21436) @@ -36,8 +36,8 @@ /*Intermediaries*/ void GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss); - IssmDouble EplSpecificStoring(Element* element); - IssmDouble SedimentStoring(Element* element); - IssmDouble GetHydrologyKMatrixTransfer(Element* element, Gauss* gauss, Input* sed_head_input, Input* epl_head_input, Input* thick_input, Input* base_input); - IssmDouble GetHydrologyPVectorTransfer(Element* element, Gauss* gauss, Input* sed_head_input, Input* epl_head_input, Input* thick_input, Input* base_input); + IssmDouble GetHydrologyKMatrixTransfer(Element* element); + IssmDouble GetHydrologyPVectorTransfer(Element* element, Gauss* gauss, Input* sed_head_input); + IssmDouble EplStoring(Element* element,Gauss* gauss, Input* epl_head_input, Input* epl_thick_input, Input* base_input); + IssmDouble EplTransmitivity(Element* element,Gauss* gauss, Input* epl_head_input, Input* epl_thick_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/analyses/HydrologyDCInefficientAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp (revision 21435) +++ /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp (revision 21436) @@ -191,5 +191,4 @@ IssmDouble D_scalar,Jdet,dt; IssmDouble sediment_transmitivity; - IssmDouble base_elev; IssmDouble transfer,sediment_storing; IssmDouble *xyz_list = NULL; @@ -213,6 +212,4 @@ Input* SedTrans_input = basalelement->GetInput(HydrologydcSedimentTransmitivityEnum); _assert_(SedTrans_input); Input* sed_head_input = basalelement->GetInput(SedimentHeadEnum); - Input* epl_head_input = basalelement->GetInput(EplHeadEnum); - Input* thick_input = basalelement->GetInput(ThicknessEnum); Input* base_input = basalelement->GetInput(BaseEnum); @@ -255,5 +252,5 @@ active_element_input->GetInputValue(&active_element); if(active_element){ - transfer=GetHydrologyKMatrixTransfer(basalelement,gauss,sed_head_input,epl_head_input,thick_input,base_input); + transfer=GetHydrologyKMatrixTransfer(basalelement); basalelement->NodalFunctions(&basis[0],gauss); D_scalar=dt*transfer*gauss->weight*Jdet; @@ -319,8 +316,9 @@ Input* sed_head_input = basalelement->GetInput(SedimentHeadEnum); Input* epl_head_input = basalelement->GetInput(EplHeadEnum); - Input* thick_input = basalelement->GetInput(ThicknessEnum); Input* base_input = basalelement->GetInput(BaseEnum); Input* water_input = basalelement->GetInput(BasalforcingsGroundediceMeltingRateEnum); _assert_(water_input); - if(dt!= 0.){old_wh_input = basalelement->GetInput(SedimentHeadOldEnum); _assert_(old_wh_input);} + if(dt!= 0.){ + old_wh_input = basalelement->GetInput(SedimentHeadOldEnum); _assert_(old_wh_input); + } /*Transfer related Inputs*/ if(isefficientlayer){ @@ -365,5 +363,5 @@ active_element_input->GetInputValue(&active_element); if(active_element){ - transfer=GetHydrologyPVectorTransfer(basalelement,gauss,sed_head_input,epl_head_input,thick_input,base_input); + transfer=GetHydrologyPVectorTransfer(basalelement,gauss,epl_head_input); } else{ @@ -514,8 +512,7 @@ }/*}}}*/ -/*Intermediaries*/ IssmDouble HydrologyDCInefficientAnalysis::SedimentStoring(Element* element,Gauss* gauss,Input* sed_head_input, Input* base_input){/*{{{*/ IssmDouble sediment_storing; - IssmDouble specific_porosity; + IssmDouble storing; IssmDouble base_elev,prestep_head,water_sheet; IssmDouble rho_freshwater = element->GetMaterialParameter(MaterialsRhoFreshwaterEnum); @@ -528,5 +525,5 @@ sed_head_input->GetInputValue(&prestep_head,gauss); water_sheet=max(0.0,(prestep_head-(base_elev-sediment_thickness))); - specific_porosity=sediment_porosity-rho_freshwater*g*sediment_porosity*sediment_thickness*(water_compressibility+(sediment_compressibility/sediment_porosity)); + storing=rho_freshwater*g*sediment_porosity*sediment_thickness*(water_compressibility+(sediment_compressibility/sediment_porosity)); //porosity for unconfined region if (water_sheet<=0.9*sediment_thickness){ @@ -535,10 +532,9 @@ //continuity ramp else if((water_sheet0.9*sediment_thickness)){ - sediment_storing=(sediment_thickness-water_sheet)*specific_porosity/(0.1*sediment_thickness)+ - rho_freshwater*g*sediment_porosity*sediment_thickness*(water_compressibility+(sediment_compressibility/sediment_porosity)); + sediment_storing=(sediment_thickness-water_sheet)*(sediment_porosity-storing)/(0.1*sediment_thickness)+storing; } //storing coefficient for confined else{ - sediment_storing=rho_freshwater*g*sediment_porosity*sediment_thickness*(water_compressibility+(sediment_compressibility/sediment_porosity)); + sediment_storing=storing; } return sediment_storing; @@ -561,39 +557,4 @@ } return sediment_transmitivity; -}/*}}}*/ - -IssmDouble HydrologyDCInefficientAnalysis::GetHydrologyDCInefficientHmax(Element* element, Gauss* gauss, Input* thick_input, Input* base_input){/*{{{*/ - int hmax_flag; - IssmDouble h_max; - IssmDouble rho_ice,rho_water; - IssmDouble thickness,bed; - /*Get the flag to the limitation method*/ - element->FindParam(&hmax_flag,HydrologydcSedimentlimitFlagEnum); - - /*Switch between the different cases*/ - switch(hmax_flag){ - case 0: - h_max=1.0e+10; - break; - case 1: - element->FindParam(&h_max,HydrologydcSedimentlimitEnum); - break; - case 2: - rho_water = element->GetMaterialParameter(MaterialsRhoFreshwaterEnum); - rho_ice = element->GetMaterialParameter(MaterialsRhoIceEnum); - _assert_(thick_input); - _assert_(base_input); - /*Compute max*/ - thick_input->GetInputValue(&thickness,gauss); - base_input->GetInputValue(&bed,gauss); - h_max=((rho_ice*thickness)/rho_water)+bed; - break; - case 3: - _error_("Using normal stress not supported yet"); - break; - default: - _error_("no case higher than 3 for SedimentlimitFlag"); - } - return h_max; }/*}}}*/ @@ -634,5 +595,5 @@ /*}}}*/ -IssmDouble HydrologyDCInefficientAnalysis::GetHydrologyKMatrixTransfer(Element* element, Gauss* gauss, Input* sed_head_input, Input* epl_head_input, Input* thick_input, Input* base_input){/*{{{*/ +IssmDouble HydrologyDCInefficientAnalysis::GetHydrologyKMatrixTransfer(Element* element){/*{{{*/ int transfermethod; @@ -656,5 +617,5 @@ }/*}}}*/ -IssmDouble HydrologyDCInefficientAnalysis::GetHydrologyPVectorTransfer(Element* element, Gauss* gauss, Input* sed_head_input, Input* epl_head_input, Input* thick_input, Input* base_input){/*{{{*/ +IssmDouble HydrologyDCInefficientAnalysis::GetHydrologyPVectorTransfer(Element* element, Gauss* gauss, Input* epl_head_input){/*{{{*/ int transfermethod; Index: /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.h (revision 21435) +++ /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.h (revision 21436) @@ -36,8 +36,7 @@ IssmDouble SedimentStoring(Element* element, Gauss* gauss, Input* sed_head_input, Input* base_input); IssmDouble SedimentTransmitivity(Element* element,Gauss* gauss,Input* sed_head_input, Input* base_input,Input* SedTrans_input); - IssmDouble GetHydrologyDCInefficientHmax(Element* element, Gauss* gauss, Input* thickness_input, Input* base_input); void GetHydrologyDCInefficientHmax(IssmDouble* ph_max,Element* element, Node* innode); - IssmDouble GetHydrologyKMatrixTransfer(Element* element, Gauss* gauss, Input* sed_head_input, Input* epl_head_input, Input* thick_input, Input* base_input); - IssmDouble GetHydrologyPVectorTransfer(Element* element, Gauss* gauss, Input* sed_head_input, Input* epl_head_input, Input* thick_input, Input* base_input); + IssmDouble GetHydrologyKMatrixTransfer(Element* element); + IssmDouble GetHydrologyPVectorTransfer(Element* element, Gauss* gauss, Input* epl_head_input); void ElementizeEplMask(FemModel* femmodel); };