Changeset 27225

Timestamp:

08/23/22 08:47:54 (3 years ago)

Author:

aleahsommers

Message:

CHG: cleaned up SHAKTI code

File:

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

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

@@ -242 +242 @@
                 element->NodalFunctions(basis,gauss);
 
-                /***/
                 base_input->GetInputValue(&bed,gauss);
                 thickness_input->GetInputValue(&thickness,gauss);
@@ -288 +287 @@
    IssmDouble  PMPheat,dissipation,dpressure_water[2],dbed[2];  
         IssmDouble* xyz_list = NULL;
-//        IssmDouble dgapxx; /***/
-//      IssmDouble dgapyy; /***/
 
         /*Fetch number of nodes and dof for this finite element*/
@@ -315 +312 @@
         Input* br_input             = element->GetInput(HydrologyBumpHeightEnum);        _assert_(br_input);
    Input* headold_input        = element->GetInput(HydrologyHeadOldEnum);           _assert_(headold_input);
-//      Input* dgapxx_input         = element->GetInput(HydrologyGapHeightXXEnum); /***/
-//      Input* dgapyy_input         = element->GetInput(HydrologyGapHeightYYEnum); /***/
 
         /*Get conductivity from inputs*/
@@ -346 +341 @@
                 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*/
@@ -365 +358 @@
                 frictionheat=alpha2*(vx*vx+vy*vy);
 
-                /* Take out frictional heat for large gap height */
-                if(gap>br)
-                 frictionheat=0.;
-
                 /*Get water and ice pressures*/
                 IssmDouble pressure_ice   = rho_ice*g*thickness;    _assert_(pressure_ice>0.); 
@@ -379 +368 @@
 
                 /*Compute change in sensible heat due to changes in pressure melting point*/
-        dpressure_water[0] = rho_water*g*(dh[0] - dbed[0]);
+                dpressure_water[0] = rho_water*g*(dh[0] - dbed[0]);
                 dpressure_water[1] = rho_water*g*(dh[1] - dbed[1]);
-                PMPheat=CT*CW*conductivity*rho_water*(dh[0]*dpressure_water[0]+dh[1]*dpressure_water[1]);
-                PMPheat=0; /*** TEST no PMPheat***/
-
-        meltrate = 1/latentheat*(G+frictionheat+rho_water*g*conductivity*(dh[0]*dh[0]+dh[1]*dh[1])-PMPheat);
+
+        meltrate = 1/latentheat*(G+frictionheat+rho_water*g*conductivity*(dh[0]*dh[0]+dh[1]*dh[1]));
 
                   for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*
@@ -396 +383 @@
                     )*basis[i];
 
-                /* Test with experimental diffusivity 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
-        //             +0.e-5*(dgapxx+dgapyy)
-         //            )*basis[i];
-
         
         }
@@ -458 +433 @@
                 values[i]=solution[doflist[i]];
 
-                /*overburden cap: make sure that p_water<p_ice ->  h<rho_i H/rho_w + zb*/
-//              if(values[i]>rho_ice*thickness[i]/rho_water+bed[i]){
-//                      values[i] = rho_ice*thickness[i]/rho_water+bed[i];
-//              }
-
-                /*Make sure that negative pressure is not allowed*/
-//      if(values[i]<bed[i]){
-//                      values[i] = bed[i];
-//              }
-
                 /*Under-relaxation*/
            values[i] = head_old[i] - relaxation*(head_old[i]-values[i]);
@@ -549 +514 @@
         IssmDouble  alpha2,frictionheat;
         IssmDouble* xyz_list = NULL;
-   IssmDouble  dpressure_water[2],dbed[2],PMPheat,dissipation;
+        IssmDouble  dpressure_water[2],dbed[2],PMPheat,dissipation;
         IssmDouble q = 0.;
-   IssmDouble channelization = 0.;
-//      IssmDouble dgapxx; /***/
-//      IssmDouble dgapyy; /***/
+        IssmDouble channelization = 0.;
 
         /*Retrieve all inputs and parameters*/
@@ -571 +534 @@
         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*/
@@ -598 +559 @@
                 lr_input->GetInputValue(&lr,gauss);
                 br_input->GetInputValue(&br,gauss);
-//              dgapxx_input->GetInputValue(&dgapxx,gauss); /***/
-//              dgapyy_input->GetInputValue(&dgapyy,gauss); /***/
 
                 /*Get ice A parameter*/
@@ -617 +576 @@
                 frictionheat=alpha2*(vx*vx+vy*vy);
 
-                /* Take out frictional heat for large gap height */
-                if(gap>br)
-                 frictionheat=0.;
-
                 /*Get water and ice pressures*/
                 IssmDouble pressure_ice   = rho_ice*g*thickness;    _assert_(pressure_ice>0.); 
@@ -629 +584 @@
            dpressure_water[0] = rho_water*g*(dh[0] - dbed[0]);
                 dpressure_water[1] = rho_water*g*(dh[1] - dbed[1]);
-                PMPheat=CT*CW*conductivity*rho_water*(dh[0]*dpressure_water[0]+dh[1]*dpressure_water[1]);
-PMPheat=0; /*** TEST no PMPheat***/
                 dissipation=rho_water*g*conductivity*(dh[0]*dh[0]+dh[1]*dh[1]);
 
-                meltrate = 1/latentheat*(G+frictionheat+rho_water*g*conductivity*(dh[0]*dh[0]+dh[1]*dh[1])-PMPheat);
+                meltrate = 1/latentheat*(G+frictionheat+rho_water*g*conductivity*(dh[0]*dh[0]+dh[1]*dh[1]));
 
                 element->AddInput(DummyEnum,&meltrate,P0Enum);
@@ -647 +600 @@
                                         ));
 
-                /* TEST with experimental gap height "diffusivity" for melting walls */
-//                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)
-//                                       +0.e-5*(dgapxx+dgapyy))
- //                                        );
-
 
                 totalweights +=gauss->weight*Jdet;
@@ -667 +612 @@
         /*Divide by connectivity*/
         newgap = newgap/totalweights;
-        IssmDouble mingap = 1e-6;
+        IssmDouble mingap = 1e-3;
         if(newgap<mingap) newgap=mingap;
 
-        /*Limit gap height to grow to surface*/
-//      if(newgap>thickness)
-//       newgap = thickness;
+        /*Limit gap height*/
         if(newgap>1)
          newgap = 1;