Index: /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp (revision 18632) +++ /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp (revision 18633) @@ -158,8 +158,8 @@ basalelement->GetVerticesCoordinates(&xyz_list); basalelement->FindParam(&dt,TimesteppingTimeStepEnum); - Input* thickness_input = basalelement->GetInput(HydrologydcEplThicknessEnum); _assert_(thickness_input); + + 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* sed_trans_input = basalelement->GetInput(HydrologydcSedimentTransmitivityEnum); Input* residual_input = basalelement->GetInput(SedimentHeadResidualEnum); @@ -172,5 +172,5 @@ gauss ->GaussPoint(ig); basalelement ->JacobianDeterminant(&Jdet,xyz_list,gauss); - thickness_input ->GetInputValue(&epl_thickness,gauss); + epl_thick_input ->GetInputValue(&epl_thickness,gauss); /*Diffusivity*/ @@ -195,5 +195,5 @@ /*Transfer EPL part*/ - transfer=GetHydrologyKMatrixTransfer(basalelement,gauss,thickness_input,sed_head_input,epl_head_input,sed_trans_input,residual_input); + transfer=GetHydrologyKMatrixTransfer(basalelement,gauss,sed_head_input,epl_head_input,residual_input); D_scalar=transfer*gauss->weight*Jdet*dt; TripleMultiply(basis,numnodes,1,0, @@ -266,9 +266,8 @@ basalelement->FindParam(&dt,TimesteppingTimeStepEnum); - Input* thickness_input = basalelement->GetInput(HydrologydcEplThicknessEnum); _assert_(thickness_input); - Input* sed_head_input = basalelement->GetInput(SedimentHeadEnum); - Input* epl_head_input = basalelement->GetInput(EplHeadEnum); - Input* sed_trans_input = basalelement->GetInput(HydrologydcSedimentTransmitivityEnum); - Input* residual_input = basalelement->GetInput(SedimentHeadResidualEnum); _assert_(residual_input); + 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* residual_input = basalelement->GetInput(SedimentHeadResidualEnum); _assert_(residual_input); if(dt!= 0.){old_wh_input = basalelement->GetInput(EplHeadOldEnum); _assert_(old_wh_input);} @@ -286,8 +285,8 @@ if(dt!=0.){ old_wh_input ->GetInputValue(&water_head,gauss); - thickness_input ->GetInputValue(&epl_thickness,gauss); + epl_thick_input ->GetInputValue(&epl_thickness,gauss); /*Dealing with the epl part of the transfer term*/ - transfer=GetHydrologyPVectorTransfer(basalelement,gauss,thickness_input,sed_head_input,epl_head_input,sed_trans_input,residual_input); + transfer=GetHydrologyPVectorTransfer(basalelement,gauss,sed_head_input,epl_head_input,residual_input); scalar = Jdet*gauss->weight*((water_head*epl_specificstoring*epl_thickness)+(transfer*dt)); for(int i=0;ivalues[i]+=scalar*basis[i]; @@ -322,4 +321,5 @@ void HydrologyDCEfficientAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/ + bool active_element; int domaintype,i; Element* basalelement=NULL; @@ -341,13 +341,28 @@ int numnodes = basalelement->GetNumberOfNodes(); + /*Fetch dof list and allocate solution vector*/ + IssmDouble* sedhead = xNew(numnodes); + IssmDouble* eplHeads = xNew(numnodes); + basalelement->GetDofList(&doflist,NoneApproximationEnum,GsetEnum); - IssmDouble* eplHeads = xNew(numnodes); - - /*Use the dof list to index into the solution vector: */ + basalelement->GetInputListOnVertices(&sedhead[0],SedimentHeadEnum); + + Input* active_element_input=basalelement->GetInput(HydrologydcMaskEplactiveEltEnum); _assert_(active_element_input); + active_element_input->GetInputValue(&active_element); + + /* if(!active_element){ */ + /* /\*Keeping thickness to initial value if EPL is not active*\/ */ + /* for(i=0;i(eplHeads[i])) _error_("NaN found in solution vector"); } + // } /*Add input to the element: */ element->AddBasalInput(EplHeadEnum,eplHeads,P1Enum); @@ -355,4 +370,5 @@ /*Free ressources:*/ xDelete(eplHeads); + xDelete(sedhead); xDelete(doflist); if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;}; @@ -381,17 +397,15 @@ return rho_freshwater*g*sediment_porosity*sediment_thickness*(water_compressibility+(sediment_compressibility/sediment_porosity)); }/*}}}*/ -IssmDouble HydrologyDCEfficientAnalysis::GetHydrologyKMatrixTransfer(Element* element, Gauss* gauss, Input* epl_thick_input, Input* sed_head_input, Input* epl_head_input, Input* sed_trans_input, Input* residual_input){/*{{{*/ + +IssmDouble HydrologyDCEfficientAnalysis::GetHydrologyKMatrixTransfer(Element* element, Gauss* gauss, Input* sed_head_input, Input* epl_head_input, Input* residual_input){/*{{{*/ + int transfermethod; - IssmDouble epl_thickness; - IssmDouble epl_head,sed_head; - IssmDouble sediment_transmitivity; - IssmDouble leakage,residual,transfer; - - IssmDouble sediment_thickness = element->GetMaterialParameter(HydrologydcSedimentThicknessEnum); - IssmDouble sediment_storing = SedimentStoring(element); - IssmDouble epl_specificstoring = EplSpecificStoring(element); + IssmDouble epl_thickness,sed_residual; + IssmDouble epl_head,sediment_head; + IssmDouble leakage,transfer; element->FindParam(&transfermethod,HydrologydcTransferFlagEnum); + /*Switch between the different transfer methods cases*/ switch(transfermethod){ @@ -401,37 +415,43 @@ break; case 1: - _assert_(epl_thick_input); - _assert_(sed_head_input); - _assert_(epl_head_input); - _assert_(sed_trans_input); - _assert_(residual_input); - /* get input */ - epl_thick_input->GetInputValue(&epl_thickness,gauss); - sed_head_input->GetInputValue(&sed_head,gauss); + _assert_(sed_head_input); + _assert_(epl_head_input); + _assert_(residual_input); + + sed_head_input->GetInputValue(&sediment_head,gauss); epl_head_input->GetInputValue(&epl_head,gauss); - sed_trans_input->GetInputValue(&sediment_transmitivity,gauss); - residual_input->GetInputValue(&residual,gauss); + residual_input->GetInputValue(&sed_residual,gauss); element->FindParam(&leakage,HydrologydcLeakageFactorEnum); - transfer=(sediment_transmitivity)/(sediment_thickness*leakage); + + if(epl_head>sediment_head){ + if(sed_residual>0.0){ + transfer=0.0; + } + else{ + transfer=(leakage); + } + } + else{ + transfer=(leakage); + } break; default: _error_("no case higher than 1 for the Transfer method"); } - + if(element->Id()==42){ + printf("Transferefficient Kmat %e, %e, %e, %i \n",transfer,sediment_head,epl_head,element->Id()); + } return transfer; }/*}}}*/ -IssmDouble HydrologyDCEfficientAnalysis::GetHydrologyPVectorTransfer(Element* element, Gauss* gauss, Input* epl_thick_input, Input* sed_head_input, Input* epl_head_input, Input* sed_trans_input, Input* residual_input){/*{{{*/ + +IssmDouble HydrologyDCEfficientAnalysis::GetHydrologyPVectorTransfer(Element* element, Gauss* gauss, Input* sed_head_input, Input* epl_head_input, Input* residual_input){/*{{{*/ int transfermethod; - IssmDouble epl_thickness; + IssmDouble sed_residual; IssmDouble epl_head,sediment_head; - IssmDouble sediment_transmitivity; - IssmDouble leakage,residual,transfer; - - IssmDouble sediment_thickness = element->GetMaterialParameter(HydrologydcSedimentThicknessEnum); - IssmDouble sediment_storing = SedimentStoring(element); - IssmDouble epl_specificstoring = EplSpecificStoring(element); + IssmDouble leakage,transfer; element->FindParam(&transfermethod,HydrologydcTransferFlagEnum); + /*Switch between the different transfer methods cases*/ switch(transfermethod){ @@ -441,21 +461,30 @@ break; case 1: - _assert_(epl_thick_input); - _assert_(sed_head_input); - _assert_(epl_head_input); - _assert_(sed_trans_input); - _assert_(residual_input); - /* get input */ - epl_thick_input->GetInputValue(&epl_thickness,gauss); + _assert_(sed_head_input); + _assert_(epl_head_input); + _assert_(residual_input); + sed_head_input->GetInputValue(&sediment_head,gauss); epl_head_input->GetInputValue(&epl_head,gauss); - sed_trans_input->GetInputValue(&sediment_transmitivity,gauss); - residual_input->GetInputValue(&residual,gauss); + residual_input->GetInputValue(&sed_residual,gauss); element->FindParam(&leakage,HydrologydcLeakageFactorEnum); - transfer=(sediment_transmitivity*sediment_head)/(sediment_thickness*leakage); + if(epl_head>sediment_head){ + if(sed_residual>0.0){ + transfer=0.0; + } + else{ + transfer=(sediment_head*leakage); + } + } + else{ + transfer=(sediment_head*leakage); + } break; default: _error_("no case higher than 1 for the Transfer method"); + } + if(element->Id()==42){ + printf("Transferefficient Pvec %e, %e, %e\n",transfer,sediment_head,epl_head); } return transfer; @@ -519,5 +548,5 @@ element->GetInputListOnVertices(&ice_thickness[0],ThicknessEnum); element->GetInputListOnVertices(&bed[0],BaseEnum); - + if(!active_element){ Index: /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.h (revision 18632) +++ /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.h (revision 18633) @@ -36,6 +36,6 @@ IssmDouble EplSpecificStoring(Element* element); IssmDouble SedimentStoring(Element* element); - IssmDouble GetHydrologyKMatrixTransfer(Element* element, Gauss* gauss, Input* epl_thick_input, Input* sed_head_input, Input* epl_head_input, Input* sed_trans_input, Input* residual_input); - IssmDouble GetHydrologyPVectorTransfer(Element* element, Gauss* gauss, Input* epl_thick_input, Input* sed_head_input, Input* epl_head_input, Input* sed_trans_input, Input* residual_input); + IssmDouble GetHydrologyKMatrixTransfer(Element* element, Gauss* gauss, Input* sed_head_input, Input* epl_head_input, Input* residual_input); + IssmDouble GetHydrologyPVectorTransfer(Element* element, Gauss* gauss, Input* sed_head_input, Input* epl_head_input, Input* residual_input); void HydrologyEPLGetMask(Vector* vec_mask,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 18632) +++ /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp (revision 18633) @@ -197,11 +197,9 @@ basalelement ->FindParam(&dt,TimesteppingTimeStepEnum); basalelement ->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum); - - Input* epl_thick_input = basalelement->GetInput(HydrologydcEplThicknessEnum); + Input* SedTrans_input = basalelement->GetInput(HydrologydcSedimentTransmitivityEnum); _assert_(SedTrans_input); Input* sed_head_input = basalelement->GetInput(SedimentHeadEnum); Input* epl_head_input = basalelement->GetInput(EplHeadEnum); - Input* thickness_input = basalelement->GetInput(ThicknessEnum); - Input* base_input = basalelement->GetInput(BaseEnum); - Input* SedTrans_input = basalelement->GetInput(HydrologydcSedimentTransmitivityEnum); _assert_(SedTrans_input); + Input* residual_input = basalelement->GetInput(SedimentHeadResidualEnum); + IssmDouble sediment_storing = SedimentStoring(basalelement); /*Transfer related Inputs*/ @@ -240,5 +238,5 @@ active_element_input->GetInputValue(&active_element); if(active_element){ - transfer=GetHydrologyKMatrixTransfer(basalelement,gauss,epl_thick_input,sed_head_input,epl_head_input,SedTrans_input,thickness_input,base_input); + transfer=GetHydrologyKMatrixTransfer(basalelement,gauss,sed_head_input,epl_head_input,residual_input); basalelement->NodalFunctions(&basis[0],gauss); D_scalar=transfer*gauss->weight*Jdet*dt; @@ -296,4 +294,5 @@ IssmDouble* basis = xNew(numnodes); + /*Retrieve all inputs and parameters*/ basalelement->GetVerticesCoordinates(&xyz_list); @@ -301,11 +300,8 @@ basalelement->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum); - Input* epl_thick_input = basalelement->GetInput(HydrologydcEplThicknessEnum); Input* sed_head_input = basalelement->GetInput(SedimentHeadEnum); Input* epl_head_input = basalelement->GetInput(EplHeadEnum); - Input* sed_trans_input = basalelement->GetInput(HydrologydcSedimentTransmitivityEnum); - Input* thickness_input = basalelement->GetInput(ThicknessEnum); - Input* base_input = basalelement->GetInput(BaseEnum); Input* water_input = basalelement->GetInput(BasalforcingsGroundediceMeltingRateEnum); _assert_(water_input); + Input* residual_input = basalelement->GetInput(SedimentHeadResidualEnum); if(dt!= 0.){old_wh_input = basalelement->GetInput(SedimentHeadOldEnum); _assert_(old_wh_input);} @@ -342,5 +338,5 @@ active_element_input->GetInputValue(&active_element); if(active_element){ - transfer=GetHydrologyPVectorTransfer(basalelement,gauss,epl_thick_input,sed_head_input,epl_head_input,sed_trans_input,thickness_input,base_input); + transfer=GetHydrologyPVectorTransfer(basalelement,gauss,sed_head_input,epl_head_input,residual_input); } else{ @@ -571,17 +567,14 @@ } /*}}}*/ -IssmDouble HydrologyDCInefficientAnalysis::GetHydrologyKMatrixTransfer(Element* element, Gauss* gauss, Input* epl_thick_input, Input* sed_head_input, Input* epl_head_input, Input* sed_trans_input, Input* thickness_input, Input* base_input){/*{{{*/ + +IssmDouble HydrologyDCInefficientAnalysis::GetHydrologyKMatrixTransfer(Element* element, Gauss* gauss, Input* sed_head_input, Input* epl_head_input, Input* residual_input){/*{{{*/ int transfermethod; - IssmDouble epl_thickness; - IssmDouble epl_head,sed_head; - IssmDouble sediment_transmitivity; - IssmDouble leakage,h_max,transfer; - - IssmDouble sediment_thickness = element->GetMaterialParameter(HydrologydcSedimentThicknessEnum); - IssmDouble sediment_storing = SedimentStoring(element); - IssmDouble epl_specificstoring = EplSpecificStoring(element); + IssmDouble sed_residual; + IssmDouble epl_head,sediment_head; + IssmDouble leakage,transfer; element->FindParam(&transfermethod,HydrologydcTransferFlagEnum); + /*Switch between the different transfer methods cases*/ switch(transfermethod){ @@ -591,40 +584,43 @@ break; case 1: - - _assert_(epl_thick_input); _assert_(sed_head_input); _assert_(epl_head_input); - _assert_(sed_trans_input); - _assert_(thickness_input); - _assert_(base_input); - - epl_thick_input->GetInputValue(&epl_thickness,gauss); - sed_head_input->GetInputValue(&sed_head,gauss); + _assert_(residual_input); + + sed_head_input->GetInputValue(&sediment_head,gauss); epl_head_input->GetInputValue(&epl_head,gauss); - sed_trans_input->GetInputValue(&sediment_transmitivity,gauss); + residual_input->GetInputValue(&sed_residual,gauss); element->FindParam(&leakage,HydrologydcLeakageFactorEnum); - transfer=(sediment_transmitivity)/(sediment_thickness*leakage); + if(epl_head>sediment_head){ + if(sediment_head>=400.0){ + transfer=0.0; + } + else{ + transfer=(leakage); + } + } + else{ + transfer=(leakage); + } break; default: _error_("no case higher than 1 for the Transfer method"); } - + if(element->Id()==42){ + printf("TransferInfficient Kmat %e, %e, %e, %i\n",transfer,sediment_head,epl_head,element->Id()); + } return transfer; }/*}}}*/ -IssmDouble HydrologyDCInefficientAnalysis::GetHydrologyPVectorTransfer(Element* element, Gauss* gauss, Input* epl_thick_input, Input* sed_head_input, Input* epl_head_input, Input* sed_trans_input, Input* thickness_input, Input* base_input){/*{{{*/ +IssmDouble HydrologyDCInefficientAnalysis::GetHydrologyPVectorTransfer(Element* element, Gauss* gauss, Input* sed_head_input, Input* epl_head_input, Input* residual_input){/*{{{*/ int transfermethod; - IssmDouble epl_thickness; + IssmDouble sed_residual; IssmDouble epl_head,sediment_head; - IssmDouble sediment_transmitivity; - IssmDouble leakage,h_max,transfer; - - IssmDouble sediment_thickness = element->GetMaterialParameter(HydrologydcSedimentThicknessEnum); - IssmDouble sediment_storing = SedimentStoring(element); - IssmDouble epl_specificstoring = EplSpecificStoring(element); + IssmDouble leakage,transfer; element->FindParam(&transfermethod,HydrologydcTransferFlagEnum); + /*Switch between the different transfer methods cases*/ switch(transfermethod){ @@ -634,22 +630,30 @@ break; case 1: - - _assert_(epl_thick_input); _assert_(sed_head_input); _assert_(epl_head_input); - _assert_(sed_trans_input); - _assert_(thickness_input); - _assert_(base_input); - - epl_thick_input->GetInputValue(&epl_thickness,gauss); + _assert_(residual_input); + sed_head_input->GetInputValue(&sediment_head,gauss); epl_head_input->GetInputValue(&epl_head,gauss); - sed_trans_input->GetInputValue(&sediment_transmitivity,gauss); + residual_input->GetInputValue(&sed_residual,gauss); element->FindParam(&leakage,HydrologydcLeakageFactorEnum); - transfer=(sediment_transmitivity*epl_head)/(sediment_thickness*leakage); + if(epl_head>sediment_head){ + if(sediment_head>=400.0){ + transfer=0.0; + } + else{ + transfer=(epl_head*leakage); + } + } + else{ + transfer=(epl_head*leakage); + } break; default: _error_("no case higher than 1 for the Transfer method"); + } + if(element->Id()==42){ + printf("TransferInfficient Pvec %e, %e, %e, %i\n",transfer,sediment_head,epl_head,element->Id()); } return transfer; Index: /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.h (revision 18632) +++ /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.h (revision 18633) @@ -38,6 +38,6 @@ 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* epl_thick_input, Input* sed_head_input, Input* epl_head_input, Input* sed_trans_input, Input* thickness_input, Input* base_input); - IssmDouble GetHydrologyPVectorTransfer(Element* element, Gauss* gauss, Input* epl_thick_input, Input* sed_head_input, Input* epl_head_input, Input* sed_trans_input, Input* thickness_input, Input* base_input); + IssmDouble GetHydrologyKMatrixTransfer(Element* element, Gauss* gauss, Input* sed_head_input, Input* epl_head_input, Input* residual_input); + IssmDouble GetHydrologyPVectorTransfer(Element* element, Gauss* gauss, Input* sed_head_input, Input* epl_head_input, Input* residual_input); void ElementizeEplMask(FemModel* femmodel); };