Changeset 17349

Timestamp:

02/25/14 17:45:51 (11 years ago)

Author:

bdef

Message:

CHG:Now computing the transfer between layer in an implicit way

Location:

issm/trunk-jpl/src/c

Files:

• analyses/HydrologyDCEfficientAnalysis.cpp (modified) (12 diffs)
• analyses/HydrologyDCEfficientAnalysis.h (modified) (1 diff)
• analyses/HydrologyDCInefficientAnalysis.cpp (modified) (12 diffs)
• analyses/HydrologyDCInefficientAnalysis.h (modified) (1 diff)
• classes/Elements/Element.h (modified) (1 diff)
• classes/Elements/Penta.h (modified) (1 diff)
• classes/Elements/Seg.h (modified) (1 diff)
• classes/Elements/Tria.cpp (modified) (4 diffs)
• classes/Elements/Tria.h (modified) (1 diff)
• solutionsequences/solutionsequence_hydro_nonlinear.cpp (modified) (11 diffs)

issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp

@@ -55 +55 @@
         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 +141 @@
         IssmDouble  D_scalar,Jdet,dt;
         IssmDouble  epl_thickness;
+        IssmDouble  transfer;
         IssmDouble *xyz_list = NULL;
 
@@ -154 +156 @@
         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();ig<gauss->end();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 +180 @@
                 /*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 +186 @@
                                                 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<IssmDouble>(xyz_list);
@@ -224 +233 @@
         }
         /*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 +252 @@
         /*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 +265 @@
                 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;i<numnodes;i++)pe->values[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;i<numnodes;i++)pe->values[i]+=scalar*basis[i];
                 }
@@ -281 +289 @@
                 residual_input->GetInputValue(&residual,gauss);
                 pe->values[iv]+=residual/connectivity;
+
         }
 
@@ -318 +327 @@
         
         IssmDouble* eplHeads    = xNew<IssmDouble>(numnodes);
-        /* IssmDouble* relaxed    = xNew<IssmDouble>(numnodes); */
-        /* IssmDouble* eplOldHeads    = xNew<IssmDouble>(numnodes); */
-        IssmDouble  Stepping;
-
+        IssmDouble* eplOldHeads    = xNew<IssmDouble>(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 +336 @@
                 eplHeads[i]=solution[doflist[i]];
                 if(xIsNan<IssmDouble>(eplHeads[i])) _error_("NaN found in solution vector");
-        }
-        
-        /* for(i=0;i<numnodes;i++) { */
-        /*      relaxed[i] = eplOldHeads[i]+0.8*(eplHeads[i]-eplOldHeads[i]); */
-        /* } */
+
+        }
         /*Add input to the element: */
         element->AddBasalInput(EplHeadEnum,eplHeads,P1Enum);
@@ -358 +361 @@
         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. 

issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.h

@@ -33 +33 @@
                 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

issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp

@@ -89 +89 @@
         iomodel->FetchDataToInput(elements,SedimentHeadEnum);
         iomodel->FetchDataToInput(elements,HydrologydcSedimentTransmitivityEnum);
+
 
         if(isefficientlayer)iomodel->FetchDataToInput(elements,HydrologydcMaskEplactiveNodeEnum);
@@ -172 +173 @@
 
         /*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 +192 @@
 
         /*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 +209 @@
                 basalelement   -> JacobianDeterminant(&Jdet,xyz_list,gauss);
                 SedTrans_input -> GetInputValue(&sediment_transmitivity,gauss);
-
+                
                 /*Diffusivity*/
                 D_scalar=sediment_transmitivity*gauss->weight*Jdet;
@@ -205 +217 @@
                 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 +226 @@
                         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<IssmDouble>(xyz_list);
@@ -248 +274 @@
 
         /*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 +292 @@
         /*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 +309 @@
         for(int ig=gauss->begin();ig<gauss->end();ig++){
                 gauss->GaussPoint(ig);
-
+        
                 basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
                 basalelement->NodalFunctions(basis,gauss);
@@ -279 +315 @@
                 /*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;i<numnodes;i++) pe->values[i]+=scalar*basis[i];
-
-                /*Transient term*/
+                for(int i=0;i<numnodes;i++){
+                        pe->values[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;i<numnodes;i++) pe->values[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;i<numnodes;i++)pe->values[i]+=scalar*basis[i];
+                        }
+                        else{
+                                scalar = Jdet*gauss->weight*(water_head*sediment_storing);
+                                for(int i=0;i<numnodes;i++)pe->values[i]+=scalar*basis[i];
+                        }
+                }
+        }
         /*Clean up and return*/
         xDelete<IssmDouble>(xyz_list);
@@ -371 +422 @@
 
                 for(int i=0;i<numnodes;i++){
-                        basalelement->GetHydrologyDCInefficientHmax(&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 +453 @@
         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){/*{{{*/
 

issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.h

@@ -33 +33 @@
                 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);
 };

issm/trunk-jpl/src/c/classes/Elements/Element.h

@@ -283 +283 @@
                 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<IssmDouble>* transfer)=0; 

issm/trunk-jpl/src/c/classes/Elements/Penta.h

@@ -243 +243 @@
                 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<IssmDouble>* transfer);
                 void           HydrologyEPLGetActive(Vector<IssmDouble>* active_vec);

issm/trunk-jpl/src/c/classes/Elements/Seg.h

@@ -151 +151 @@
 
                 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<IssmDouble>* transfer){_error_("not implemented yet");};
                 void    HydrologyEPLGetActive(Vector<IssmDouble>* active_vec){_error_("not implemented yet");};

issm/trunk-jpl/src/c/classes/Elements/Tria.cpp

@@ -4465 +4465 @@
 /*}}}*/
 /*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 +4500 @@
 /*}}}*/
 /*FUNCTION Tria::GetHydrologyTransfer{{{*/
-void  Tria::GetHydrologyTransfer(Vector<IssmDouble>* 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<IssmDouble>* 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;i<numdof;i++){
-                                        this->GetHydrologyDCInefficientHmax(&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;i<numdof;i++){ */
+/*                                      this->GetHydrologyDCInefficientHmax(&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<int>(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<int>(doflist); */
+/* } */
 /*}}}*/
 /*FUNCTION Tria::HydrologyEPLGetActive {{{*/
@@ -4690 +4659 @@
                         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;j<numdof;j++){
@@ -4772 +4742 @@
                                 thickness[i] = old_thickness[i]*(1+((rho_water*gravity*dt)/(rho_ice*latentheat))*epl_conductivity*EPLgrad2-2.0*(A*dt/(pow(n,n)))*(pow(EPL_N,n)));
                                         
-                                /*Relaxation stuff*/
-                                if(thickness[i]<10.0*init_thick){
-                                        thickness[i] = preceding_thickness[i]+0.8*(thickness[i]-preceding_thickness[i]);
-                                }
-                                else{
+                                /*Take care of otherthikening*/
+                                if(thickness[i]>10.0*init_thick){
                                         thickness[i] = 10.0*init_thick;
                                 }

issm/trunk-jpl/src/c/classes/Elements/Tria.h

@@ -254 +254 @@
                 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<IssmDouble>* transfer);

issm/trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp

@@ -44 +44 @@
         IssmDouble ndu_sed,nu_sed;
         IssmDouble ndu_epl,nu_epl;
-        IssmDouble SedConv,EplConv;
+        IssmDouble EplConv;
 
         /*Recover parameters: */
@@ -76 +76 @@
         /*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 +97 @@
                 /*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 +106 @@
                         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 +129 @@
                                 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 +149 @@
                         if (xIsNan<IssmDouble>(ndu_sed) || xIsNan<IssmDouble>(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)<eps_hyd){
+                        if(VerboseConvergence()) _printf0_(setw(50) << left << "   Inner Sediment Convergence criterion:" << ndu_sed/nu_sed*100 << "%, aiming lower than " << eps_hyd*10*100 << " %\n");
+                        if((ndu_sed/nu_sed)<eps_hyd*10.){
                                 if(VerboseConvergence()) _printf0_("   # Inner sediment convergence achieve \n");
                                 sedconverged=true;
                         }
                         delete uf_sed_sub_iter;
-                        if(sedconverged && transfered){
+                        if(sedconverged){
                                 femmodel->parameters->SetParam(sediment_kmax,HydrologySedimentKmaxEnum);
                                 InputUpdateFromConstantx(femmodel,sedconverged,ConvergedEnum);
@@ -179 +173 @@
                         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 +193 @@
                                         //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 +222 @@
                                 if (xIsNan<IssmDouble>(ndu_epl) || xIsNan<IssmDouble>(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)<eps_hyd) eplconverged=true;
+                                if(VerboseConvergence()) _printf0_(setw(50) << left << "   Inner EPL Convergence criterion:" << ndu_epl/nu_epl*100 << "%, aiming lower than " << eps_hyd*10*100 << " %\n");
+                                EplConv=ndu_epl/nu_epl;
+                                if((ndu_epl/nu_epl)<eps_hyd*10.) eplconverged=true;
                                 if (eplcount>=hydro_maxiter){
                                         _error_("   maximum number of EPL iterations (" << hydro_maxiter << ") exceeded");
@@ -246 +232 @@
                                 
                                 delete ug_epl_sub_iter;
-                                if(eplconverged && transfered){
+                                if(eplconverged){
                                         if(VerboseSolution()) _printf0_("eplconverged...\n");
                                         InputUpdateFromConstantx(femmodel,eplconverged,ConvergedEnum);
@@ -283 +269 @@
                                 }
                                 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;
                                 }