Changeset 18592

Timestamp:

10/08/14 09:13:19 (10 years ago)

Author:

bdef

Message:

BUG: fixing of python hydrology

Location:

issm/trunk-jpl/src/m

Files:

• classes/hydrologydc.py (modified) (1 diff)
• classes/initialization.m (modified) (1 diff)
• classes/initialization.py (modified) (6 diffs)
• classes/model.py (modified) (3 diffs)
• exp/expdisp.m (modified) (2 diffs)

issm/trunk-jpl/src/m/classes/hydrologydc.py ¶

@@ -125 +125 @@
                         return md
 
-                md = checkfield(md,'fieldname','hydrology.water_compressibility','>',0,'numel',1)
+                md = checkfield(md,'fieldname','hydrology.water_compressibility','>',0,'numel',1)
                 md = checkfield(md,'fieldname','hydrology.isefficientlayer','numel',[1],'values',[0,1])
                 md = checkfield(md,'fieldname','hydrology.penalty_factor','>',0,'numel',1)

issm/trunk-jpl/src/m/classes/initialization.m ¶

@@ -91 +91 @@
                                 if isa(md.hydrology,'hydrologydc'),
                                         if md.hydrology.isefficientlayer==1,
-                                                md = checkfield(md,'fieldname', ...
-                                                                                                                'initialization.epl_head','NaN',1,'size',[md.mesh.numberofvertices 1]);
+                                                md = checkfield(md,'fieldname','initialization.epl_head','NaN',1,'size',[md.mesh.numberofvertices 1]);
                                                 md = checkfield(md,'fieldname','initialization.epl_thickness','NaN',1,'size',[md.mesh.numberofvertices 1]);
                                         end

issm/trunk-jpl/src/m/classes/initialization.py ¶

@@ -9 +9 @@
         """
         INITIALIZATION class definition
-
-           Usage:
-              initialization=initialization();
+        
+        Usage:
+        initialization=initialization();
         """
 
         def __init__(self): # {{{
-                
+                                        
                 self.vx            = float('NaN')
                 self.vy            = float('NaN')
@@ -24 +24 @@
                 self.waterfraction = float('NaN')
                 self.watercolumn   = float('NaN')
-                self.sediment_head  = float('NaN')
+                self.sediment_head = float('NaN')
+                self.epl_head      = float('NaN')
+                self.epl_thickness = float('NaN')
 
                 #set defaults
@@ -32 +34 @@
         def __repr__(self): # {{{
                 string='   initial field values:'
-
                 string="%s\n%s"%(string,fielddisplay(self,'vx','x component of velocity [m/yr]'))
                 string="%s\n%s"%(string,fielddisplay(self,'vy','y component of velocity [m/yr]'))
@@ -42 +43 @@
                 string="%s\n%s"%(string,fielddisplay(self,'watercolumn','thickness of subglacial water [m]'))
                 string="%s\n%s"%(string,fielddisplay(self,'sediment_head','sediment water head of subglacial system [m]'))
+                string="%s\n%s"%(string,fielddisplay(self,'epl_head','epl water head of subglacial system [m]'))
+                string="%s\n%s"%(string,fielddisplay(self,'epl_thickness','thickness of the epl [m]'))
 
                 return string
@@ -69 +72 @@
                                 md = checkfield(md,'fieldname','initialization.vz','NaN',1,'size',[md.mesh.numberofvertices])
                         md = checkfield(md,'fieldname','initialization.pressure','NaN',1,'size',[md.mesh.numberofvertices])
-                if (EnthalpyAnalysisEnum() in analyses and md.thermal.isenthalpy):
+                        if (EnthalpyAnalysisEnum() in analyses and md.thermal.isenthalpy):
                         md = checkfield(md,'fieldname','initialization.waterfraction','>=',0,'size',[md.mesh.numberofvertices])
                         md = checkfield(md,'fieldname','initialization.watercolumn'  ,'>=',0,'size',[md.mesh.numberofvertices])
                 if HydrologyShreveAnalysisEnum() in analyses:
                         md = checkfield(md,'fieldname','initialization.watercolumn','NaN',1,'size',[md.mesh.numberofvertices])
+                if HydrologyDCInefficientAnalysisEnum() in analyses:
+                        if hasattr(md.hydrology,'hydrologydc'):
+                                md = checkfield(md,'fieldname','initialization.sediment_head','NaN',1,'size',[md.mesh.numberofvertices 1]);
+                if HydrologyDCEfficientAnalysisEnum() in analyses:
+                        if hasattr(md.hydrology,'hydrologydc'):
+                                if md.hydrology.isefficientlayer==1,
+                                        md = checkfield(md,'fieldname','initialization.epl_head','NaN',1,'size',[md.mesh.numberofvertices 1]);
+                                        md = checkfield(md,'fieldname','initialization.epl_thickness','NaN',1,'size',[md.mesh.numberofvertices 1]);
+                                
+
+
                 return md
         # }}}
@@ -88 +102 @@
                 WriteData(fid,'data',self.watercolumn,'format','DoubleMat','mattype',1,'enum',WatercolumnEnum())
                 WriteData(fid,'data',self.sediment_head,'format','DoubleMat','mattype',1,'enum',SedimentHeadEnum())
+                WriteData(fid,'data',self.epl_head,'format','DoubleMat','mattype',1,'enum',EplHeadEnum);
+                WriteData(fid,'data',self.epl_thickness,'format','DoubleMat','mattype',1,'enum',HydrologydcEplThicknessEnum);
 
+                
                 if md.thermal.isenthalpy:
                         tpmp = md.materials.meltingpoint - md.materials.beta*md.initialization.pressure;

issm/trunk-jpl/src/m/classes/model.py ¶

@@ -667 +667 @@
                         md.hydrology.sediment_transmitivity=project3d(md,'vector',md.hydrology.sediment_transmitivity,'type','node','layer',1)
                         md.hydrology.basal_moulin_input=project3d(md,'vector',md.hydrology.basal_moulin_input,'type','node','layer',1)
+                if md.hydrology.isefficientlayer==1 :
+                        md.hydrology.spcepl_head=project3d(md,'vector',md.hydrology.spcepl_head,'type','node','layer',1)
 
                 #connectivity
@@ -760 +762 @@
                 if not numpy.isnan(md.initialization.vel).all(): md.initialization.vel=DepthAverage(md,md.initialization.vel)
                 if not numpy.isnan(md.initialization.temperature).all(): md.initialization.temperature=DepthAverage(md,md.initialization.temperature)
+                if not numpy.isnan(md.initialization.pressure).all(): md.initialization.pressure=project2d(md,md.initialization.pressure,1)
+                if not numpy.isnan(md.initialization.sediment_head).all(): md.initialization.sediment_head=project2d(md,md.initialization.sediment_head,1)
+                if not numpy.isnan(md.initialization.epl_head).all(): md.initialization.epl_head=project2d(md,md.initialization.epl_head,1)
+                if not numpy.isnan(md.initialization.epl_thickness).all(): md.initialization.epl_thickness=project2d(md,md.initialization.epl_thickness,1)
 
                 #gia
@@ -772 +778 @@
                         md.flowequation.borderHO=project2d(md,md.flowequation.borderHO,1)
                         md.flowequation.borderFS=project2d(md,md.flowequation.borderFS,1)
+
+
+                # Hydrologydc variables
+                if hasattr(md.hydrology,'hydrologydc'):
+                        md.hydrology.spcsediment_head=project2d(md,md.hydrology.spcsediment_head,1)
+                        md.hydrology.mask_eplactive_node=project2d(md,md.hydrology.mask_eplactive_node,1)
+                        md.hydrology.sediment_transmitivity=project2d(md,md.hydrology.sediment_transmitivity,1)
+                        md.hydrology.basal_moulin_input=project2d(md,md.hydrology.basal_moulin_input,1)
+                if md.hydrology.isefficientlayer == 1:
+                        md.hydrology.spcepl_head=project2d(md,md.hydrology.spcepl_head,1)
 
                 #boundary conditions

issm/trunk-jpl/src/m/exp/expdisp.m ¶

@@ -47 +47 @@
         if domain(i).nods==1
                 plot(domain(i).x*unitmultiplier,domain(i).y*unitmultiplier,'o','MarkerEdgeColor','k','MarkerFaceColor','r','MarkerSize',10);
+                text(domain(i).x*unitmultiplier,domain(i).y*unitmultiplier, ...
+                                 domain(i).name,'BackgroundColor',[1. .0 .0]);
         else
                 if (isnumeric(linestyle))
@@ -53 +55 @@
                         plot(domain(i).x*unitmultiplier,domain(i).y*unitmultiplier,linestyle,'linewidth',linewidth);
           end
+                text(domain(i).x(1)*unitmultiplier,domain(i).y(1)*unitmultiplier,domain(i).name,'BackgroundColor',[.7 .9 .7]);
   end
 end