Changeset 24385 for issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp

Timestamp:

11/22/19 02:51:08 (5 years ago)

Author:

bdef

Message:

BUG: fix in smb GradComp, and some shuffling in Hydro

File:

• issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp (modified) (5 diffs)

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

@@ -383 +383 @@
         return pe;
 }/*}}}*/
+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.
+         * For node i, Bi can be expressed in the actual coordinate system
+         * by:
+         *       Bi=[ dN/dx ]
+         *          [ dN/dy ]
+         * where N is the finiteelement function for node i.
+         *
+         * We assume B has been allocated already, of size: 3x(NDOF2*numnodes)
+         */
+
+        /*Fetch number of nodes for this finite element*/
+        int numnodes = element->GetNumberOfNodes();
+
+        /*Get nodal functions derivatives*/
+        IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
+        element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
+
+        /*Build B: */
+        for(int i=0;i<numnodes;i++){
+                B[numnodes*0+i] = dbasis[0*numnodes+i];
+                B[numnodes*1+i] = dbasis[1*numnodes+i];
+        }
+
+        /*Clean-up*/
+        xDelete<IssmDouble>(dbasis);
+}/*}}}*/
 void HydrologyDCEfficientAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
         element->GetSolutionFromInputsOneDof(solution,EplHeadSubstepEnum);
@@ -437 +464 @@
 
 /*Intermediaries*/
+IssmDouble HydrologyDCEfficientAnalysis::EplStoring(Element* element,Gauss* gauss, Input2* epl_thick_input, Input2* epl_head_input, Input2* base_input){/*{{{*/
+        IssmDouble epl_storing;
+        IssmDouble water_sheet,storing;
+        IssmDouble epl_thickness,prestep_head,base_elev;
+        IssmDouble rho_freshwater        = element->FindParam(MaterialsRhoFreshwaterEnum);
+        IssmDouble g                     = element->FindParam(ConstantsGEnum);
+        IssmDouble epl_porosity                                  = element->FindParam(HydrologydcEplPorosityEnum);
+        IssmDouble epl_compressibility   = element->FindParam(HydrologydcEplCompressibilityEnum);
+        IssmDouble water_compressibility = element->FindParam(HydrologydcWaterCompressibilityEnum);
+
+        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 epl_storing; */
+        return storing;
+}/*}}}*/
+IssmDouble HydrologyDCEfficientAnalysis::EplTransmitivity(Element* element,Gauss* gauss, Input2* epl_thick_input, Input2* epl_head_input, Input2* base_input){/*{{{*/
+        IssmDouble epl_transmitivity;
+        IssmDouble water_sheet;
+        IssmDouble epl_thickness,base_elev,prestep_head;
+        IssmDouble epl_conductivity      = element->FindParam(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_transmitivity=epl_conductivity*epl_thickness;
+        //epl_transmitivity=max(1.0e-6,(epl_conductivity*min(water_sheet,epl_thickness)));
+        return epl_transmitivity;
+}/*}}}*/
+void HydrologyDCEfficientAnalysis::GetHydrologyDCInefficientHmax(IssmDouble* ph_max,Element* element, Node* innode){/*{{{*/
+
+        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:
+                /*Compute max*/
+                rho_water = element->FindParam(MaterialsRhoFreshwaterEnum);
+                rho_ice   = element->FindParam(MaterialsRhoIceEnum);
+                element-> GetInputValue(&thickness,innode,ThicknessEnum);
+                element-> GetInputValue(&bed,innode,BaseEnum);
+                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;
+}
+/*}}}*/
 IssmDouble HydrologyDCEfficientAnalysis::GetHydrologyKMatrixTransfer(Element* element){/*{{{*/
 
@@ -582 +689 @@
         }
 }/*}}}*/
-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.
-         * For node i, Bi can be expressed in the actual coordinate system
-         * by:
-         *       Bi=[ dN/dx ]
-         *          [ dN/dy ]
-         * where N is the finiteelement function for node i.
-         *
-         * We assume B has been allocated already, of size: 3x(NDOF2*numnodes)
-         */
-
-        /*Fetch number of nodes for this finite element*/
-        int numnodes = element->GetNumberOfNodes();
-
-        /*Get nodal functions derivatives*/
-        IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
-        element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
-
-        /*Build B: */
-        for(int i=0;i<numnodes;i++){
-                B[numnodes*0+i] = dbasis[0*numnodes+i];
-                B[numnodes*1+i] = dbasis[1*numnodes+i];
-        }
-
-        /*Clean-up*/
-        xDelete<IssmDouble>(dbasis);
-}/*}}}*/
 void  HydrologyDCEfficientAnalysis::HydrologyEPLGetMask(Vector<IssmDouble>* vec_mask, Vector<IssmDouble>* recurence, Element* element){/*{{{*/
 
@@ -703 +783 @@
 }
 /*}}}*/
-IssmDouble HydrologyDCEfficientAnalysis::EplStoring(Element* element,Gauss* gauss, Input2* epl_thick_input, Input2* epl_head_input, Input2* base_input){/*{{{*/
-        IssmDouble epl_storing;
-        IssmDouble water_sheet,storing;
-        IssmDouble epl_thickness,prestep_head,base_elev;
-        IssmDouble rho_freshwater        = element->FindParam(MaterialsRhoFreshwaterEnum);
-        IssmDouble g                     = element->FindParam(ConstantsGEnum);
-        IssmDouble epl_porosity                                  = element->FindParam(HydrologydcEplPorosityEnum);
-        IssmDouble epl_compressibility   = element->FindParam(HydrologydcEplCompressibilityEnum);
-        IssmDouble water_compressibility = element->FindParam(HydrologydcWaterCompressibilityEnum);
-
-        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 epl_storing; */
-        return storing;
-}/*}}}*/
-IssmDouble HydrologyDCEfficientAnalysis::EplTransmitivity(Element* element,Gauss* gauss, Input2* epl_thick_input, Input2* epl_head_input, Input2* base_input){/*{{{*/
-        IssmDouble epl_transmitivity;
-        IssmDouble water_sheet;
-        IssmDouble epl_thickness,base_elev,prestep_head;
-        IssmDouble epl_conductivity      = element->FindParam(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_transmitivity=epl_conductivity*epl_thickness;
-        //epl_transmitivity=max(1.0e-6,(epl_conductivity*min(water_sheet,epl_thickness)));
-        return epl_transmitivity;
-}/*}}}*/
 void HydrologyDCEfficientAnalysis::HydrologyEPLGetActive(Vector<IssmDouble>* active_vec, Element* element){/*{{{*/
         /*Constants*/
@@ -796 +831 @@
 }
 /*}}}*/
-void HydrologyDCEfficientAnalysis::GetHydrologyDCInefficientHmax(IssmDouble* ph_max,Element* element, Node* innode){/*{{{*/
-
-        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:
-                /*Compute max*/
-                rho_water = element->FindParam(MaterialsRhoFreshwaterEnum);
-                rho_ice   = element->FindParam(MaterialsRhoIceEnum);
-                element-> GetInputValue(&thickness,innode,ThicknessEnum);
-                element-> GetInputValue(&bed,innode,BaseEnum);
-                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;
-}
-/*}}}*/