Changeset 17882

Timestamp:

04/29/14 14:14:29 (11 years ago)

Author:

bdef

Message:

CHG:reverting hydroshreve from Johnson to LeBrocq

Location:

issm/trunk-jpl/src/c

Files:

• analyses/HydrologyShreveAnalysis.cpp (modified) (4 diffs)
• classes/Materials/Matpar.cpp (modified) (1 diff)

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

@@ -331 +331 @@
         IssmDouble  rho_water = element->GetMaterialParameter(MaterialsRhoWaterEnum);
         IssmDouble  g         = element->GetMaterialParameter(ConstantsGEnum);
-        IssmDouble  CR        = element->GetMaterialParameter(HydrologyshreveCREnum);
-        IssmDouble  n_man     = element->GetMaterialParameter(HydrologyshreveNEnum);
         IssmDouble  mu_water  = element->GetMaterialParameter(MaterialsMuWaterEnum);
         Input* surfaceslopex_input = element->GetInput(SurfaceSlopeXEnum); _assert_(surfaceslopex_input);
@@ -339 +337 @@
         Input* bedslopey_input     = element->GetInput(BedSlopeYEnum);     _assert_(bedslopey_input);
         Input* watercolumn_input   = element->GetInput(WatercolumnEnum);   _assert_(watercolumn_input);
-
-        /* compute VelocityFactor */
-        IssmDouble VelocityFactor = n_man*CR*CR*rho_water*g/mu_water;
 
         /*Fetch number of vertices and allocate output*/
@@ -358 +353 @@
 
                 /* Water velocity x and y components */
-        //      vx[iv]= - w*w/(12 * mu_water)*(rho_ice*g*dsdx+(rho_water-rho_ice)*g*dbdx);
-        //      vy[iv]= - w*w/(12 * mu_water)*(rho_ice*g*dsdy+(rho_water-rho_ice)*g*dbdy);
-                vx[iv]= - w*w/(VelocityFactor* mu_water)*(rho_ice*g*dsdx+(rho_water-rho_ice)*g*dbdx);
-                vy[iv]= - w*w/(VelocityFactor* mu_water)*(rho_ice*g*dsdy+(rho_water-rho_ice)*g*dbdy);
+                vx[iv]= - w*w/(12 * mu_water)*(rho_ice*g*dsdx+(rho_water-rho_ice)*g*dbdx);
+                vy[iv]= - w*w/(12 * mu_water)*(rho_ice*g*dsdy+(rho_water-rho_ice)*g*dbdy);
         }
 
@@ -373 +366 @@
         xDelete<IssmDouble>(vy);
 }/*}}}*/
+
+
+
+/*Needed changes to switch to the Johnson formulation*//*{{{*/
+/*All the changes are to be done in the velocity computation.
+        The new velocity needs some new parameter that should be introduce in the hydrologyshreve class:
+        'p' and 'q' which are the exponent of the Manning formula for laminar (p=2,q=1) or turbulent (p=2/3,q=1/2) flow
+        'R' the hydraulic radius
+        'n' the manning roughness coeficient
+
+        With these, the velocity reads ;
+
+        v= - (1/n)* pow(R,p)*pow((grad phi(rho_water*g)),q)
+
+        you should also redefine the water pressure potential 'phi' with respect to the effective pressure deffinition given in Johson:
+        phi=(rho_ice*g*( surface + ((rho_water/rho_ice)-1)*base - k_n*((thickness* grad(base))/omega) )
+
+        where 
+        'omega' is the fractional area of the base occupied by the water film
+        'k_n' is a parameter
+        This last equation derives from the effective pressure definition developped in Alley 1989
+*/
+/*}}}*/

issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp ¶

@@ -60 +60 @@
                         _error_("Surface mass balance model "<<EnumToStringx(smb_model)<<" not supported yet");
         }
-
-        if(hydrology_model==HydrologyshreveEnum){
-                iomodel->Constant(&this->hydro_CR,HydrologyshreveCREnum);
-                iomodel->Constant(&this->hydro_kn,HydrologyshreveKnEnum);
-                iomodel->Constant(&this->hydro_n,HydrologyshreveNEnum);
-                iomodel->Constant(&this->hydro_p,HydrologyshrevePEnum);
-                iomodel->Constant(&this->hydro_q,HydrologyshreveQEnum);
-        }
-        else if(hydrology_model==HydrologydcEnum){
+        if(hydrology_model==HydrologydcEnum){
                 iomodel->Constant(&this->sediment_compressibility,HydrologydcSedimentCompressibilityEnum);
                 iomodel->Constant(&this->sediment_porosity,HydrologydcSedimentPorosityEnum);