Index: /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp (revision 17348) +++ /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp (revision 17349) @@ -55,7 +55,8 @@ iomodel->FetchDataToInput(elements,MeshElementonbedEnum); iomodel->FetchDataToInput(elements,MeshElementonsurfaceEnum); - iomodel->FetchDataToInput(elements,BasalforcingsMeltingRateEnum); iomodel->FetchDataToInput(elements,EplHeadEnum); + iomodel->FetchDataToInput(elements,SedimentHeadEnum); iomodel->FetchDataToInput(elements,HydrologydcEplInitialThicknessEnum); + iomodel->FetchDataToInput(elements,HydrologydcSedimentTransmitivityEnum); // iomodel->FetchDataToInput(elements,HydrologydcMaskEplactiveNodeEnum); @@ -140,4 +141,5 @@ IssmDouble D_scalar,Jdet,dt; IssmDouble epl_thickness; + IssmDouble transfer; IssmDouble *xyz_list = NULL; @@ -154,15 +156,14 @@ basalelement->GetVerticesCoordinates(&xyz_list); basalelement->FindParam(&dt,TimesteppingTimeStepEnum); - Input* thickness_input=basalelement->GetInput(HydrologydcEplThicknessEnum); _assert_(thickness_input); - IssmDouble epl_specificstoring = EplSpecificStoring(basalelement); - IssmDouble epl_conductivity = basalelement->GetMaterialParameter(HydrologydcEplConductivityEnum); + Input* thickness_input = basalelement->GetInput(HydrologydcEplThicknessEnum); _assert_(thickness_input); + IssmDouble epl_specificstoring = EplSpecificStoring(basalelement); + IssmDouble epl_conductivity = basalelement->GetMaterialParameter(HydrologydcEplConductivityEnum); /* Start looping on the number of gaussian points: */ Gauss* gauss=basalelement->NewGauss(2); for(int ig=gauss->begin();igend();ig++){ - gauss->GaussPoint(ig); - - basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss); - thickness_input->GetInputValue(&epl_thickness,gauss); + gauss ->GaussPoint(ig); + basalelement ->JacobianDeterminant(&Jdet,xyz_list,gauss); + thickness_input ->GetInputValue(&epl_thickness,gauss); /*Diffusivity*/ @@ -179,5 +180,5 @@ /*Transient*/ if(dt!=0.){ - basalelement->NodalFunctions(basis,gauss); + basalelement->NodalFunctions(&basis[0],gauss); D_scalar=epl_specificstoring*epl_thickness*gauss->weight*Jdet; TripleMultiply(basis,numnodes,1,0, @@ -185,7 +186,15 @@ basis,1,numnodes,0, &Ke->values[0],1); - } - } - + + /*Transfer EPL part*/ + transfer=GetHydrologyKMatrixTransfer(basalelement,gauss); + basalelement->NodalFunctions(&basis[0],gauss); + D_scalar=transfer*gauss->weight*Jdet*dt; + TripleMultiply(basis,numnodes,1,0, + &D_scalar,1,1,0, + basis,1,numnodes,0, + &Ke->values[0],1); + } + } /*Clean up and return*/ xDelete(xyz_list); @@ -224,7 +233,10 @@ } /*Intermediaries */ - IssmDouble dt,scalar,water_head,connectivity; - IssmDouble transfer,residual,epl_thickness; + IssmDouble dt,scalar,water_head; + IssmDouble transfer; + IssmDouble epl_thickness; IssmDouble Jdet; + IssmDouble residual,connectivity; + IssmDouble *xyz_list = NULL; Input* old_wh_input = NULL; @@ -240,10 +252,11 @@ /*Retrieve all inputs and parameters*/ basalelement->GetVerticesCoordinates(&xyz_list); + basalelement->FindParam(&dt,TimesteppingTimeStepEnum); + + Input* residual_input = basalelement->GetInput(SedimentHeadResidualEnum); _assert_(residual_input); + Input* thickness_input = basalelement->GetInput(HydrologydcEplThicknessEnum); _assert_(thickness_input); + if(dt!= 0.){old_wh_input = basalelement->GetInput(EplHeadOldEnum); _assert_(old_wh_input);} + IssmDouble epl_specificstoring = EplSpecificStoring(basalelement); - basalelement->FindParam(&dt,TimesteppingTimeStepEnum); - Input* residual_input = basalelement->GetInput(SedimentHeadResidualEnum); _assert_(residual_input); - Input* transfer_input = basalelement->GetInput(WaterTransferEnum); _assert_(transfer_input); - Input* thickness_input = basalelement->GetInput(HydrologydcEplThicknessEnum); _assert_(thickness_input); - if(dt!= 0.){old_wh_input = basalelement->GetInput(EplHeadOldEnum); _assert_(old_wh_input);} /* Start looping on the number of gaussian points: */ @@ -252,20 +265,15 @@ gauss->GaussPoint(ig); - basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss); - basalelement->NodalFunctions(basis,gauss); - - /*Loading term*/ - transfer_input->GetInputValue(&transfer,gauss); - scalar = Jdet*gauss->weight*(-transfer); - if(dt!=0.) scalar = scalar*dt; - - for(int i=0;ivalues[i]+=scalar*basis[i]; - - /*Transient term*/ + basalelement ->JacobianDeterminant(&Jdet,xyz_list,gauss); + basalelement ->NodalFunctions(basis,gauss); + + /*Transient and transfer terms*/ if(dt!=0.){ - thickness_input->GetInputValue(&epl_thickness,gauss); - old_wh_input->GetInputValue(&water_head,gauss); - scalar = Jdet*gauss->weight*water_head*epl_specificstoring*epl_thickness; + old_wh_input ->GetInputValue(&water_head,gauss); + thickness_input ->GetInputValue(&epl_thickness,gauss); + /*Dealing with the sediment part of the transfer term*/ + transfer=GetHydrologyPVectorTransfer(basalelement,gauss); + scalar = Jdet*gauss->weight*((water_head*epl_specificstoring*epl_thickness)+(transfer*dt)); for(int i=0;ivalues[i]+=scalar*basis[i]; } @@ -281,4 +289,5 @@ residual_input->GetInputValue(&residual,gauss); pe->values[iv]+=residual/connectivity; + } @@ -318,11 +327,8 @@ IssmDouble* eplHeads = xNew(numnodes); - /* IssmDouble* relaxed = xNew(numnodes); */ - /* IssmDouble* eplOldHeads = xNew(numnodes); */ - IssmDouble Stepping; - + IssmDouble* eplOldHeads = xNew(numnodes); /*Get previous water head*/ - //basalelement->GetInputListOnNodes(&eplOldHeads[0],EplHeadEnum); + basalelement->GetInputListOnNodes(&eplOldHeads[0],EplHeadEnum); /*Use the dof list to index into the solution vector: */ @@ -330,9 +336,6 @@ eplHeads[i]=solution[doflist[i]]; if(xIsNan(eplHeads[i])) _error_("NaN found in solution vector"); - } - - /* for(i=0;iAddBasalInput(EplHeadEnum,eplHeads,P1Enum); @@ -358,4 +361,135 @@ 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::GetHydrologyDCInefficientHmax(Element* element, Gauss* gauss){/*{{{*/ + 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); + element->GetInputValue(&thickness,gauss,ThicknessEnum); + element->GetInputValue(&bed,gauss,BedEnum); + 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; +}/*}}}*/ +IssmDouble HydrologyDCEfficientAnalysis::GetHydrologyKMatrixTransfer(Element* element, Gauss* gauss){/*{{{*/ + + 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); + + element->FindParam(&transfermethod,HydrologydcTransferFlagEnum); + /*Switch between the different transfer methods cases*/ + switch(transfermethod){ + case 0: + /*Just keepping the transfer to zero, should be OK with the initial value of transfer*/ + break; + case 1: + + element->GetInputValue(&epl_thickness,gauss,HydrologydcEplThicknessEnum); + element->GetInputValue(&sed_head,gauss,SedimentHeadEnum); + element->GetInputValue(&epl_head,gauss,EplHeadEnum); + element->GetInputValue(&sediment_transmitivity,gauss,HydrologydcSedimentTransmitivityEnum); + element->GetInputValue(&residual,gauss,SedimentHeadResidualEnum); + element->FindParam(&leakage,HydrologydcLeakageFactorEnum); + + if(epl_head>sed_head){ + if(residual>0.0){ + transfer=0.0; + } + else{ + transfer=(epl_specificstoring*epl_thickness*sediment_transmitivity)/(sediment_thickness*leakage); + } + } + else{ + transfer=(sediment_storing*sediment_transmitivity)/(sediment_thickness*leakage); + } + break; + default: + _error_("no case higher than 1 for the Transfer method"); + } + + return transfer; +}/*}}}*/ + +IssmDouble HydrologyDCEfficientAnalysis::GetHydrologyPVectorTransfer(Element* element, Gauss* gauss){/*{{{*/ + + int transfermethod; + IssmDouble epl_thickness; + 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); + + element->FindParam(&transfermethod,HydrologydcTransferFlagEnum); + /*Switch between the different transfer methods cases*/ + switch(transfermethod){ + case 0: + /*Just keepping the transfer to zero, should be OK with the initial value of transfer*/ + break; + case 1: + + element->GetInputValue(&epl_thickness,gauss,HydrologydcEplThicknessEnum); + element->GetInputValue(&sediment_head,gauss,SedimentHeadEnum); + element->GetInputValue(&epl_head,gauss,EplHeadEnum); + element->GetInputValue(&sediment_transmitivity,gauss,HydrologydcSedimentTransmitivityEnum); + element->GetInputValue(&residual,gauss,SedimentHeadResidualEnum); + element->FindParam(&leakage,HydrologydcLeakageFactorEnum); + + if(epl_head>sediment_head){ + if(residual>0.0){ + transfer=0.0; + } + else{ + transfer=(epl_specificstoring*epl_thickness*sediment_transmitivity*sediment_head)/(sediment_thickness*leakage); + } + } + else{ + transfer=(sediment_storing*sediment_transmitivity*sediment_head)/(sediment_thickness*leakage); + } + break; + default: + _error_("no case higher than 1 for the Transfer method"); + } + return transfer; +}/*}}}*/ + void HydrologyDCEfficientAnalysis::GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ /*Compute B matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2. Index: /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.h (revision 17348) +++ /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.h (revision 17349) @@ -33,4 +33,8 @@ void GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss); IssmDouble EplSpecificStoring(Element* element); + IssmDouble SedimentStoring(Element* element); + IssmDouble GetHydrologyDCInefficientHmax(Element* element, Gauss* gauss); + IssmDouble GetHydrologyKMatrixTransfer(Element* element, Gauss* gauss); + IssmDouble GetHydrologyPVectorTransfer(Element* element, Gauss* gauss); }; #endif Index: /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp (revision 17348) +++ /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp (revision 17349) @@ -89,4 +89,5 @@ iomodel->FetchDataToInput(elements,SedimentHeadEnum); iomodel->FetchDataToInput(elements,HydrologydcSedimentTransmitivityEnum); + if(isefficientlayer)iomodel->FetchDataToInput(elements,HydrologydcMaskEplactiveNodeEnum); @@ -172,7 +173,12 @@ /*Intermediaries */ + bool active_element,isefficientlayer; IssmDouble D_scalar,Jdet,dt; IssmDouble sediment_transmitivity; + IssmDouble transfer; IssmDouble *xyz_list = NULL; + + /*Define transfer related variables*/ + Input* active_element_input =NULL; /*Fetch number of nodes and dof for this finite element*/ @@ -186,9 +192,15 @@ /*Retrieve all inputs and parameters*/ - basalelement ->GetVerticesCoordinates(&xyz_list); - IssmDouble sediment_storing = SedimentStoring(basalelement); - Input* SedTrans_input=basalelement ->GetInput(HydrologydcSedimentTransmitivityEnum); _assert_(SedTrans_input); - basalelement ->FindParam(&dt,TimesteppingTimeStepEnum); - + basalelement ->GetVerticesCoordinates(&xyz_list); + basalelement ->FindParam(&dt,TimesteppingTimeStepEnum); + basalelement ->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum); + + Input* SedTrans_input = basalelement->GetInput(HydrologydcSedimentTransmitivityEnum); _assert_(SedTrans_input); + IssmDouble sediment_storing = SedimentStoring(basalelement); + /*Transfer related Inputs*/ + if(isefficientlayer){ + active_element_input = basalelement->GetInput(HydrologydcMaskEplactiveEltEnum); _assert_(active_element_input); + } + /* Start looping on the number of gaussian points: */ Gauss* gauss=basalelement->NewGauss(2); @@ -197,5 +209,5 @@ basalelement -> JacobianDeterminant(&Jdet,xyz_list,gauss); SedTrans_input -> GetInputValue(&sediment_transmitivity,gauss); - + /*Diffusivity*/ D_scalar=sediment_transmitivity*gauss->weight*Jdet; @@ -205,7 +217,7 @@ GetB(B,basalelement,xyz_list,gauss); TripleMultiply(B,2,numnodes,1, - &D[0][0],2,2,0, - B,2,numnodes,0, - &Ke->values[0],1); + &D[0][0],2,2,0, + B,2,numnodes,0, + &Ke->values[0],1); /*Transient*/ @@ -214,10 +226,24 @@ D_scalar=sediment_storing*gauss->weight*Jdet; TripleMultiply(basis,numnodes,1,0, - &D_scalar,1,1,0, - basis,1,numnodes,0, - &Ke->values[0],1); - } - } - + &D_scalar,1,1,0, + basis,1,numnodes,0, + &Ke->values[0],1); + + /*Transfer EPL part*/ + if(isefficientlayer){ + active_element_input->GetInputValue(&active_element); + if(active_element){ + transfer=GetHydrologyKMatrixTransfer(basalelement,gauss); + + basalelement->NodalFunctions(&basis[0],gauss); + D_scalar=transfer*gauss->weight*Jdet*dt; + TripleMultiply(basis,numnodes,1,0, + &D_scalar,1,1,0, + basis,1,numnodes,0, + &Ke->values[0],1); + } + } + } + } /*Clean up and return*/ xDelete(xyz_list); @@ -248,9 +274,12 @@ /*Intermediaries */ + bool active_element,isefficientlayer; IssmDouble dt,scalar,water_head; IssmDouble water_load,transfer; IssmDouble Jdet; - IssmDouble *xyz_list = NULL; - Input* old_wh_input = NULL; + + IssmDouble *xyz_list = NULL; + Input* old_wh_input = NULL; + Input* active_element_input = NULL; /*Fetch number of nodes and dof for this finite element*/ @@ -263,9 +292,16 @@ /*Retrieve all inputs and parameters*/ basalelement->GetVerticesCoordinates(&xyz_list); - IssmDouble sediment_storing = SedimentStoring(basalelement); basalelement->FindParam(&dt,TimesteppingTimeStepEnum); + basalelement->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum); + Input* water_input = basalelement->GetInput(BasalforcingsMeltingRateEnum); _assert_(water_input); - Input* transfer_input = basalelement->GetInput(WaterTransferEnum); _assert_(transfer_input); if(dt!= 0.){old_wh_input = basalelement->GetInput(SedimentHeadOldEnum); _assert_(old_wh_input);} + + IssmDouble sediment_storing = SedimentStoring(basalelement); + + /*Transfer related Inputs*/ + if(isefficientlayer){ + active_element_input = basalelement->GetInput(HydrologydcMaskEplactiveEltEnum); _assert_(active_element_input); + } /* Start looping on the number of gaussian points: */ @@ -273,5 +309,5 @@ for(int ig=gauss->begin();igend();ig++){ gauss->GaussPoint(ig); - + basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss); basalelement->NodalFunctions(basis,gauss); @@ -279,17 +315,32 @@ /*Loading term*/ water_input->GetInputValue(&water_load,gauss); - transfer_input->GetInputValue(&transfer,gauss); - scalar = Jdet*gauss->weight*(water_load+transfer); + + scalar = Jdet*gauss->weight*(water_load); if(dt!=0.) scalar = scalar*dt; - for(int i=0;ivalues[i]+=scalar*basis[i]; - - /*Transient term*/ + for(int i=0;ivalues[i]+=scalar*basis[i]; + } + + /*Transient and transfer terms*/ if(dt!=0.){ - old_wh_input->GetInputValue(&water_head,gauss); - scalar = Jdet*gauss->weight*water_head*sediment_storing; - for(int i=0;ivalues[i]+=scalar*basis[i]; - } - } - + old_wh_input ->GetInputValue(&water_head,gauss); + if(isefficientlayer){ + /*Dealing with the sediment part of the transfer term*/ + active_element_input->GetInputValue(&active_element); + if(active_element){ + transfer=GetHydrologyPVectorTransfer(basalelement,gauss); + } + else{ + transfer=0.0; + } + scalar = Jdet*gauss->weight*((water_head*sediment_storing)+(dt*transfer)); + for(int i=0;ivalues[i]+=scalar*basis[i]; + } + else{ + scalar = Jdet*gauss->weight*(water_head*sediment_storing); + for(int i=0;ivalues[i]+=scalar*basis[i]; + } + } + } /*Clean up and return*/ xDelete(xyz_list); @@ -371,5 +422,5 @@ for(int i=0;iGetHydrologyDCInefficientHmax(&h_max,i); + basalelement->GetHydrologyDCInefficientHmax(&h_max,basalelement->GetNode(i)); if(values[i]>h_max) residual[i] = kappa*(values[i]-h_max); else residual[i] = 0.; @@ -402,4 +453,143 @@ return rho_freshwater*g*sediment_porosity*sediment_thickness*(water_compressibility+(sediment_compressibility/sediment_porosity)); }/*}}}*/ + +IssmDouble HydrologyDCInefficientAnalysis::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 HydrologyDCInefficientAnalysis::GetHydrologyDCInefficientHmax(Element* element, Gauss* gauss){/*{{{*/ + 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); + + /*Compute max*/ + /* thick_input->GetInputValue(&thickness,gauss); */ + /* bed_input->GetInputValue(&bed,gauss); */ + element->GetInputValue(&thickness,gauss,ThicknessEnum); + element->GetInputValue(&bed,gauss,BedEnum); + 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; +}/*}}}*/ + +IssmDouble HydrologyDCInefficientAnalysis::GetHydrologyKMatrixTransfer(Element* element, Gauss* gauss){/*{{{*/ + + 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); + + element->FindParam(&transfermethod,HydrologydcTransferFlagEnum); + /*Switch between the different transfer methods cases*/ + switch(transfermethod){ + case 0: + /*Just keepping the transfer to zero, should be OK with the initial value of transfer*/ + break; + case 1: + + element->GetInputValue(&epl_thickness,gauss,HydrologydcEplThicknessEnum); + element->GetInputValue(&sed_head,gauss,SedimentHeadEnum); + element->GetInputValue(&epl_head,gauss,EplHeadEnum); + element->GetInputValue(&sediment_transmitivity,gauss,HydrologydcSedimentTransmitivityEnum); + + element->FindParam(&leakage,HydrologydcLeakageFactorEnum); + + if(epl_head>sed_head){ + h_max=GetHydrologyDCInefficientHmax(element,gauss); + if(sed_head>=h_max){ + transfer=0.0; + } + else{ + transfer=(epl_specificstoring*epl_thickness*sediment_transmitivity)/(sediment_thickness*leakage); + } + } + else{ + transfer=(sediment_storing*sediment_transmitivity)/(sediment_thickness*leakage); + } + break; + default: + _error_("no case higher than 1 for the Transfer method"); + } + + return transfer; +}/*}}}*/ + +IssmDouble HydrologyDCInefficientAnalysis::GetHydrologyPVectorTransfer(Element* element, Gauss* gauss){/*{{{*/ + + int transfermethod; + IssmDouble epl_thickness; + 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); + + element->FindParam(&transfermethod,HydrologydcTransferFlagEnum); + /*Switch between the different transfer methods cases*/ + switch(transfermethod){ + case 0: + /*Just keepping the transfer to zero, should be OK with the initial value of transfer*/ + break; + case 1: + + element->GetInputValue(&epl_thickness,gauss,HydrologydcEplThicknessEnum); + element->GetInputValue(&sediment_head,gauss,SedimentHeadEnum); + element->GetInputValue(&epl_head,gauss,EplHeadEnum); + element->GetInputValue(&sediment_transmitivity,gauss,HydrologydcSedimentTransmitivityEnum); + + element->FindParam(&leakage,HydrologydcLeakageFactorEnum); + if(epl_head>sediment_head){ + h_max=GetHydrologyDCInefficientHmax(element,gauss); + if(sediment_head>=h_max){ + transfer=0.0; + } + else{ + transfer=(epl_specificstoring*epl_thickness*sediment_transmitivity*epl_head)/(sediment_thickness*leakage); + } + } + else{ + transfer=(sediment_storing*sediment_transmitivity*epl_head)/(sediment_thickness*leakage); + } + break; + default: + _error_("no case higher than 1 for the Transfer method"); + } + return transfer; +}/*}}}*/ + void HydrologyDCInefficientAnalysis::ElementizeEplMask(FemModel* femmodel){/*{{{*/ Index: /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.h (revision 17348) +++ /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.h (revision 17349) @@ -33,4 +33,8 @@ void GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss); IssmDouble SedimentStoring(Element* element); + IssmDouble EplSpecificStoring(Element* element); + IssmDouble GetHydrologyDCInefficientHmax(Element* element, Gauss* gauss); + IssmDouble GetHydrologyKMatrixTransfer(Element* element, Gauss* gauss); + IssmDouble GetHydrologyPVectorTransfer(Element* element, Gauss* gauss); void ElementizeEplMask(FemModel* femmodel); }; Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 17348) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 17349) @@ -283,5 +283,5 @@ virtual void UpdateConstraintsExtrudeFromTop(void)=0; - virtual void GetHydrologyDCInefficientHmax(IssmDouble* ph_max,int index)=0; + // virtual void GetHydrologyDCInefficientHmax(IssmDouble* ph_max,int index)=0; virtual void GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode)=0; virtual void GetHydrologyTransfer(Vector* transfer)=0; Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 17348) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 17349) @@ -243,5 +243,5 @@ ElementMatrix* CreateEPLDomainMassMatrix(void); void GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode); - void GetHydrologyDCInefficientHmax(IssmDouble* ph_max,int index){_error_("not implemented yet");}; + // void GetHydrologyDCInefficientHmax(IssmDouble* ph_max,int index){_error_("not implemented yet");}; void GetHydrologyTransfer(Vector* transfer); void HydrologyEPLGetActive(Vector* active_vec); Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 17348) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 17349) @@ -151,5 +151,5 @@ void GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode){_error_("not implemented yet");}; - void GetHydrologyDCInefficientHmax(IssmDouble* ph_max,int index){_error_("not implemented yet");}; + // void GetHydrologyDCInefficientHmax(IssmDouble* ph_max,int index){_error_("not implemented yet");}; void GetHydrologyTransfer(Vector* transfer){_error_("not implemented yet");}; void HydrologyEPLGetActive(Vector* active_vec){_error_("not implemented yet");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 17348) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 17349) @@ -4465,39 +4465,4 @@ /*}}}*/ /*FUNCTION Tria::GetHydrologyDCInefficientHmax{{{*/ -void Tria::GetHydrologyDCInefficientHmax(IssmDouble* ph_max,int index){ - - int hmax_flag; - IssmDouble h_max; - IssmDouble rho_ice,rho_water; - IssmDouble thickness,bed; - /*Get the flag to the limitation method*/ - this->parameters->FindParam(&hmax_flag,HydrologydcSedimentlimitFlagEnum); - - /*Switch between the different cases*/ - switch(hmax_flag){ - case 0: - h_max=1.0e+10; - break; - case 1: - parameters->FindParam(&h_max,HydrologydcSedimentlimitEnum); - break; - case 2: - rho_ice=matpar->GetRhoIce(); - rho_water=matpar->GetRhoFreshwater(); - this->GetInputValue(&thickness,this->nodes[index],ThicknessEnum); - this->GetInputValue(&bed,this->nodes[index],BedEnum); - 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"); - } - /*Assign output pointer*/ - *ph_max=h_max; -} -/*}}}*/ -/*FUNCTION Tria::GetHydrologyDCInefficientHmax{{{*/ void Tria::GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode){ @@ -4535,91 +4500,95 @@ /*}}}*/ /*FUNCTION Tria::GetHydrologyTransfer{{{*/ -void Tria::GetHydrologyTransfer(Vector* transfer){ - - const int numdof = NDOF1 *NUMVERTICES; - int *doflist = NULL; - int analysis_type; - bool isefficientlayer; - bool active_element; - int transfermethod; - IssmDouble leakage,h_max; - IssmDouble wh_trans,sed_thick; - IssmDouble epl_specificstoring,sedstoring; - IssmDouble activeEpl[numdof],epl_thickness[numdof],old_epl_thickness[numdof]; - IssmDouble epl_head[numdof],sed_head[numdof]; - IssmDouble preceding_transfer[numdof],sed_trans[numdof]; - - Input* active_element_input=NULL; - - parameters->FindParam(&analysis_type,AnalysisTypeEnum); +/*This thing should be useless and deleted soon*/ + +/* void Tria::GetHydrologyTransfer(Vector* transfer){ */ + +/* const int numdof = NDOF1 *NUMVERTICES; */ +/* int *doflist = NULL; */ +/* int analysis_type; */ +/* bool isefficientlayer; */ +/* bool active_element; */ +/* int transfermethod; */ +/* IssmDouble leakage,h_max,dt,test; */ +/* IssmDouble wh_trans,sed_thick,relaxed; */ +/* IssmDouble epl_specificstoring,sedstoring; */ +/* IssmDouble activeEpl[numdof],epl_thickness[numdof],old_epl_thickness[numdof]; */ +/* IssmDouble epl_head[numdof],sed_head[numdof]; */ +/* IssmDouble epl_head_old[numdof],sed_head_old[numdof]; */ +/* IssmDouble preceding_transfer[numdof],sed_trans[numdof]; */ + +/* Input* active_element_input=NULL; */ + + +/* parameters->FindParam(&analysis_type,AnalysisTypeEnum); */ - GetDofList(&doflist,NoneApproximationEnum,GsetEnum); - - /*Get the flag to know if the efficient layer is present*/ - this->parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum); - - if(isefficientlayer){ - /*Also get the flag to the transfer method*/ - this->parameters->FindParam(&transfermethod,HydrologydcTransferFlagEnum); - - /*Switch between the different transfer methods cases*/ - switch(transfermethod){ - case 0: - /*Just keepping the transfer to zero, should be OK with the initial value of transfer*/ - break; - case 1: +/* GetDofList(&doflist,NoneApproximationEnum,GsetEnum); */ + +/* /\*Get the flag to know if the efficient layer is present*\/ */ +/* this->parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum); */ + +/* if(isefficientlayer){ */ +/* /\*Also get the flag to the transfer method*\/ */ +/* this->parameters->FindParam(&transfermethod,HydrologydcTransferFlagEnum); */ +/* this->parameters->FindParam(&dt,TimesteppingTimeStepEnum); */ +/* /\*Switch between the different transfer methods cases*\/ */ +/* switch(transfermethod){ */ +/* case 0: */ +/* /\*Just keepping the transfer to zero, should be OK with the initial value of transfer*\/ */ +/* break; */ +/* case 1: */ - active_element_input=inputs->GetInput(HydrologydcMaskEplactiveEltEnum); _assert_(active_element_input); - active_element_input->GetInputValue(&active_element); - - GetInputListOnVertices(&sed_head[0],SedimentHeadEnum); - GetInputListOnVertices(&sed_trans[0],HydrologydcSedimentTransmitivityEnum); - GetInputListOnVertices(&epl_head[0],EplHeadEnum); - GetInputListOnVertices(&epl_thickness[0],HydrologydcEplThicknessEnum); - - this->parameters->FindParam(&leakage,HydrologydcLeakageFactorEnum); - - sed_thick = matpar->GetSedimentThickness(); - - if(!active_element){ - /*No transfer if the EPL is not active*/ - } - else{ - GetInputListOnVertices(&preceding_transfer[0],WaterTransferEnum); - sedstoring=matpar->GetSedimentStoring(); - epl_specificstoring=matpar->GetEplSpecificStoring(); +/* active_element_input=inputs->GetInput(HydrologydcMaskEplactiveEltEnum); _assert_(active_element_input); */ +/* active_element_input->GetInputValue(&active_element); */ + +/* GetInputListOnVertices(&sed_head[0],SedimentHeadEnum); */ +/* GetInputListOnVertices(&sed_head_old[0],SedimentHeadOldEnum); */ +/* GetInputListOnVertices(&sed_trans[0],HydrologydcSedimentTransmitivityEnum); */ +/* GetInputListOnVertices(&epl_head[0],EplHeadEnum); */ +/* GetInputListOnVertices(&epl_head_old[0],EplHeadOldEnum); */ +/* GetInputListOnVertices(&epl_thickness[0],HydrologydcEplThicknessEnum); */ + +/* this->parameters->FindParam(&leakage,HydrologydcLeakageFactorEnum); */ + +/* sed_thick = matpar->GetSedimentThickness(); */ + +/* if(!active_element){ */ + + +/* /\*No transfer if the EPL is not active*\/ */ +/* } */ +/* else{ */ +/* GetInputListOnVertices(&preceding_transfer[0],WaterTransferEnum); */ +/* sedstoring=matpar->GetSedimentStoring(); */ +/* epl_specificstoring=matpar->GetEplSpecificStoring(); */ - for(int i=0;iGetHydrologyDCInefficientHmax(&h_max,i); - - /*EPL head higher than sediment head, transfer from the epl to the sediment*/ - if(epl_head[i]>sed_head[i]){ - wh_trans=epl_specificstoring*epl_thickness[i]*sed_trans[i]*(epl_head[i]-sed_head[i])/(leakage*sed_thick); - - /*No transfer if the sediment head is allready at the maximum*/ - if(sed_head[i]>=h_max){ - wh_trans=0.0; - } - } - /* EPL head lower than sediment head, transfer from the sediment to the epl */ - else if(epl_head[i]<=sed_head[i]){ - wh_trans=sedstoring*sed_trans[i]*(epl_head[i]-sed_head[i])/(leakage*sed_thick); - } +/* for(int i=0;iGetHydrologyDCInefficientHmax(&h_max,this->nodes[i]); */ +/* /\*EPL head higher than sediment head, transfer from the epl to the sediment*\/ */ +/* if(epl_head[i]>sed_head[i]){ */ +/* if(sed_head[i]>=h_max){ */ +/* wh_trans=0.0; */ +/* } */ +/* else{ */ +/* wh_trans=epl_specificstoring*epl_thickness[i]*sed_trans[i]*(epl_head[i]-sed_head[i])/(leakage*sed_thick); */ +/* } */ +/* } */ +/* /\* EPL head lower than sediment head, transfer from the sediment to the epl *\/ */ +/* else if(epl_head[i]<=sed_head[i]){ */ +/* wh_trans=sedstoring*sed_trans[i]*(epl_head[i]-sed_head[i])/(leakage*sed_thick); */ +/* } */ - /*Relaxation stuff*/ - wh_trans=preceding_transfer[i]+0.8*(wh_trans-preceding_transfer[i]); - - /*Assign output pointer*/ - transfer->SetValue(doflist[i],wh_trans,INS_VAL); - } - } - break; - default: - _error_("no case higher than 1 for the Transfer method"); - } - } - /*Free ressources:*/ - xDelete(doflist); -} +/* /\*Assign output pointer*\/ */ +/* transfer->SetValue(doflist[i],wh_trans,INS_VAL); */ +/* } */ +/* } */ +/* break; */ +/* default: */ +/* _error_("no case higher than 1 for the Transfer method"); */ +/* } */ +/* } */ +/* /\*Free ressources:*\/ */ +/* xDelete(doflist); */ +/* } */ /*}}}*/ /*FUNCTION Tria::HydrologyEPLGetActive {{{*/ @@ -4690,8 +4659,9 @@ if(epl_thickness[i]<0.001*init_thick){ vec_mask->SetValue(nodes[i]->Sid(),0.,INS_VAL); + epl_thickness[i]=init_thick; } } /*Increase of the efficient system is needed if the epl head reach the maximum value (sediment max value for now)*/ - this->GetHydrologyDCInefficientHmax(&h_max,i); + this->GetHydrologyDCInefficientHmax(&h_max,this->nodes[i]); if(eplhead[i]>=h_max && active_element){ for(j=0;j10.0*init_thick){ thickness[i] = 10.0*init_thick; } Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 17348) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 17349) @@ -254,5 +254,5 @@ ElementMatrix* CreateEPLDomainMassMatrix(void); void CreateHydrologyWaterVelocityInput(void); - void GetHydrologyDCInefficientHmax(IssmDouble* ph_max,int index); + // void GetHydrologyDCInefficientHmax(IssmDouble* ph_max,int index); void GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode); void GetHydrologyTransfer(Vector* transfer); Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp (revision 17348) +++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp (revision 17349) @@ -44,5 +44,5 @@ IssmDouble ndu_sed,nu_sed; IssmDouble ndu_epl,nu_epl; - IssmDouble SedConv,EplConv; + IssmDouble EplConv; /*Recover parameters: */ @@ -76,8 +76,8 @@ /*The real computation starts here, outermost loop is on the two layer system*/ for(;;){ + + EplConv=1.0; sedcount=1; eplcount=1; - SedConv=1.0; - EplConv=1.0; /*If there is two layers we need an outer loop value to compute convergence*/ @@ -97,7 +97,8 @@ /*Reset constraint on the ZigZag Lock, this thing doesn't work, it have to disapear*/ ResetConstraintsx(femmodel); - sedconverged=false; + /* {{{ *//*Treating the sediment*/ for(;;){ + sedconverged=false; uf_sed_sub_iter=uf_sed->Duplicate();_assert_(uf_sed_sub_iter); uf_sed->Copy(uf_sed_sub_iter); @@ -105,4 +106,5 @@ for(;;){ /* {{{ *//*Core of the computation*/ + if(VerboseSolution()) _printf0_("Building Sediment Matrix...\n"); SystemMatricesx(&Kff,&Kfs,&pf,&df,&sediment_kmax,femmodel); CreateNodalConstraintsx(&ys,femmodel->nodes,HydrologyDCInefficientAnalysisEnum); @@ -127,22 +129,15 @@ if(sedconverged)break; } + /* }}} *//*End of the sediment penalization loop*/ + /*Update EPL mask*/ + if(isefficientlayer){ + HydrologyDCInefficientAnalysis* analysis = new HydrologyDCInefficientAnalysis(); + analysis->ElementizeEplMask(femmodel); + delete analysis; + } + sedconverged=false; - /*Update Elemental Mask and compute transfer*/ - if(isefficientlayer){ - if(SedConv<0.3){ - HydrologyDCInefficientAnalysis* analysis = new HydrologyDCInefficientAnalysis(); - analysis->ElementizeEplMask(femmodel); - delete analysis; - femmodel->HydrologyTransferx(); - transfered=true; - } - else{ - transfered=false; - } - } - else{ - transfered=true; - } + /*Checking convegence on the value of the sediment head*/ duf=uf_sed_sub_iter->Duplicate();_assert_(duf); @@ -154,12 +149,11 @@ if (xIsNan(ndu_sed) || xIsNan(nu_sed)) _error_("convergence criterion is NaN!"); if (ndu_sed==0.0 && nu_sed==0.0) nu_sed=1.0e-6; /*Hacking the case where the layer is empty*/ - if(VerboseConvergence()) _printf0_(setw(50) << left << " Inner Sediment Convergence criterion:" << ndu_sed/nu_sed*100 << " aiming lower than " << eps_hyd*100 << " %\n"); - SedConv = ndu_sed/nu_sed*100; - if((ndu_sed/nu_sed)parameters->SetParam(sediment_kmax,HydrologySedimentKmaxEnum); InputUpdateFromConstantx(femmodel,sedconverged,ConvergedEnum); @@ -179,14 +173,12 @@ InputUpdateFromConstantx(femmodel,false,ConvergedEnum); femmodel->parameters->SetParam(HydrologyEfficientEnum,HydrologyLayerEnum); - - eplconverged=false; + EplConv=1.0; for(;;){ - + eplconverged=false; ug_epl_sub_iter=ug_epl->Duplicate();_assert_(ug_epl_sub_iter); ug_epl->Copy(ug_epl_sub_iter); - if(EplConv<0.3){ - /* {{{ *//*Retriev the EPL head slopes and compute EPL Thickness*/ + /* {{{ *//*Retriev the EPL head slopes and compute EPL Thickness*/ if(VerboseSolution()) _printf0_("computing EPL Head slope...\n"); femmodel->SetCurrentConfiguration(L2ProjectionEPLAnalysisEnum); @@ -201,17 +193,11 @@ //updating mask after the computation of the epl thickness (Allow to close too thin EPL) femmodel->HydrologyEPLupdateDomainx(); - /* }}} */ - - /*Update Elemental Mask and compute transfer*/ + HydrologyDCInefficientAnalysis* analysis = new HydrologyDCInefficientAnalysis(); analysis->ElementizeEplMask(femmodel); delete analysis; - femmodel->HydrologyTransferx(); - transfered=true; - } - else{ - transfered=false; - } - + /* }}} */ + + if(VerboseSolution()) _printf0_("Building EPL Matrix...\n"); SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel); CreateNodalConstraintsx(&ys,femmodel->nodes,HydrologyDCEfficientAnalysisEnum); @@ -236,7 +222,7 @@ if (xIsNan(ndu_epl) || xIsNan(nu_epl)) _error_("convergence criterion is NaN!"); if (ndu_epl==0.0 && nu_epl==0.0) nu_epl=1.0e-6; /*Hacking the case where the EPL is used but empty*/ - if(VerboseConvergence()) _printf0_(setw(50) << left << " Inner EPL Convergence criterion:" << ndu_epl/nu_epl*100 << " aiming lower than " << eps_hyd*100 << " %\n"); - EplConv=ndu_epl/nu_epl*100; - if((ndu_epl/nu_epl)=hydro_maxiter){ _error_(" maximum number of EPL iterations (" << hydro_maxiter << ") exceeded"); @@ -246,5 +232,5 @@ delete ug_epl_sub_iter; - if(eplconverged && transfered){ + if(eplconverged){ if(VerboseSolution()) _printf0_("eplconverged...\n"); InputUpdateFromConstantx(femmodel,eplconverged,ConvergedEnum); @@ -283,6 +269,6 @@ } else{ - if(VerboseConvergence()) _printf0_(setw(50) << left << " Sediment Convergence criterion:" << ndu_sed/nu_sed*100 << " aiming lower than " << eps_hyd*100 << " %\n"); - if(VerboseConvergence()) _printf0_(setw(50) << left << " EPL Convergence criterion:" << ndu_epl/nu_epl*100 << " aiming lower than " << eps_hyd*100 << " %\n"); + if(VerboseConvergence()) _printf0_(setw(50) << left << " Sediment Convergence criterion:" << ndu_sed/nu_sed*100 << "%, aiming lower than " << eps_hyd*100 << " %\n"); + if(VerboseConvergence()) _printf0_(setw(50) << left << " EPL Convergence criterion:" << ndu_epl/nu_epl*100 << "%, aiming lower than " << eps_hyd*100 << " %\n"); hydroconverged=false; }