Changeset 25252 for issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp

Timestamp:

07/10/20 10:50:37 (5 years ago)

Author:

jdquinn

Message:

BUG: Reverting recent changes until Basile can review

File:

• issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp (modified) (16 diffs)

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

@@ -238 +238 @@
         Input2* sed_head_input = basalelement->GetInput2(SedimentHeadSubstepEnum);
         Input2* base_input     = basalelement->GetInput2(BaseEnum);
+        Input2* old_wh_input   = basalelement->GetInput2(SedimentHeadOldEnum);                  _assert_(old_wh_input);
 
         /*Transfer related Inputs*/
@@ -243 +244 @@
                 basalelement->GetInput2Value(&active_element,HydrologydcMaskEplactiveEltEnum);
         }
+
         /* 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);
+                gauss          -> GaussPoint(ig);
+                basalelement   -> JacobianDeterminant(&Jdet,xyz_list,gauss);
                 basalelement->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
                 basalelement->NodalFunctions(basis,gauss);
@@ -257 +259 @@
                 /*Diffusivity*/
                 D_scalar=sediment_transmitivity*gauss->weight*Jdet;
+                //D_scalar=gauss->weight*Jdet;
                 if(dt!=0.) D_scalar=D_scalar*dt;
                 for(int i=0;i<numnodes;i++){
@@ -263 +266 @@
                         }
                 }
+
                 /*Transient*/
                 if(dt!=0.){
                         D_scalar=sediment_storing*gauss->weight*Jdet;
+                        //D_scalar=(sediment_storing/sediment_transmitivity)*gauss->weight*Jdet;
                         for(int i=0;i<numnodes;i++) for(int j=0;j<numnodes;j++) Ke->values[i*numnodes+j] += D_scalar*basis[j]*basis[i];
+
                         /*Transfer EPL part*/
                         if(isefficientlayer){
@@ -272 +278 @@
                                         transfer=GetHydrologyKMatrixTransfer(basalelement);
                                         D_scalar=dt*transfer*gauss->weight*Jdet;
+                                        //D_scalar=dt*(transfer/sediment_transmitivity)*gauss->weight*Jdet;
                                         for(int i=0;i<numnodes;i++) for(int j=0;j<numnodes;j++) Ke->values[i*numnodes+j] += D_scalar*basis[j]*basis[i];
                                 }
@@ -291 +298 @@
 
         /*Intermediaries*/
-        bool            thawed_element;
-        int             domaintype;
+        bool             thawed_element;
+        int                      domaintype;
         Element* basalelement;
 
@@ -321 +328 @@
         /*Intermediaries */
         bool       active_element,isefficientlayer;
-        int        smb_model,smbsubstepping;
-        int        hydrologysubstepping,smb_averaging;
+        int        smb_model;
+        int        smbsubstepping;
+        int        hydrologysubstepping;
+        int        smb_averaging;
         IssmDouble dt,scalar,sediment_storing;
         IssmDouble water_head,sediment_transmitivity;
         IssmDouble water_load,runoff_value,transfer;
-        IssmDouble Jdet,time;
+        IssmDouble Jdet;
 
         IssmDouble *xyz_list             = NULL;
@@ -346 +355 @@
         basalelement->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum);
         basalelement->FindParam(&smb_model,SmbEnum);
+        basalelement->FindParam(&smb_averaging,SmbAveragingEnum);
 
         Input2* sed_head_input   = basalelement->GetInput2(SedimentHeadSubstepEnum);
@@ -353 +363 @@
         Input2* SedTrans_input   = basalelement->GetInput2(HydrologydcSedimentTransmitivityEnum); _assert_(SedTrans_input);
 
+        IssmDouble time;
+        basalelement->FindParam(&time,TimeEnum);
+
         if(dt!= 0.){
                 old_wh_input = basalelement->GetInput2(SedimentHeadOldEnum); _assert_(old_wh_input);
         }
         if(smb_model==SMBgradientscomponentsEnum){
-                basalelement->FindParam(&time,TimeEnum);
                 basalelement->FindParam(&smbsubstepping,SmbStepsPerStepEnum);
                 basalelement->FindParam(&hydrologysubstepping,HydrologyStepsPerStepEnum);
@@ -371 +383 @@
                 else{
                         //finer stepping in smb, we average the runoff from transient input
-                        basalelement->FindParam(&smb_averaging,SmbAveragingEnum);
                         dummy_input = basalelement->GetInput2(SmbRunoffTransientEnum,time-dt,time,smb_averaging); _assert_(dummy_input);
                 }
@@ -384 +395 @@
         /* Start  looping on the number of gaussian points: */
         Gauss* gauss=basalelement->NewGauss(2);
+
+        IssmDouble yts;
+        basalelement->FindParam(&yts,ConstantsYtsEnum);
 
         for(int ig=gauss->begin();ig<gauss->end();ig++){
@@ -397 +411 @@
                         else                     runoff_value = 0.;
                         scalar = Jdet*gauss->weight*(water_load+runoff_value);
+                        //scalar = Jdet*gauss->weight*(water_load)/sediment_transmitivity;
                         if(dt!=0.) scalar = scalar*dt;
                         for(int i=0;i<numnodes;i++){
@@ -409 +424 @@
                                 else runoff_value = 0.;
                                 scalar = Jdet*gauss->weight*(water_load+runoff_value);
+                                //scalar = Jdet*gauss->weight*(water_load)/sediment_transmitivity;
                                 if(dt!=0.) scalar = scalar*dt;
                                 for(int i=0;i<numnodes;i++){
@@ -429 +445 @@
                                 }
                                 scalar = Jdet*gauss->weight*((water_head*sediment_storing)+(dt*transfer));
+                                //scalar = Jdet*gauss->weight*((water_head*sediment_storing)+(dt*transfer))/sediment_transmitivity;
                                 for(int i=0;i<numnodes;i++)pe->values[i]+=scalar*basis[i];
                         }
                         else{
                                 scalar = Jdet*gauss->weight*(water_head*sediment_storing);
+                                //scalar = Jdet*gauss->weight*(water_head*sediment_storing)/sediment_transmitivity;
                                 for(int i=0;i<numnodes;i++)pe->values[i]+=scalar*basis[i];
                         }
@@ -553 +571 @@
         IssmDouble storing,yield;
         IssmDouble base_elev,prestep_head,water_sheet;
-        IssmDouble porewater_mass           = element->FindParam(HydrologydcSedimentPoreWaterMassEnum);
-        IssmDouble layer_compressibility    = element->FindParam(HydrologydcSedimentLayerCompressibilityEnum);
+        IssmDouble rho_freshwater           = element->FindParam(MaterialsRhoFreshwaterEnum);
+        IssmDouble g                        = element->FindParam(ConstantsGEnum);
+        IssmDouble sediment_porosity        = element->FindParam(HydrologydcSedimentPorosityEnum);
         IssmDouble sediment_thickness       = element->FindParam(HydrologydcSedimentThicknessEnum);
+        IssmDouble sediment_compressibility = element->FindParam(HydrologydcSedimentCompressibilityEnum);
+        IssmDouble water_compressibility    = element->FindParam(HydrologydcWaterCompressibilityEnum);
         element->FindParam(&unconf_scheme,HydrologydcUnconfinedFlagEnum);
         switch(unconf_scheme){
         case 0:
-                sediment_storing=porewater_mass*sediment_thickness*layer_compressibility;
+                sediment_storing=rho_freshwater*g*sediment_porosity*sediment_thickness*(water_compressibility+(sediment_compressibility/sediment_porosity));
                 break;
         case 1:
-                yield = element->FindParam(HydrologydcSedimentPorosityEnum);
                 base_input->GetInputValue(&base_elev,gauss);
                 sed_head_input->GetInputValue(&prestep_head,gauss);
-
                 water_sheet=max(0.0,(prestep_head-(base_elev-sediment_thickness)));
-                storing=porewater_mass*sediment_thickness*layer_compressibility;
+
+                /* if (water_sheet<sediment_thickness){ */
+                /*      sediment_storing=rho_freshwater*g*sediment_porosity*sediment_thickness*(water_compressibility+(sediment_compressibility/sediment_porosity)); */
+                /* } */
+                /* else{ */
+                /*      sediment_storing=sediment_porosity; */
+                /* } */
+                storing=rho_freshwater*g*sediment_porosity*sediment_thickness*(water_compressibility+(sediment_compressibility/sediment_porosity));
                 //using logistic function for heavyside approximation
                 expfac=10.;
+                yield=sediment_porosity;
                 sediment_storing=yield+(storing-yield)/(1+exp(-2*expfac*(water_sheet-0.99*sediment_thickness)));
                 break;
@@ -579 +606 @@
 IssmDouble HydrologyDCInefficientAnalysis::SedimentTransmitivity(Element* element,Gauss* gauss,Input2* sed_head_input, Input2* base_input,Input2* SedTrans_input){/*{{{*/
         int unconf_scheme;
+        IssmDouble ratio,expfac;
         IssmDouble sediment_transmitivity;
         IssmDouble FullLayer_transmitivity;
+        IssmDouble meltingrate;
+        IssmDouble groundedice;
         IssmDouble base_elev,prestep_head,water_sheet;
         IssmDouble sediment_thickness       = element->FindParam(HydrologydcSedimentThicknessEnum);