![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 23833 for issm/trunk-jpl/src/m/classes/initialization.py
- Timestamp:
- 04/08/19 16:35:05 (6 years ago)
- File:
-
- 1 edited
-
issm/trunk-jpl/src/m/classes/initialization.py (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/src/m/classes/initialization.py
r23812 r23833 7 7 8 8 class initialization(object): 9 10 11 12 13 14 9 """ 10 INITIALIZATION class definition 11 12 Usage: 13 initialization=initialization(); 14 """ 15 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 16 def __init__(self): # {{{ 17 18 self.vx = float('NaN') 19 self.vy = float('NaN') 20 self.vz = float('NaN') 21 self.vel = float('NaN') 22 self.enthalpy = float('NaN') 23 self.pressure = float('NaN') 24 self.temperature = float('NaN') 25 self.waterfraction = float('NaN') 26 self.watercolumn = float('NaN') 27 self.sediment_head = float('NaN') 28 self.epl_head = float('NaN') 29 self.epl_thickness = float('NaN') 30 30 31 32 31 #set defaults 32 self.setdefaultparameters() 33 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 34 #}}} 35 def __repr__(self): # {{{ 36 string=' initial field values:' 37 string="%s\n%s"%(string,fielddisplay(self,'vx','x component of velocity [m/yr]')) 38 string="%s\n%s"%(string,fielddisplay(self,'vy','y component of velocity [m/yr]')) 39 string="%s\n%s"%(string,fielddisplay(self,'vz','z component of velocity [m/yr]')) 40 string="%s\n%s"%(string,fielddisplay(self,'vel','velocity norm [m/yr]')) 41 string="%s\n%s"%(string,fielddisplay(self,'pressure','pressure [Pa]')) 42 string="%s\n%s"%(string,fielddisplay(self,'temperature','temperature [K]')) 43 string="%s\n%s"%(string,fielddisplay(self,'enthalpy','enthalpy [J]')) 44 string="%s\n%s"%(string,fielddisplay(self,'waterfraction','fraction of water in the ice')) 45 string="%s\n%s"%(string,fielddisplay(self,'watercolumn','thickness of subglacial water [m]')) 46 string="%s\n%s"%(string,fielddisplay(self,'sediment_head','sediment water head of subglacial system [m]')) 47 string="%s\n%s"%(string,fielddisplay(self,'epl_head','epl water head of subglacial system [m]')) 48 string="%s\n%s"%(string,fielddisplay(self,'epl_thickness','thickness of the epl [m]')) 49 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 50 return string 51 #}}} 52 def extrude(self,md): # {{{ 53 self.vx=project3d(md,'vector',self.vx,'type','node') 54 self.vy=project3d(md,'vector',self.vy,'type','node') 55 self.vz=project3d(md,'vector',self.vz,'type','node') 56 self.vel=project3d(md,'vector',self.vel,'type','node') 57 self.temperature=project3d(md,'vector',self.temperature,'type','node') 58 self.enthalpy=project3d(md,'vector',self.enthalpy,'type','node') 59 self.waterfraction=project3d(md,'vector',self.waterfraction,'type','node') 60 self.watercolumn=project3d(md,'vector',self.watercolumn,'type','node') 61 self.sediment_head=project3d(md,'vector',self.sediment_head,'type','node','layer',1) 62 self.epl_head=project3d(md,'vector',self.epl_head,'type','node','layer',1) 63 self.epl_thickness=project3d(md,'vector',self.epl_thickness,'type','node','layer',1) 64 64 65 66 #self.pressure=md.constants.g*md.materials.rho_ice*(md.geometry.surface[:,0]-md.mesh.z)67 65 #Lithostatic pressure by default 66 # self.pressure=md.constants.g*md.materials.rho_ice*(md.geometry.surface[:,0]-md.mesh.z) 67 #self.pressure=md.constants.g*md.materials.rho_ice*(md.geometry.surface-md.mesh.z.reshape(-1,)) 68 68 69 70 71 72 73 69 if np.ndim(md.geometry.surface)==2: 70 print('Reshaping md.geometry.surface for you convenience but you should fix it in you files') 71 self.pressure=md.constants.g*md.materials.rho_ice*(md.geometry.surface.reshape(-1,)-md.mesh.z) 72 else: 73 self.pressure=md.constants.g*md.materials.rho_ice*(md.geometry.surface-md.mesh.z) 74 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 md = checkfield(md,'fieldname', 'delta Tpmp', 'field', np.absolute(md.initialization.temperature[pos]-(md.materials.meltingpoint-md.materials.beta*md.initialization.pressure[pos])),'<',1e-11,'message','set temperature to pressure melting point at locations with waterfraction>0');108 109 110 111 112 113 114 115 116 117 118 75 return self 76 #}}} 77 def setdefaultparameters(self): # {{{ 78 return self 79 #}}} 80 def checkconsistency(self,md,solution,analyses): # {{{ 81 if 'StressbalanceAnalysis' in analyses: 82 if not np.any(np.logical_or(np.isnan(md.initialization.vx),np.isnan(md.initialization.vy))): 83 md = checkfield(md,'fieldname','initialization.vx','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices]) 84 md = checkfield(md,'fieldname','initialization.vy','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices]) 85 if 'MasstransportAnalysis' in analyses: 86 md = checkfield(md,'fieldname','initialization.vx','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices]) 87 md = checkfield(md,'fieldname','initialization.vy','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices]) 88 if 'BalancethicknessAnalysis' in analyses: 89 md = checkfield(md,'fieldname','initialization.vx','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices]) 90 md = checkfield(md,'fieldname','initialization.vy','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices]) 91 #Triangle with zero velocity 92 if np.any(np.logical_and(np.sum(np.abs(md.initialization.vx[md.mesh.elements-1]),axis=1)==0,\ 93 np.sum(np.abs(md.initialization.vy[md.mesh.elements-1]),axis=1)==0)): 94 md.checkmessage("at least one triangle has all its vertices with a zero velocity") 95 if 'ThermalAnalysis' in analyses: 96 md = checkfield(md,'fieldname','initialization.vx','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices]) 97 md = checkfield(md,'fieldname','initialization.vy','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices]) 98 md = checkfield(md,'fieldname','initialization.temperature','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices]) 99 if md.mesh.dimension()==3: 100 md = checkfield(md,'fieldname','initialization.vz','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices]) 101 md = checkfield(md,'fieldname','initialization.pressure','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices]) 102 if ('EnthalpyAnalysis' in analyses and md.thermal.isenthalpy): 103 md = checkfield(md,'fieldname','initialization.waterfraction','>=',0,'size',[md.mesh.numberofvertices]) 104 md = checkfield(md,'fieldname','initialization.watercolumn' ,'>=',0,'size',[md.mesh.numberofvertices]) 105 pos = np.nonzero(md.initialization.waterfraction > 0.)[0] 106 if(pos.size): 107 md = checkfield(md,'fieldname', 'delta Tpmp', 'field', np.absolute(md.initialization.temperature[pos]-(md.materials.meltingpoint-md.materials.beta*md.initialization.pressure[pos])),'<',1e-11, 'message','set temperature to pressure melting point at locations with waterfraction>0'); 108 if 'HydrologyShreveAnalysis' in analyses: 109 if hasattr(md.hydrology,'hydrologyshreve'): 110 md = checkfield(md,'fieldname','initialization.watercolumn','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices]) 111 if 'HydrologyDCInefficientAnalysis' in analyses: 112 if hasattr(md.hydrology,'hydrologydc'): 113 md = checkfield(md,'fieldname','initialization.sediment_head','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices]) 114 if 'HydrologyDCEfficientAnalysis' in analyses: 115 if hasattr(md.hydrology,'hydrologydc'): 116 if md.hydrology.isefficientlayer==1: 117 md = checkfield(md,'fieldname','initialization.epl_head','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices]) 118 md = checkfield(md,'fieldname','initialization.epl_thickness','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices]) 119 119 120 121 122 120 return md 121 # }}} 122 def marshall(self,prefix,md,fid): # {{{ 123 123 124 124 yts=md.constants.yts 125 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 126 WriteData(fid,prefix,'object',self,'fieldname','vx','format','DoubleMat','mattype',1,'scale',1./yts) 127 WriteData(fid,prefix,'object',self,'fieldname','vy','format','DoubleMat','mattype',1,'scale',1./yts) 128 WriteData(fid,prefix,'object',self,'fieldname','vz','format','DoubleMat','mattype',1,'scale',1./yts) 129 WriteData(fid,prefix,'object',self,'fieldname','pressure','format','DoubleMat','mattype',1) 130 WriteData(fid,prefix,'object',self,'fieldname','temperature','format','DoubleMat','mattype',1) 131 WriteData(fid,prefix,'object',self,'fieldname','waterfraction','format','DoubleMat','mattype',1) 132 WriteData(fid,prefix,'object',self,'fieldname','sediment_head','format','DoubleMat','mattype',1) 133 WriteData(fid,prefix,'object',self,'fieldname','epl_head','format','DoubleMat','mattype',1) 134 WriteData(fid,prefix,'object',self,'fieldname','epl_thickness','format','DoubleMat','mattype',1) 135 WriteData(fid,prefix,'object',self,'fieldname','watercolumn','format','DoubleMat','mattype',1) 136 137 if md.thermal.isenthalpy: 138 if (np.size(self.enthalpy)<=1): 139 tpmp = md.materials.meltingpoint - md.materials.beta*md.initialization.pressure; 140 pos = np.nonzero(md.initialization.waterfraction > 0.)[0] 141 self.enthalpy = md.materials.heatcapacity*(md.initialization.temperature-md.constants.referencetemperature); 142 self.enthalpy[pos] = md.materials.heatcapacity*(tpmp[pos].reshape(-1,) - md.constants.referencetemperature) + md.materials.latentheat*md.initialization.waterfraction[pos].reshape(-1,) 143 143 144 144 WriteData(fid,prefix,'data',self.enthalpy,'format','DoubleMat','mattype',1,'name','md.initialization.enthalpy'); 145 145 146 146 # }}}
Note:
See TracChangeset
for help on using the changeset viewer.
