Changeset 21438

Timestamp:

12/12/16 01:17:08 (8 years ago)

Author:

bdef

Message:

CHG:reverting to confined for EPL but with colapse if head gets under base elevation

Location:

issm/trunk-jpl/src/c/analyses

Files:

• HydrologyDCEfficientAnalysis.cpp (modified) (7 diffs)
• HydrologyDCEfficientAnalysis.h (modified) (1 diff)

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

@@ -235 +235 @@
                 basalelement    ->JacobianDeterminant(&Jdet,xyz_list,gauss);
 
-                epl_transmitivity = EplTransmitivity(basalelement,gauss,epl_head_input,epl_thick_input,base_input);
-                epl_storing                             = EplStoring(basalelement,gauss,epl_head_input,epl_thick_input,base_input);
+                epl_transmitivity = EplTransmitivity(basalelement,gauss,epl_thick_input);
+                epl_storing                             = EplStoring(basalelement,gauss,epl_thick_input);
 
                 /*Diffusivity*/
@@ -348 +348 @@
                 basalelement ->NodalFunctions(basis,gauss);
 
-                epl_storing     = EplStoring(basalelement,gauss,epl_head_input,epl_thick_input,base_input);
+                epl_storing     = EplStoring(basalelement,gauss,epl_thick_input);
                 /*Loading term*/
                 water_input->GetInputValue(&water_load,gauss);
@@ -655 +655 @@
         IssmDouble* eplhead       =xNew<IssmDouble>(numnodes);
         IssmDouble* residual      =xNew<IssmDouble>(numnodes);
+        IssmDouble* base          =xNew<IssmDouble>(numnodes);
         
         IssmDouble init_thick    =basalelement->GetMaterialParameter(HydrologydcEplInitialThicknessEnum);
@@ -667 +668 @@
         basalelement-> GetInputListOnVertices(&eplhead[0],EplHeadEnum);
         basalelement-> GetInputListOnVertices(&residual[0],SedimentHeadResidualEnum);
+        basalelement-> GetInputListOnVertices(&base[0],BaseEnum);
 
         /*Get minimum sediment head of the element*/
@@ -688 +690 @@
                         /*If epl thickness gets under colapse thickness, close the layer*/
                         if(epl_thickness[i]<colapse_thick){
+                                vec_mask->SetValue(basalelement->nodes[i]->Sid(),0.,INS_VAL);
+                                recurence->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL);
+                        }
+                        /*If epl head gets under base elevation, close the layer*/
+                        else if(eplhead[i]<base[i]){
                                 vec_mask->SetValue(basalelement->nodes[i]->Sid(),0.,INS_VAL);
                                 recurence->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL);
@@ -714 +721 @@
         xDelete<IssmDouble>(eplhead);
         xDelete<IssmDouble>(residual);
+        xDelete<IssmDouble>(base);
 }
 /*}}}*/
-IssmDouble HydrologyDCEfficientAnalysis::EplStoring(Element* element,Gauss* gauss, Input* epl_head_input, Input* epl_thick_input, Input* base_input){/*{{{*/
+IssmDouble HydrologyDCEfficientAnalysis::EplStoring(Element* element,Gauss* gauss, Input* epl_thick_input){/*{{{*/
         IssmDouble epl_storing;
-        IssmDouble storing,water_sheet;
-        IssmDouble base_elev,prestep_head,epl_thick;
+        IssmDouble epl_thick;
         IssmDouble rho_freshwater        = element->GetMaterialParameter(MaterialsRhoFreshwaterEnum);
         IssmDouble g                     = element->GetMaterialParameter(ConstantsGEnum);
@@ -726 +733 @@
         IssmDouble water_compressibility = element->GetMaterialParameter(HydrologydcWaterCompressibilityEnum);
 
-        base_input->GetInputValue(&base_elev,gauss);
-        epl_head_input->GetInputValue(&prestep_head,gauss);
         epl_thick_input->GetInputValue(&epl_thick,gauss);
         
-        water_sheet=max(0.0,(prestep_head-base_elev));
-        
-        storing=rho_freshwater*g*epl_porosity*epl_thick*(water_compressibility+(epl_compressibility/epl_porosity));
-        //porosity for unconfined region
-        if (water_sheet<=0.9*epl_thick){
-                epl_storing=epl_porosity;
-        }
-        //continuity ramp
-        else if((water_sheet<epl_thick) && (water_sheet>0.9*epl_thick)){
-                epl_storing=(epl_thick-water_sheet)*(epl_porosity-storing)/(0.1*epl_thick)+storing;
-        }
-        //storing coefficient for confined
-        else{
-                epl_storing=storing;
-        }
+        epl_storing=rho_freshwater*g*epl_porosity*epl_thick*(water_compressibility+(epl_compressibility/epl_porosity));
         return epl_storing;
 }/*}}}*/
 
-IssmDouble HydrologyDCEfficientAnalysis::EplTransmitivity(Element* element,Gauss* gauss, Input* epl_head_input, Input* epl_thick_input, Input* base_input){/*{{{*/
+IssmDouble HydrologyDCEfficientAnalysis::EplTransmitivity(Element* element,Gauss* gauss, Input* epl_thick_input){/*{{{*/
         IssmDouble epl_transmitivity;
-        IssmDouble base_elev,prestep_head,epl_thick;
-        IssmDouble water_sheet;
+        IssmDouble epl_thick;
         IssmDouble epl_conductivity      = element->GetMaterialParameter(HydrologydcEplConductivityEnum);
-        base_input->GetInputValue(&base_elev,gauss);
-        epl_head_input->GetInputValue(&prestep_head,gauss);
         epl_thick_input->GetInputValue(&epl_thick,gauss);
 
-        water_sheet=max(0.0,(prestep_head-base_elev));
-        epl_transmitivity=epl_conductivity*min(water_sheet,epl_thick);
+        epl_transmitivity=epl_conductivity*epl_thick;
 
         return epl_transmitivity;

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

@@ -38 +38 @@
                 IssmDouble GetHydrologyKMatrixTransfer(Element* element);
                 IssmDouble GetHydrologyPVectorTransfer(Element* element, Gauss* gauss, Input* sed_head_input);
-                IssmDouble EplStoring(Element* element,Gauss* gauss, Input* epl_head_input, Input* epl_thick_input, Input* base_input);
-                IssmDouble EplTransmitivity(Element* element,Gauss* gauss, Input* epl_head_input, Input* epl_thick_input, Input* base_input);
+                IssmDouble EplStoring(Element* element,Gauss* gauss, Input* epl_thick_input);
+                IssmDouble EplTransmitivity(Element* element,Gauss* gauss, Input* epl_thick_input);
                 void HydrologyEPLGetMask(Vector<IssmDouble>* vec_mask, Vector<IssmDouble>* recurence, int* eplzigzag_counter, Element* element);
                 void HydrologyEPLGetActive(Vector<IssmDouble>* active_vec, Element* element);