Changeset 17884

Timestamp:

04/29/14 15:09:47 (11 years ago)

Author:

Mathieu Morlighem

Message:

BUG: fixed Shreve's new implementation

Location:

issm/trunk-jpl/src/c/classes

Files:

• Elements/Tria.cpp (modified) (1 diff)
• Elements/Tria.h (modified) (1 diff)
• Materials/Matpar.cpp (modified) (3 diffs)
• Materials/Matpar.h (modified) (2 diffs)

issm/trunk-jpl/src/c/classes/Elements/Tria.cpp

@@ -4244 +4244 @@
 }
 /*}}}*/
-/*FUNCTION Tria::CreateHydrologyWaterVelocityInput {{{*/
-void Tria::CreateHydrologyWaterVelocityInput(void){
-
-        /*material parameters: */
-        IssmDouble mu_water;
-        IssmDouble VelocityFactor;  // This factor represents the number 12 in laminar flow velocity which can vary by differnt hydrology.CR
-        IssmDouble n_man,CR;
-        IssmDouble w;
-        IssmDouble rho_ice, rho_water, g;
-        IssmDouble dsdx,dsdy,dbdx,dbdy;
-        IssmDouble vx[NUMVERTICES];
-        IssmDouble vy[NUMVERTICES];
-        GaussTria *gauss = NULL;
-
-        /*Retrieve all inputs and parameters*/
-        rho_ice=matpar->GetRhoIce();
-        rho_water=matpar->GetRhoWater();
-        g=matpar->GetG();
-        CR=matpar->GetHydrologyCR(); // To have Lebrocq equavalent equation: CR=0.01,n_man=0.02
-        n_man=matpar->GetHydrologyN(); 
-        mu_water=matpar->GetMuWater();
-        Input* surfaceslopex_input=inputs->GetInput(SurfaceSlopeXEnum); _assert_(surfaceslopex_input);
-        Input* surfaceslopey_input=inputs->GetInput(SurfaceSlopeYEnum); _assert_(surfaceslopey_input);
-        Input* bedslopex_input=inputs->GetInput(BedSlopeXEnum);         _assert_(bedslopex_input);
-        Input* bedslopey_input=inputs->GetInput(BedSlopeYEnum);         _assert_(bedslopey_input);
-        Input* watercolumn_input=inputs->GetInput(WatercolumnEnum);     _assert_(watercolumn_input);
-
-        /* compute VelocityFactor */
-        VelocityFactor= n_man*CR*CR*rho_water*g/mu_water;
-
-        gauss=new GaussTria();
-        for(int iv=0;iv<NUMVERTICES;iv++){
-                gauss->GaussVertex(iv);
-                surfaceslopex_input->GetInputValue(&dsdx,gauss);
-                surfaceslopey_input->GetInputValue(&dsdy,gauss);
-                bedslopex_input->GetInputValue(&dbdx,gauss);
-                bedslopey_input->GetInputValue(&dbdy,gauss);
-                watercolumn_input->GetInputValue(&w,gauss);
-
-                /* 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);
-        }
-
-        /*clean-up*/
-        delete gauss;
-
-        /*Add to inputs*/
-        this->inputs->AddInput(new TriaInput(HydrologyWaterVxEnum,vx,P1Enum));
-        this->inputs->AddInput(new TriaInput(HydrologyWaterVyEnum,vy,P1Enum));
-}
-/*}}}*/
 /*FUNCTION Tria::GetSolutionFromInputsOneDof{{{*/
 void  Tria::GetSolutionFromInputsOneDof(Vector<IssmDouble>* solution, int enum_type){

issm/trunk-jpl/src/c/classes/Elements/Tria.h

@@ -226 +226 @@
                 void UpdateConstraintsExtrudeFromBase(void);
                 void UpdateConstraintsExtrudeFromTop(void);
-
-                void           CreateHydrologyWaterVelocityInput(void);
                 /*}}}*/
 

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

@@ -74 +74 @@
                                 iomodel->Constant(&this->epl_conductivity,HydrologydcEplConductivityEnum);
                 }
+        }
+        else if(hydrology_model==HydrologyshreveEnum){
+                /*Nothing to add*/
         }
         else{
@@ -258 +261 @@
                 case HydrologydcEplMaxThicknessEnum:         return this->epl_max_thickness;
                 case HydrologydcWaterCompressibilityEnum:    return this->water_compressibility;
-                case HydrologyshreveCREnum:                  return this->hydro_CR;
-                case HydrologyshreveKnEnum:                  return this->hydro_kn;
-                case HydrologyshreveNEnum:                   return this->hydro_n;
-                case HydrologyshrevePEnum:                   return this->hydro_p;
-                case HydrologyshreveQEnum:                   return this->hydro_q;
                 case ConstantsGEnum:                         return this->g;
                 default: _error_("Enum "<<EnumToStringx(enum_in)<<" not supported yet");
@@ -350 +348 @@
 }
 /*}}}*/
-/*FUNCTION Matpar::GetHydrologyKn {{{*/
-IssmDouble Matpar::GetHydrologyKn(){
-        return hydro_kn;                 
-}                
-/*}}}*/                  
-/*FUNCTION Matpar::GetHydrologyP {{{*/                   
-IssmDouble Matpar::GetHydrologyP(){              
-        return hydro_p;                  
-}                
-/*}}}*/                  
-/*FUNCTION Matqar::GetHydrologyQ {{{*/                   
-IssmDouble Matpar::GetHydrologyQ(){              
-        return hydro_q;                  
-}                
-/*}}}*/                  
-/*FUNCTION Matpar::GetHydrologyCR {{{*/          
-IssmDouble Matpar::GetHydrologyCR(){             
-        return hydro_CR;                 
-}                
-/*}}}*/                  
-/*FUNCTION Matpar::GetHydrologyN {{{*/                   
-IssmDouble Matpar::GetHydrologyN(){              
-        return hydro_n;                  
-}                
-/*}}}*/ 
 /*FUNCTION Matpar::GetSedimentStoring {{{*/
 IssmDouble Matpar::GetSedimentStoring(){

issm/trunk-jpl/src/c/classes/Materials/Matpar.h

@@ -32 +32 @@
                 IssmDouble  desfac;
                 IssmDouble  s0p;
-
-                /*hydrology Shreve: */   
-                IssmDouble  hydro_kn;                    
-                IssmDouble  hydro_p;             
-                IssmDouble  hydro_q;             
-                IssmDouble  hydro_CR;                    
-                IssmDouble  hydro_n; 
 
                 /*hydrology Dual Porous Continuum: */    
@@ -113 +106 @@
                 IssmDouble GetMeltingPoint();
                 IssmDouble GetReferenceTemperature();
-                IssmDouble GetHydrologyKn();
-                IssmDouble GetHydrologyP();
-                IssmDouble GetHydrologyQ();
-                IssmDouble GetHydrologyCR();
-                IssmDouble GetHydrologyN();
                 IssmDouble GetSedimentStoring();
                 IssmDouble GetEplSpecificStoring();