Changeset 21468

Timestamp:

01/09/17 08:32:15 (8 years ago)

Author:

bdef

Message:

CHG:changes in efficient layer computations

Location:

issm/trunk-jpl/src/c

Files:

• analyses/HydrologyDCEfficientAnalysis.cpp (modified) (9 diffs)
• analyses/HydrologyDCEfficientAnalysis.h (modified) (1 diff)
• solutionsequences/solutionsequence_hydro_nonlinear.cpp (modified) (2 diffs)

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

@@ -235 +235 @@
                 basalelement    ->JacobianDeterminant(&Jdet,xyz_list,gauss);
 
-                epl_transmitivity = EplTransmitivity(basalelement,gauss,epl_thick_input);
-                epl_storing                             = EplStoring(basalelement,gauss,epl_thick_input);
+                epl_transmitivity = EplTransmitivity(basalelement,gauss,epl_thick_input,epl_head_input,base_input);
+                epl_storing                             = EplStoring(basalelement,gauss,epl_thick_input,epl_head_input,base_input);
 
                 /*Diffusivity*/
@@ -348 +348 @@
                 basalelement ->NodalFunctions(basis,gauss);
 
-                epl_storing     = EplStoring(basalelement,gauss,epl_thick_input);
+                epl_storing     = EplStoring(basalelement,gauss,epl_thick_input,epl_head_input,base_input);
                 /*Loading term*/
                 water_input->GetInputValue(&water_load,gauss);
                 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];
-                }
+                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);
-                        
+                        old_wh_input->GetInputValue(&water_head,gauss);
                         /*Dealing with the epl part of the transfer term*/
                         transfer=GetHydrologyPVectorTransfer(basalelement,gauss,sed_head_input);
@@ -425 +422 @@
 
         /*Use the dof list to index into the solution vector: */
-        /*If the EPL is not active we revert to the initialisation vallue*/
-        /*       For now we keep OldValue but it is probably not the best*/
+        /*If the EPL is not active we revert to the bedrock elevation*/
         if(active_element){
                 for(int i=0;i<numnodes;i++){
@@ -435 +431 @@
         }
         else{
-                basalelement->GetInputListOnVertices(&eplHeads[0],EplHeadOldEnum);
+                basalelement->GetInputListOnVertices(&eplHeads[0],BaseEnum);
                 for(int i=0;i<numnodes;i++){
                         if(xIsNan<IssmDouble>(eplHeads[i])) _error_("NaN found in solution vector");
@@ -519 +515 @@
                 
                 Element* element=(Element*)femmodel->elements->GetObjectByOffset(j);
-                /* element->parameters->FindParam(&iseplthickcomp,HydrologydcEplThickCompEnum); */
-                /* if(iseplthickcomp==0) return; */
                 
                 switch(domaintype){
@@ -576 +570 @@
                                 /*Compute first the effective pressure in the EPL*/
                                 EPL_N=gravity*((rho_ice*ice_thickness[i])-(rho_water*(eplhead[i]-bed[i])));
-                                if(EPL_N<0.0)EPL_N=0.0;
+                                //                              if(EPL_N<0.0)EPL_N=0.0;
                                 /*Get then the square of the gradient of EPL heads*/
                                 EPLgrad2 = (epl_slopeX[i]*epl_slopeX[i])+(epl_slopeY[i]*epl_slopeY[i]);
@@ -688 +682 @@
                 else if(old_active[i]>0.){
                         vec_mask->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL);
-                        /*If epl thickness gets under colapse thickness, close the layer*/
+                        /* 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); */
-                        /* } */
+                        /*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);
+                        }
                 }
                 /*Increase of the efficient system is needed if the epl head reach the maximum value (sediment max value for now)*/
@@ -724 +718 @@
 }
 /*}}}*/
-IssmDouble HydrologyDCEfficientAnalysis::EplStoring(Element* element,Gauss* gauss, Input* epl_thick_input){/*{{{*/
+IssmDouble HydrologyDCEfficientAnalysis::EplStoring(Element* element,Gauss* gauss, Input* epl_thick_input, Input* epl_head_input, Input* base_input){/*{{{*/
         IssmDouble epl_storing;
-        IssmDouble epl_thick;
+        IssmDouble water_sheet,storing;
+        IssmDouble epl_thickness,prestep_head,base_elev;
         IssmDouble rho_freshwater        = element->GetMaterialParameter(MaterialsRhoFreshwaterEnum);
         IssmDouble g                     = element->GetMaterialParameter(ConstantsGEnum);
@@ -733 +728 @@
         IssmDouble water_compressibility = element->GetMaterialParameter(HydrologydcWaterCompressibilityEnum);
 
-        epl_thick_input->GetInputValue(&epl_thick,gauss);
+        epl_thick_input->GetInputValue(&epl_thickness,gauss);
+        epl_head_input->GetInputValue(&prestep_head,gauss);
+        base_input->GetInputValue(&base_elev,gauss);
+        water_sheet=max(0.0,(prestep_head-base_elev));
+
+        storing=rho_freshwater*g*epl_porosity*epl_thickness*(water_compressibility+(epl_compressibility/epl_porosity));
+
+        //porosity for unconfined region
+        if (water_sheet<=0.9*epl_thickness){
+                epl_storing=epl_porosity;
+        }
+        //continuity ramp
+        else if((water_sheet<epl_thickness) && (water_sheet>0.9*epl_thickness)){
+                epl_storing=(epl_thickness-water_sheet)*(epl_porosity-storing)/(0.1*epl_thickness)+storing;
+        }
+        //storing coefficient for confined
+        else{
+                epl_storing=storing;
+        }
+        //return storing;
+        return epl_storing;
+}/*}}}*/
+
+IssmDouble HydrologyDCEfficientAnalysis::EplTransmitivity(Element* element,Gauss* gauss, Input* epl_thick_input, Input* epl_head_input, Input* base_input){/*{{{*/
+        IssmDouble epl_transmitivity;
+        IssmDouble water_sheet;
+        IssmDouble epl_thickness,base_elev,prestep_head;
+        IssmDouble epl_conductivity      = element->GetMaterialParameter(HydrologydcEplConductivityEnum);
+        epl_thick_input->GetInputValue(&epl_thickness,gauss);
+        epl_head_input->GetInputValue(&prestep_head,gauss);
+        base_input->GetInputValue(&base_elev,gauss);
+
+        water_sheet=max(0.0,(prestep_head-base_elev));
         
-        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_thick_input){/*{{{*/
-        IssmDouble epl_transmitivity;
-        IssmDouble epl_thick;
-        IssmDouble epl_conductivity      = element->GetMaterialParameter(HydrologydcEplConductivityEnum);
-        epl_thick_input->GetInputValue(&epl_thick,gauss);
-
-        epl_transmitivity=epl_conductivity*epl_thick;
-
+        //epl_transmitivity=epl_conductivity*epl_thickness;
+        epl_transmitivity=epl_conductivity*min(water_sheet,epl_thickness);
         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_thick_input);
-                IssmDouble EplTransmitivity(Element* element,Gauss* gauss, Input* epl_thick_input);
+                IssmDouble EplStoring(Element* element,Gauss* gauss, Input* epl_thick_input, Input* epl_head_input, Input* base_input);
+                IssmDouble EplTransmitivity(Element* element,Gauss* gauss, Input* epl_thick_input, Input* epl_head_input, Input* base_input);
                 void HydrologyEPLGetMask(Vector<IssmDouble>* vec_mask, Vector<IssmDouble>* recurence, int* eplzigzag_counter, Element* element);
                 void HydrologyEPLGetActive(Vector<IssmDouble>* active_vec, Element* element);

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

@@ -170 +170 @@
                                 /*{{{*//*Retrieve the EPL head slopes and compute EPL Thickness*/
                                 if(VerboseSolution()) _printf0_("computing EPL Head slope...\n");
-                                //inefanalysis->ElementizeEplMask(femmodel);
                                 femmodel->SetCurrentConfiguration(L2ProjectionEPLAnalysisEnum);
                                 femmodel->UpdateConstraintsL2ProjectionEPLx(&L2Count);
@@ -182 +181 @@
                                 effanalysis->ComputeEPLThickness(femmodel);
                                 //updating mask after the computation of the epl thickness (Allow to close too thin EPL)
-                                femmodel->HydrologyEPLupdateDomainx(&ThickCount);
-                                inefanalysis->ElementizeEplMask(femmodel);
+                                /* femmodel->HydrologyEPLupdateDomainx(&ThickCount); */
+                                /* inefanalysis->ElementizeEplMask(femmodel); */
                                 /*}}}*/