Source code for issm.initialization

import numpy as np
from issm.project3d import project3d
from issm.fielddisplay import fielddisplay
from issm.checkfield import checkfield
from issm.WriteData import WriteData
import issm.MatlabFuncs as m

[docs]class initialization(object): """ INITIALIZATION class definition Usage: initialization=initialization(); """ def __init__(self): # {{{ self.vx = float('NaN') self.vy = float('NaN') self.vz = float('NaN') self.vel = float('NaN') self.pressure = float('NaN') self.temperature = float('NaN') self.waterfraction = float('NaN') self.watercolumn = float('NaN') self.sediment_head = float('NaN') self.epl_head = float('NaN') self.epl_thickness = float('NaN') #set defaults self.setdefaultparameters() #}}} 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]')) string="%s\n%s"%(string,fielddisplay(self,'vz','z component of velocity [m/yr]')) string="%s\n%s"%(string,fielddisplay(self,'vel','velocity norm [m/yr]')) string="%s\n%s"%(string,fielddisplay(self,'pressure','pressure [Pa]')) string="%s\n%s"%(string,fielddisplay(self,'temperature','temperature [K]')) string="%s\n%s"%(string,fielddisplay(self,'waterfraction','fraction of water in the ice')) 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 #}}}
[docs] def extrude(self,md): # {{{ self.vx=project3d(md,'vector',self.vx,'type','node') self.vy=project3d(md,'vector',self.vy,'type','node') self.vz=project3d(md,'vector',self.vz,'type','node') self.vel=project3d(md,'vector',self.vel,'type','node') self.temperature=project3d(md,'vector',self.temperature,'type','node') self.waterfraction=project3d(md,'vector',self.waterfraction,'type','node') self.watercolumn=project3d(md,'vector',self.watercolumn,'type','node') self.sediment_head=project3d(md,'vector',self.sediment_head,'type','node','layer',1) self.epl_head=project3d(md,'vector',self.epl_head,'type','node','layer',1) self.epl_thickness=project3d(md,'vector',self.epl_thickness,'type','node','layer',1) #Lithostatic pressure by default # self.pressure=md.constants.g*md.materials.rho_ice*(md.geometry.surface[:,0]-md.mesh.z) #self.pressure=md.constants.g*md.materials.rho_ice*(md.geometry.surface-md.mesh.z.reshape(-1,)) if np.ndim(md.geometry.surface)==2: print('Reshaping md.geometry.surface for you convenience but you should fix it in you files') self.pressure=md.constants.g*md.materials.rho_ice*(md.geometry.surface.reshape(-1,)-md.mesh.z) else: self.pressure=md.constants.g*md.materials.rho_ice*(md.geometry.surface-md.mesh.z) return self
#}}}
[docs] def setdefaultparameters(self): # {{{ return self
#}}}
[docs] def checkconsistency(self,md,solution,analyses): # {{{ if 'StressbalanceAnalysis' in analyses: if not np.any(np.logical_or(np.isnan(md.initialization.vx),np.isnan(md.initialization.vy))): md = checkfield(md,'fieldname','initialization.vx','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices]) md = checkfield(md,'fieldname','initialization.vy','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices]) if 'MasstransportAnalysis' in analyses: md = checkfield(md,'fieldname','initialization.vx','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices]) md = checkfield(md,'fieldname','initialization.vy','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices]) if 'BalancethicknessAnalysis' in analyses: md = checkfield(md,'fieldname','initialization.vx','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices]) md = checkfield(md,'fieldname','initialization.vy','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices]) #Triangle with zero velocity if np.any(np.logical_and(np.sum(np.abs(md.initialization.vx[md.mesh.elements-1]),axis=1)==0,\ np.sum(np.abs(md.initialization.vy[md.mesh.elements-1]),axis=1)==0)): md.checkmessage("at least one triangle has all its vertices with a zero velocity") if 'ThermalAnalysis' in analyses: md = checkfield(md,'fieldname','initialization.vx','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices]) md = checkfield(md,'fieldname','initialization.vy','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices]) md = checkfield(md,'fieldname','initialization.temperature','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices]) if md.mesh.dimension()==3: md = checkfield(md,'fieldname','initialization.vz','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices]) md = checkfield(md,'fieldname','initialization.pressure','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices]) if ('EnthalpyAnalysis' 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]) pos = np.nonzero(md.initialization.waterfraction > 0.)[0] if(pos.size): 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'); if 'HydrologyShreveAnalysis' in analyses: if hasattr(md.hydrology,'hydrologyshreve'): md = checkfield(md,'fieldname','initialization.watercolumn','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices]) if 'HydrologyDCInefficientAnalysis' in analyses: if hasattr(md.hydrology,'hydrologydc'): md = checkfield(md,'fieldname','initialization.sediment_head','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices]) if 'HydrologyDCEfficientAnalysis' in analyses: if hasattr(md.hydrology,'hydrologydc'): if md.hydrology.isefficientlayer==1: md = checkfield(md,'fieldname','initialization.epl_head','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices]) md = checkfield(md,'fieldname','initialization.epl_thickness','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices]) return md
# }}}
[docs] def marshall(self,prefix,md,fid): # {{{ yts=md.constants.yts WriteData(fid,prefix,'object',self,'fieldname','vx','format','DoubleMat','mattype',1,'scale',1./yts) WriteData(fid,prefix,'object',self,'fieldname','vy','format','DoubleMat','mattype',1,'scale',1./yts) WriteData(fid,prefix,'object',self,'fieldname','vz','format','DoubleMat','mattype',1,'scale',1./yts) WriteData(fid,prefix,'object',self,'fieldname','pressure','format','DoubleMat','mattype',1) WriteData(fid,prefix,'object',self,'fieldname','temperature','format','DoubleMat','mattype',1) WriteData(fid,prefix,'object',self,'fieldname','waterfraction','format','DoubleMat','mattype',1) WriteData(fid,prefix,'object',self,'fieldname','sediment_head','format','DoubleMat','mattype',1) WriteData(fid,prefix,'object',self,'fieldname','epl_head','format','DoubleMat','mattype',1) WriteData(fid,prefix,'object',self,'fieldname','epl_thickness','format','DoubleMat','mattype',1) WriteData(fid,prefix,'object',self,'fieldname','watercolumn','format','DoubleMat','mattype',1) if md.thermal.isenthalpy: tpmp = md.materials.meltingpoint - md.materials.beta*md.initialization.pressure; pos = np.nonzero(md.initialization.waterfraction > 0.)[0] enthalpy = md.materials.heatcapacity*(md.initialization.temperature-md.constants.referencetemperature); enthalpy[pos] = md.materials.heatcapacity*(tpmp[pos].reshape(-1,) - md.constants.referencetemperature) + md.materials.latentheat*md.initialization.waterfraction[pos].reshape(-1,) WriteData(fid,prefix,'data',enthalpy,'format','DoubleMat','mattype',1,'name','md.initialization.enthalpy');
# }}}