Changeset 26491

Timestamp:

10/19/21 14:15:16 (3 years ago)

Author:

aleahsommers

Message:

CHG: includes diffusivity-like term for gap height - this is an experiment

File:

• issm/trunk-jpl/src/c/analyses/HydrologyShaktiAnalysis.cpp (modified) (8 diffs)

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

@@ -272 +272 @@
    IssmDouble  PMPheat,dpressure_water[2],dbed[2];      
         IssmDouble* xyz_list = NULL;
-        IssmDouble dgap[2]; /***/
+        IssmDouble dgapxx; /***/
+        IssmDouble dgapyy; /***/
 
         /*Fetch number of nodes and dof for this finite element*/
@@ -298 +299 @@
         Input* br_input             = element->GetInput(HydrologyBumpHeightEnum);        _assert_(br_input);
    Input* headold_input        = element->GetInput(HydrologyHeadOldEnum);           _assert_(headold_input);
-//      Input* dgapxx_input         = element->GetInput(HydrologyGapXXEnum); /***/
-//      Input* dgapyy_input         = element->GetInput(HydrologyGapYYEnum); /***/
+        Input* dgapxx_input         = element->GetInput(HydrologyGapHeightXXEnum); /***/
+        Input* dgapyy_input         = element->GetInput(HydrologyGapHeightYYEnum); /***/
 
         /*Get conductivity from inputs*/
@@ -329 +330 @@
                 br_input->GetInputValue(&br,gauss);
                 headold_input->GetInputValue(&head_old,gauss);
+                dgapxx_input->GetInputValue(&dgapxx,gauss); /***/
+                dgapyy_input->GetInputValue(&dgapyy,gauss); /***/
 
                 /*Get ice A parameter*/
@@ -382 +385 @@
 
                  /* With weighted creep term*/
-                 for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*
-                  (
-                   meltrate*(1/rho_water-1/rho_ice)
-                   +A*pow(fabs(pressure_ice - pressure_water),n-1)*(pressure_ice + rho_water*g*bed)*gap
-                   +(n-1)*A*pow(fabs(pressure_ice - pressure_water),n-1)*(rho_water*g*head)*gap
-                   -beta*sqrt(vx*vx+vy*vy)
-                   +ieb
-                   +storage*head_old/dt
-                   )*basis[i];
+//                 for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*
+ //                 (
+  //                 meltrate*(1/rho_water-1/rho_ice)
+   //                +A*pow(fabs(pressure_ice - pressure_water),n-1)*(pressure_ice + rho_water*g*bed)*gap
+    //               +(n-1)*A*pow(fabs(pressure_ice - pressure_water),n-1)*(rho_water*g*head)*gap
+     //              -beta*sqrt(vx*vx+vy*vy)
+      //             +ieb
+       //            +storage*head_old/dt
+        //           )*basis[i];
+
+                 /* With weighted creep term*/
+                  for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*
+                   (
+                    meltrate*(1/rho_water-1/rho_ice)
+                    +A*pow(fabs(pressure_ice - pressure_water),n-1)*(pressure_ice + rho_water*g*bed)*gap
+                    +(n-1)*A*pow(fabs(pressure_ice - pressure_water),n-1)*(rho_water*g*head)*gap
+                    -beta*sqrt(vx*vx+vy*vy)
+                    +ieb
+                    +storage*head_old/dt
+                    +meltrate*gap/rho_ice*(dgapxx+dgapyy)
+                    )*basis[i];
 
         
@@ -548 +563 @@
         IssmDouble q = 0.;
    IssmDouble channelization = 0.;
+        IssmDouble dgapxx; /***/
+        IssmDouble dgapyy; /***/
 
         /*Retrieve all inputs and parameters*/
@@ -565 +582 @@
         Input* lr_input             = element->GetInput(HydrologyBumpSpacingEnum);       _assert_(lr_input);
         Input* br_input             = element->GetInput(HydrologyBumpHeightEnum);        _assert_(br_input);
+        Input* dgapxx_input         = element->GetInput(HydrologyGapHeightXXEnum); /***/
+        Input* dgapyy_input         = element->GetInput(HydrologyGapHeightYYEnum); /***/
 
         /*Get conductivity from inputs*/
@@ -590 +609 @@
                 lr_input->GetInputValue(&lr,gauss);
                 br_input->GetInputValue(&br,gauss);
+                dgapxx_input->GetInputValue(&dgapxx,gauss); /***/
+                dgapyy_input->GetInputValue(&dgapyy,gauss); /***/
 
                 /*Get ice A parameter*/
@@ -628 +649 @@
                 _assert_(meltrate>0.);
 
-                newgap += gauss->weight*Jdet*(gap+dt*(
-                                        meltrate/rho_ice
-                                        -A*pow(fabs(pressure_ice-pressure_water),n-1)*(pressure_ice-pressure_water)*gap
-                                        +beta*sqrt(vx*vx+vy*vy)
-                                        ));
+//              newgap += gauss->weight*Jdet*(gap+dt*(
+//                                      meltrate/rho_ice
+//                                      -A*pow(fabs(pressure_ice-pressure_water),n-1)*(pressure_ice-pressure_water)*gap
+//                                      +beta*sqrt(vx*vx+vy*vy)
+//                                      ));
 
                 /* TEST with gap height "diffusivity" for melting walls - not yet completed*/
-//                newgap += gauss->weight*Jdet*(gap+dt*(
-//                                         meltrate/rho_ice
-//                                         -A*pow(fabs(pressure_ice-pressure_water),n-1)*(pressure_ice-pressure_water)*lc
-//                                         +beta*sqrt(vx*vx+vy*vy)
-//                                       +(meltrate*gap/rho_ice*(dgapxx+dgapyy))
-//                                         ));
+                newgap += gauss->weight*Jdet*(gap+dt*(
+                                         meltrate/rho_ice
+                                         -A*pow(fabs(pressure_ice-pressure_water),n-1)*(pressure_ice-pressure_water)*lc
+                                         +beta*sqrt(vx*vx+vy*vy)
+                                         +(meltrate*gap/rho_ice*(dgapxx+dgapyy))
+                                         ));