Changeset 23951

Timestamp:

05/29/19 15:56:20 (6 years ago)

Author:

Mathieu Morlighem

Message:

CHG: added transient term

File:

• issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp (modified) (7 diffs)

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

@@ -112 +112 @@
         parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.englacial_void_ratio",HydrologyEnglacialVoidRatioEnum));
 
+        /*Deal with friction parameters*/
+        int frictionlaw;
+        iomodel->FindConstant(&frictionlaw,"md.friction.law");
+        if(frictionlaw==4 || frictionlaw==6){
+                parameters->AddObject(iomodel->CopyConstantObject("md.friction.gamma",FrictionGammaEnum));
+        }
+        if(frictionlaw==3 || frictionlaw==1 || frictionlaw==7){
+                parameters->AddObject(iomodel->CopyConstantObject("md.friction.coupling",FrictionCouplingEnum));
+        }
+        if(frictionlaw==9){
+                parameters->AddObject(iomodel->CopyConstantObject("md.friction.gamma",FrictionGammaEnum));
+                parameters->AddObject(new IntParam(FrictionCouplingEnum,0));
+        }
+
   /*Requested outputs*/
   iomodel->FindConstant(&requestedoutputs,&numoutputs,"md.hydrology.requested_outputs");
@@ -153 +167 @@
         /*Get all inputs and parameters*/
         IssmDouble dt        = element->FindParam(TimesteppingTimeStepEnum);
-        IssmDouble c_t       = element->FindParam(HydrologyPressureMeltCoefficientEnum);
         IssmDouble rho_water = element->FindParam(MaterialsRhoFreshwaterEnum);
         IssmDouble g         = element->FindParam(ConstantsGEnum);
+        IssmDouble e_v       = element->FindParam(HydrologyEnglacialVoidRatioEnum);
         Input* k_input   = element->GetInput(HydrologySheetConductivityEnum);_assert_(k_input);
         Input* phi_input = element->GetInput(HydraulicPotentialEnum);      _assert_(phi_input);
@@ -179 +193 @@
                 IssmDouble coeff = k*pow(h,alpha)*pow(normgradphi,beta-2.);
 
+                /*Diffusive term*/
                 for(int i=0;i<numnodes;i++){
                         for(int j=0;j<numnodes;j++){
-
-                                /*Diffusive term*/
                                 Ke->values[i*numnodes+j] += gauss->weight*Jdet*(
                                                         coeff*dbasis[0*numnodes+i]*dbasis[0*numnodes+j]
@@ -189 +202 @@
                 }
 
-                /*Transient term if dt*/
+                /*Transient term if dt>0*/
                 if(dt>0.){
+                        /*Diffusive term*/
+                        for(int i=0;i<numnodes;i++){
+                                for(int j=0;j<numnodes;j++){
+                                        Ke->values[i*numnodes+j] += gauss->weight*Jdet*e_v/(rho_water*g*dt)*basis[i]*basis[j];
+                                }
+                        }
                 }
+
         }
 
@@ -209 +229 @@
         IssmDouble  Jdet,w,v,vx,vy,ub,h,N,h_r;
         IssmDouble  G,m,frictionheat,alpha2;
-        IssmDouble  A,B,n;
+        IssmDouble  A,B,n,phi_old;
         IssmDouble* xyz_list = NULL;
 
@@ -221 +241 @@
         /*Retrieve all inputs and parameters*/
         element->GetVerticesCoordinates(&xyz_list);
-        IssmDouble  L       = element->FindParam(MaterialsLatentheatEnum);
-        IssmDouble  rho_ice = element->FindParam(MaterialsRhoIceEnum);
-        IssmDouble l_r      = element->FindParam(HydrologyCavitySpacingEnum);
-        Input* hr_input = element->GetInput(HydrologyBumpHeightEnum);_assert_(hr_input);
-        Input* vx_input = element->GetInput(VxEnum);_assert_(vx_input);
-        Input* vy_input = element->GetInput(VyEnum);_assert_(vy_input);
-        Input*  N_input = element->GetInput(EffectivePressureEnum); _assert_(N_input);
-        Input*  h_input = element->GetInput(HydrologySheetThicknessEnum);_assert_(h_input);
-        Input*  G_input = element->GetInput(BasalforcingsGeothermalfluxEnum);_assert_(G_input);
-        Input* B_input  = element->GetInput(MaterialsRheologyBEnum);         _assert_(B_input);
-        Input* n_input  = element->GetInput(MaterialsRheologyNEnum);         _assert_(n_input);
+        IssmDouble L         = element->FindParam(MaterialsLatentheatEnum);
+        IssmDouble rho_ice   = element->FindParam(MaterialsRhoIceEnum);
+        IssmDouble rho_water = element->FindParam(MaterialsRhoFreshwaterEnum);
+        IssmDouble l_r       = element->FindParam(HydrologyCavitySpacingEnum);
+        IssmDouble dt        = element->FindParam(TimesteppingTimeStepEnum);
+        IssmDouble g         = element->FindParam(ConstantsGEnum);
+        IssmDouble e_v       = element->FindParam(HydrologyEnglacialVoidRatioEnum);
+        Input* hr_input     = element->GetInput(HydrologyBumpHeightEnum);_assert_(hr_input);
+        Input* vx_input     = element->GetInput(VxEnum);_assert_(vx_input);
+        Input* vy_input     = element->GetInput(VyEnum);_assert_(vy_input);
+        Input* N_input      = element->GetInput(EffectivePressureEnum); _assert_(N_input);
+        Input* h_input      = element->GetInput(HydrologySheetThicknessEnum);_assert_(h_input);
+        Input* G_input      = element->GetInput(BasalforcingsGeothermalfluxEnum);_assert_(G_input);
+        Input* B_input      = element->GetInput(MaterialsRheologyBEnum);         _assert_(B_input);
+        Input* n_input      = element->GetInput(MaterialsRheologyNEnum);         _assert_(n_input);
+        Input* phiold_input = element->GetInput(HydraulicPotentialOldEnum);      _assert_(phiold_input);
 
         /*Build friction element, needed later: */
@@ -273 +298 @@
 
                 for(int i=0;i<numnodes;i++) pe->values[i]+= - Jdet*gauss->weight*(w-v-m)*basis[i];
+
+                /*Transient term if dt>0*/
+                if(dt>0.){
+                        phiold_input->GetInputValue(&phi_old,gauss);
+                        for(int i=0;i<numnodes;i++) pe->values[i] += gauss->weight*Jdet*e_v/(rho_water*g*dt)*phi_old*basis[i];
+                }
         }