Changeset 22849

Timestamp:

06/18/18 13:51:28 (7 years ago)

Author:

Mathieu Morlighem

Message:

NEW: added support for different effective pressure, only works for friction right now, to be extended later

Location:

issm/trunk-jpl/src

Files:

• c/classes/Loads/Friction.cpp (modified) (2 diffs)
• m/classes/friction.m (modified) (4 diffs)
• m/classes/friction.py (modified) (4 diffs)

issm/trunk-jpl/src/c/classes/Loads/Friction.cpp

@@ -670 +670 @@
         *palpha2=alpha2;
 }/*}}}*/
+
 IssmDouble Friction::EffectivePressure(Gauss* gauss){/*{{{*/
-        /*Get effective pressure as a function of 
-         *  - coupling
-         *  - type
-         * Neff=rho_ice*g*thickness+rho_ice*g*base
-         */
-
+        /*Get effective pressure as a function of  flag */
 
         /*diverse: */
         int         coupled_flag;
+        IssmDouble  thickness,base,sealevel;
+        IssmDouble  p_ice,p_water;
         IssmDouble  Neff;
 
@@ -688 +686 @@
         switch(coupled_flag){
                 case 0:{
-                   IssmDouble  thickness,base,sealevel;
-                        element->GetInputValue(&thickness, gauss,ThicknessEnum);
-                        element->GetInputValue(&base, gauss,BaseEnum);
-                        element->GetInputValue(&sealevel, gauss,SealevelEnum);
-                        IssmDouble rho_water = element->GetMaterialParameter(MaterialsRhoSeawaterEnum);
-                        IssmDouble rho_ice   = element->GetMaterialParameter(MaterialsRhoIceEnum);
-                        IssmDouble gravity   = element->GetMaterialParameter(ConstantsGEnum);
-                        Neff=gravity*(rho_ice*thickness+rho_water*(base-sealevel));
+                                         element->GetInputValue(&thickness, gauss,ThicknessEnum);
+                                         IssmDouble rho_ice   = element->GetMaterialParameter(MaterialsRhoIceEnum);
+                                         IssmDouble gravity   = element->GetMaterialParameter(ConstantsGEnum);
+                                         p_ice   = gravity*rho_ice*thickness;
+                                         p_water = 0.;
+                                         Neff = p_ice - p_water;
                                  }
-                        break;
-                case 1:
+                          break;
+                case 1:{
+                                         element->GetInputValue(&thickness, gauss,ThicknessEnum);
+                                         element->GetInputValue(&base, gauss,BaseEnum);
+                                         element->GetInputValue(&sealevel, gauss,SealevelEnum);
+                                         IssmDouble rho_water = element->GetMaterialParameter(MaterialsRhoSeawaterEnum);
+                                         IssmDouble rho_ice   = element->GetMaterialParameter(MaterialsRhoIceEnum);
+                                         IssmDouble gravity   = element->GetMaterialParameter(ConstantsGEnum);
+                                         p_ice   = gravity*rho_ice*thickness;
+                                         p_water = max(0.,rho_water*gravity*(sealevel-base));
+                                         Neff = p_ice - p_water;
+                                 }
+                          break;
+                case 2:{
+                                         element->GetInputValue(&thickness, gauss,ThicknessEnum);
+                                         element->GetInputValue(&base, gauss,BaseEnum);
+                                         element->GetInputValue(&sealevel, gauss,SealevelEnum);
+                                         IssmDouble rho_water = element->GetMaterialParameter(MaterialsRhoSeawaterEnum);
+                                         IssmDouble rho_ice   = element->GetMaterialParameter(MaterialsRhoIceEnum);
+                                         IssmDouble gravity   = element->GetMaterialParameter(ConstantsGEnum);
+                                         p_ice   = gravity*rho_ice*thickness;
+                                         p_water = rho_water*gravity*(sealevel-base);
+                                         Neff = p_ice - p_water;
+                                 }
+                          break;
+                case 3:
                         element->GetInputValue(&Neff,gauss,FrictionEffectivePressureEnum);
                         break;
-                case 2:
+                case 4:
                         element->GetInputValue(&Neff,gauss,EffectivePressureTimeAverageEnum);
                         break;

issm/trunk-jpl/src/m/classes/friction.m

@@ -20 +20 @@
                                 case 0
                                 case 1
+                                case 2
+                                case 3
                                         self.effective_pressure=project3d(md,'vector',self.effective_pressure,'type','node','layer',1);
-                                case 2
+                                case 4
                                         error('not implemented yet');
                                 otherwise
@@ -39 +41 @@
                 function self = setdefaultparameters(self) % {{{
 
+                        self.coupling = 0;
+
                 end % }}}
                 function md = checkconsistency(self,md,solution,analyses) % {{{
@@ -49 +53 @@
                         md = checkfield(md,'fieldname','friction.q','NaN',1,'Inf',1,'size',[md.mesh.numberofelements 1]);
                         md = checkfield(md,'fieldname','friction.p','NaN',1,'Inf',1,'size',[md.mesh.numberofelements 1]);
-                        md = checkfield(md,'fieldname','friction.coupling','numel',[1],'values',[0 1 2]);
+                        md = checkfield(md,'fieldname','friction.coupling','numel',[1],'values',[0:4]);
                         switch self.coupling
                                 case 0
                                 case 1
+                                case 2
+                                case 3
                                         md = checkfield(md,'fieldname','friction.effective_pressure','NaN',1,'Inf',1,'timeseries',1);
-                                case 2
+                                case 4
                                         error('not implemented yet');
                                 otherwise
@@ -66 +72 @@
                         fielddisplay(self,'q','q exponent');
                         fielddisplay(self,'effective_pressure','Effective Pressure for the forcing if not coupled [Pa]');
-                        fielddisplay(self,'coupling','Coupling flag: 0 for default, 1 for forcing(provide md.friction.effective_pressure)  and 2 for coupled(not implemented yet)');
+                        fielddisplay(self,'coupling','Coupling flag 0: uniform sheet (negative pressure ok, default), 1: ice pressure only, 2: water pressure assuming uniform sheet (no negative pressure), 3: use provided effective_pressure, 4: used coupled model (not implemented yet)');
                 end % }}}
                 function marshall(self,prefix,md,fid) % {{{

issm/trunk-jpl/src/m/classes/friction.py

@@ -29 +29 @@
                 string="%s\n%s"%(string,fielddisplay(self,"p","p exponent"))
                 string="%s\n%s"%(string,fielddisplay(self,"q","q exponent"))
-                string="%s\n%s"%(string,fielddisplay(self,'coupling','Coupling flag: 0 for default, 1 for forcing(provide md.friction.effective_pressure)  and 2 for coupled(not implemented yet)'))
+                string="%s\n%s"%(string,fielddisplay(self,'coupling','Coupling flag 0: uniform sheet (negative pressure ok, default), 1: ice pressure only, 2: water pressure assuming uniform sheet (no negative pressure), 3: use provided effective_pressure, 4: used coupled model (not implemented yet)'))
                 string="%s\n%s"%(string,fielddisplay(self,'effective_pressure','Effective Pressure for the forcing if not coupled [Pa]'))
                 return string
@@ -38 +38 @@
                 self.q=project3d(md,'vector',self.q,'type','element')
                 #if self.coupling==0: #doesnt work with empty loop, so just skip it?
-                if self.coupling==1:
+                if self.coupling==3:
                         self.effective_pressure=project3d(md,'vector',self.effective_pressure,'type','node','layer',1)
-                elif self.coupling==2:
+                elif self.coupling==4:
                         raise ValueError('coupling not supported yet')
-                elif self.coupling > 2:
-                        raise ValueError('md.friction.coupling larger than 2, not supported yet')
+                elif self.coupling > 4:
+                        raise ValueError('md.friction.coupling larger than 4, not supported yet')
                 return self
         #}}}
@@ -58 +58 @@
                 md = checkfield(md,'fieldname','friction.q','NaN',1,'Inf',1,'size',[md.mesh.numberofelements])
                 md = checkfield(md,'fieldname','friction.p','NaN',1,'Inf',1,'size',[md.mesh.numberofelements])
-                md = checkfield(md,'fieldname','friction.coupling','numel',[1],'values',[0,1,2])
-                if self.coupling==1:
+                md = checkfield(md,'fieldname','friction.coupling','numel',[1],'values',[0,1,2,3,4])
+                if self.coupling==3:
                         md = checkfield(md,'fieldname','friction.effective_pressure','NaN',1,'Inf',1,'timeseries',1)
-                elif self.coupling > 2:
-                        raise ValueError('md.friction.coupling larger than 2, not supported yet')
+                elif self.coupling > 4:
+                        raise ValueError('md.friction.coupling larger than 4, not supported yet')
                 return md
         # }}}
@@ -71 +71 @@
                 WriteData(fid,prefix,'object',self,'fieldname','q','format','DoubleMat','mattype',2)
                 WriteData(fid,prefix,'class','friction','object',self,'fieldname','coupling','format','Integer')
-                if self.coupling==1:
+                if self.coupling==3:
                         WriteData(fid,prefix,'class','friction','object',self,'fieldname','effective_pressure','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts)
-                elif self.coupling > 2:
+                elif self.coupling > 4:
                         raise ValueError('md.friction.coupling larger than 2, not supported yet')
         # }}}