Index: /issm/trunk-jpl/src/m/classes/SMBgemb.m =================================================================== --- /issm/trunk-jpl/src/m/classes/SMBgemb.m (revision 22266) +++ /issm/trunk-jpl/src/m/classes/SMBgemb.m (revision 22267) @@ -358,5 +358,5 @@ WriteData(fid,prefix,'data',dt,'name','md.smb.dt','format','Double','scale',yts); - + % Check if smb_dt goes evenly into transient core time step if (mod(md.timestepping.time_step,dt) >= 1e-10) Index: /issm/trunk-jpl/src/m/classes/SMBgemb.py =================================================================== --- /issm/trunk-jpl/src/m/classes/SMBgemb.py (revision 22267) +++ /issm/trunk-jpl/src/m/classes/SMBgemb.py (revision 22267) @@ -0,0 +1,368 @@ +import numpy as np +from fielddisplay import fielddisplay +from checkfield import checkfield +from WriteData import WriteData +from project3d import project3d + +class SMBgemb(object): + """ + SMBgemb Class definition + + Usage: + SMB = SMBgemb() + """ + + def __init__(self): # {{{ + #each one of these properties is a transient forcing to the GEMB model, loaded from meteorological data derived + #from an automatic weather stations (AWS). Each property is therefore a matrix, of size (numberofvertices x number + #of time steps. ) + + #solution choices + #check these: + #isgraingrowth + #isalbedo + #isshortwave + #isthermal + #isaccumulation + #ismelt + #isdensification + #isturbulentflux + + #inputs: + Ta = float('NaN') #2 m air temperature, in Kelvin + V = float('NaN') #wind speed (m/s-1) + dswrf = float('NaN') #downward shortwave radiation flux [W/m^2] + dlwrf = float('NaN') #downward longwave radiation flux [W/m^2] + P = float('NaN') #precipitation [mm w.e. / m^2] + eAir = float('NaN') #screen level vapor pressure [Pa] + pAir = float('NaN') #surface pressure [Pa] + + Tmean = float('NaN') #mean annual temperature [K] + C = float('NaN') #mean annual snow accumulation [kg m-2 yr-1] + Tz = float('NaN') #height above ground at which temperature (T) was sampled [m] + Vz = float('NaN') #height above ground at which wind (V) eas sampled [m] + + # Initialization of snow properties + Dzini = float('NaN') #cell depth (m) + Dini = float('NaN') #snow density (kg m-3) + Reini = float('NaN') #effective grain size (mm) + Gdnini = float('NaN') #grain dricity (0-1) + Gspini = float('NaN') #grain sphericity (0-1) + ECini = float('NaN') #evaporation/condensation (kg m-2) + Wini = float('NaN') #Water content (kg m-2) + Aini = float('NaN') #albedo (0-1) + Tini = float('NaN') #snow temperature (K) + Sizeini = float('NaN') #Number of layers + + #settings: + aIdx = float('NaN') #method for calculating albedo and subsurface absorption (default is 1) + # 1: effective grain radius [Gardner & Sharp, 2009] + # 2: effective grain radius [Brun et al., 2009] + # 3: density and cloud amount [Greuell & Konzelmann, 1994] + # 4: exponential time decay & wetness [Bougamont & Bamber, 2005] + swIdx = float('NaN') # apply all SW to top grid cell (0) or allow SW to penetrate surface (1) (default 1) + + denIdx = float('NaN') #densification model to use (default is 2): + # 1 = emperical model of Herron and Langway (1980) + # 2 = semi-emerical model of Anthern et al. (2010) + # 3 = DO NOT USE: physical model from Appix B of Anthern et al. (2010) + # 4 = DO NOT USE: emperical model of Li and Zwally (2004) + # 5 = DO NOT USE: modified emperical model (4) by Helsen et al. (2008) + + zTop = float('NaN') # depth over which grid length is constant at the top of the snopack (default 10) [m] + dzTop = float('NaN') # initial top vertical grid spacing (default .05) [m] + dzMin = float('NaN') # initial min vertical allowable grid spacing (default dzMin/2) [m] + zY = float('NaN') # strech grid cells bellow top_z by a [top_dz * y ^ (cells bellow top_z)] + zMax = float('NaN') #initial max model depth (default is min(thickness,500)) [m] + zMin = float('NaN') #initial min model depth (default is min(thickness,30)) [m] + outputFreq = float('NaN') #output frequency in days (default is monthly, 30) + + #specific albedo parameters: + #Method 1 and 2: + aSnow = float('NaN') # new snow albedo (0.64 - 0.89) + aIce = float('NaN') # range 0.27-0.58 for old snow + #Method 3: Radiation Correction Factors -> only used for met station data and Greuell & Konzelmann, 1994 albedo + cldFrac = float('NaN') # average cloud amount + #Method 4: additonal tuning parameters albedo as a funtion of age and water content (Bougamont et al., 2005) + t0wet = float('NaN') # time scale for wet snow (15-21.9) + t0dry = float('NaN') # warm snow timescale (30) + K = float('NaN') # time scale temperature coef. (7) + + #densities: + InitDensityScaling = float('NaN') #initial scaling factor multiplying the density of ice, which describes the density of the snowpack. + + requested_outputs = [] + + #Several fields are missing from the standard GEMB model, which are capture intrinsically by ISSM. + #dateN: that's the last row of the above fields. + #dt: included in dateN. Not an input. + #elev: this is taken from the ISSM surface itself. + + #}}} + def __repr__(self): # {{{ + #string = " surface forcings parameters:" + #string = "#s\n#s"%(string,fielddisplay(self,'mass_balance','surface mass balance [m/yr ice eq]')) + #string = "#s\n#s"%(string,fielddisplay(self,'requested_outputs','additional outputs requested')) + string = ' surface forcings for SMB GEMB model :' + + string="%s\n%s"%(string,fielddisplay(self,'issmbgradients','is smb gradients method activated (0 or 1, default is 0)')) + string = "%s\n%s"%(string,fielddisplay(self,'isgraingrowth','run grain growth module (default true)')) + string = "%s\n%s"%(string,fielddisplay(self,'isalbedo','run albedo module (default true)')) + string = "%s\n%s"%(string,fielddisplay(self,'isshortwave','run short wave module (default true)')) + string = "%s\n%s"%(string,fielddisplay(self,'isthermal','run thermal module (default true)')) + string = "%s\n%s"%(string,fielddisplay(self,'isaccumulation','run accumulation module (default true)')) + string = "%s\n%s"%(string,fielddisplay(self,'ismelt','run melting module (default true)')) + string = "%s\n%s"%(string,fielddisplay(self,'isdensification','run densification module (default true)')) + string = "%s\n%s"%(string,fielddisplay(self,'isturbulentflux','run turbulant heat fluxes module (default true)')) + string = "%s\n%s"%(string,fielddisplay(self,'Ta','2 m air temperature, in Kelvin')) + string = "%s\n%s"%(string,fielddisplay(self,'V','wind speed (m/s-1)')) + string = "%s\n%s"%(string,fielddisplay(self,'dlwrf','downward shortwave radiation flux [W/m^2]')) + string = "%s\n%s"%(string,fielddisplay(self,'dswrf','downward longwave radiation flux [W/m^2]')) + string = "%s\n%s"%(string,fielddisplay(self,'P','precipitation [mm w.e. / m^2]')) + string = "%s\n%s"%(string,fielddisplay(self,'eAir','screen level vapor pressure [Pa]')) + string = "%s\n%s"%(string,fielddisplay(self,'pAir','surface pressure [Pa]')) + string = "%s\n%s"%(string,fielddisplay(self,'Tmean','mean annual temperature [K]')) + string = "%s\n%s"%(string,fielddisplay(self,'C','mean annual snow accumulation [kg m-2 yr-1]')) + string = "%s\n%s"%(string,fielddisplay(self,'Tz','height above ground at which temperature (T) was sampled [m]')) + string = "%s\n%s"%(string,fielddisplay(self,'Vz','height above ground at which wind (V) eas sampled [m]')) + string = "%s\n%s"%(string,fielddisplay(self,'zTop','depth over which grid length is constant at the top of the snopack (default 10) [m]')) + string = "%s\n%s"%(string,fielddisplay(self,'dzTop','initial top vertical grid spacing (default .05) [m] ')) + string = "%s\n%s"%(string,fielddisplay(self,'dzMin','initial min vertical allowable grid spacing (default dzMin/2) [m] ')) + string = "%s\n%s"%(string,fielddisplay(self,'zMax','initial max model depth (default is min(thickness,500)) [m]')) + string = "%s\n%s"%(string,fielddisplay(self,'zMin','initial min model depth (default is min(thickness,30)) [m]')) + string = "%s\n%s"%(string,fielddisplay(self,'zY','strech grid cells bellow top_z by a [top_dz * y ^ (cells bellow top_z)]')) + string = "%s\n%s"%(string,fielddisplay(self,'InitDensityScaling',['initial scaling factor multiplying the density of ice','which describes the density of the snowpack.'])) + string = "%s\n%s"%(string,fielddisplay(self,'outputFreq','output frequency in days (default is monthly, 30)')) + string = "%s\n%s"%(string,fielddisplay(self,'aIdx',['method for calculating albedo and subsurface absorption (default is 1)', + '1: effective grain radius [Gardner & Sharp, 2009]', + '2: effective grain radius [Brun et al., 2009]', + '3: density and cloud amount [Greuell & Konzelmann, 1994]', + '4: exponential time decay & wetness [Bougamont & Bamber, 2005]'])) + + #snow properties init + string = "%s\n%s"%(string,fielddisplay(self,'Dzini','Initial cel depth when restart [m]')) + string = "%s\n%s"%(string,fielddisplay(self,'Dini','Initial snow density when restart [kg m-3]')) + string = "%s\n%s"%(string,fielddisplay(self,'Reini','Initial grain size when restart [mm]')) + string = "%s\n%s"%(string,fielddisplay(self,'Gdnini','Initial grain dricity when restart [-]')) + string = "%s\n%s"%(string,fielddisplay(self,'Gspini','Initial grain sphericity when restart [-]')) + string = "%s\n%s"%(string,fielddisplay(self,'ECini','Initial evaporation/condensation when restart [kg m-2]')) + string = "%s\n%s"%(string,fielddisplay(self,'Wini','Initial snow water content when restart [kg m-2]')) + string = "%s\n%s"%(string,fielddisplay(self,'Aini','Initial albedo when restart [-]')) + string = "%s\n%s"%(string,fielddisplay(self,'Tini','Initial snow temperature when restart [K]')) + string = "%s\n%s"%(string,fielddisplay(self,'Sizeini','Initial number of layers when restart [K]')) + + #additional albedo parameters: + if type(self.aIdx) == list or type(self.aIdx) == type(np.array([1,2])) and (self.aIdx == [1,2] or (1 in self.aIdx and 2 in self.aIdx)): + string = "%s\n%s"%(string,fielddisplay(self,'aSnow','new snow albedo (0.64 - 0.89)')) + string = "%s\n%s"%(string,fielddisplay(self,'aIce','albedo of ice (0.27-0.58)')) + elif self.aIdx == 3: + string = "%s\n%s"%(string,fielddisplay(self,'cldFrac','average cloud amount')) + elif self.aIdx == 4: + string = "%s\n%s"%(string,fielddisplay(self,'t0wet','time scale for wet snow (15-21.9) [d]')) + string = "%s\n%s"%(string,fielddisplay(self,'t0dry','warm snow timescale (30) [d]')) + string = "%s\n%s"%(string,fielddisplay(self,'K','time scale temperature coef. (7) [d]')) + + string = "%s\n%s"%(string,fielddisplay(self,'swIdx','apply all SW to top grid cell (0) or allow SW to penetrate surface (1) [default 1]')) + string = "%s\n%s"%(string,fielddisplay(self,'denIdx',['densification model to use (default is 2):', + '1 = emperical model of Herron and Langway (1980)', + '2 = semi-emerical model of Anthern et al. (2010)', + '3 = DO NOT USE: physical model from Appix B of Anthern et al. (2010)', + '4 = DO NOT USE: emperical model of Li and Zwally (2004)', + '5 = DO NOT USE: modified emperical model (4) by Helsen et al. (2008)'])) + string = "%s\n%s"%(string,fielddisplay(self,'requested_outputs','additional outputs requested')) + return string + #}}} + + def extrude(self,md): # {{{ + + self.Ta = project3d(md,'vector',self.Ta,'type','node') + self.V = project3d(md,'vector',self.V,'type','node') + self.dswrf = project3d(md,'vector',self.dswrf,'type','node') + self.dswrf = project3d(md,'vector',self.dswrf,'type','node') + self.P = project3d(md,'vector',self.P,'type','node') + self.eAir = project3d(md,'vector',self.eAir,'type','node') + self.pAir = project3d(md,'vector',self.pAir,'type','node') + return self + #}}} + def defaultoutputs(self,md): # {{{ + return ['SmbMassBalance'] + #}}} + + def setdefaultparameters(self,mesh,geometry): # {{{ + self.isgraingrowth = 1 + self.isalbedo = 1 + self.isshortwave = 1 + self.isthermal = 1 + self.isaccumulation = 1 + self.ismelt = 1 + self.isdensification = 1 + self.isturbulentflux = 1 + + self.aIdx = 1 + self.swIdx = 1 + self.denIdx = 2 + self.zTop = 10*np.ones((mesh.numberofelements,)) + self.dzTop = .05* np.ones((mesh.numberofelements,)) + self.dzMin = self.dzTop/2 + self.InitDensityScaling = 1.0 + + self.zMax = 250*np.ones((mesh.numberofelements,)) + self.zMin = 130*np.ones((mesh.numberofelements,)) + self.zY = 1.10*np.ones((mesh.numberofelements,)) + self.outputFreq = 30 + + #additional albedo parameters + self.aSnow = 0.85 + self.aIce = 0.48 + self.cldFrac = 0.1 + self.t0wet = 15 + self.t0dry = 30 + self.K = 7 + + self.Dzini = 0.05*np.ones((mesh.numberofelements,2)) + self.Dini = 910.0*np.ones((mesh.numberofelements,2)) + self.Reini = 2.5*np.ones((mesh.numberofelements,2)) + self.Gdnini = 0.0*np.ones((mesh.numberofelements,2)) + self.Gspini = 0.0*np.ones((mesh.numberofelements,2)) + self.ECini = 0.0*np.ones((mesh.numberofelements,)) + self.Wini = 0.0*np.ones((mesh.numberofelements,2)) + self.Aini = self.aSnow*np.ones((mesh.numberofelements,2)) + self.Tini = 273.15*np.ones((mesh.numberofelements,2)) +# /!\ Default value of Tini must be equal to Tmean but don't know Tmean yet (computed when atmospheric forcings are interpolated on mesh). +# If initialization without restart, this value will be overwritten when snow parameters are retrieved in Element.cpp + self.Sizeini = 2*np.ones((mesh.numberofelements,)) + #}}} + + def checkconsistency(self,md,solution,analyses): # {{{ + + md = checkfield(md,'fieldname','smb.isgraingrowth','values',[0,1]) + md = checkfield(md,'fieldname','smb.isalbedo','values',[0,1]) + md = checkfield(md,'fieldname','smb.isshortwave','values',[0,1]) + md = checkfield(md,'fieldname','smb.isthermal','values',[0,1]) + md = checkfield(md,'fieldname','smb.isaccumulation','values',[0,1]) + md = checkfield(md,'fieldname','smb.ismelt','values',[0,1]) + md = checkfield(md,'fieldname','smb.isdensification','values',[0,1]) + md = checkfield(md,'fieldname','smb.isturbulentflux','values',[0,1]) + + md = checkfield(md,'fieldname','smb.Ta','timeseries',1,'NaN',1,'Inf',1,'>',273-100,'<',273+100) #-100/100 celsius min/max value + md = checkfield(md,'fieldname','smb.V','timeseries',1,'NaN',1,'Inf',1,'> = ',0,'<',45) #max 500 km/h + md = checkfield(md,'fieldname','smb.dswrf','timeseries',1,'NaN',1,'Inf',1,'> = ',0,'< = ',1400) + md = checkfield(md,'fieldname','smb.dlwrf','timeseries',1,'NaN',1,'Inf',1,'> = ',0) + md = checkfield(md,'fieldname','smb.P','timeseries',1,'NaN',1,'Inf',1,'> = ',0,'< = ',100) + md = checkfield(md,'fieldname','smb.eAir','timeseries',1,'NaN',1,'Inf',1) + + md = checkfield(md,'fieldname','smb.Tmean','size',[md.mesh.numberofelements],'NaN',1,'Inf',1,'>',273-100,'<',273+100) #-100/100 celsius min/max value + md = checkfield(md,'fieldname','smb.C','size',[md.mesh.numberofelements],'NaN',1,'Inf',1,'> = ',0) + md = checkfield(md,'fieldname','smb.Tz','size',[md.mesh.numberofelements],'NaN',1,'Inf',1,'> = ',0,'< = ',5000) + md = checkfield(md,'fieldname','smb.Vz','size',[md.mesh.numberofelements],'NaN',1,'Inf',1,'> = ',0,'< = ',5000) + + md = checkfield(md,'fieldname','smb.aIdx','NaN',1,'Inf',1,'values',[1,2,3,4]) + md = checkfield(md,'fieldname','smb.swIdx','NaN',1,'Inf',1,'values',[0,1]) + md = checkfield(md,'fieldname','smb.denIdx','NaN',1,'Inf',1,'values',[1,2,3,4,5]) + + md = checkfield(md,'fieldname','smb.zTop','NaN',1,'Inf',1,'> = ',0) + md = checkfield(md,'fieldname','smb.dzTop','NaN',1,'Inf',1,'>',0) + md = checkfield(md,'fieldname','smb.dzMin','NaN',1,'Inf',1,'>',0) + md = checkfield(md,'fieldname','smb.zY','NaN',1,'Inf',1,'> = ',1) + md = checkfield(md,'fieldname','smb.outputFreq','NaN',1,'Inf',1,'>',0,'<',10*365) #10 years max + md = checkfield(md,'fieldname','smb.InitDensityScaling','NaN',1,'Inf',1,'> = ',0,'< = ',1) + + if type(self.aIdx) == list or type(self.aIdx) == type(np.array([1,2])) and (self.aIdx == [1,2] or (1 in self.aIdx and 2 in self.aIdx)): + md = checkfield(md,'fieldname','smb.aSnow','NaN',1,'Inf',1,'> = ',.64,'< = ',.89) + md = checkfield(md,'fieldname','smb.aIce','NaN',1,'Inf',1,'> = ',.27,'< = ',.58) + elif self.aIdx == 3: + md = checkfield(md,'fieldname','smb.cldFrac','NaN',1,'Inf',1,'> = ',0,'< = ',1) + elif self.aIdx == 4: + md = checkfield(md,'fieldname','smb.t0wet','NaN',1,'Inf',1,'> = ',15,'< = ',21.9) + md = checkfield(md,'fieldname','smb.t0dry','NaN',1,'Inf',1,'> = ',30,'< = ',30) + md = checkfield(md,'fieldname','smb.K','NaN',1,'Inf',1,'> = ',7,'< = ',7) + + + #check zTop is < local thickness: + he = np.sum(md.geometry.thickness[md.mesh.elements-1],axis=1)/np.size(md.mesh.elements,1) + if np.any(he= 1e-10): + error('smb_dt/dt = #f. The number of SMB time steps in one transient core time step has to be an an integer',md.timestepping.time_step/dt) + + #process requested outputs + outputs = self.requested_outputs + indices = [i for i, x in enumerate(outputs) if x == 'default'] + if len(indices) > 0: + outputscopy=outputs[0:max(0,indices[0]-1)]+self.defaultoutputs(md)+outputs[indices[0]+1:] + outputs =outputscopy + + WriteData(fid,prefix,'data',outputs,'name','md.smb.requested_outputs','format','StringArray') + # }}} + Index: /issm/trunk-jpl/src/m/classes/SMBgradientsela.py =================================================================== --- /issm/trunk-jpl/src/m/classes/SMBgradientsela.py (revision 22266) +++ /issm/trunk-jpl/src/m/classes/SMBgradientsela.py (revision 22267) @@ -9,5 +9,5 @@ Usage: - SMBgradientsela=SMBgradientsela(); + SMBgradientsela=SMBgradientsela() """ @@ -19,9 +19,10 @@ self.b_min = float('NaN') self.requested_outputs = [] + self.setdefaultparameters() #}}} def __repr__(self): # {{{ - string=" surface forcings parameters:" + string = " surface forcings parameters:" + string+= '\n SMB gradients ela parameters:' - string="%s\n%s"%(string,fielddisplay(self,'issmbgradientsela','is smb gradients ela method activated (0 or 1, default is 0)')) string="%s\n%s"%(string,fielddisplay(self,'ela',' equilibrium line altitude from which deviation is used to calculate smb using the smb gradients ela method [m a.s.l.]')) string="%s\n%s"%(string,fielddisplay(self,'b_pos',' vertical smb gradient (dB/dz) above ela')) @@ -47,4 +48,9 @@ return self #}}} + def setdefaultparameters(self): # {{{ + self.b_max=9999. + self.b_min=-9999. + return self + #}}} def checkconsistency(self,md,solution,analyses): # {{{ @@ -56,5 +62,5 @@ md = checkfield(md,'fieldname','smb.b_min','timeseries',1,'NaN',1,'Inf',1) - md = checkfield(md,'fieldname','masstransport.requested_outputs','stringrow',1) + md = checkfield(md,'fieldname','smb.requested_outputs','stringrow',1) return md # }}} Index: /issm/trunk-jpl/src/m/classes/amr.py =================================================================== --- /issm/trunk-jpl/src/m/classes/amr.py (revision 22266) +++ /issm/trunk-jpl/src/m/classes/amr.py (revision 22267) @@ -28,4 +28,5 @@ self.deviatoricerror_threshold = 0. self.deviatoricerror_maximum = 0. + self.restart=0. #set defaults self.setdefaultparameters() @@ -48,4 +49,5 @@ string="%s\n%s"%(string,fielddisplay(self,"deviatoricerror_threshold","maximum threshold deviatoricstress error permitted")) string="%s\n%s"%(string,fielddisplay(self,"deviatoricerror_maximum","maximum deviatoricstress error permitted")) + string="%s\n%s"%(string,fielddisplay(self,'restart',['indicates if ReMesh() will call before first time step'])) return string #}}} @@ -67,4 +69,5 @@ self.deviatoricerror_threshold = 0 self.deviatoricerror_maximum = 0 + self.restart = 0. return self #}}} @@ -84,4 +87,5 @@ md = checkfield(md,'fieldname','amr.deviatoricerror_threshold','numel',[1],'>=',0,'<=',1,'NaN',1); md = checkfield(md,'fieldname','amr.deviatoricerror_maximum','numel',[1],'>=',0,'NaN',1,'Inf',1); + md = checkfield(md,'fieldname','amr.restart','numel',[1],'>=',0,'<=',1,'NaN',1) return md # }}} @@ -104,3 +108,4 @@ WriteData(fid,prefix,'object',self,'fieldname','deviatoricerror_threshold','format','Double'); WriteData(fid,prefix,'object',self,'fieldname','deviatoricerror_maximum','format','Double'); + WriteData(fid,prefix,'object',self,'class','amr','fieldname','restart','format','Integer') # }}} Index: /issm/trunk-jpl/src/m/classes/basalforcings.py =================================================================== --- /issm/trunk-jpl/src/m/classes/basalforcings.py (revision 22266) +++ /issm/trunk-jpl/src/m/classes/basalforcings.py (revision 22267) @@ -45,4 +45,7 @@ self.floatingice_melting_rate=np.zeros((md.mesh.numberofvertices)) print " no basalforcings.floatingice_melting_rate specified: values set as zero" + #if np.all(np.isnan(self.geothermalflux)): + #self.geothermalflux=np.zeros((md.mesh.numberofvertices)) + #print " no basalforcings.geothermalflux specified: values set as zero" return self Index: /issm/trunk-jpl/src/m/classes/calvingdev.py =================================================================== --- /issm/trunk-jpl/src/m/classes/calvingdev.py (revision 22267) +++ /issm/trunk-jpl/src/m/classes/calvingdev.py (revision 22267) @@ -0,0 +1,60 @@ +from fielddisplay import fielddisplay +from project3d import project3d +from checkfield import checkfield +from WriteData import WriteData + +class calvingdev(object): + """ + CALVINGDEV class definition + + Usage: + calvingdev=calvingdev(); + """ + + def __init__(self): # {{{ + + self.stress_threshold_groundedice = 0. + self.stress_threshold_floatingice = 0. + self.meltingrate = float('NaN') + + #set defaults + self.setdefaultparameters() + + #}}} + def __repr__(self): # {{{ + string=' Calving Pi parameters:' + string="%s\n%s"%(string,fielddisplay(self,'stress_threshold_groundedice','sigma_max applied to grounded ice only [Pa]')) + string="%s\n%s"%(string,fielddisplay(self,'stress_threshold_floatingice','sigma_max applied to floating ice only [Pa]')) + + string="%s\n%s"%(string,fielddisplay(self,'meltingrate','melting rate at given location [m/a]')) + return string + #}}} + def extrude(self,md): # {{{ + self.meltingrate=project3d(md,'vector',self.meltingrate,'type','node') + return self + #}}} + def setdefaultparameters(self): # {{{ + #Default sigma max + self.stress_threshold_groundedice = 1e6 + self.stress_threshold_floatingice = 150e3 + return self + #}}} + def checkconsistency(self,md,solution,analyses): # {{{ + #Early return + if solution == 'TransientSolution' or md.transient.ismovingfront == 0: + return + + md = checkfield(md,'fieldname','calving.stress_threshold_groundedice','>',0,'nan',1,'Inf',1) + md = checkfield(md,'fieldname','calving.stress_threshold_floatingice','>',0,'nan',1,'Inf',1) + md = checkfield(md,'fieldname','calving.meltingrate','NaN',1,'Inf',1,'timeseries',1,'>=',0) + + return md + # }}} + def marshall(self,prefix,md,fid): # {{{ + yts=md.constants.yts + + WriteData(fid,prefix,'name','md.calving.law','data',2,'format','Integer') + WriteData(fid,prefix,'object',self,'fieldname','stress_threshold_groundedice','format','DoubleMat','mattype',1) + WriteData(fid,prefix,'object',self,'fieldname','stress_threshold_floatingice','format','DoubleMat','mattype',1) + WriteData(fid,prefix,'object',self,'fieldname','meltingrate','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts,'scale',1./yts) + # }}} Index: /issm/trunk-jpl/src/m/classes/calvingminthickness.py =================================================================== --- /issm/trunk-jpl/src/m/classes/calvingminthickness.py (revision 22267) +++ /issm/trunk-jpl/src/m/classes/calvingminthickness.py (revision 22267) @@ -0,0 +1,52 @@ +from fielddisplay import fielddisplay +from checkfield import checkfield +from WriteData import WriteData + +class calvingminthickness(object): + """ + CALVINGMINTHICKNESS class definition + + Usage: + calvingminthickness=calvingminthickness() + """ + + def __init__(self): # {{{ + + self.min_thickness = 0. + self.meltingrate = float('NaN') + + #set defaults + self.setdefaultparameters() + + #}}} + def __repr__(self): # {{{ + string=' Calving Minimum thickness:' + string="%s\n%s"%(string,fielddisplay(self,'min_thickness','minimum thickness below which no ice is allowed')) + string="%s\n%s"%(string,fielddisplay(self,'meltingrate','melting rate at given location [m/a]')) + return string + #}}} + def extrude(self,md): # {{{ + self.meltingrate=project3d(md,'vector',self.meltingrate,'type','node') + return self + #}}} + def setdefaultparameters(self): # {{{ + + #minimum thickness is 100 m by default + self.min_thickness = 100. + #}}} + def checkconsistency(self,md,solution,analyses): # {{{ + + #Early return + if solution == 'TransientSolution' or md.transient.ismovingfront == 0: + return + + md = checkfield(md,'fieldname','calving.min_thickness','>',0,'NaN',1,'Inf',1) + md = checkfield(md,'fieldname','calving.meltingrate','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices],'>=',0) + return md + # }}} + def marshall(self,prefix,md,fid): # {{{ + yts=md.constants.yts + WriteData(fid,prefix,'name','md.calving.law','data',4,'format','Integer') + WriteData(fid,prefix,'object',self,'fieldname','min_thickness','format','Double') + WriteData(fid,prefix,'object',self,'fieldname','meltingrate','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'scale',1./yts) + # }}} Index: /issm/trunk-jpl/src/m/classes/calvingvonmises.py =================================================================== --- /issm/trunk-jpl/src/m/classes/calvingvonmises.py (revision 22267) +++ /issm/trunk-jpl/src/m/classes/calvingvonmises.py (revision 22267) @@ -0,0 +1,60 @@ +from fielddisplay import fielddisplay +from project3d import project3d +from checkfield import checkfield +from WriteData import WriteData + +class calvingvonmises(object): + """ + CALVINGVONMISES class definition + + Usage: + calvingvonmises=calvingvonmises() + """ + + def __init__(self): # {{{ + + self.stress_threshold_groundedice = 0. + self.stress_threshold_floatingice = 0. + self.meltingrate = float('NaN') + + #set defaults + self.setdefaultparameters() + + #}}} + def __repr__(self): # {{{ + string=' Calving VonMises parameters:' + string="%s\n%s"%(string,fielddisplay(self,'stress_threshold_groundedice','sigma_max applied to grounded ice only [Pa]')) + string="%s\n%s"%(string,fielddisplay(self,'stress_threshold_floatingice','sigma_max applied to floating ice only [Pa]')) + + string="%s\n%s"%(string,fielddisplay(self,'meltingrate','melting rate at given location [m/a]')) + return string + #}}} + def extrude(self,md): # {{{ + self.meltingrate=project3d(md,'vector',self.meltingrate,'type','node') + return self + #}}} + def setdefaultparameters(self): # {{{ + #Default sigma max + self.stress_threshold_groundedice = 1e6 + self.stress_threshold_floatingice = 150e3 + return self + #}}} + def checkconsistency(self,md,solution,analyses): # {{{ + #Early return + if solution == 'TransientSolution' or md.transient.ismovingfront == 0: + return + + md = checkfield(md,'fieldname','calving.stress_threshold_groundedice','>',0,'nan',1,'Inf',1) + md = checkfield(md,'fieldname','calving.stress_threshold_floatingice','>',0,'nan',1,'Inf',1) + md = checkfield(md,'fieldname','calving.meltingrate','NaN',1,'Inf',1,'timeseries',1,'>=',0) + + return md + # }}} + def marshall(self,prefix,md,fid): # {{{ + yts=md.constants.yts + + WriteData(fid,prefix,'name','md.calving.law','data',2,'format','Integer') + WriteData(fid,prefix,'object',self,'fieldname','stress_threshold_groundedice','format','DoubleMat','mattype',1) + WriteData(fid,prefix,'object',self,'fieldname','stress_threshold_floatingice','format','DoubleMat','mattype',1) + WriteData(fid,prefix,'object',self,'fieldname','meltingrate','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts,'scale',1./yts) + # }}} Index: /issm/trunk-jpl/src/m/classes/constants.py =================================================================== --- /issm/trunk-jpl/src/m/classes/constants.py (revision 22266) +++ /issm/trunk-jpl/src/m/classes/constants.py (revision 22267) @@ -12,8 +12,8 @@ def __init__(self): # {{{ - self.g = 0 - self.omega = 0 - self.yts = 0 - self.referencetemperature = 0 + self.g = 0. + self.omega = 0. + self.yts = 0. + self.referencetemperature = 0. #set defaults @@ -49,8 +49,8 @@ def checkconsistency(self,md,solution,analyses): # {{{ - md = checkfield(md,'fieldname','constants.g','>',0,'size',[1]) - md = checkfield(md,'fieldname','constants.omega','>=',0,'size',[1]) - md = checkfield(md,'fieldname','constants.yts','>',0,'size',[1]) - md = checkfield(md,'fieldname','constants.referencetemperature','size',[1]) + md = checkfield(md,'fieldname','constants.g','>=',0,'size',[1,1]) + md = checkfield(md,'fieldname','constants.omega','>=',0,'size',[1,1]) + md = checkfield(md,'fieldname','constants.yts','>',0,'size',[1,1]) + md = checkfield(md,'fieldname','constants.referencetemperature','size',[1,1]) return md Index: /issm/trunk-jpl/src/m/classes/flowequation.py =================================================================== --- /issm/trunk-jpl/src/m/classes/flowequation.py (revision 22266) +++ /issm/trunk-jpl/src/m/classes/flowequation.py (revision 22267) @@ -92,5 +92,5 @@ md = checkfield(md,'fieldname','flowequation.fe_SSA','values',['P1','P1bubble','P1bubblecondensed','P2','P2bubble']) md = checkfield(md,'fieldname','flowequation.fe_HO' ,'values',['P1','P1bubble','P1bubblecondensed','P1xP2','P2xP1','P2','P2bubble','P1xP3','P2xP4']) - md = checkfield(md,'fieldname','flowequation.fe_FS' ,'values',['P1P1','P1P1GLS','MINIcondensed','MINI','TaylorHood','XTaylorHood','OneLayerP4z','CrouzeixRaviart']) + md = checkfield(md,'fieldname','flowequation.fe_FS' ,'values',['P1P1','P1P1GLS','MINIcondensed','MINI','TaylorHood','LATaylorHood','XTaylorHood','OneLayerP4z','CrouzeixRaviart','LACrouzeixRaviart']) md = checkfield(md,'fieldname','flowequation.borderSSA','size',[md.mesh.numberofvertices],'values',[0,1]) md = checkfield(md,'fieldname','flowequation.borderHO','size',[md.mesh.numberofvertices],'values',[0,1]) @@ -104,4 +104,7 @@ md = checkfield(md,'fieldname','flowequation.vertex_equation','size',[md.mesh.numberofvertices],'values',[1,2]) md = checkfield(md,'fieldname','flowequation.element_equation','size',[md.mesh.numberofelements],'values',[1,2]) + elif m.strcmp(md.mesh.domaintype(),'2Dvertical'): + md = checkfield(md,'fieldname','flowequation.vertex_equation','size',[md.mesh.numberofvertices],'values',[2,4,5]) + md = checkfield(md,'fieldname','flowequation.element_equation','size',[md.mesh.numberofelements],'values',[2,4,5]) elif m.strcmp(md.mesh.domaintype(),'3D'): md = checkfield(md,'fieldname','flowequation.vertex_equation','size',[md.mesh.numberofvertices],'values',np.arange(0,8+1)) Index: /issm/trunk-jpl/src/m/classes/frictionsommers.py =================================================================== --- /issm/trunk-jpl/src/m/classes/frictionsommers.py (revision 22267) +++ /issm/trunk-jpl/src/m/classes/frictionsommers.py (revision 22267) @@ -0,0 +1,47 @@ +from fielddisplay import fielddisplay +from project3d import project3d +from checkfield import checkfield +from WriteData import WriteData + +class frictionsommers(object): + """ + FRICTIONSOMMERS class definition + + Usage: + friction=frictionsommers() + """ + + def __init__(self,md): # {{{ + self.coefficient = md.friction.coefficient + #set defaults + self.setdefaultparameters() + + #}}} + def __repr__(self): # {{{ + string="Basal shear stress parameters: Sigma_b = coefficient^2 * Neff * u_b\n(effective stress Neff=rho_ice*g*thickness+rho_water*g*(head-b))" + + string="%s\n%s"%(string,fielddisplay(self,"coefficient","friction coefficient [SI]")) + return string + #}}} + def extrude(self,md): # {{{ + self.coefficient=project3d(md,'vector',self.coefficient,'type','node','layer',1) + return self + #}}} + def setdefaultparameters(self): # {{{ + return self + #}}} + def checkconsistency(self,md,solution,analyses): # {{{ + + #Early return + if 'StressbalanceAnalysis' not in analyses and 'ThermalAnalysis' not in analyses: + return md + + md = checkfield(md,'fieldname','friction.coefficient','timeseries',1,'NaN',1,'Inf',1) + return md + + # }}} + def marshall(self,prefix,md,fid): # {{{ + yts=md.constants.yts + WriteData(fid,prefix,'name','md.friction.law','data',8,'format','Integer') + WriteData(fid,prefix,'object',self,'fieldname','coefficient','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts) + # }}} Index: /issm/trunk-jpl/src/m/classes/hydrologysommers.py =================================================================== --- /issm/trunk-jpl/src/m/classes/hydrologysommers.py (revision 22267) +++ /issm/trunk-jpl/src/m/classes/hydrologysommers.py (revision 22267) @@ -0,0 +1,90 @@ +from fielddisplay import fielddisplay +from checkfield import checkfield +from WriteData import WriteData + +class hydrologysommers(object): + """ + HYDROLOGYSOMMERS class definition + + Usage: + hydrologysommers=hydrologysommers() + """ + + def __init__(self): # {{{ + self.head = float('NaN') + self.gap_height = float('NaN') + self.bump_spacing = float('NaN') + self.bump_height = float('NaN') + self.englacial_input = float('NaN') + self.moulin_input = float('NaN') + self.reynolds = float('NaN') + self.spchead = float('NaN') + self.neumannflux = float('NaN') + self.relaxation = 0 + self.storage = 0 + + #set defaults + self.setdefaultparameters() + + #}}} + def __repr__(self): # {{{ + + string=' hydrologysommers solution parameters:' + string="%s\n%s"%(string,fielddisplay(self,'head','subglacial hydrology water head (m)')) + string="%s\n%s"%(string,fielddisplay(self,'gap_height','height of gap separating ice to bed (m)')) + string="%s\n%s"%(string,fielddisplay(self,'bump_spacing','characteristic bedrock bump spacing (m)')) + string="%s\n%s"%(string,fielddisplay(self,'bump_height','characteristic bedrock bump height (m)')) + string="%s\n%s"%(string,fielddisplay(self,'englacial_input','liquid water input from englacial to subglacial system (m/yr)')) + string="%s\n%s"%(string,fielddisplay(self,'moulin_input','liquid water input from moulins (at the vertices) to subglacial system (m^3/s)')) + string="%s\n%s"%(string,fielddisplay(self,'reynolds','Reynolds'' number')) + string="%s\n%s"%(string,fielddisplay(self,'neumannflux','water flux applied along the model boundary (m^2/s)')) + string="%s\n%s"%(string,fielddisplay(self,'spchead','water head constraints (NaN means no constraint) (m)')) + string="%s\n%s"%(string,fielddisplay(self,'relaxation','under-relaxation coefficient for nonlinear iteration')) + string="%s\n%s"%(string,fielddisplay(self,'storage','englacial storage coefficient (void ratio)')) + return string + #}}} + def extrude(self,md): # {{{ + return self + #}}} + def setdefaultparameters(self): # {{{ + # Set under-relaxation parameter to be 1 (no under-relaxation of nonlinear iteration) + self.relaxation=1 + self.storage=0 + return self + #}}} + def checkconsistency(self,md,solution,analyses): # {{{ + + #Early return + if 'HydrologySommersAnalysis' not in analyses: + return md + + md = checkfield(md,'fieldname','hydrology.head','size',[md.mesh.numberofvertices],'NaN',1,'Inf',1) + md = checkfield(md,'fieldname','hydrology.gap_height','>=',0,'size',[md.mesh.numberofelements],'NaN',1,'Inf',1) + md = checkfield(md,'fieldname','hydrology.bump_spacing','>',0,'size',[md.mesh.numberofelements],'NaN',1,'Inf',1) + md = checkfield(md,'fieldname','hydrology.bump_height','>=',0,'size',[md.mesh.numberofelements],'NaN',1,'Inf',1) + md = checkfield(md,'fieldname','hydrology.englacial_input','>=',0,'NaN',1,'Inf',1,'timeseries',1) + md = checkfield(md,'fieldname','hydrology.moulin_input','>=',0,'NaN',1,'Inf',1,'timeseries',1) + md = checkfield(md,'fieldname','hydrology.reynolds','>',0,'size',[md.mesh.numberofelements],'NaN',1,'Inf',1) + md = checkfield(md,'fieldname','hydrology.neumannflux','timeseries',1,'NaN',1,'Inf',1) + md = checkfield(md,'fieldname','hydrology.spchead','size',[md.mesh.numberofvertices]) + md = checkfield(md,'fieldname','hydrology.relaxation','>=',0) + md = checkfield(md,'fieldname','hydrology.storage','>=',0) + + return md + # }}} + def marshall(self,prefix,md,fid): # {{{ + yts=md.constants.yts + + WriteData(fid,prefix,'name','md.hydrology.model','data',3,'format','Integer') + WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','head','format','DoubleMat','mattype',1) + WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','gap_height','format','DoubleMat','mattype',2) + WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','bump_spacing','format','DoubleMat','mattype',2) + WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','bump_height','format','DoubleMat','mattype',2) + WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','englacial_input','format','DoubleMat','mattype',1,'scale',1./yts,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts) + WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','moulin_input','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts) + WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','reynolds','format','DoubleMat','mattype',2) + WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','neumannflux','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts) + WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','spchead','format','DoubleMat','mattype',1) + WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','relaxation','format','Double') + WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','storage','format','Double') + # }}} Index: /issm/trunk-jpl/src/m/classes/matenhancedice.py =================================================================== --- /issm/trunk-jpl/src/m/classes/matenhancedice.py (revision 22267) +++ /issm/trunk-jpl/src/m/classes/matenhancedice.py (revision 22267) @@ -0,0 +1,171 @@ +from fielddisplay import fielddisplay +from project3d import project3d +from checkfield import checkfield +from WriteData import WriteData + +class matenhancedice(object): + """ + MATICE class definition + + Usage: + matenhancedice=matenhancedice(); + """ + + def __init__(self): # {{{ + self.rho_ice = 0. + self.rho_water = 0. + self.rho_freshwater = 0. + self.mu_water = 0. + self.heatcapacity = 0. + self.latentheat = 0. + self.thermalconductivity = 0. + self.temperateiceconductivity = 0. + self.meltingpoint = 0. + self.beta = 0. + self.mixed_layer_capacity = 0. + self.thermal_exchange_velocity = 0. + self.rheology_E = float('NaN') + self.rheology_B = float('NaN') + self.rheology_n = float('NaN') + self.rheology_law = '' + + #giaivins: + self.lithosphere_shear_modulus = 0. + self.lithosphere_density = 0. + self.mantle_shear_modulus = 0. + self.mantle_density = 0. + + #SLR + self.earth_density= 0 # average density of the Earth, (kg/m^3) + + self.setdefaultparameters() + #}}} + def __repr__(self): # {{{ + string=" Materials:" + + string="%s\n%s"%(string,fielddisplay(self,"rho_ice","ice density [kg/m^3]")) + string="%s\n%s"%(string,fielddisplay(self,"rho_water","water density [kg/m^3]")) + string="%s\n%s"%(string,fielddisplay(self,"rho_freshwater","fresh water density [kg/m^3]")) + string="%s\n%s"%(string,fielddisplay(self,"mu_water","water viscosity [N s/m^2]")) + string="%s\n%s"%(string,fielddisplay(self,"heatcapacity","heat capacity [J/kg/K]")) + string="%s\n%s"%(string,fielddisplay(self,"thermalconductivity","ice thermal conductivity [W/m/K]")) + string="%s\n%s"%(string,fielddisplay(self,"temperateiceconductivity","temperate ice thermal conductivity [W/m/K]")) + string="%s\n%s"%(string,fielddisplay(self,"meltingpoint","melting point of ice at 1atm in K")) + string="%s\n%s"%(string,fielddisplay(self,"latentheat","latent heat of fusion [J/m^3]")) + string="%s\n%s"%(string,fielddisplay(self,"beta","rate of change of melting point with pressure [K/Pa]")) + string="%s\n%s"%(string,fielddisplay(self,"mixed_layer_capacity","mixed layer capacity [W/kg/K]")) + string="%s\n%s"%(string,fielddisplay(self,"thermal_exchange_velocity","thermal exchange velocity [m/s]")) + string="%s\n%s"%(string,fielddisplay(self,"rheology_E","enhancement factor")) + string="%s\n%s"%(string,fielddisplay(self,"rheology_B","flow law parameter [Pa/s^(1/n)]")) + string="%s\n%s"%(string,fielddisplay(self,"rheology_n","Glen's flow law exponent")) + string="%s\n%s"%(string,fielddisplay(self,"rheology_law","law for the temperature dependance of the rheology: 'None', 'BuddJacka', 'Cuffey', 'CuffeyTemperate', 'Paterson', 'Arrhenius' or 'LliboutryDuval'")) + string="%s\n%s"%(string,fielddisplay(self,"lithosphere_shear_modulus","Lithosphere shear modulus [Pa]")) + string="%s\n%s"%(string,fielddisplay(self,"lithosphere_density","Lithosphere density [g/cm^-3]")) + string="%s\n%s"%(string,fielddisplay(self,"mantle_shear_modulus","Mantle shear modulus [Pa]")) + string="%s\n%s"%(string,fielddisplay(self,"mantle_density","Mantle density [g/cm^-3]")) + string="%s\n%s"%(string,fielddisplay(self,"earth_density","Mantle density [kg/m^-3]")) + + return string + #}}} + def extrude(self,md): # {{{ + self.rheology_E=project3d(md,'vector',self.rheology_E,'type','node') + self.rheology_B=project3d(md,'vector',self.rheology_B,'type','node') + self.rheology_n=project3d(md,'vector',self.rheology_n,'type','element') + return self + #}}} + def setdefaultparameters(self): # {{{ + #ice density (kg/m^3) + self.rho_ice=917. + + #ocean water density (kg/m^3) + self.rho_water=1023. + + #fresh water density (kg/m^3) + self.rho_freshwater=1000. + + #water viscosity (N.s/m^2) + self.mu_water=0.001787 + + #ice heat capacity cp (J/kg/K) + self.heatcapacity=2093. + + #ice latent heat of fusion L (J/kg) + self.latentheat=3.34*10**5 + + #ice thermal conductivity (W/m/K) + self.thermalconductivity=2.4 + + #temperate ice thermal conductivity (W/m/K) + self.temperateiceconductivity=0.24 + + #the melting point of ice at 1 atmosphere of pressure in K + self.meltingpoint=273.15 + + #rate of change of melting point with pressure (K/Pa) + self.beta=9.8*10**-8 + + #mixed layer (ice-water interface) heat capacity (J/kg/K) + self.mixed_layer_capacity=3974. + + #thermal exchange velocity (ice-water interface) (m/s) + self.thermal_exchange_velocity=1.00*10**-4 + + #Rheology law: what is the temperature dependence of B with T + #available: none, paterson and arrhenius + self.rheology_law='Paterson' + + # GIA: + self.lithosphere_shear_modulus = 6.7*10**10 # (Pa) + self.lithosphere_density = 3.32 # (g/cm^-3) + self.mantle_shear_modulus = 1.45*10**11 # (Pa) + self.mantle_density = 3.34 # (g/cm^-3) + + #SLR + self.earth_density= 5512 #average density of the Earth, (kg/m^3) + + return self + #}}} + def checkconsistency(self,md,solution,analyses): # {{{ + md = checkfield(md,'fieldname','materials.rho_ice','>',0) + md = checkfield(md,'fieldname','materials.rho_water','>',0) + md = checkfield(md,'fieldname','materials.rho_freshwater','>',0) + md = checkfield(md,'fieldname','materials.mu_water','>',0) + md = checkfield(md,'fieldname','materials.rheology_E','>',0,'timeseries',1,'NaN',1,'Inf',1) + md = checkfield(md,'fieldname','materials.rheology_B','>',0,'timeseries',1,'NaN',1,'Inf',1) + md = checkfield(md,'fieldname','materials.rheology_n','>',0,'size',[md.mesh.numberofelements]) + md = checkfield(md,'fieldname','materials.rheology_law','values',['None','BuddJacka','Cuffey','CuffeyTemperate','Paterson','Arrhenius','LliboutryDuval']) + + if 'GiaAnalysis' in analyses: + md = checkfield(md,'fieldname','materials.lithosphere_shear_modulus','>',0,'numel',1) + md = checkfield(md,'fieldname','materials.lithosphere_density','>',0,'numel',1) + md = checkfield(md,'fieldname','materials.mantle_shear_modulus','>',0,'numel',1) + md = checkfield(md,'fieldname','materials.mantle_density','>',0,'numel',1) + if 'SealevelriseAnalysis' in analyses: + md = checkfield(md,'fieldname','materials.earth_density','>',0,'numel',1) + return md + # }}} + def marshall(self,prefix,md,fid): # {{{ + WriteData(fid,prefix,'name','md.materials.type','data',4,'format','Integer') + WriteData(fid,prefix,'object',self,'class','materials','fieldname','rho_ice','format','Double') + WriteData(fid,prefix,'object',self,'class','materials','fieldname','rho_water','format','Double') + WriteData(fid,prefix,'object',self,'class','materials','fieldname','rho_freshwater','format','Double') + WriteData(fid,prefix,'object',self,'class','materials','fieldname','mu_water','format','Double') + WriteData(fid,prefix,'object',self,'class','materials','fieldname','heatcapacity','format','Double') + WriteData(fid,prefix,'object',self,'class','materials','fieldname','latentheat','format','Double') + WriteData(fid,prefix,'object',self,'class','materials','fieldname','thermalconductivity','format','Double') + WriteData(fid,prefix,'object',self,'class','materials','fieldname','temperateiceconductivity','format','Double') + WriteData(fid,prefix,'object',self,'class','materials','fieldname','meltingpoint','format','Double') + WriteData(fid,prefix,'object',self,'class','materials','fieldname','beta','format','Double') + WriteData(fid,prefix,'object',self,'class','materials','fieldname','mixed_layer_capacity','format','Double') + WriteData(fid,prefix,'object',self,'class','materials','fieldname','thermal_exchange_velocity','format','Double') + WriteData(fid,prefix,'object',self,'class','materials','fieldname','rheology_E','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts) + WriteData(fid,prefix,'object',self,'class','materials','fieldname','rheology_B','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts) + WriteData(fid,prefix,'object',self,'class','materials','fieldname','rheology_n','format','DoubleMat','mattype',2) + WriteData(fid,prefix,'data',self.rheology_law,'name','md.materials.rheology_law','format','String') + + WriteData(fid,prefix,'object',self,'class','materials','fieldname','lithosphere_shear_modulus','format','Double') + WriteData(fid,prefix,'object',self,'class','materials','fieldname','lithosphere_density','format','Double','scale',10^3) + WriteData(fid,prefix,'object',self,'class','materials','fieldname','mantle_shear_modulus','format','Double') + WriteData(fid,prefix,'object',self,'class','materials','fieldname','mantle_density','format','Double','scale',10^3) + WriteData(fid,prefix,'object',self,'class','materials','fieldname','earth_density','format','Double') + # }}} Index: /issm/trunk-jpl/src/m/classes/matestar.py =================================================================== --- /issm/trunk-jpl/src/m/classes/matestar.py (revision 22267) +++ /issm/trunk-jpl/src/m/classes/matestar.py (revision 22267) @@ -0,0 +1,175 @@ +import numpy as np +from fielddisplay import fielddisplay +from project3d import project3d +from checkfield import checkfield +from WriteData import WriteData + +class matestar(object): + """ + matestar class definition + + Usage: + matestar=matestar() + """ + + def __init__(self): # {{{ + + rho_ice = 0. + rho_water = 0. + rho_freshwater = 0. + mu_water = 0. + heatcapacity = 0. + latentheat = 0. + thermalconductivity = 0. + temperateiceconductivity = 0. + meltingpoint = 0. + beta = 0. + mixed_layer_capacity = 0. + thermal_exchange_velocity = 0. + rheology_B = float('NaN') + rheology_Ec = float('NaN') + rheology_Es = float('NaN') + rheology_law = '' + + #giaivins: + lithosphere_shear_modulus = 0. + lithosphere_density = 0. + mantle_shear_modulus = 0. + mantle_density = 0. + + #slr + earth_density = 0 + + #set default parameters: + self.setdefaultparameters() + #}}} + def __repr__(self): # {{{ + string=" Materials:" + + string="%s\n%s"%(string,fielddisplay(self,'rho_ice','ice density [kg/m^3]')) + string="%s\n%s"%(string,fielddisplay(self,'rho_water','ocean water density [kg/m^3]')) + string="%s\n%s"%(string,fielddisplay(self,'rho_freshwater','fresh water density [kg/m^3]')) + string="%s\n%s"%(string,fielddisplay(self,'mu_water','water viscosity [N s/m^2]')) + string="%s\n%s"%(string,fielddisplay(self,'heatcapacity','heat capacity [J/kg/K]')) + string="%s\n%s"%(string,fielddisplay(self,'thermalconductivity',['ice thermal conductivity [W/m/K]'])) + string="%s\n%s"%(string,fielddisplay(self,'temperateiceconductivity','temperate ice thermal conductivity [W/m/K]')) + string="%s\n%s"%(string,fielddisplay(self,'meltingpoint','melting point of ice at 1atm in K')) + string="%s\n%s"%(string,fielddisplay(self,'latentheat','latent heat of fusion [J/kg]')) + string="%s\n%s"%(string,fielddisplay(self,'beta','rate of change of melting point with pressure [K/Pa]')) + string="%s\n%s"%(string,fielddisplay(self,'mixed_layer_capacity','mixed layer capacity [W/kg/K]')) + string="%s\n%s"%(string,fielddisplay(self,'thermal_exchange_velocity','thermal exchange velocity [m/s]')) + string="%s\n%s"%(string,fielddisplay(self,'rheology_B','flow law parameter [Pa/s^(1/3)]')) + string="%s\n%s"%(string,fielddisplay(self,'rheology_Ec','compressive enhancement factor')) + string="%s\n%s"%(string,fielddisplay(self,'rheology_Es','shear enhancement factor')) + string="%s\n%s"%(string,fielddisplay(self,'rheology_law',['law for the temperature dependance of the rheology: ''None'', ''BuddJacka'', ''Cuffey'', ''CuffeyTemperate'', ''Paterson'', ''Arrhenius'' or ''LliboutryDuval'''])) + string="%s\n%s"%(string,fielddisplay(self,'lithosphere_shear_modulus','Lithosphere shear modulus [Pa]')) + string="%s\n%s"%(string,fielddisplay(self,'lithosphere_density','Lithosphere density [g/cm^-3]')) + string="%s\n%s"%(string,fielddisplay(self,'mantle_shear_modulus','Mantle shear modulus [Pa]')) + string="%s\n%s"%(string,fielddisplay(self,'mantle_density','Mantle density [g/cm^-3]')) + string="%s\n%s"%(string,fielddisplay(self,'earth_density','Mantle density [kg/m^-3]')) + + return string + #}}} + def extrude(self,md): # {{{ + self.rheology_B=project3d(md,'vector',self.rheology_B,'type','node') + self.rheology_Ec=project3d(md,'vector',self.rheology_Ec,'type','node') + self.rheology_Es=project3d(md,'vector',self.rheology_Es,'type','node') + return self + #}}} + def setdefaultparameters(self): # {{{ + #ice density (kg/m^3) + self.rho_ice=917. + + #ocean water density (kg/m^3) + self.rho_water=1023. + + #fresh water density (kg/m^3) + self.rho_freshwater=1000. + + #water viscosity (N.s/m^2) + self.mu_water=0.001787 + + #ice heat capacity cp (J/kg/K) + self.heatcapacity=2093. + + #ice latent heat of fusion L (J/kg) + self.latentheat=3.34*10**5 + + #ice thermal conductivity (W/m/K) + self.thermalconductivity=2.4 + + #wet ice thermal conductivity (W/m/K) + self.temperateiceconductivity=.24 + + #the melting point of ice at 1 atmosphere of pressure in K + self.meltingpoint=273.15 + + #rate of change of melting point with pressure (K/Pa) + self.beta=9.8*10**-8 + + #mixed layer (ice-water interface) heat capacity (J/kg/K) + self.mixed_layer_capacity=3974. + + #thermal exchange velocity (ice-water interface) (m/s) + self.thermal_exchange_velocity=1.00*10**-4 + + #Rheology law: what is the temperature dependence of B with T + #available: none, paterson and arrhenius + self.rheology_law='Paterson' + + # GIA: + self.lithosphere_shear_modulus = 6.7*10**10 # (Pa) + self.lithosphere_density = 3.32 # (g/cm^-3) + self.mantle_shear_modulus = 1.45*10**11 # (Pa) + self.mantle_density = 3.34 # (g/cm^-3) + + #SLR + self.earth_density= 5512 # average density of the Earth, (kg/m^3) + + return self + #}}} + def checkconsistency(self,md,solution,analyses): # {{{ + md = checkfield(md,'fieldname','materials.rho_ice','>',0) + md = checkfield(md,'fieldname','materials.rho_water','>',0) + md = checkfield(md,'fieldname','materials.rho_freshwater','>',0) + md = checkfield(md,'fieldname','materials.mu_water','>',0) + md = checkfield(md,'fieldname','materials.rheology_B','>',0,'size',[md.mesh.numberofvertices],'NaN',1,'Inf',1) + md = checkfield(md,'fieldname','materials.rheology_Ec','>',0,'size',[md.mesh.numberofvertices],'NaN',1,'Inf',1) + md = checkfield(md,'fieldname','materials.rheology_Es','>',0,'size',[md.mesh.numberofvertices],'NaN',1,'Inf',1) + md = checkfield(md,'fieldname','materials.rheology_law','values',['None','BuddJacka', 'Cuffey','CuffeyTemperate','Paterson','Arrhenius','LliboutryDuval']) + + if 'GiaAnalysis' in analyses: + md = checkfield(md,'fieldname','materials.lithosphere_shear_modulus','>',0,'numel',1) + md = checkfield(md,'fieldname','materials.lithosphere_density','>',0,'numel',1) + md = checkfield(md,'fieldname','materials.mantle_shear_modulus','>',0,'numel',1) + md = checkfield(md,'fieldname','materials.mantle_density','>',0,'numel',1) + if 'SealevelriseAnalysis' in analyses: + md = checkfield(md,'fieldname','materials.earth_density','>',0,'numel',1) + + return md + # }}} + def marshall(self,prefix,md,fid): # {{{ + WriteData(fid,prefix,'name','md.materials.type','data',2,'format','Integer') + WriteData(fid,prefix,'object',self,'class','materials','fieldname','rho_ice','format','Double') + WriteData(fid,prefix,'object',self,'class','materials','fieldname','rho_water','format','Double') + WriteData(fid,prefix,'object',self,'class','materials','fieldname','rho_freshwater','format','Double') + WriteData(fid,prefix,'object',self,'class','materials','fieldname','mu_water','format','Double') + WriteData(fid,prefix,'object',self,'class','materials','fieldname','heatcapacity','format','Double') + WriteData(fid,prefix,'object',self,'class','materials','fieldname','latentheat','format','Double') + WriteData(fid,prefix,'object',self,'class','materials','fieldname','thermalconductivity','format','Double') + WriteData(fid,prefix,'object',self,'class','materials','fieldname','temperateiceconductivity','format','Double') + WriteData(fid,prefix,'object',self,'class','materials','fieldname','meltingpoint','format','Double') + WriteData(fid,prefix,'object',self,'class','materials','fieldname','beta','format','Double') + WriteData(fid,prefix,'object',self,'class','materials','fieldname','mixed_layer_capacity','format','Double') + WriteData(fid,prefix,'object',self,'class','materials','fieldname','thermal_exchange_velocity','format','Double') + WriteData(fid,prefix,'object',self,'class','materials','fieldname','rheology_B','format','DoubleMat','mattype',1) + WriteData(fid,prefix,'object',self,'class','materials','fieldname','rheology_Ec','format','DoubleMat','mattype',1) + WriteData(fid,prefix,'object',self,'class','materials','fieldname','rheology_Es','format','DoubleMat','mattype',1) + WriteData(fid,prefix,'data',self.rheology_law,'name','md.materials.rheology_law','format','String') + + WriteData(fid,prefix,'object',self,'class','materials','fieldname','lithosphere_shear_modulus','format','Double') + WriteData(fid,prefix,'object',self,'class','materials','fieldname','lithosphere_density','format','Double','scale',10^3) + WriteData(fid,prefix,'object',self,'class','materials','fieldname','mantle_shear_modulus','format','Double') + WriteData(fid,prefix,'object',self,'class','materials','fieldname','mantle_density','format','Double','scale',10**3) + WriteData(fid,prefix,'object',self,'class','materials','fieldname','earth_density','format','Double') + # }}} Index: /issm/trunk-jpl/src/m/classes/mesh2dvertical.py =================================================================== --- /issm/trunk-jpl/src/m/classes/mesh2dvertical.py (revision 22267) +++ /issm/trunk-jpl/src/m/classes/mesh2dvertical.py (revision 22267) @@ -0,0 +1,130 @@ +import numpy as np +from fielddisplay import fielddisplay +from checkfield import checkfield +import MatlabFuncs as m +from WriteData import WriteData + +class mesh2dvertical(object): + """ + MESH2DVERTICAL class definition + + Usage: + mesh2dvertical=mesh2dvertical(); + """ + + def __init__(self): # {{{ + self.x = float('NaN') + self.y = float('NaN') + self.elements = float('NaN') + self.numberofelements = 0 + self.numberofvertices = 0 + self.numberofedges = 0 + + self.lat = float('NaN') + self.long = float('NaN') + self.epsg = float('NaN') + + self.vertexonboundary = float('NaN') + self.vertexonbase = float('NaN') + self.vertexonsurface = float('NaN') + + self.edges = float('NaN') + self.segments = float('NaN') + self.segmentmarkers = float('NaN') + self.vertexconnectivity = float('NaN') + self.elementconnectivity = float('NaN') + self.average_vertex_connectivity = 0 + + #set defaults + self.setdefaultparameters() + + #}}} + def __repr__(self): # {{{ + string=" 2D tria Mesh (vertical):" + + string="%s\n%s"%(string,"\n Elements and vertices:") + string="%s\n%s"%(string,fielddisplay(self,"numberofelements","number of elements")) + string="%s\n%s"%(string,fielddisplay(self,"numberofvertices","number of vertices")) + string="%s\n%s"%(string,fielddisplay(self,"elements","vertex indices of the mesh elements")) + string="%s\n%s"%(string,fielddisplay(self,"x","vertices x coordinate [m]")) + string="%s\n%s"%(string,fielddisplay(self,"y","vertices y coordinate [m]")) + string="%s\n%s"%(string,fielddisplay(self,"edges","edges of the 2d mesh (vertex1 vertex2 element1 element2)")) + string="%s\n%s"%(string,fielddisplay(self,"numberofedges","number of edges of the 2d mesh")) + + string="%s%s"%(string,"\n\n Properties:") + string="%s\n%s"%(string,fielddisplay(self,"vertexonboundary","vertices on the boundary of the domain flag list")) + string="%s\n%s"%(string,fielddisplay(self,'vertexonbase','vertices on the bed of the domain flag list')) + string="%s\n%s"%(string,fielddisplay(self,'vertexonsurface','vertices on the surface of the domain flag list')) + string="%s\n%s"%(string,fielddisplay(self,"segments","edges on domain boundary (vertex1 vertex2 element)")) + string="%s\n%s"%(string,fielddisplay(self,"segmentmarkers","number associated to each segment")) + string="%s\n%s"%(string,fielddisplay(self,"vertexconnectivity","list of vertices connected to vertex_i")) + string="%s\n%s"%(string,fielddisplay(self,"elementconnectivity","list of vertices connected to element_i")) + string="%s\n%s"%(string,fielddisplay(self,"average_vertex_connectivity","average number of vertices connected to one vertex")) + + string="%s%s"%(string,"\n\n Projection:") + string="%s\n%s"%(string,fielddisplay(self,"lat","vertices latitude [degrees]")) + string="%s\n%s"%(string,fielddisplay(self,"long","vertices longitude [degrees]")) + string="%s\n%s"%(string,fielddisplay(self,"epsg","EPSG code (ex: 3413 for UPS Greenland, 3031 for UPS Antarctica)")) + return string + #}}} + def setdefaultparameters(self): # {{{ + + #the connectivity is the averaged number of nodes linked to a + #given node through an edge. This connectivity is used to initially + #allocate memory to the stiffness matrix. A value of 16 seems to + #give a good memory/time ration. This value can be checked in + #trunk/test/Miscellaneous/runme.m + self.average_vertex_connectivity=25. + + return self + #}}} + def checkconsistency(self,md,solution,analyses): # {{{ + if(solution=='LoveSolution'): + return + + md = checkfield(md,'fieldname','mesh.x','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices]) + md = checkfield(md,'fieldname','mesh.y','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices]) + md = checkfield(md,'fieldname','mesh.elements','NaN',1,'Inf',1,'>',0,'values',np.arange(1,md.mesh.numberofvertices+1)) + md = checkfield(md,'fieldname','mesh.elements','size',[md.mesh.numberofelements,3]) + if np.any(np.logical_not(m.ismember(np.arange(1,md.mesh.numberofvertices+1),md.mesh.elements))): + md.checkmessage("orphan nodes have been found. Check the mesh outline") + md = checkfield(md,'fieldname','mesh.numberofelements','>',0) + md = checkfield(md,'fieldname','mesh.numberofvertices','>',0) + md = checkfield(md,'fieldname','mesh.vertexonbase','size',[md.mesh.numberofvertices],'values',[0,1]) + md = checkfield(md,'fieldname','mesh.vertexonsurface','size',[md.mesh.numberofvertices],'values',[0,1]) + md = checkfield(md,'fieldname','mesh.average_vertex_connectivity','>=',9,'message',"'mesh.average_vertex_connectivity' should be at least 9 in 2d") + + if solution=='ThermalSolution': + md.checkmessage("thermal not supported for 2d mesh") + + return md + # }}} + def domaintype(self): # {{{ + return "2Dvertical" + #}}} + def dimension(self): # {{{ + return 2 + #}}} + def elementtype(self): # {{{ + return "Tria" + #}}} + def vertexflags(self,value): # {{{ + flags = np.zeros((self.numberofvertices,)) + pos = self.segments[np.where(self.segmentmarkers==value),0:2]-1 + flags[pos] = 1 + return flags + #}}} + def marshall(self,prefix,md,fid): # {{{ + WriteData(fid,prefix,'name','md.mesh.domain_type','data',"Domain"+self.domaintype(),'format','String'); + WriteData(fid,prefix,'name','md.mesh.domain_dimension','data',self.dimension(),'format','Integer'); + WriteData(fid,prefix,'name','md.mesh.elementtype','data',self.elementtype(),'format','String'); + WriteData(fid,prefix,'object',self,'class','mesh','fieldname','x','format','DoubleMat','mattype',1) + WriteData(fid,prefix,'object',self,'class','mesh','fieldname','y','format','DoubleMat','mattype',1) + WriteData(fid,prefix,'name','md.mesh.z','data',np.zeros(self.numberofvertices),'format','DoubleMat','mattype',1); + WriteData(fid,prefix,'object',self,'class','mesh','fieldname','elements','format','DoubleMat','mattype',2) + WriteData(fid,prefix,'object',self,'class','mesh','fieldname','numberofelements','format','Integer') + WriteData(fid,prefix,'object',self,'class','mesh','fieldname','numberofvertices','format','Integer') + WriteData(fid,prefix,'object',self,'class','mesh','fieldname','vertexonbase','format','BooleanMat','mattype',1) + WriteData(fid,prefix,'object',self,'class','mesh','fieldname','vertexonsurface','format','BooleanMat','mattype',1) + WriteData(fid,prefix,'object',self,'class','mesh','fieldname','average_vertex_connectivity','format','Integer') + # }}} Index: /issm/trunk-jpl/src/m/classes/mismipbasalforcings.py =================================================================== --- /issm/trunk-jpl/src/m/classes/mismipbasalforcings.py (revision 22266) +++ /issm/trunk-jpl/src/m/classes/mismipbasalforcings.py (revision 22267) @@ -40,6 +40,9 @@ def initialize(self,md): # {{{ if np.all(np.isnan(self.groundedice_melting_rate)): - self.groundedice_melting_rate=np.zeros(md.mesh.numberofvertices) + self.groundedice_melting_rate=np.zeros((md.mesh.numberofvertices)) print ' no basalforcings.groundedice_melting_rate specified: values set as zero' + if np.all(np.isnan(self.geothermalflux)): + self.geothermalflux=np.zeros((md.mesh.numberofvertices)) + print " no basalforcings.geothermalflux specified: values set as zero" return self #}}} @@ -84,5 +87,5 @@ floatingice_melting_rate = np.zeros((md.mesh.numberofvertices)) - floatingice_melting_rate = md.basalforcings.meltrate_factor*np.tanh((md.geometry.base-md.geometry.bed)/md.basalforcings.threshold_thickness)*np.amax(md.basalforcings.upperdepth_melt-md.geometry.base,0) + floatingice_melting_rate = md.basalforcings.meltrate_factor*np.tanh((md.geometry.base-md.geometry.bed)/md.basalforcings.threshold_thickness)*(md.basalforcings.upperdepth_melt-md.geometry.base) WriteData(fid,prefix,'name','md.basalforcings.model','data',3,'format','Integer') Index: /issm/trunk-jpl/src/m/classes/plumebasalforcings.py =================================================================== --- /issm/trunk-jpl/src/m/classes/plumebasalforcings.py (revision 22267) +++ /issm/trunk-jpl/src/m/classes/plumebasalforcings.py (revision 22267) @@ -0,0 +1,126 @@ +from fielddisplay import fielddisplay +from checkfield import checkfield +from WriteData import WriteData +from project3d import project3d + +class plumebasalforcings(object): + ''' + PLUME BASAL FORCINGS class definition + + Usage: + plumebasalforcings=plumebasalforcings() + ''' + + def __init__(self): # {{{ + floatingice_melting_rate = float('NaN') + groundedice_melting_rate = float('NaN') + mantleconductivity = float('NaN') + nusselt = float('NaN') + dtbg = float('NaN') + plumeradius = float('NaN') + topplumedepth = float('NaN') + bottomplumedepth = float('NaN') + plumex = float('NaN') + plumey = float('NaN') + crustthickness = float('NaN') + uppercrustthickness = float('NaN') + uppercrustheat = float('NaN') + lowercrustheat = float('NaN') + + self.setdefaultparameters() + #}}} + + def __repr__(self): # {{{ + print ' mantle plume basal melt parameterization:' + + string="%s\n%s"%(string,fielddisplay(self,'groundedice_melting_rate','basal melting rate (positive if melting) [m/yr]')) + string="%s\n%s"%(string,fielddisplay(self,'floatingice_melting_rate','basal melting rate (positive if melting) [m/yr]')) + string="%s\n%s"%(string,fielddisplay(self,'mantleconductivity','mantle heat conductivity [W/m^3]')) + string="%s\n%s"%(string,fielddisplay(self,'nusselt','nusselt number, ratio of mantle to plume [1]')) + string="%s\n%s"%(string,fielddisplay(self,'dtbg','background temperature gradient [degree/m]')) + string="%s\n%s"%(string,fielddisplay(self,'plumeradius','radius of the mantle plume [m]')) + string="%s\n%s"%(string,fielddisplay(self,'topplumedepth','depth of the mantle plume top below the crust [m]')) + string="%s\n%s"%(string,fielddisplay(self,'bottomplumedepth','depth of the mantle plume base below the crust [m]')) + string="%s\n%s"%(string,fielddisplay(self,'plumex','x coordinate of the center of the plume [m]')) + string="%s\n%s"%(string,fielddisplay(self,'plumey','y coordinate of the center of the plume [m]')) + string="%s\n%s"%(string,fielddisplay(self,'crustthickness','thickness of the crust [m]')) + string="%s\n%s"%(string,fielddisplay(self,'uppercrustthickness','thickness of the upper crust [m]')) + string="%s\n%s"%(string,fielddisplay(self,'uppercrustheat','volumic heat of the upper crust [w/m^3]')) + string="%s\n%s"%(string,fielddisplay(self,'lowercrustheat','volumic heat of the lowercrust [w/m^3]')) + + return string + #}}} + + def initialize(self,md): #{{{ + if np.all(np.isnan(self.groundedice_melting_rate)): + self.groundedice_melting_rate=np.zeros((md.mesh.numberofvertices,)) + print ' no basalforcings.groundedice_melting_rate specified: values set as zero' + if np.all(np.isnan(self.floatingice_melting_rate)): + self.floatingice_melting_rate=np.zeros((md.mesh.numberofvertices,)) + print ' no basalforcings.floatingice_melting_rate specified: values set as zero' + return + #}}} + + def extrude(self,md): # {{{ + self.groundedice_melting_rate=project3d(md,'vector',self.groundedice_melting_rate,'type','node','layer',1); + self.floatingice_melting_rate=project3d(md,'vector',self.floatingice_melting_rate,'type','node','layer',1); + return self + #}}} + + def setdefaultparameters(self): # {{{ + #default values for melting parameterization + self.mantleconductivity = 2.2 + self.nusselt = 300 + self.dtbg = 11/1000. + self.plumeradius = 100000 + self.topplumedepth = 10000 + self.bottomplumedepth = 1050000 + self.crustthickness = 30000 + self.uppercrustthickness = 14000 + self.uppercrustheat = 1.7*10**-6 + self.lowercrustheat = 0.4*10**-6 + return self + #}}} + + def checkconsistency(self,md,solution,analyses): # {{{ + if 'MasstransportAnalysis' in analyses and not (solution == 'TransientSolution' and md.transient.ismasstransport==0): + md = checkfield(md,'fieldname','basalforcings.groundedice_melting_rate','NaN',1,'timeseries',1) + md = checkfield(md,'fieldname','basalforcings.floatingice_melting_rate','NaN',1,'timeseries',1) + if 'BalancethicknessAnalysis' in analyses: + md = checkfield(md,'fieldname','basalforcings.groundedice_melting_rate','NaN',1,'size',[md.mesh.numberofvertices]) + md = checkfield(md,'fieldname','basalforcings.floatingice_melting_rate','NaN',1,'size',[md.mesh.numberofvertices]) + if 'ThermalAnalysis' in analyses and not (solution == 'TransientSolution' and md.transient.isthermal==0): + md = checkfield(md,'fieldname','basalforcings.groundedice_melting_rate','NaN',1,'timeseries',1) + md = checkfield(md,'fieldname','basalforcings.floatingice_melting_rate','NaN',1,'timeseries',1) + md = checkfield(md,'fieldname','basalforcings.mantleconductivity','>=',0,'numel',1) + md = checkfield(md,'fieldname','basalforcings.nusselt','>',0,'numel',1) + md = checkfield(md,'fieldname','basalforcings.dtbg','>',0,'numel',1) + md = checkfield(md,'fieldname','basalforcings.topplumedepth','>',0,'numel',1) + md = checkfield(md,'fieldname','basalforcings.bottomplumedepth','>',0,'numel',1) + md = checkfield(md,'fieldname','basalforcings.plumex','numel',1) + md = checkfield(md,'fieldname','basalforcings.plumey','numel',1) + md = checkfield(md,'fieldname','basalforcings.crustthickness','>',0,'numel',1) + md = checkfield(md,'fieldname','basalforcings.uppercrustthickness','>',0,'numel',1) + md = checkfield(md,'fieldname','basalforcings.uppercrustheat','>',0,'numel',1) + md = checkfield(md,'fieldname','basalforcings.lowercrustheat','>',0,'numel',1) + return md + # }}} + def marshall(self,prefix,md,fid): # {{{ + yts=md.constants.yts + + WriteData(fid,prefix,'name','md.basalforcings.model','data',4,'format','Integer') + WriteData(fid,prefix,'object',self,'fieldname','floatingice_melting_rate','format','DoubleMat','name','md.basalforcings.floatingice_melting_rate','mattype',1,'scale',1./yts,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts) + WriteData(fid,prefix,'object',self,'fieldname','groundedice_melting_rate','format','DoubleMat','name','md.basalforcings.groundedice_melting_rate','mattype',1,'scale',1./yts,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts) + WriteData(fid,prefix,'object',self,'fieldname','mantleconductivity','format','Double') + WriteData(fid,prefix,'object',self,'fieldname','nusselt','format','Double') + WriteData(fid,prefix,'object',self,'fieldname','dtbg','format','Double') + WriteData(fid,prefix,'object',self,'fieldname','plumeradius','format','Double') + WriteData(fid,prefix,'object',self,'fieldname','topplumedepth','format','Double') + WriteData(fid,prefix,'object',self,'fieldname','bottomplumedepth','format','Double') + WriteData(fid,prefix,'object',self,'fieldname','plumex','format','Double') + WriteData(fid,prefix,'object',self,'fieldname','plumey','format','Double') + WriteData(fid,prefix,'object',self,'fieldname','crustthickness','format','Double') + WriteData(fid,prefix,'object',self,'fieldname','uppercrustthickness','format','Double') + WriteData(fid,prefix,'object',self,'fieldname','uppercrustheat','format','Double') + WriteData(fid,prefix,'object',self,'fieldname','lowercrustheat','format','Double') + # }}} Index: /issm/trunk-jpl/src/m/classes/taoinversion.py =================================================================== --- /issm/trunk-jpl/src/m/classes/taoinversion.py (revision 22266) +++ /issm/trunk-jpl/src/m/classes/taoinversion.py (revision 22267) @@ -6,29 +6,31 @@ from IssmConfig import IssmConfig from marshallcostfunctions import marshallcostfunctions +from supportedcontrols import * +from supportedcostfunctions import * - -class taoinversion: +class taoinversion(object): def __init__(self): - iscontrol = 0 - incomplete_adjoint = 0 - control_parameters = float('NaN') - maxsteps = 0 - maxiter = 0 - fatol = 0 - frtol = 0 - gatol = 0 - grtol = 0 - gttol = 0 - algorithm = '' - cost_functions = float('NaN') - cost_functions_coefficients = float('NaN') - min_parameters = float('NaN') - max_parameters = float('NaN') - vx_obs = float('NaN') - vy_obs = float('NaN') - vz_obs = float('NaN') - vel_obs = float('NaN') - thickness_obs = float('NaN') - surface_obs = float('NaN') + self.iscontrol = 0 + self.incomplete_adjoint = 0 + self.control_parameters = float('NaN') + self.maxsteps = 0 + self.maxiter = 0 + self.fatol = 0 + self.frtol = 0 + self.gatol = 0 + self.grtol = 0 + self.gttol = 0 + self.algorithm = '' + self.cost_functions = float('NaN') + self.cost_functions_coefficients = float('NaN') + self.min_parameters = float('NaN') + self.max_parameters = float('NaN') + self.vx_obs = float('NaN') + self.vy_obs = float('NaN') + self.vz_obs = float('NaN') + self.vel_obs = float('NaN') + self.thickness_obs = float('NaN') + self.surface_obs = float('NaN') + self.setdefaultparameters() def __repr__(self): @@ -98,5 +100,4 @@ #several responses can be used: self.cost_functions=101; - return self @@ -119,5 +120,5 @@ def checkconsistency(self,md,solution,analyses): - if not self.control: + if not self.iscontrol: return md if not IssmConfig('_HAVE_TAO_')[0]: @@ -125,9 +126,9 @@ - num_controls= np.numel(md.inversion.control_parameters) - num_costfunc= np.size(md.inversion.cost_functions,2) + num_controls= np.size(md.inversion.control_parameters) + num_costfunc= np.size(md.inversion.cost_functions) md = checkfield(md,'fieldname','inversion.iscontrol','values',[0, 1]) - md = checkfield(md,'fieldname','inversion.incomplete_adjoint','values',[0, 1]) + md = checkfield(md,'fieldname','inversion.incomplete_adjoint','values',[0,1]) md = checkfield(md,'fieldname','inversion.control_parameters','cell',1,'values',supportedcontrols()) md = checkfield(md,'fieldname','inversion.maxsteps','numel',1,'>=',0) @@ -143,10 +144,10 @@ PETSCMINOR = IssmConfig('_PETSC_MINOR_')[0] if(PETSCMAJOR>3 or (PETSCMAJOR==3 and PETSCMINOR>=5)): - md = checkfield(md,'fieldname','inversion.algorithm','values',{'blmvm','cg','lmvm'}) + md = checkfield(md,'fieldname','inversion.algorithm','values',['blmvm','cg','lmvm']) else: - md = checkfield(md,'fieldname','inversion.algorithm','values',{'tao_blmvm','tao_cg','tao_lmvm'}) + md = checkfield(md,'fieldname','inversion.algorithm','values',['tao_blmvm','tao_cg','tao_lmvm']) - md = checkfield(md,'fieldname','inversion.cost_functions','size',[1, num_costfunc],'values',supportedcostfunctions()) + md = checkfield(md,'fieldname','inversion.cost_functions','size', [num_costfunc],'values',supportedcostfunctions()) md = checkfield(md,'fieldname','inversion.cost_functions_coefficients','size',[md.mesh.numberofvertices, num_costfunc],'>=',0) md = checkfield(md,'fieldname','inversion.min_parameters','size',[md.mesh.numberofvertices, num_controls]) @@ -162,37 +163,37 @@ md = checkfield(md,'fieldname','inversion.vy_obs','size',[md.mesh.numberofvertices],'NaN',1,'Inf',1) - def marshall(self, md, fid): + def marshall(self,prefix,md,fid): - yts=md.constants.yts; - WriteData(fid,prefix,'object',self,'class','inversion','fieldname','iscontrol','format','Boolean') - WriteData(fid,prefix,'name','md.inversion.type','data',1,'format','Integer') - if not self.iscontrol: - return - WriteData(fid,prefix,'object',self,'class','inversion','fieldname','incomplete_adjoint','format','Boolean') - WriteData(fid,prefix,'object',self,'class','inversion','fieldname','maxsteps','format','Integer') - WriteData(fid,prefix,'object',self,'class','inversion','fieldname','maxiter','format','Integer') - WriteData(fid,prefix,'object',self,'class','inversion','fieldname','fatol','format','Double') - WriteData(fid,prefix,'object',self,'class','inversion','fieldname','frtol','format','Double') - WriteData(fid,prefix,'object',self,'class','inversion','fieldname','gatol','format','Double') - WriteData(fid,prefix,'object',self,'class','inversion','fieldname','grtol','format','Double') - WriteData(fid,prefix,'object',self,'class','inversion','fieldname','gttol','format','Double') - WriteData(fid,prefix,'object',self,'class','inversion','fieldname','algorithm','format','String') - WriteData(fid,prefix,'object',self,'class','inversion','fieldname','cost_functions_coefficients','format','DoubleMat','mattype',1) - WriteData(fid,prefix,'object',self,'class','inversion','fieldname','min_parameters','format','DoubleMat','mattype',3) - WriteData(fid,prefix,'object',self,'class','inversion','fieldname','max_parameters','format','DoubleMat','mattype',3) - WriteData(fid,prefix,'object',self,'class','inversion','fieldname','vx_obs','format','DoubleMat','mattype',1,'scale',1./yts) - WriteData(fid,prefix,'object',self,'class','inversion','fieldname','vy_obs','format','DoubleMat','mattype',1,'scale',1./yts) - WriteData(fid,prefix,'object',self,'class','inversion','fieldname','vz_obs','format','DoubleMat','mattype',1,'scale',1./yts) - WriteData(fid,prefix,'object',self,'class','inversion','fieldname','thickness_obs','format','DoubleMat','mattype',1) - WriteData(fid,prefix,'object',self,'class','inversion','fieldname','surface_obs','format','DoubleMat','mattype',1) + yts=md.constants.yts; + WriteData(fid,prefix,'object',self,'class','inversion','fieldname','iscontrol','format','Boolean') + WriteData(fid,prefix,'name','md.inversion.type','data',1,'format','Integer') + if not self.iscontrol: + return + WriteData(fid,prefix,'object',self,'class','inversion','fieldname','incomplete_adjoint','format','Boolean') + WriteData(fid,prefix,'object',self,'class','inversion','fieldname','maxsteps','format','Integer') + WriteData(fid,prefix,'object',self,'class','inversion','fieldname','maxiter','format','Integer') + WriteData(fid,prefix,'object',self,'class','inversion','fieldname','fatol','format','Double') + WriteData(fid,prefix,'object',self,'class','inversion','fieldname','frtol','format','Double') + WriteData(fid,prefix,'object',self,'class','inversion','fieldname','gatol','format','Double') + WriteData(fid,prefix,'object',self,'class','inversion','fieldname','grtol','format','Double') + WriteData(fid,prefix,'object',self,'class','inversion','fieldname','gttol','format','Double') + WriteData(fid,prefix,'object',self,'class','inversion','fieldname','algorithm','format','String') + WriteData(fid,prefix,'object',self,'class','inversion','fieldname','cost_functions_coefficients','format','DoubleMat','mattype',1) + WriteData(fid,prefix,'object',self,'class','inversion','fieldname','min_parameters','format','DoubleMat','mattype',3) + WriteData(fid,prefix,'object',self,'class','inversion','fieldname','max_parameters','format','DoubleMat','mattype',3) + WriteData(fid,prefix,'object',self,'class','inversion','fieldname','vx_obs','format','DoubleMat','mattype',1,'scale',1./yts) + WriteData(fid,prefix,'object',self,'class','inversion','fieldname','vy_obs','format','DoubleMat','mattype',1,'scale',1./yts) + WriteData(fid,prefix,'object',self,'class','inversion','fieldname','vz_obs','format','DoubleMat','mattype',1,'scale',1./yts) + WriteData(fid,prefix,'object',self,'class','inversion','fieldname','thickness_obs','format','DoubleMat','mattype',1) + WriteData(fid,prefix,'object',self,'class','inversion','fieldname','surface_obs','format','DoubleMat','mattype',1) - #process control parameters - num_control_parameters = np.numel(self.control_parameters) - WriteData(fid,prefix,'object',self,'fieldname','control_parameters','format','StringArray') - WriteData(fid,prefix,'data',num_control_parameters,'name','md.inversion.num_control_parameters','format','Integer') + #process control parameters + num_control_parameters = np.size(self.control_parameters) + WriteData(fid,prefix,'object',self,'fieldname','control_parameters','format','StringArray') + WriteData(fid,prefix,'data',num_control_parameters,'name','md.inversion.num_control_parameters','format','Integer') - #process cost functions - num_cost_functions = np.size(self.cost_functions,2) - data= marshallcostfunctions(self.cost_functions) - WriteData(fid,prefix,'data',data,'name','md.inversion.cost_functions','format','StringArray') - WriteData(fid,prefix,'data',num_cost_functions,'name','md.inversion.num_cost_functions','format','Integer') + #process cost functions + num_cost_functions = np.size(self.cost_functions) + data= marshallcostfunctions(self.cost_functions) + WriteData(fid,prefix,'data',data,'name','md.inversion.cost_functions','format','StringArray') + WriteData(fid,prefix,'data',num_cost_functions,'name','md.inversion.num_cost_functions','format','Integer') Index: /issm/trunk-jpl/src/m/classes/thermal.py =================================================================== --- /issm/trunk-jpl/src/m/classes/thermal.py (revision 22266) +++ /issm/trunk-jpl/src/m/classes/thermal.py (revision 22267) @@ -101,12 +101,16 @@ if 'EnthalpyAnalysis' in analyses and md.thermal.isenthalpy and md.mesh.dimension()==3: - pos=np.where(~np.isnan(md.thermal.spctemperature[0:md.mesh.numberofvertices])) + TEMP = md.thermal.spctemperature[:-1].flatten(-1) + pos=np.where(~np.isnan(TEMP)) try: spccol=np.size(md.thermal.spctemperature,1) except IndexError: spccol=1 - replicate=np.tile(md.geometry.surface-md.mesh.z,(spccol)) - control=md.materials.meltingpoint-md.materials.beta*md.materials.rho_ice*md.constants.g*replicate - md = checkfield(md,'fieldname','thermal.spctemperature','field',md.thermal.spctemperature[pos],'<=',control[pos],'message',"spctemperature should be below the adjusted melting point") + + replicate=np.tile(md.geometry.surface-md.mesh.z,(spccol)).flatten(-1) + + control=md.materials.meltingpoint-md.materials.beta*md.materials.rho_ice*md.constants.g*replicate+10**-5 + + md = checkfield(md,'fieldname','thermal.spctemperature','field',md.thermal.spctemperature.flatten(-1)[pos],'<=',control[pos],'message',"spctemperature should be below the adjusted melting point") md = checkfield(md,'fieldname','thermal.isenthalpy','numel',[1],'values',[0,1]) md = checkfield(md,'fieldname','thermal.isdynamicbasalspc','numel',[1],'values',[0,1]); Index: /issm/trunk-jpl/src/m/classes/transient.py =================================================================== --- /issm/trunk-jpl/src/m/classes/transient.py (revision 22266) +++ /issm/trunk-jpl/src/m/classes/transient.py (revision 22267) @@ -12,5 +12,5 @@ def __init__(self): # {{{ - self.issmb = False + self.issmb = False self.ismasstransport = False self.isstressbalance = False @@ -23,7 +23,7 @@ self.ishydrology = False self.isslr = False + self.iscoupler = False + self.amr_frequency = 0 self.isoceancoupling = False - self.iscoupler = False - amr_frequency = 0 self.requested_outputs = [] @@ -76,4 +76,24 @@ self.iscoupler = False self.amr_frequency = 0 + + #default output + self.requested_outputs=[] + return self + #}}} + def deactivateall(self):#{{{ + self.issmb = False + self.ismasstransport = False + self.isstressbalance = False + self.isthermal = False + self.isgroundingline = False + self.isgia = False + self.isesa = False + self.isdamageevolution = False + self.ismovingfront = False + self.ishydrology = False + self.isslr = False + self.isoceancoupling = False + self.iscoupler = False + self.amr_frequency = 0 #default output Index: /issm/trunk-jpl/src/m/consistency/checkfield.py =================================================================== --- /issm/trunk-jpl/src/m/consistency/checkfield.py (revision 22266) +++ /issm/trunk-jpl/src/m/consistency/checkfield.py (revision 22267) @@ -84,5 +84,5 @@ if options.exist('numel'): fieldnumel=options.getfieldvalue('numel') - if np.size(field) not in fieldnumel: + if (type(fieldnumel) == int and np.size(field) != fieldnumel) or (type(fieldnumel) == list and np.size(field) not in fieldnumel): if len(fieldnumel)==1: md = md.checkmessage(options.getfieldvalue('message',\ @@ -101,4 +101,5 @@ "NaN values found in field '%s'" % fieldname)) + #check Inf if options.getfieldvalue('Inf',0): @@ -106,4 +107,5 @@ md = md.checkmessage(options.getfieldvalue('message',\ "Inf values found in field '%s'" % fieldname)) + #check cell @@ -129,11 +131,26 @@ #check greater if options.exist('>='): - lowerbound=options.getfieldvalue('>=') - if np.any(field=') + field2 = np.reshape(field,(np.prod(np.shape(field)),1),order='F').copy() + if options.getfieldvalue('timeseries',0): + field2 = np.reshape(field[:-1],np.prod(np.shape(field[:-1])),1) + + if options.getfieldvalue('singletimeseries',0): + field2 = np.reshape(field[0],np.prod(np.shape(field[0])),1) + + if np.any(field2'): lowerbound=options.getfieldvalue('>') - if np.any(field<=lowerbound): + field2 = np.reshape(field,(np.prod(np.shape(field)),1),order='F').copy() + if options.getfieldvalue('timeseries',0): + field2 = np.reshape(field[:-1],np.prod(np.shape(field[:-1])),1) + + if options.getfieldvalue('singletimeseries',0): + field2 = np.reshape(field[0],np.prod(np.shape(field[0])),1) + + if np.any(field2<=lowerbound): md = md.checkmessage(options.getfieldvalue('message',\ "field '%s' should have values above %d" % (fieldname,lowerbound))) @@ -142,10 +159,24 @@ if options.exist('<='): upperbound=options.getfieldvalue('<=') - if np.any(field>upperbound): + field2 = np.reshape(field,(np.prod(np.shape(field)),1),order='F').copy() + if options.getfieldvalue('timeseries',0): + field2 = np.reshape(field[:-1],np.prod(np.shape(field[:-1])),1) + + if options.getfieldvalue('singletimeseries',0): + field2 = np.reshape(field[0],np.prod(np.shape(field[0])),1) + + if np.any(field2>upperbound): md = md.checkmessage(options.getfieldvalue('message',\ "field '%s' should have values below %d" % (fieldname,upperbound))) if options.exist('<'): upperbound=options.getfieldvalue('<') - if np.any(field>=upperbound): + field2 = np.reshape(field,(np.prod(np.shape(field)),1),order='F').copy() + if options.getfieldvalue('timeseries',0): + field2 = np.reshape(field[:-1],np.prod(np.shape(field[:-1])),1) + + if options.getfieldvalue('singletimeseries',0): + field2 = np.reshape(field[0],np.prod(np.shape(field[0])),1) + + if np.any(field2>=upperbound): md = md.checkmessage(options.getfieldvalue('message',\ "field '%s' should have values below %d" % (fieldname,upperbound))) @@ -164,13 +195,13 @@ #Check forcings (size and times) if options.getfieldvalue('timeseries',0): - if np.size(field,0)==md.mesh.numberofvertices: + if np.size(field,0)==md.mesh.numberofvertices or np.size(field,0)==md.mesh.numberofelements: if np.ndim(field)>1 and not np.size(field,1)==1: md = md.checkmessage(options.getfieldvalue('message',\ "field '%s' should have only one column as there are md.mesh.numberofvertices lines" % fieldname)) - elif np.size(field,0)==md.mesh.numberofvertices+1 or np.size(field,0)==2: - if not all(field[-1,:]==np.sort(field[-1,:])): + elif np.size(field,0)==md.mesh.numberofvertices+1 or np.size(field,0)==md.mesh.numberofelements+1: + if np.ndim(field) > 1 and not all(field[-1,:]==np.sort(field[-1,:])): md = md.checkmessage(options.getfieldvalue('message',\ "field '%s' columns should be sorted chronologically" % fieldname)) - if any(field[-1,0:-1]==field[-1,1:]): + if np.ndim(field) > 1 and any(field[-1,0:-1]==field[-1,1:]): md = md.checkmessage(options.getfieldvalue('message',\ "field '%s' columns must not contain duplicate timesteps" % fieldname)) @@ -198,2 +229,3 @@ return md + Index: /issm/trunk-jpl/src/m/dev/devpath.py =================================================================== --- /issm/trunk-jpl/src/m/dev/devpath.py (revision 22266) +++ /issm/trunk-jpl/src/m/dev/devpath.py (revision 22267) @@ -31,4 +31,5 @@ #Manual imports for commonly used functions from plotmodel import plotmodel +from runme import runme #c = get_ipython().config Index: /issm/trunk-jpl/src/m/geometry/NowickiProfile.py =================================================================== --- /issm/trunk-jpl/src/m/geometry/NowickiProfile.py (revision 22267) +++ /issm/trunk-jpl/src/m/geometry/NowickiProfile.py (revision 22267) @@ -0,0 +1,44 @@ +import numpy as np + +def NowickiProfile(x): + """ + NOWICKIPROFILE - Create profile at the transition zone based on Sophie Nowicki's thesis + + Usage: + [b h] = NowickiProfile(x) + + - h = ice thickness + - b = ice base + - x = along flow coordinate + """ + #Constant for theoretical profile + delta = 0.1 #ratio of water density and ice density -1 + hg = 1. #ice thickness at grounding line + sea = hg / (1+delta) #sea level + lamda = 0.1 #ration of deviatoric stress and water pressure + beta = 5. #friction coefficient + ms = 0.005 #surface accumulation rat + mu = 5. #viscosity + q = 0.801 #ice mass flux + + #mesh parameters + b = np.zeros((np.size(x),)) + h = np.zeros((np.size(x),)) + s = np.zeros((np.size(x),)) + + #upstream of the GL + for i in range(np.size(x) / 2): + ss = np.roots([1, 4 * lamda * beta, 0, 0, 6 * lamda * ms * x[i]**2 + + 12 * lamda * q * x[i] - hg**4 - 4 * lamda * beta * hg**3]) + for j in range(4): + if (np.isreal(ss[j]) > 0) and (np.imag(ss[j]) == 0): + s[i] = ss[j] + h[i] = s[i] + b[i] = 0. + + #downstream of the GL + for i in range(np.size(x) / 2, np.size(x)): + h[i] = (x[i] / (4. * (delta+1) * q) + hg**(-2))**(-0.5) # ice thickness for ice shelf from (3.1) + b[i] = sea - h[i] * (1. / (1+delta)) + + return [b, h, sea] Index: /issm/trunk-jpl/src/m/mech/newforcing.py =================================================================== --- /issm/trunk-jpl/src/m/mech/newforcing.py (revision 22267) +++ /issm/trunk-jpl/src/m/mech/newforcing.py (revision 22267) @@ -0,0 +1,34 @@ +import numpy as np + +def newforcing(t0,t1,deltaT,f0,f1,nodes): + ''' +FUNCTION NEWFORCING Build forcing that extends temporally from t0 to t1, and in magnitude from f0 to f1. Equal time + and magnitude spacing. + + Usage: forcing=newforcing(t0,t1,deltaT,f0,f1,nodes); + Where: + t0:t1: time interval. + deltaT: time step + f0:f1: magnitude interval. + nodes: number of vertices where we have a temporal forcing + + Example: + md.smb.mass_balance=newforcing(md.timestepping.start_time,md.timestepping.final_time, + md.timestepping.time_step,-1,+2,md.mesh.numberofvertices) + ''' + #Number of time steps: + nsteps = (t1 - t0) / deltaT + 1 + + #delta forcing: + deltaf = (f1 - f0) / (nsteps - 1) + + #creates times: + times = np.arange(t0,t1+deltaT,deltaT) #Add deltaT to fix python/matlab discrepency + + #create forcing: + forcing = np.arange(f0,f1+deltaf,deltaf)#Add deltaf to fix python/matlab discrepency + + #replicate for all nodes + forcing = np.tile(forcing, (nodes+1,1)) + forcing[-1,:] = times + return forcing Index: /issm/trunk-jpl/src/m/mesh/bamg.py =================================================================== --- /issm/trunk-jpl/src/m/mesh/bamg.py (revision 22266) +++ /issm/trunk-jpl/src/m/mesh/bamg.py (revision 22267) @@ -1,5 +1,5 @@ import os.path import numpy as np -from mesh2d import mesh2d +from mesh2dvertical import * from collections import OrderedDict from pairoptions import pairoptions @@ -80,7 +80,10 @@ #Check that file exists domainfile=options.getfieldvalue('domain') - if not os.path.exists(domainfile): - raise IOError("bamg error message: file '%s' not found" % domainfile) - domain=expread(domainfile) + if type(domainfile) == str: + if not os.path.exists(domainfile): + raise IOError("bamg error message: file '%s' not found" % domainfile) + domain=expread(domainfile) + else: + domain=domainfile #Build geometry @@ -89,16 +92,16 @@ #Check that the domain is closed - if (domaini['x'][0]!=domaini['x'][-1] or domaini['y'][0]!=domaini['y'][-1]): + if (domaini.x[0]!=domaini.x[-1] or domaini.y[0]!=domaini.y[-1]): raise RuntimeError("bamg error message: all contours provided in ''domain'' should be closed") #Checks that all holes are INSIDE the principle domain outline if i: - flags=ContourToNodes(domaini['x'],domaini['y'],domainfile,0)[0] + flags=ContourToNodes(domaini.x,domaini.y,domainfile,0)[0] if np.any(np.logical_not(flags)): raise RuntimeError("bamg error message: All holes should be strictly inside the principal domain") #Add all points to bamg_geometry - nods=domaini['nods']-1 #the domain are closed 0=end - bamg_geometry.Vertices=np.vstack((bamg_geometry.Vertices,np.vstack((domaini['x'][0:nods],domaini['y'][0:nods],np.ones((nods)))).T)) + nods=domaini.nods-1 #the domain are closed 0=end + bamg_geometry.Vertices=np.vstack((bamg_geometry.Vertices,np.vstack((domaini.x[0:nods],domaini.y[0:nods],np.ones((nods)))).T)) bamg_geometry.Edges =np.vstack((bamg_geometry.Edges,np.vstack((np.arange(count+1,count+nods+1),np.hstack((np.arange(count+2,count+nods+1),count+1)),1.*np.ones((nods)))).T)) if i: @@ -115,4 +118,10 @@ #update counter count+=nods + + if options.getfieldvalue('vertical',0): + if np.size(options.getfieldvalue('Markers',[]))!=np.size(bamg_geometry.Edges,0): + raise RuntimeError('for 2d vertical mesh, ''Markers'' option is required, and should be of size ' + str(np.size(bamg_geometry.Edges,0))) + if np.size(options.getfieldvalue('Markers',[]))==np.size(bamg_geometry.Edges,0): + bamg_geometry.Edges[:,2]=options.getfieldvalue('Markers') #take care of rifts @@ -323,21 +332,62 @@ # plug results onto model + if options.getfieldvalue('vertical',0): + md.mesh=mesh2dvertical() + md.mesh.x=bamgmesh_out['Vertices'][:,0].copy() + md.mesh.y=bamgmesh_out['Vertices'][:,1].copy() + md.mesh.elements=bamgmesh_out['Triangles'][:,0:3].astype(int) + md.mesh.edges=bamgmesh_out['IssmEdges'].astype(int) + md.mesh.segments=bamgmesh_out['IssmSegments'][:,0:3].astype(int) + md.mesh.segmentmarkers=bamgmesh_out['IssmSegments'][:,3].astype(int) + + #Fill in rest of fields: + md.mesh.numberofelements=np.size(md.mesh.elements,axis=0) + md.mesh.numberofvertices=np.size(md.mesh.x) + md.mesh.numberofedges=np.size(md.mesh.edges,axis=0) + md.mesh.vertexonboundary=np.zeros(md.mesh.numberofvertices,bool) + md.mesh.vertexonboundary[md.mesh.segments[:,0:2]-1]=True + + #Now, build the connectivity tables for this mesh. Doubled in matlab for some reason + md.mesh.vertexonboundary=np.zeros(md.mesh.numberofvertices,) + md.mesh.vertexonboundary[md.mesh.segments[:,0:2]-1]=1 + + elif options.getfieldvalue('3dsurface',0): + md.mesh=mesh3dsurface() + md.mesh.x=bamgmesh_out['Vertices'][:,0].copy() + md.mesh.y=bamgmesh_out['Vertices'][:,1].copy() + md.mesh.z=md.mesh.x + md.mesh.z[:]=0 + md.mesh.elements=bamgmesh_out['Triangles'][:,0:3].astype(int) + md.mesh.edges=bamgmesh_out['IssmEdges'].astype(int) + md.mesh.segments=bamgmesh_out['IssmSegments'][:,0:3].astype(int) + md.mesh.segmentmarkers=bamgmesh_out['IssmSegments'][:,3].astype(int) + + #Fill in rest of fields: + md.mesh.numberofelements=np.size(md.mesh.elements,axis=0) + md.mesh.numberofvertices=np.size(md.mesh.x) + md.mesh.numberofedges=np.size(md.mesh.edges,axis=0) + md.mesh.vertexonboundary=np.zeros(md.mesh.numberofvertices,bool) + md.mesh.vertexonboundary[md.mesh.segments[:,0:2]-1]=True + + else: + md.mesh=mesh2d() + md.mesh.x=bamgmesh_out['Vertices'][:,0].copy() + md.mesh.y=bamgmesh_out['Vertices'][:,1].copy() + md.mesh.elements=bamgmesh_out['Triangles'][:,0:3].astype(int) + md.mesh.edges=bamgmesh_out['IssmEdges'].astype(int) + md.mesh.segments=bamgmesh_out['IssmSegments'][:,0:3].astype(int) + md.mesh.segmentmarkers=bamgmesh_out['IssmSegments'][:,3].astype(int) + + #Fill in rest of fields: + md.mesh.numberofelements=np.size(md.mesh.elements,axis=0) + md.mesh.numberofvertices=np.size(md.mesh.x) + md.mesh.numberofedges=np.size(md.mesh.edges,axis=0) + md.mesh.vertexonboundary=np.zeros(md.mesh.numberofvertices,bool) + md.mesh.vertexonboundary[md.mesh.segments[:,0:2]-1]=True + + #Bamg private fields md.private.bamg=OrderedDict() md.private.bamg['mesh']=bamgmesh(bamgmesh_out) md.private.bamg['geometry']=bamggeom(bamggeom_out) - md.mesh = mesh2d() - md.mesh.x=bamgmesh_out['Vertices'][:,0].copy() - md.mesh.y=bamgmesh_out['Vertices'][:,1].copy() - md.mesh.elements=bamgmesh_out['Triangles'][:,0:3].astype(int) - md.mesh.edges=bamgmesh_out['IssmEdges'].astype(int) - md.mesh.segments=bamgmesh_out['IssmSegments'][:,0:3].astype(int) - md.mesh.segmentmarkers=bamgmesh_out['IssmSegments'][:,3].astype(int) - - #Fill in rest of fields: - md.mesh.numberofelements=np.size(md.mesh.elements,axis=0) - md.mesh.numberofvertices=np.size(md.mesh.x) - md.mesh.numberofedges=np.size(md.mesh.edges,axis=0) - md.mesh.vertexonboundary=np.zeros(md.mesh.numberofvertices,bool) - md.mesh.vertexonboundary[md.mesh.segments[:,0:2]-1]=True md.mesh.elementconnectivity=md.private.bamg['mesh'].ElementConnectivity md.mesh.elementconnectivity[np.nonzero(np.isnan(md.mesh.elementconnectivity))]=0 Index: /issm/trunk-jpl/src/m/mesh/bamgflowband.py =================================================================== --- /issm/trunk-jpl/src/m/mesh/bamgflowband.py (revision 22267) +++ /issm/trunk-jpl/src/m/mesh/bamgflowband.py (revision 22267) @@ -0,0 +1,47 @@ +import numpy as np +from model import * +from collections import OrderedDict +from bamg import * +from mesh2dvertical import * + +def bamgflowband(md,x,surf,base,*args): + """ + BAMGFLOWBAND - create flowband mesh with bamg + + Usage: + md=bamgflowband(md,x,surf,base,OPTIONS) + + surf and bed are the surface elevation and base for each x provided + x must be increasing + OPTIONS are bamg options + + Example: + x =np.arrange(1,3001,100) + h=linspace(1000,300,numel(x)) + b=-917/1023*h + md=bamgflowband(model,b+h,b,'hmax',80,'vertical',1,'Markers',m) + """ + + #Write expfile with domain outline + A = OrderedDict() + A.x = np.concatenate((x,np.flipud(x),[x[0]])) + A.y = np.concatenate((base,np.flipud(surf),[base[0]])) + A.nods = np.size(A.x) + + #markers: + m = np.ones((np.size(A.x)-1,)) # base = 1 + m[np.size(x) - 1] = 2 # right side = 2 + m[np.size(x):2 * np.size(x) - 1] = 3 # top surface = 3 + m[2 * np.size(x) - 1] = 4 # left side = 4 + + #mesh domain + md = bamg(model(),'domain',[A],'vertical',1,'Markers',m,*args) + #print md.mesh.numberofvertices + + #Deal with vertices on bed + md.mesh.vertexonbase = np.zeros((md.mesh.numberofvertices,)) + md.mesh.vertexonbase[np.where(md.mesh.vertexflags(1))] = 1 + md.mesh.vertexonsurface = np.zeros((md.mesh.numberofvertices,)) + md.mesh.vertexonsurface[np.where(md.mesh.vertexflags(3))] = 1 + + return md Index: /issm/trunk-jpl/src/m/print/printmodel.py =================================================================== --- /issm/trunk-jpl/src/m/print/printmodel.py (revision 22267) +++ /issm/trunk-jpl/src/m/print/printmodel.py (revision 22267) @@ -0,0 +1,109 @@ + +import numpy as np +#from model import * +#from socket import gethostname +#from bamg import * +#from setmask import * +#from parameterize import * +#from setflowequation import * +#from solve import * +from pairoptions import * + +def printmodel(filename,format,*args): +''' PRINTMODEL - save an image of current figure + + filename: output name of image file (no extension) + format: image format (ex: 'tiff','jpg','pdf') + + List of options to printfmodel: + + figure: number of figure to print (default: current figure) + resolution: use higher resolution to anti-alias (default 2) + margin: add margin around final image + marginsize: size of margin around final image (default 5) + frame: add frame around final image + framesize: size of frame around final image (default 5) + framecolor: color of frame around final image (default 'black') + trim: trim empty space around image (default 'off') + hardcopy: 'off' to impose MATLAB to use the same colors (default 'off') + + Usage: + printmodel(filename,format,varargin); + + Examples: + printmodel('image','tiff') + printmodel('image','eps','margin','on','frame','on','hardcopy','on') +''' + +#get options: +options = pairoptions(*args) + +#set defaults +options = addfielddefault(options,'figure','gcf') +options = addfielddefault(options,'format','tiff') +options = addfielddefault(options,'resolution',1) +options = addfielddefault(options,'margin','on') +options = addfielddefault(options,'marginsize',25) +options = addfielddefault(options,'frame','on') +options = addfielddefault(options,'framesize',3) +options = addfielddefault(options,'framecolor','black') +options = addfielddefault(options,'trim','on') +options = addfielddefault(options,'hardcopy','off') + +#get fig: +fig = getfieldvalue(options,'figure') +if len(fig) == 1: + fig=gcf +else: + figure(fig) + fig=gcf + +#In auto mode, MATLAB prints the figure the same size as it appears on the computer screen, centered on the page +set(fig, 'PaperPositionMode', 'auto'); + +#InvertHardcopy off imposes MATLAB to use the same colors +set(fig, 'InvertHardcopy', getfieldvalue(options,'hardcopy')); + +#we could have several formats, as a cell array of strings. +formats=format; +if ~iscell(formats), + formats={formats}; +end + +#loop on formats: +for i=1:length(formats), + format=formats{i}; + + #Use higher resolution to anti-alias and use zbuffer to have smooth colors + print(fig, '-zbuffer','-dtiff',['-r' num2str(get(0,'ScreenPixelsPerInch')*getfieldvalue(options,'resolution'))],filename); + + #some trimming involved? + if ~strcmpi(format,'pdf'), + if strcmpi(getfieldvalue(options,'trim'),'on'), + system(['convert -trim ' filename '.tif ' filename '.tif']); + end + end + + #margin? + if ~strcmpi(format,'pdf'), + if strcmpi(getfieldvalue(options,'margin'),'on'), + marginsize=getfieldvalue(options,'marginsize'); + system(['convert -border ' num2str(marginsize) 'x' num2str(marginsize) ' -bordercolor "white" ' filename '.tif ' filename '.tif']); + end + end + + #frame? + if ~strcmpi(format,'pdf'), + if strcmpi(getfieldvalue(options,'frame'),'on'), + framesize=getfieldvalue(options,'framesize'); + framecolor=getfieldvalue(options,'framecolor'); + system(['convert -border ' num2str(framesize) 'x' num2str(framesize) ' -bordercolor "' framecolor '" ' filename '.tif ' filename '.tif']); + end + end + + #convert image to correct format + if ~strcmpi(format,'tiff') & ~strcmpi(format,'tif'), + system(['convert ' filename '.tif ' filename '.' format]); + delete([ filename '.tif']); + end +end Index: /issm/trunk-jpl/src/m/solve/WriteData.py =================================================================== --- /issm/trunk-jpl/src/m/solve/WriteData.py (revision 22266) +++ /issm/trunk-jpl/src/m/solve/WriteData.py (revision 22267) @@ -51,7 +51,7 @@ yts = options.getfieldvalue('yts') if np.ndim(data) > 1: - data[-1,:] = data[-1,:]*yts - else: - data[-1] = data[-1]*yts + data[-1,:] = yts*data[-1,:] + else: + data[-1] = yts*data[-1] #Step 1: write the enum to identify this record uniquely Index: /issm/trunk-jpl/src/m/solve/parseresultsfromdisk.py =================================================================== --- /issm/trunk-jpl/src/m/solve/parseresultsfromdisk.py (revision 22266) +++ /issm/trunk-jpl/src/m/solve/parseresultsfromdisk.py (revision 22267) @@ -175,6 +175,4 @@ if fieldname=='BalancethicknessThickeningRate': field = field*yts - elif fieldname=='Time': - field = field/yts elif fieldname=='HydrologyWaterVx': field = field*yts @@ -191,14 +189,33 @@ elif fieldname=='BasalforcingsGroundediceMeltingRate': field = field*yts + elif fieldname=='BasalforcingsFloatingiceMeltingRate': + field = field*yts elif fieldname=='TotalFloatingBmb': - field = field/10.**12.*yts #(GigaTon/year) + field = field/10.**12*yts #(GigaTon/year) elif fieldname=='TotalGroundedBmb': - field = field/10.**12.*yts #(GigaTon/year) + field = field/10.**12*yts #(GigaTon/year) elif fieldname=='TotalSmb': - field = field/10.**12.*yts #(GigaTon/year) + field = field/10.**12*yts #(GigaTon/year) elif fieldname=='SmbMassBalance': field = field*yts + elif fieldname=='SmbPrecipitation': + field = field*yts + elif fieldname=='SmbRunoff': + field = field*yts + elif fieldname=='SmbEC': + field = field*yts + elif fieldname=='SmbAccumulation': + field = field*yts + elif fieldname=='SmbMelt': + field = field*yts + elif fieldname=='SmbDz_add': + field = field*yts + elif fieldname=='SmbM_add': + field = field*yts elif fieldname=='CalvingCalvingrate': field = field*yts + + if time !=-9999: + time = time/yts saveres=OrderedDict() Index: /issm/trunk-jpl/test/NightlyRun/test2424.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test2424.py (revision 22267) +++ /issm/trunk-jpl/test/NightlyRun/test2424.py (revision 22267) @@ -0,0 +1,51 @@ +#Test Name: SquareSheetShelfGroundingLine2dAggressive. From test424, with sea level increasing. +import numpy as np +from model import * +from socket import gethostname +from triangle import * +from setmask import * +from parameterize import * +from setflowequation import * +from solve import * +from newforcing import * + +md = triangle(model(),'../Exp/Square.exp',150000.) +md = setmask(md,'../Exp/SquareShelf.exp','') +md = parameterize(md,'../Par/SquareSheetShelf.py') +md = setflowequation(md,'SSA','all') +md.initialization.vx[:] = 0. +md.initialization.vy[:] = 0. +md.smb.mass_balance[:] = 0. + +md.geometry.base = -700. - np.abs(md.mesh.y-500000.) / 1000. +md.geometry.bed = -700. - np.abs(md.mesh.y-500000.) / 1000. +md.geometry.thickness[:] = 1000. +md.geometry.surface = md.geometry.base + md.geometry.thickness + +md.transient.isstressbalance = 0 +md.transient.isgroundingline = 1 +md.transient.isthermal = 0 +md.groundingline.migration = 'AggressiveMigration' +md.transient.requested_outputs = ['IceVolume','IceVolumeAboveFloatation','Sealevel'] + +md.timestepping.time_step = .1 +md.slr.sealevel = newforcing(md.timestepping.start_time, md.timestepping.final_time, + md.timestepping.time_step, -200., 200., md.mesh.numberofvertices) + +md.cluster = generic('name',gethostname(),'np',3) +md = solve(md,'Transient') + +#we are checking that the grounding line position is near the theorical one, which is the 0 contour level +#of surface - sealevel - (1-di)* thickness + +nsteps = len(md.results.TransientSolution) +field_names = [] +field_tolerances = [] +field_values = [] +#time is off by the year constant +for i in range(nsteps): + field_names.append('Time-' + str(md.results.TransientSolution[i].time) + + '-yr-ice_levelset-S-sl-(1-di)*H') + field_tolerances.append(1e-12) + field_values.append(md.results.TransientSolution[i].MaskGroundediceLevelset.reshape(-1,) - (md.geometry.surface - md.results.TransientSolution[i].Sealevel.reshape(-1,) - (1 - md.materials.rho_ice / md.materials.rho_water) * md.geometry.thickness)) + Index: /issm/trunk-jpl/test/NightlyRun/test2425.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test2425.py (revision 22267) +++ /issm/trunk-jpl/test/NightlyRun/test2425.py (revision 22267) @@ -0,0 +1,50 @@ +#Test Name: SquareSheetShelfGroundingLine2dSoft +import numpy as np +from model import * +from socket import gethostname +from triangle import * +from setmask import * +from parameterize import * +from setflowequation import * +from solve import * + +md = triangle(model(),'../Exp/Square.exp',150000.) +md = setmask(md,'../Exp/SquareShelf.exp','') +md = parameterize(md,'../Par/SquareSheetShelf.py') +md = setflowequation(md,'SSA','all') +md.initialization.vx[:] = 0. +md.initialization.vy[:] = 0. +md.geometry.base = -700. - (md.mesh.y - 500000.) / 1000. +md.geometry.bed = -700. - (md.mesh.y - 500000.) / 1000. +md.geometry.thickness[:] = 1300. +md.geometry.surface = md.geometry.base + md.geometry.thickness +md.transient.isstressbalance = 1 +md.transient.isgroundingline = 1 +md.groundingline.migration = 'AggressiveMigration' + +md.timestepping.time_step = .1 +md.timestepping.final_time = 1 + +md.cluster = generic('name',gethostname(),'np',3) +md = solve(md,'Transient') +vel1 = md.results.TransientSolution[-1].Vel + +#get same results with offset in bed and sea level: +md.geometry.base = -700. - (md.mesh.y - 500000.) / 1000. +md.geometry.bed = -700. - (md.mesh.y - 500000.) / 1000. +md.geometry.thickness[:] = 1300. +md.geometry.surface = md.geometry.base + md.geometry.thickness + +md.geometry.base = md.geometry.base + 1000 +md.geometry.bed = md.geometry.bed + 1000 +md.geometry.surface = md.geometry.surface + 1000 +md.slr.sealevel = 1000 * np.ones((md.mesh.numberofvertices,)) + +md = solve(md,'Transient','checkconsistency','no') +vel2 = md.results.TransientSolution[-1].Vel + +#Fields and tolerances to track changes +field_names = ['Vel','Veloffset'] +field_tolerances = [1e-13,1e-13] +field_values = [vel1,vel2] + Index: /issm/trunk-jpl/test/NightlyRun/test243.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test243.py (revision 22267) +++ /issm/trunk-jpl/test/NightlyRun/test243.py (revision 22267) @@ -0,0 +1,72 @@ +#Test Name: SquareShelfSMBGemb +import numpy as np +import scipy.io as spio +from model import * +from socket import gethostname +from triangle import * +from setmask import * +from parameterize import * +from setflowequation import * +from solve import * +from SMBgemb import * + +md = triangle(model(),'../Exp/Square.exp',200000.) +md = setmask(md,'all','') +md = parameterize(md,'../Par/SquareShelf.py') +md = setflowequation(md,'SSA','all') +md.materials.rho_ice = 910 +md.cluster = generic('name',gethostname(),'np',3) + +#Use of Gemb method for SMB computation +md.smb = SMBgemb() +md.smb.setdefaultparameters(md.mesh,md.geometry) + +#load hourly surface forcing date from 1979 to 2009: +inputs = spio.loadmat('../Data/gemb_input.mat',squeeze_me=True) + +#setup the inputs: +md.smb.Ta = np.append(np.tile(np.conjugate(inputs['Ta0']),(md.mesh.numberofelements,1)),np.conjugate([inputs['dateN']]),axis=0) +md.smb.V = np.append(np.tile(np.conjugate(inputs['V0']),(md.mesh.numberofelements,1)),np.conjugate([inputs['dateN']]),axis=0) +md.smb.dswrf = np.append(np.tile(np.conjugate(inputs['dsw0']),(md.mesh.numberofelements,1)),np.conjugate([inputs['dateN']]),axis=0) +md.smb.dlwrf = np.append(np.tile(np.conjugate(inputs['dlw0']),(md.mesh.numberofelements,1)),np.conjugate([inputs['dateN']]),axis=0) +md.smb.P = np.append(np.tile(np.conjugate(inputs['P0']),(md.mesh.numberofelements,1)),np.conjugate([inputs['dateN']]),axis=0) +md.smb.eAir = np.append(np.tile(np.conjugate(inputs['eAir0']),(md.mesh.numberofelements,1)),np.conjugate([inputs['dateN']]),axis=0) +md.smb.pAir = np.append(np.tile(np.conjugate(inputs['pAir0']),(md.mesh.numberofelements,1)),np.conjugate([inputs['dateN']]),axis=0) +md.smb.pAir = np.append(np.tile(np.conjugate(inputs['pAir0']),(md.mesh.numberofelements,1)),np.conjugate([inputs['dateN']]),axis=0) +md.smb.Vz = np.tile(np.conjugate(inputs['LP']['Vz']),(md.mesh.numberofelements,1)).flatten() +md.smb.Tz = np.tile(np.conjugate(inputs['LP']['Tz']),(md.mesh.numberofelements,1)).flatten() +md.smb.Tmean = np.tile(np.conjugate(inputs['LP']['Tmean']),(md.mesh.numberofelements,1)).flatten() +md.smb.C = np.tile(np.conjugate(inputs['LP']['C']),(md.mesh.numberofelements,1)).flatten() + +#smb settings +md.smb.requested_outputs = ['SmbDz','SmbT','SmbD','SmbRe','SmbGdn','SmbGsp','SmbEC','SmbA','SmbMassBalance'] + +#only run smb core: +md.transient.isstressbalance = 0 +md.transient.ismasstransport = 0 +md.transient.isthermal = 0 + +#time stepping: +md.timestepping.start_time = 1965. +md.timestepping.final_time = 1966. +md.timestepping.time_step = 1. / 365.0 +md.timestepping.interp_forcings = 0. + +#Run transient +md = solve(md,'Transient') + +#Fields and tolerances to track changes +field_names = ['SmbDz','SmbT' ,'SmbD' ,'SmbRe','SmbGdn','SmbGsp','SmbA' ,'SmbEC','SmbMassBalance'] +field_tolerances = [5e-4,5e-5,0.0006,0.0002,1e-5,0.0003,2e-5,2e-7,1e-7] +#shape is different in python solution (fixed using reshape) which can cause test failure: +field_values = [ + md.results.TransientSolution[-1].SmbDz[0,0:240].reshape(1,-1), + md.results.TransientSolution[-1].SmbT[0,0:240].reshape(1,-1), + md.results.TransientSolution[-1].SmbD[0,0:240].reshape(1,-1), + md.results.TransientSolution[-1].SmbRe[0,0:240].reshape(1,-1), + md.results.TransientSolution[-1].SmbGdn[0,0:240].reshape(1,-1), + md.results.TransientSolution[-1].SmbGsp[0,0:240].reshape(1,-1), + md.results.TransientSolution[-1].SmbA[0,0:240].reshape(1,-1), + md.results.TransientSolution[-1].SmbEC[0], + md.results.TransientSolution[-1].SmbMassBalance[0], + ] Index: /issm/trunk-jpl/test/NightlyRun/test260.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test260.py (revision 22267) +++ /issm/trunk-jpl/test/NightlyRun/test260.py (revision 22267) @@ -0,0 +1,31 @@ +#Test Name: SquareShelfStressSSA2dEnhanced +import numpy as np +from model import * +from socket import gethostname +from triangle import * +from setmask import * +from parameterize import * +from setflowequation import * +from solve import * +from matenhancedice import * + +md = triangle(model(),'../Exp/Square.exp',150000.) +md = setmask(md,'all','') +md.materials = matenhancedice() +md.materials.rheology_B = 3.15e8 * np.ones(md.mesh.numberofvertices,) +md.materials.rheology_n = 3 * np.ones(md.mesh.numberofelements,) +md.materials.rheology_E = np.ones(md.mesh.numberofvertices,) +md = parameterize(md,'../Par/SquareShelf.py') +md = setflowequation(md,'SSA','all') +md.cluster = generic('name',gethostname(),'np',3) +md = solve(md,'Stressbalance') + +#Fields and tolerances to track changes +field_names = ['Vx','Vy','Vel','Pressure'] +field_tolerances = [1e-13,1e-13,1e-13,1e-13] +field_values = [ + md.results.StressbalanceSolution.Vx, + md.results.StressbalanceSolution.Vy, + md.results.StressbalanceSolution.Vel, + md.results.StressbalanceSolution.Pressure, + ] Index: /issm/trunk-jpl/test/NightlyRun/test261.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test261.py (revision 22267) +++ /issm/trunk-jpl/test/NightlyRun/test261.py (revision 22267) @@ -0,0 +1,38 @@ +#Test Name: SquareShelfConstrainedTranEnhanced +import numpy as np +from model import * +from socket import gethostname +from triangle import * +from setmask import * +from parameterize import * +from setflowequation import * +from solve import * +from matenhancedice import * + +md = triangle(model(),'../Exp/Square.exp',180000.) +md = setmask(md,'all','') +md = parameterize(md,'../Par/SquareShelfConstrained.py') +md = md.extrude(3,1.) +md = setflowequation(md,'SSA','all') +md.cluster = generic('name',gethostname(),'np',3) +md.materials = matenhancedice() +md.materials.rheology_B = 3.15e8 * np.ones(md.mesh.numberofvertices,) +md.materials.rheology_n = 3 * np.ones(md.mesh.numberofelements,) +md.materials.rheology_E = np.ones(md.mesh.numberofvertices,) +md.transient.isstressbalance = 1 +md.transient.ismasstransport = 0 +md.transient.issmb = 1 +md.transient.isthermal = 1 +md.transient.isgroundingline = 0 +md = solve(md,'Transient') + +#Fields and tolerances to track changes +field_names = ['Vx','Vy','Vel','Temperature','BasalforcingsGroundediceMeltingRate'] +field_tolerances = [1e-13,1e-13,1e-13,1e-13,1e-13] +field_values = [ + md.results.TransientSolution[0].Vx, + md.results.TransientSolution[0].Vy, + md.results.TransientSolution[0].Vel, + md.results.TransientSolution[0].Temperature, + md.results.TransientSolution[0].BasalforcingsGroundediceMeltingRate, + ] Index: /issm/trunk-jpl/test/NightlyRun/test293.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test293.py (revision 22267) +++ /issm/trunk-jpl/test/NightlyRun/test293.py (revision 22267) @@ -0,0 +1,60 @@ +#Test Name: SquareShelfTranSSA2dMismipFloatingMeltParam +import numpy as np +from model import * +from socket import gethostname +from triangle import * +from setmask import * +from parameterize import * +from setflowequation import * +from solve import * +from mismipbasalforcings import * + +md = triangle(model(),'../Exp/Square.exp',150000.) +md = setmask(md,'all','') +md = parameterize(md,'../Par/SquareShelf.py') +md = setflowequation(md,'SSA','all') +md.cluster = generic('name',gethostname(),'np',3) +md.constants.yts = 365.2422 * 24. * 3600. +md.basalforcings = mismipbasalforcings() +md.basalforcings = md.basalforcings.initialize(md) +md.transient.isgroundingline = 1 +md.geometry.bed = min(md.geometry.base) * np.ones(md.mesh.numberofvertices,) +md.transient.requested_outputs = ['default','BasalforcingsFloatingiceMeltingRate'] +print md.constants.yts + +md = solve(md,'Transient') + +#Fields and tolerances to track changes +field_names =['Vx1','Vy1','Vel1','Pressure1','Bed1','Surface1','Thickness1','BasalforcingsFloatingiceMeltingRate1', + 'Vx2','Vy2','Vel2','Pressure2','Bed2','Surface2','Thickness2','BasalforcingsFloatingiceMeltingRate2', + 'Vx3','Vy3','Vel3','Pressure3','Bed3','Surface3','Thickness3','BasalforcingsFloatingiceMeltingRate3'] +field_tolerances = [ + 1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,2e-13, + 1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,2e-13, + 1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,2e-13] +field_values = [ + md.results.TransientSolution[0].Vx, + md.results.TransientSolution[0].Vy, + md.results.TransientSolution[0].Vel, + md.results.TransientSolution[0].Pressure, + md.results.TransientSolution[0].Base, + md.results.TransientSolution[0].Surface, + md.results.TransientSolution[0].Thickness, + md.results.TransientSolution[0].BasalforcingsFloatingiceMeltingRate, + md.results.TransientSolution[1].Vx, + md.results.TransientSolution[1].Vy, + md.results.TransientSolution[1].Vel, + md.results.TransientSolution[1].Pressure, + md.results.TransientSolution[1].Base, + md.results.TransientSolution[1].Surface, + md.results.TransientSolution[1].Thickness, + md.results.TransientSolution[1].BasalforcingsFloatingiceMeltingRate, + md.results.TransientSolution[2].Vx, + md.results.TransientSolution[2].Vy, + md.results.TransientSolution[2].Vel, + md.results.TransientSolution[2].Pressure, + md.results.TransientSolution[2].Base, + md.results.TransientSolution[2].Surface, + md.results.TransientSolution[2].Thickness, + md.results.TransientSolution[2].BasalforcingsFloatingiceMeltingRate, + ] Index: /issm/trunk-jpl/test/NightlyRun/test340.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test340.py (revision 22267) +++ /issm/trunk-jpl/test/NightlyRun/test340.py (revision 22267) @@ -0,0 +1,46 @@ +#Test Name: SquareSheetConstrainedSmbGradientsEla3d +import numpy as np +from model import * +from socket import gethostname +from triangle import * +from setmask import * +from parameterize import * +from setflowequation import * +from solve import * +from taoinversion import * + +md = triangle(model(),'../Exp/Square.exp',200000.) +md = setmask(md,'','') +md = parameterize(md,'../Par/SquareSheetConstrained.py') +md.extrude(3,1.) +md = setflowequation(md,'HO','all') + +#control parameters +md.inversion = taoinversion() +md.inversion.iscontrol = 1 +md.inversion.control_parameters = ['FrictionCoefficient'] +md.inversion.min_parameters = 1. * np.ones((md.mesh.numberofvertices,)) +md.inversion.max_parameters = 200. * np.ones((md.mesh.numberofvertices,)) +md.inversion.maxsteps = 2 +md.inversion.maxiter = 6 +md.inversion.cost_functions = [102,501] +md.inversion.cost_functions_coefficients = np.ones((md.mesh.numberofvertices,2)) +md.inversion.cost_functions_coefficients[:,1] = 2. * 1e-7 +md.inversion.vx_obs = md.initialization.vx +md.inversion.vy_obs = md.initialization.vy + +md.cluster = generic('name',gethostname(),'np',3) +md = solve(md,'Stressbalance') + +#Fields and tolerances to track changes +field_names = ['Gradient','Misfits','FrictionCoefficient','Pressure','Vel','Vx','Vy'] +field_tolerances = [3e-08,1e-07,5e-10,1e-10,1e-09,1e-09,1e-09] +field_values = [ + md.results.StressbalanceSolution.Gradient1, + md.results.StressbalanceSolution.J, + md.results.StressbalanceSolution.FrictionCoefficient, + md.results.StressbalanceSolution.Pressure, + md.results.StressbalanceSolution.Vel, + md.results.StressbalanceSolution.Vx, + md.results.StressbalanceSolution.Vy +] Index: /issm/trunk-jpl/test/NightlyRun/test342.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test342.py (revision 22267) +++ /issm/trunk-jpl/test/NightlyRun/test342.py (revision 22267) @@ -0,0 +1,35 @@ +#Test Name: SquareSheetTherSteaPlume +import numpy as np +from model import * +from socket import gethostname +from triangle import * +from setmask import * +from parameterize import * +from setflowequation import * +from solve import * +from plumebasalforcings import * + +md = triangle(model(),'../Exp/Square.exp',180000.) +md = setmask(md,'','') +md = parameterize(md,'../Par/SquareSheetConstrained.py') +md.basalforcings = plumebasalforcings() +md.basalforcings = md.basalforcings.setdefaultparameters() +md.basalforcings.floatingice_melting_rate = np.zeros((md.mesh.numberofvertices,)) +md.basalforcings.groundedice_melting_rate = np.zeros((md.mesh.numberofvertices,)) +md.basalforcings.plumex = 500000 +md.basalforcings.plumey = 500000 +md.extrude(3,1.) +md = setflowequation(md,'SSA','all') +md.timestepping.time_step = 0. +md.thermal.requested_outputs = ['default','BasalforcingsGeothermalflux'] +md.cluster = generic('name',gethostname(),'np',3) +md = solve(md,'Thermal') + +#Fields and tolerances to track changes +field_names = ['Temperature','BasalforcingsGroundediceMeltingRate','BasalforcingsGeothermalflux'] +field_tolerances = [1e-13,1e-8,1e-13] +field_values = [ + md.results.ThermalSolution.Temperature, + md.results.ThermalSolution.BasalforcingsGroundediceMeltingRate, + md.results.ThermalSolution.BasalforcingsGeothermalflux, + ] Index: /issm/trunk-jpl/test/NightlyRun/test343.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test343.py (revision 22267) +++ /issm/trunk-jpl/test/NightlyRun/test343.py (revision 22267) @@ -0,0 +1,66 @@ +#Test Name: SquareSheetConstrainedSmbGradientsEla2d +import numpy as np +from model import * +from socket import gethostname +from triangle import * +from setmask import * +from parameterize import * +from setflowequation import * +from solve import * +from SMBgradientsela import * + +md = triangle(model(),'../Exp/Square.exp',150000.) +md = setmask(md,'','') +md = parameterize(md,'../Par/SquareSheetConstrained.py') +md = setflowequation(md,'SSA','all') +md.smb = SMBgradientsela() +md.smb.ela = 1500. * np.ones((md.mesh.numberofvertices+1,)) +md.smb.b_pos = 0.002 * np.ones((md.mesh.numberofvertices+1,)) +md.smb.b_neg = 0.005 * np.ones((md.mesh.numberofvertices+1,)) +md.smb.b_max = 4. * (md.materials.rho_freshwater / md.materials.rho_ice) * np.ones((md.mesh.numberofvertices+1,)) +md.smb.b_min = -4. * (md.materials.rho_freshwater / md.materials.rho_ice) * np.ones((md.mesh.numberofvertices+1,)) + +#Change geometry +md.geometry.thickness = md.geometry.surface * 30. +md.geometry.surface = md.geometry.base + md.geometry.thickness + +#Transient options +md.transient.requested_outputs = ['default','TotalSmb'] +md.cluster = generic('name',gethostname(),'np',3) +md = solve(md,'Transient') + +#Fields and tolerances to track changes +field_names = [ + 'Vx1','Vy1','Vel1','Bed1','Surface1','Thickness1','SMB1','TotalSmb1', + 'Vx2','Vy2','Vel2','Bed2','Surface2','Thickness2','SMB2','TotalSmb2', + 'Vx3','Vy3','Vel3','Bed3','Surface3','Thickness3','SMB3','TotalSmb3'] +field_tolerances = [ + 1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13, + 1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13, + 1e-12,1e-12,1e-12,1e-13,1e-13,1e-13,1.5e-13,1e-13] +field_values = [ + md.results.TransientSolution[0].Vx, + md.results.TransientSolution[0].Vy, + md.results.TransientSolution[0].Vel, + md.results.TransientSolution[0].Base, + md.results.TransientSolution[0].Surface, + md.results.TransientSolution[0].Thickness, + md.results.TransientSolution[0].SmbMassBalance, + md.results.TransientSolution[0].TotalSmb, + md.results.TransientSolution[1].Vx, + md.results.TransientSolution[1].Vy, + md.results.TransientSolution[1].Vel, + md.results.TransientSolution[1].Base, + md.results.TransientSolution[1].Surface, + md.results.TransientSolution[1].Thickness, + md.results.TransientSolution[1].TotalSmb, + md.results.TransientSolution[1].SmbMassBalance, + md.results.TransientSolution[2].Vx, + md.results.TransientSolution[2].Vy, + md.results.TransientSolution[2].Vel, + md.results.TransientSolution[2].Base, + md.results.TransientSolution[2].Surface, + md.results.TransientSolution[2].Thickness, + md.results.TransientSolution[2].SmbMassBalance, + md.results.TransientSolution[2].TotalSmb + ] Index: /issm/trunk-jpl/test/NightlyRun/test344.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test344.py (revision 22267) +++ /issm/trunk-jpl/test/NightlyRun/test344.py (revision 22267) @@ -0,0 +1,73 @@ +#Test Name: SquareSheetConstrainedSmbGradientsEla3d +import numpy as np +from model import * +from socket import gethostname +from triangle import * +from setmask import * +from parameterize import * +from setflowequation import * +from solve import * +from SMBgradientsela import * + +md = triangle(model(),'../Exp/Square.exp',150000.) +md = setmask(md,'','') +md = parameterize(md,'../Par/SquareSheetConstrained.py') + +#Change geometry +md.geometry.thickness = md.geometry.surface * 30. +md.geometry.surface = md.geometry.base + md.geometry.thickness + +md = md.extrude(3,1.) +md = setflowequation(md,'HO','all') +md.smb = SMBgradientsela() +md.smb.ela = 1500. * np.ones((md.mesh.numberofvertices+1,)) +md.smb.b_pos = 0.002 * np.ones((md.mesh.numberofvertices+1,)) +md.smb.b_neg = 0.005 * np.ones((md.mesh.numberofvertices+1,)) +md.smb.b_max = 4. * (md.materials.rho_freshwater / md.materials.rho_ice) * np.ones((md.mesh.numberofvertices+1,)) +md.smb.b_min = -4. * (md.materials.rho_freshwater / md.materials.rho_ice) * np.ones((md.mesh.numberofvertices+1,)) + + +#Transient options +md.transient.requested_outputs = ['default','TotalSmb'] +md.cluster = generic('name',gethostname(),'np',3) +md = solve(md,'Transient') + +#Fields and tolerances to track changes +field_names = ['Vx1','Vy1','Vz1','Vel1','Bed1','Surface1','Thickness1','Temperature1','SMB1','TotalSmb1', + 'Vx2','Vy2','Vz2','Vel2','Bed2','Surface2','Thickness2','Temperature2','SMB2','TotalSmb2', + 'Vx3','Vy3','Vz3','Vel3','Bed3','Surface3','Thickness3','Temperature3','SMB3','TotalSmb3'] +field_tolerances = [1e-09,1e-09,1e-09,1e-09,1e-10,1e-10,1e-10,1e-10,1e-10,1e-10, + 1e-09,1e-09,1e-10,1e-09,1e-10,1e-10,1e-10,1e-10,1e-10,1e-10, + 1e-09,5e-09,1e-09,1e-09,1e-10,1e-10,1e-10,1e-10,1e-10,1e-10] +field_values = [ + md.results.TransientSolution[0].Vx, + md.results.TransientSolution[0].Vy, + md.results.TransientSolution[0].Vz, + md.results.TransientSolution[0].Vel, + md.results.TransientSolution[0].Base, + md.results.TransientSolution[0].Surface, + md.results.TransientSolution[0].Thickness, + md.results.TransientSolution[0].Temperature, + md.results.TransientSolution[0].SmbMassBalance, + md.results.TransientSolution[0].TotalSmb, + md.results.TransientSolution[1].Vx, + md.results.TransientSolution[1].Vy, + md.results.TransientSolution[1].Vz, + md.results.TransientSolution[1].Vel, + md.results.TransientSolution[1].Base, + md.results.TransientSolution[1].Surface, + md.results.TransientSolution[1].Thickness, + md.results.TransientSolution[1].Temperature, + md.results.TransientSolution[1].SmbMassBalance, + md.results.TransientSolution[1].TotalSmb, + md.results.TransientSolution[2].Vx, + md.results.TransientSolution[2].Vy, + md.results.TransientSolution[2].Vz, + md.results.TransientSolution[2].Vel, + md.results.TransientSolution[2].Base, + md.results.TransientSolution[2].Surface, + md.results.TransientSolution[2].Thickness, + md.results.TransientSolution[2].Temperature, + md.results.TransientSolution[2].SmbMassBalance, + md.results.TransientSolution[2].TotalSmb, + ] Index: /issm/trunk-jpl/test/NightlyRun/test350.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test350.py (revision 22267) +++ /issm/trunk-jpl/test/NightlyRun/test350.py (revision 22267) @@ -0,0 +1,94 @@ +#Test Name: SquareSheetHydrologySommers +from operator import itemgetter +import numpy as np +from model import * +from socket import gethostname +from triangle import * +from setmask import * +from parameterize import * +from setflowequation import * +from solve import * +from frictionsommers import * +from hydrologysommers import * +from transient import * + +md = triangle(model(),'../Exp/Square.exp',50000.) +md.mesh.x = md.mesh.x / 1000 +md.mesh.y = md.mesh.y / 1000 +md = setmask(md,'','') +md = parameterize(md,'../Par/SquareSheetConstrained.py') +md.transient = transient().deactivateall() +md.transient.ishydrology = 1 +md = setflowequation(md,'SSA','all') +md.cluster = generic('name',gethostname(),'np',2) + +#Use hydroogy coupled friciton law +md.friction = frictionsommers(md) + +#Change hydrology class to Sommers' model +md.hydrology = hydrologysommers() + +#Change geometry +md.geometry.base = -.02 * md.mesh.x + 20. +md.geometry.thickness = 300. * np.ones((md.mesh.numberofvertices,)) +md.geometry.bed = md.geometry.base +md.geometry.surface = md.geometry.bed + md.geometry.thickness + +#define the initial water head as being such that the water pressure is 50% of the ice overburden pressure +md.hydrology.head = 0.5 * md.materials.rho_ice / md.materials.rho_freshwater * md.geometry.thickness + md.geometry.base +md.hydrology.gap_height = 0.01 * np.ones((md.mesh.numberofelements,)) +md.hydrology.bump_spacing = 2 * np.ones((md.mesh.numberofelements,)) +md.hydrology.bump_height = 0.05 * np.ones((md.mesh.numberofelements,)) +md.hydrology.englacial_input = 0.5 * np.ones((md.mesh.numberofvertices,)) +md.hydrology.reynolds= 1000. * np.ones((md.mesh.numberofelements,)) +md.hydrology.spchead = float('NaN') * np.ones((md.mesh.numberofvertices,)) +pos = np.intersect1d(np.array(np.where(md.mesh.vertexonboundary)), np.array(np.where(md.mesh.x == 1000))) +md.hydrology.spchead[pos] = md.geometry.base[pos] + +#Define velocity +md.initialization.vx = 1e-6 * md.constants.yts * np.ones((md.mesh.numberofvertices,)) +md.initialization.vy = np.zeros((md.mesh.numberofvertices,)) + +md.timestepping.time_step = 3. * 3600. / md.constants.yts +md.timestepping.final_time = .5 / 365. +md.materials.rheology_B = (5e-25)**(-1. / 3.) * np.ones((md.mesh.numberofvertices,)) + +#Add one moulin and Neumann BC, varying in time +a = np.sqrt((md.mesh.x - 500.)**2 + (md.mesh.y - 500.)**2) +pos = min(enumerate(a), key=itemgetter(1))[0] +time = np.arange(0,md.timestepping.final_time+1,md.timestepping.time_step) +md.hydrology.moulin_input = np.zeros((md.mesh.numberofvertices+1,np.size(time))) +md.hydrology.moulin_input[-1,:] = time +md.hydrology.moulin_input[pos,:] = 5. * (1. - np.sin(2. * np.pi / (1. / 365.) * time)) +md.hydrology.neumannflux = np.zeros((md.mesh.numberofelements+1,np.size(time))) +md.hydrology.neumannflux[-1,:] = time +segx = md.mesh.x[md.mesh.segments[:,0]-1] +segy = md.mesh.y[md.mesh.segments[:,0]-1] +posA = np.intersect1d(np.intersect1d(np.array(np.where(segx < 1.)),np.array(np.where(segy > 400.))), np.array(np.where(segy < 600.))) +pos = (md.mesh.segments[posA]-1)[:,2] +md.hydrology.neumannflux[pos,:] = np.tile(0.05*(1.-np.sin(2.*np.pi/(1./365.)*time)),(len(pos),1)) + +md = solve(md,'Transient') + +#Fields and tolerances to track changes +field_names = [ + 'HydrologyHead1','HydrologyGapHeight1', + 'HydrologyHead2','HydrologyGapHeight2', + 'HydrologyHead3','HydrologyGapHeight3', + 'HydrologyHead4','HydrologyGapHeight4'] +field_tolerances = [ + 1e-13, 1e-13, + 1e-13, 1e-13, + 1e-13, 1e-13, + 1e-13, 1e-12] +field_values = [ + md.results.TransientSolution[0].HydrologyHead, + md.results.TransientSolution[0].HydrologyGapHeight, + md.results.TransientSolution[1].HydrologyHead, + md.results.TransientSolution[1].HydrologyGapHeight, + md.results.TransientSolution[2].HydrologyHead, + md.results.TransientSolution[2].HydrologyGapHeight, + md.results.TransientSolution[3].HydrologyHead, + md.results.TransientSolution[3].HydrologyGapHeight + ] + Index: /issm/trunk-jpl/test/NightlyRun/test430.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test430.py (revision 22267) +++ /issm/trunk-jpl/test/NightlyRun/test430.py (revision 22267) @@ -0,0 +1,95 @@ +#Test Name: MISMIP3D +import numpy as np +from model import * +from socket import gethostname +from triangle import * +from setmask import * +from parameterize import * +from setflowequation import * +from solve import * + +md = triangle(model(),'../Exp/Square.exp',100000.) +md = setmask(md,'../Exp/SquareShelf.exp','') +md = parameterize(md,'../Par/SquareSheetShelf.py') +md.initialization.vx[:] = 1. +md.initialization.vy[:] = 1. +md.geometry.thickness[:] = 500. - md.mesh.x / 10000. +md.geometry.bed = -100. - md.mesh.x / 1000. +md.geometry.base = -md.geometry.thickness * md.materials.rho_ice / md.materials.rho_water +md.mask.groundedice_levelset = md.geometry.thickness + md.materials.rho_water / md.materials.rho_ice * md.geometry.bed +pos = np.where(md.mask.groundedice_levelset >= 0.) +md.geometry.base[pos] = md.geometry.bed[pos] +md.geometry.surface = md.geometry.base + md.geometry.thickness +md = setflowequation(md,'SSA','all') + +#Boundary conditions: +md.mask.ice_levelset = -np.ones((md.mesh.numberofvertices,)) +md.mask.ice_levelset[np.where(md.mesh.x == max(md.mesh.x))] = 0. +md.stressbalance.spcvx[:] = float('NaN') +md.stressbalance.spcvy[:] = float('NaN') +md.stressbalance.spcvz[:] = float('NaN') +posA = np.intersect1d(np.array(np.where(md.mesh.y < 1000000.1)),np.array(np.where(md.mesh.y > 999999.9))) +posB = np.intersect1d(np.array(np.where(md.mesh.y < 0.1)),np.array(np.where(md.mesh.y > -0.1))) +pos = np.unique(np.concatenate((posA,posB))) +md.stressbalance.spcvy[pos] = 0. +pos2 = np.intersect1d(np.array(np.where(md.mesh.x < 0.1)), np.array(np.where(md.mesh.x > -0.1))) +md.stressbalance.spcvx[pos2] = 0. +md.stressbalance.spcvy[pos2] = 0. + +md.materials.rheology_B = 1. / ((10**-25)**(1. / 3.)) * np.ones((md.mesh.numberofvertices,)) +md.materials.rheology_law = 'None' +md.friction.coefficient[:] = np.sqrt(10**7) * np.ones((md.mesh.numberofvertices,)) +md.friction.p = 3. * np.ones((md.mesh.numberofelements,)) +md.smb.mass_balance[:] = 1. +md.basalforcings.groundedice_melting_rate[:] = 0. +md.basalforcings.floatingice_melting_rate[:] = 30. +md.transient.isthermal = 0 +md.transient.isstressbalance = 1 +md.transient.isgroundingline = 1 +md.transient.ismasstransport = 1 +md.transient.issmb = 1 +md.transient.requested_outputs = ['default','BasalforcingsFloatingiceMeltingRate'] +md.groundingline.migration = 'SubelementMigration' +md.timestepping.final_time = 30 +md.timestepping.time_step = 10 + +md.cluster = generic('name',gethostname(),'np',3) +md = solve(md,'Transient') +#print md.results.TransientSolution[0].BasalforcingsFloatingiceMeltingRate +#print md.results.TransientSolution[1].BasalforcingsFloatingiceMeltingRate +#print md.results.TransientSolution[2].BasalforcingsFloatingiceMeltingRate + +#Fields and tolerances to track changes +field_names = [ + 'Bed1','Surface1','Thickness1','Floatingice1','Vx1','Vy1','Pressure1','FloatingiceMeltingrate1', + 'Bed2','Surface2','Thickness2','Floatingice2','Vx2','Vy2','Pressure2','FloatingiceMeltingrate2', + 'Bed3','Surface3','Thickness3','Floatingice3','Vx3','Vy3','Pressure3','FloatingiceMeltingrate3'] +field_tolerances = [2e-11,5e-12,2e-11,1e-11,5e-10,1e-08,1e-13,1e-13, + 3e-11,3e-11,9e-10,7e-11,1e-09,5e-08,1e-10,1e-13, + 1e-10,3e-11,1e-10,7e-11,1e-09,5e-08,1e-10,1e-13] +field_values = [ + md.results.TransientSolution[0].Base, + md.results.TransientSolution[0].Surface, + md.results.TransientSolution[0].Thickness, + md.results.TransientSolution[0].MaskGroundediceLevelset, + md.results.TransientSolution[0].Vx, + md.results.TransientSolution[0].Vy, + md.results.TransientSolution[0].Pressure, + md.results.TransientSolution[0].BasalforcingsFloatingiceMeltingRate, + md.results.TransientSolution[1].Base, + md.results.TransientSolution[1].Surface, + md.results.TransientSolution[1].Thickness, + md.results.TransientSolution[1].MaskGroundediceLevelset, + md.results.TransientSolution[1].Vx, + md.results.TransientSolution[1].Vy, + md.results.TransientSolution[1].Pressure, + md.results.TransientSolution[1].BasalforcingsFloatingiceMeltingRate, + md.results.TransientSolution[2].Base, + md.results.TransientSolution[2].Surface, + md.results.TransientSolution[2].Thickness, + md.results.TransientSolution[2].MaskGroundediceLevelset, + md.results.TransientSolution[2].Vx, + md.results.TransientSolution[2].Vy, + md.results.TransientSolution[2].Pressure, + md.results.TransientSolution[2].BasalforcingsFloatingiceMeltingRate, + ] Index: /issm/trunk-jpl/test/NightlyRun/test435.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test435.py (revision 22267) +++ /issm/trunk-jpl/test/NightlyRun/test435.py (revision 22267) @@ -0,0 +1,97 @@ +#Test Name: MISMIP3DHO +import numpy as np +from model import * +from socket import gethostname +from triangle import * +from setmask import * +from parameterize import * +from setflowequation import * +from solve import * + +md = triangle(model(),'../Exp/Square.exp',100000.) +md = setmask(md,'../Exp/SquareShelf.exp','') +md = parameterize(md,'../Par/SquareSheetShelf.py') +md.initialization.vx[:] = 1. +md.initialization.vy[:] = 1. +md.geometry.thickness[:] = 500. - md.mesh.x / 10000. +md.geometry.bed = -100. - md.mesh.x / 1000. +md.geometry.base = -md.geometry.thickness * md.materials.rho_ice / md.materials.rho_water +md.mask.groundedice_levelset = md.geometry.thickness + md.materials.rho_water / md.materials.rho_ice * md.geometry.bed +pos = np.where(md.mask.groundedice_levelset >= 0) +md.geometry.base[pos] = md.geometry.bed[pos] +md.geometry.surface = md.geometry.base + md.geometry.thickness +md = md.extrude(4,1.) +md = setflowequation(md,'HO','all') + +#Boundary conditions: +md.mask.ice_levelset = -np.ones((md.mesh.numberofvertices,)) +md.mask.ice_levelset[np.where(md.mesh.x == max(md.mesh.x))] = 0. +md.stressbalance.spcvx[:] = float('Nan') +md.stressbalance.spcvy[:] = float('Nan') +md.stressbalance.spcvz[:] = float('Nan') +posA = np.intersect1d(np.array(np.where(md.mesh.y < 1000000.1)),np.array(np.where(md.mesh.y > 999999.9))) +posB = np.intersect1d(np.array(np.where(md.mesh.y < 0.1)),np.array(np.where(md.mesh.y > -0.1))) +pos = np.unique(np.concatenate((posA,posB))) +md.stressbalance.spcvy[pos] = 0. +pos2 = np.intersect1d(np.array(np.where(md.mesh.x < 0.1)), np.array(np.where(md.mesh.x > -0.1))) +md.stressbalance.spcvx[pos2] = 0. +md.stressbalance.spcvy[pos2] = 0. + +md.materials.rheology_B = 1. / ((10**-25)**(1./3.)) * np.ones((md.mesh.numberofvertices,)) +md.materials.rheology_law = 'None' +md.friction.coefficient[:] = np.sqrt(1e7) * np.ones((md.mesh.numberofvertices,)) +md.friction.p = 3. * np.ones((md.mesh.numberofelements,)) +md.smb.mass_balance[:] = 1. +md.basalforcings.groundedice_melting_rate[:] = 0. +md.basalforcings.floatingice_melting_rate[:] = 30. +md.transient.isthermal = 0 +md.transient.isstressbalance = 1 +md.transient.isgroundingline = 1 +md.transient.ismasstransport = 1 +md.transient.issmb = 1 +md.transient.requested_outputs = ['default','BasalforcingsFloatingiceMeltingRate'] +md.groundingline.migration = 'SubelementMigration' +md.timestepping.final_time = 30 +md.timestepping.time_step = 10 + +md.cluster = generic('name',gethostname(),'np',3) +md = solve(md,'Transient') + +#Fields and tolerances to track changes +field_names = [ + 'Bed1','Surface1','Thickness1','Floatingice1','Vx1','Vy1','Vz1','Pressure1','FloatingiceMeltingrate1', + 'Bed2','Surface2','Thickness2','Floatingice2','Vx2','Vy2','Vz2','Pressure2','FloatingiceMeltingrate2', + 'Bed3','Surface3','Thickness3','Floatingice3','Vx3','Vy3','Vz3','Pressure3','FloatingiceMeltingrate3'] +field_tolerances = [ + 2e-11,5e-12,2e-11,1e-11,5e-10,3e-08,6e-10,1e-13,1e-13, + 3e-11,3e-11,9e-10,7e-11,6e-09,1e-07,1e-09,1e-10,1e-13, + 1e-9,2e-08,7e-09,2e-7 ,1e-03,8e-04,2e-09,1e-10,1e-13] +field_values = [ + md.results.TransientSolution[0].Base, + md.results.TransientSolution[0].Surface, + md.results.TransientSolution[0].Thickness, + md.results.TransientSolution[0].MaskGroundediceLevelset, + md.results.TransientSolution[0].Vx, + md.results.TransientSolution[0].Vy, + md.results.TransientSolution[0].Vz, + md.results.TransientSolution[0].Pressure, + md.results.TransientSolution[0].BasalforcingsFloatingiceMeltingRate, + md.results.TransientSolution[1].Base, + md.results.TransientSolution[1].Surface, + md.results.TransientSolution[1].Thickness, + md.results.TransientSolution[1].MaskGroundediceLevelset, + md.results.TransientSolution[1].Vx, + md.results.TransientSolution[1].Vy, + md.results.TransientSolution[1].Vz, + md.results.TransientSolution[1].Pressure, + md.results.TransientSolution[1].BasalforcingsFloatingiceMeltingRate, + md.results.TransientSolution[2].Base, + md.results.TransientSolution[2].Surface, + md.results.TransientSolution[2].Thickness, + md.results.TransientSolution[2].MaskGroundediceLevelset, + md.results.TransientSolution[2].Vx, + md.results.TransientSolution[2].Vy, + md.results.TransientSolution[2].Vz, + md.results.TransientSolution[2].Pressure, + md.results.TransientSolution[2].BasalforcingsFloatingiceMeltingRate, + ] Index: /issm/trunk-jpl/test/NightlyRun/test436.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test436.py (revision 22267) +++ /issm/trunk-jpl/test/NightlyRun/test436.py (revision 22267) @@ -0,0 +1,45 @@ +#Test Name: SquareSheetShelfSteaEnthalpyRheologiesHO +import numpy as np +from model import * +from socket import gethostname +from triangle import * +from setmask import * +from parameterize import * +from setflowequation import * +from solve import * + +md = triangle(model(),'../Exp/Square.exp',150000.) +md = setmask(md,'../Exp/SquareShelf.exp','') +md = parameterize(md,'../Par/SquareSheetShelf.py') +md = md.extrude(3,2.) +md = setflowequation(md,'HO','all') +md.cluster = generic('name',gethostname(),'np',3) +md.timestepping.time_step = 0. +md.thermal.isenthalpy = 1 +md.initialization.waterfraction = np.zeros((md.mesh.numberofvertices,)) +md.initialization.watercolumn = np.zeros((md.mesh.numberofvertices,)) + +#Go solve +field_names = [] +field_tolerances = [] +field_values = [] +for i in ['LliboutryDuval', 'CuffeyTemperate']: + print ' ' + print '====== Testing rheology law: ' + i + ' =====' + + md.materials.rheology_law = i + md = solve(md,'Steadystate') + field_names += ['Vx'+i,'Vy'+i,'Vz'+i,'Vel'+i,'Pressure'+i, + 'Temperature'+i,'Waterfraction'+i,'Enthalpy'+i] + field_tolerances += [2e-09,1e-09,1e-09,1e-09,1e-13,1e-10,6e-10,5e-10] + field_values += [ + md.results.SteadystateSolution.Vx, + md.results.SteadystateSolution.Vy, + md.results.SteadystateSolution.Vz, + md.results.SteadystateSolution.Vel, + md.results.SteadystateSolution.Pressure, + md.results.SteadystateSolution.Temperature, + md.results.SteadystateSolution.Waterfraction, + md.results.SteadystateSolution.Enthalpy, + ] + Index: /issm/trunk-jpl/test/NightlyRun/test437.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test437.py (revision 22267) +++ /issm/trunk-jpl/test/NightlyRun/test437.py (revision 22267) @@ -0,0 +1,96 @@ +#Test Name: ThermalEnthBasalcondsTrans +import numpy as np +from model import * +from socket import gethostname +from triangle import * +from setmask import * +from parameterize import * +from setflowequation import * +from solve import * + +md = triangle(model(),'../Exp/Square.exp',300000.) +md = setmask(md,'','') +md = parameterize(md,'../Par/SquareThermal.py') + +h = 100. +md.geometry.thickness = h * np.ones((md.mesh.numberofvertices,)) +md.geometry.base = -h * np.ones((md.mesh.numberofvertices,)) +md.geometry.surface = md.geometry.base + md.geometry.thickness + +md.extrude(41,2.) +md = setflowequation(md,'HO','all') +md.thermal.isenthalpy = True +md.thermal.isdynamicbasalspc = True + +#Basal forcing +Ts = 273.15-3. +Tb = 273.15-1. +Tsw = Tb +qgeo = md.materials.thermalconductivity / max(md.geometry.thickness) * (Tb - Ts) #qgeo=kappa*(Tb-Ts)/H +md.basalforcings.geothermalflux[np.where(md.mesh.vertexonbase)] = qgeo +md.initialization.temperature = qgeo / md.materials.thermalconductivity * (md.geometry.surface - md.mesh.z) + Ts + +#Surface forcing +pos = np.where(md.mesh.vertexonsurface) +SPC_cold = float('NaN') * np.ones((md.mesh.numberofvertices,)) +SPC_warm = float('NaN') * np.ones((md.mesh.numberofvertices,)) +SPC_cold[pos] = Ts +SPC_warm[pos] = Tsw +md.thermal.spctemperature = SPC_cold +md.timestepping.time_step = 5. +t0 = 0. +t1 = 10. +t2 = 100. +md.timestepping.final_time = 400. +md.thermal.spctemperature = np.array([SPC_cold,SPC_cold,SPC_warm,SPC_warm,SPC_cold]).T +md.thermal.spctemperature = np.vstack((md.thermal.spctemperature,[t0, t1-1, t1, t2-1, t2])) +#print np.shape(md.thermal.spctemperature) +#print md.thermal.spctemperature + +#Additional settings +md.transient.ismasstransport = False +md.transient.isstressbalance = False +md.transient.issmb = True +md.transient.isthermal = True +md.thermal.stabilization = False + +#Go solve +md.verbose = verbose(0) +md.cluster = generic('name',gethostname(),'np',3) +md = solve(md,'Transient') + +#Fields and tolerances to track changes +field_names = ['Enthalpy1','Temperature1','Waterfraction1','BasalMeltingRate1','Watercolumn1', + 'Enthalpy2','Temperature2','Waterfraction2','BasalMeltingRate2','Watercolumn2', + 'Enthalpy3','Temperature3','Waterfraction3','BasalMeltingRate3','Watercolumn3', + 'Enthalpy4','Temperature4','Waterfraction4','BasalMeltingRate4','Watercolumn4'] +field_tolerances = [1.e-10,1.e-10,1.e-10,1.e-9,1.e-10, + 1.e-10,1.e-10,1.e-10,2.e-9,2.e-10, + 1.e-10,1.e-10,1.e-10,2.e-9,1.e-10, + 1.e-10,1.e-10,1.e-10,2.e-9,1.e-10] +i1 = 0 +i2 = int(np.ceil(t2 / md.timestepping.time_step) + 2)-1 +i3 = int(np.ceil(md.timestepping.final_time / (2. * md.timestepping.time_step)))-1 +i4 = np.shape(md.results.TransientSolution)[0]-1 +field_values = [ + md.results.TransientSolution[i1].Enthalpy, + md.results.TransientSolution[i1].Temperature, + md.results.TransientSolution[i1].Waterfraction, + md.results.TransientSolution[i1].BasalforcingsGroundediceMeltingRate, + md.results.TransientSolution[i1].Watercolumn, + md.results.TransientSolution[i2].Enthalpy, + md.results.TransientSolution[i2].Temperature, + md.results.TransientSolution[i2].Waterfraction, + md.results.TransientSolution[i2].BasalforcingsGroundediceMeltingRate, + md.results.TransientSolution[i2].Watercolumn, + md.results.TransientSolution[i3].Enthalpy, + md.results.TransientSolution[i3].Temperature, + md.results.TransientSolution[i3].Waterfraction, + md.results.TransientSolution[i3].BasalforcingsGroundediceMeltingRate, + md.results.TransientSolution[i3].Watercolumn, + md.results.TransientSolution[i4].Enthalpy, + md.results.TransientSolution[i4].Temperature, + md.results.TransientSolution[i4].Waterfraction, + md.results.TransientSolution[i4].BasalforcingsGroundediceMeltingRate, + md.results.TransientSolution[i4].Watercolumn, + ] Index: /issm/trunk-jpl/test/NightlyRun/test438.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test438.py (revision 22267) +++ /issm/trunk-jpl/test/NightlyRun/test438.py (revision 22267) @@ -0,0 +1,53 @@ +#Test Name: TransientFrictionWaterlayer2D +import numpy as np +from model import * +from socket import gethostname +from triangle import * +from setmask import * +from parameterize import * +from setflowequation import * +from solve import * +from frictionwaterlayer import * + +md = triangle(model(),'../Exp/Square.exp',150000.) +md = setmask(md,'../Exp/SquareShelf.exp','') +md = parameterize(md,'../Par/SquareSheetShelf.py') +md = setflowequation(md,'SSA','all') +md.friction = frictionwaterlayer(md) +md.friction.water_layer = np.zeros((md.mesh.numberofvertices+1,2)) +md.friction.water_layer[:,1] = 1. +md.friction.water_layer[md.mesh.numberofvertices,:] = [1.,2.] +md.friction.f = 0.8 +md.cluster = generic('name',gethostname(),'np',3) +md = solve(md,'Transient') + +#Fields and tolerances to track changes +field_names =['Vx1','Vy1','Vel1','Pressure1','Bed1','Surface1','Thickness1', + 'Vx2','Vy2','Vel2','Pressure2','Bed2','Surface2','Thickness2', + 'Vx3','Vy3','Vel3','Pressure3','Bed3','Surface3','Thickness3'] +field_tolerances=[1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13, + 1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13, + 1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13] +field_values = [ + md.results.TransientSolution[0].Vx, + md.results.TransientSolution[0].Vy, + md.results.TransientSolution[0].Vel, + md.results.TransientSolution[0].Pressure, + md.results.TransientSolution[0].Base, + md.results.TransientSolution[0].Surface, + md.results.TransientSolution[0].Thickness, + md.results.TransientSolution[1].Vx, + md.results.TransientSolution[1].Vy, + md.results.TransientSolution[1].Vel, + md.results.TransientSolution[1].Pressure, + md.results.TransientSolution[1].Base, + md.results.TransientSolution[1].Surface, + md.results.TransientSolution[1].Thickness, + md.results.TransientSolution[2].Vx, + md.results.TransientSolution[2].Vy, + md.results.TransientSolution[2].Vel, + md.results.TransientSolution[2].Pressure, + md.results.TransientSolution[2].Base, + md.results.TransientSolution[2].Surface, + md.results.TransientSolution[2].Thickness + ] Index: /issm/trunk-jpl/test/NightlyRun/test439.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test439.py (revision 22267) +++ /issm/trunk-jpl/test/NightlyRun/test439.py (revision 22267) @@ -0,0 +1,54 @@ +#Test Name: TransientFrictionWaterlayer3D +import numpy as np +from model import * +from socket import gethostname +from triangle import * +from setmask import * +from parameterize import * +from setflowequation import * +from solve import * +from frictionwaterlayer import * + +md = triangle(model(),'../Exp/Square.exp',150000.) +md = setmask(md,'../Exp/SquareShelf.exp','') +md = parameterize(md,'../Par/SquareSheetShelf.py') +md = md.extrude(4,1.) +md = setflowequation(md,'HO','all') +md.friction = frictionwaterlayer(md) +md.friction.water_layer = np.zeros((md.mesh.numberofvertices+1,2)) +md.friction.water_layer[:,1] = 1. +md.friction.water_layer[md.mesh.numberofvertices,:] = [1.,2.] +md.friction.f = 0.8 +md.cluster = generic('name',gethostname(),'np',3) +md = solve(md,'Transient') + +#Fields and tolerances to track changes +field_names =['Vx1','Vy1','Vel1','Pressure1','Bed1','Surface1','Thickness1', + 'Vx2','Vy2','Vel2','Pressure2','Bed2','Surface2','Thickness2', + 'Vx3','Vy3','Vel3','Pressure3','Bed3','Surface3','Thickness3'] +field_tolerances = [1e-09,1e-09,1e-09,1e-09,1e-09,1e-09,1e-09, + 1e-09,1e-09,1e-09,1e-09,1e-09,1e-09,1e-09, + 1e-09,1e-09,1e-09,1e-09,1e-09,1e-09,1e-09] +field_values = [ + md.results.TransientSolution[0].Vx, + md.results.TransientSolution[0].Vy, + md.results.TransientSolution[0].Vel, + md.results.TransientSolution[0].Pressure, + md.results.TransientSolution[0].Base, + md.results.TransientSolution[0].Surface, + md.results.TransientSolution[0].Thickness, + md.results.TransientSolution[1].Vx, + md.results.TransientSolution[1].Vy, + md.results.TransientSolution[1].Vel, + md.results.TransientSolution[1].Pressure, + md.results.TransientSolution[1].Base, + md.results.TransientSolution[1].Surface, + md.results.TransientSolution[1].Thickness, + md.results.TransientSolution[2].Vx, + md.results.TransientSolution[2].Vy, + md.results.TransientSolution[2].Vel, + md.results.TransientSolution[2].Pressure, + md.results.TransientSolution[2].Base, + md.results.TransientSolution[2].Surface, + md.results.TransientSolution[2].Thickness + ] Index: /issm/trunk-jpl/test/NightlyRun/test441.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test441.py (revision 22267) +++ /issm/trunk-jpl/test/NightlyRun/test441.py (revision 22267) @@ -0,0 +1,93 @@ +#Test Name: MISMIP3D +import numpy as np +from model import * +from socket import gethostname +from triangle import * +from setmask import * +from parameterize import * +from setflowequation import * +from solve import * + +md = triangle(model(),'../Exp/Square.exp',100000.) +md = setmask(md,'../Exp/SquareShelf.exp','') +md = parameterize(md,'../Par/SquareSheetShelf.py') +md.initialization.vx[:] = 1. +md.initialization.vy[:] = 1. +md.geometry.thickness[:] = 500. - md.mesh.x / 10000. +md.geometry.bed = -100. - md.mesh.x / 1000. +md.geometry.base = -md.geometry.thickness * md.materials.rho_ice / md.materials.rho_water +md.mask.groundedice_levelset = md.geometry.thickness + md.materials.rho_water / md.materials.rho_ice * md.geometry.bed +pos = np.array(np.where(md.mask.groundedice_levelset >= 0.)) +md.geometry.base[pos] = md.geometry.bed[pos] +md.geometry.surface = md.geometry.base + md.geometry.thickness +md = setflowequation(md,'SSA','all') + +#Boundary conditions: +md.mask.ice_levelset = -np.ones((md.mesh.numberofvertices,)) +md.mask.ice_levelset[np.where(md.mesh.x == max(md.mesh.x))] = 0. +md.stressbalance.spcvx[:] = float('Nan') +md.stressbalance.spcvy[:] = float('Nan') +md.stressbalance.spcvz[:] = float('Nan') +posA = np.intersect1d(np.array(np.where(md.mesh.y < 1000000.1)),np.array(np.where(md.mesh.y > 999999.9))) +posB = np.intersect1d(np.array(np.where(md.mesh.y < 0.1)),np.array(np.where(md.mesh.y > -0.1))) +pos = np.unique(np.concatenate((posA,posB))) +md.stressbalance.spcvy[pos] = 0. +pos2 = np.intersect1d(np.array(np.where(md.mesh.x < 0.1)), np.array(np.where(md.mesh.x > -0.1))) +md.stressbalance.spcvx[pos2] = 0. +md.stressbalance.spcvy[pos2] = 0. + +md.materials.rheology_B = 1. / ((10**-25)**(1./3.)) * np.ones((md.mesh.numberofvertices,)) +md.materials.rheology_law = 'None' +md.friction.coefficient[:] = np.sqrt(1e7) * np.ones((md.mesh.numberofvertices,)) +md.friction.p = 3. * np.ones((md.mesh.numberofelements,)) +md.smb.mass_balance[:] = 1. +md.basalforcings.groundedice_melting_rate[:] = 0. +md.basalforcings.floatingice_melting_rate[:] = 30. +md.transient.isthermal = 0 +md.transient.isstressbalance = 1 +md.transient.isgroundingline = 1 +md.transient.ismasstransport = 1 +md.transient.issmb = 1 +md.transient.requested_outputs = ['default','BasalforcingsFloatingiceMeltingRate'] +md.groundingline.migration = 'SubelementMigration2' +md.timestepping.final_time = 30. +md.timestepping.time_step = 10. + +md.cluster = generic('name',gethostname(),'np',3) +md = solve(md,'Transient') + +#Fields and tolerances to track changes +field_names = [ + 'Bed1','Surface1','Thickness1','Floatingice1','Vx1','Vy1','Pressure1','FloatingiceMeltingrate1', + 'Bed2','Surface2','Thickness2','Floatingice2','Vx2','Vy2','Pressure2','FloatingiceMeltingrate2', + 'Bed3','Surface3','Thickness3','Floatingice3','Vx3','Vy3','Pressure3','FloatingiceMeltingrate3'] +field_tolerances = [ + 2e-11,5e-12,2e-11,1e-11,5e-10,1e-08,1e-13,1e-13, + 3e-11,3e-11,9e-10,7e-11,1e-09,5e-08,1e-10,1e-13, + 1e-10,3e-11,1e-10,7e-11,1e-09,5e-08,1e-10,1e-13] +field_values = [ + md.results.TransientSolution[0].Base, + md.results.TransientSolution[0].Surface, + md.results.TransientSolution[0].Thickness, + md.results.TransientSolution[0].MaskGroundediceLevelset, + md.results.TransientSolution[0].Vx, + md.results.TransientSolution[0].Vy, + md.results.TransientSolution[0].Pressure, + md.results.TransientSolution[0].BasalforcingsFloatingiceMeltingRate, + md.results.TransientSolution[1].Base, + md.results.TransientSolution[1].Surface, + md.results.TransientSolution[1].Thickness, + md.results.TransientSolution[1].MaskGroundediceLevelset, + md.results.TransientSolution[1].Vx, + md.results.TransientSolution[1].Vy, + md.results.TransientSolution[1].Pressure, + md.results.TransientSolution[1].BasalforcingsFloatingiceMeltingRate, + md.results.TransientSolution[2].Base, + md.results.TransientSolution[2].Surface, + md.results.TransientSolution[2].Thickness, + md.results.TransientSolution[2].MaskGroundediceLevelset, + md.results.TransientSolution[2].Vx, + md.results.TransientSolution[2].Vy, + md.results.TransientSolution[2].Pressure, + md.results.TransientSolution[2].BasalforcingsFloatingiceMeltingRate, + ] Index: /issm/trunk-jpl/test/NightlyRun/test442.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test442.py (revision 22267) +++ /issm/trunk-jpl/test/NightlyRun/test442.py (revision 22267) @@ -0,0 +1,97 @@ +#Test Name: MISMIP3DHO +import numpy as np +from model import * +from socket import gethostname +from triangle import * +from setmask import * +from parameterize import * +from setflowequation import * +from solve import * + +md = triangle(model(),'../Exp/Square.exp',100000.) +md = setmask(md,'../Exp/SquareShelf.exp','') +md = parameterize(md,'../Par/SquareSheetShelf.py') +md.initialization.vx[:] = 1. +md.initialization.vy[:] = 1. +md.geometry.thickness[:] = 500. - md.mesh.x / 10000. +md.geometry.bed = -100. - md.mesh.x / 1000. +md.geometry.base = -md.geometry.thickness * md.materials.rho_ice / md.materials.rho_water +md.mask.groundedice_levelset = md.geometry.thickness + md.materials.rho_water / md.materials.rho_ice * md.geometry.bed +pos = np.where(md.mask.groundedice_levelset >= 0.) +md.geometry.base[pos] = md.geometry.bed[pos] +md.geometry.surface = md.geometry.base + md.geometry.thickness +md = md.extrude(4,1.) +md = setflowequation(md,'HO','all') + +#Boundary conditions: +md.mask.ice_levelset = -np.ones((md.mesh.numberofvertices,)) +md.mask.ice_levelset[np.where(md.mesh.x == max(md.mesh.x))] = 0. +md.stressbalance.spcvx[:] = float('Nan') +md.stressbalance.spcvy[:] = float('Nan') +md.stressbalance.spcvz[:] = float('Nan') +posA = np.intersect1d(np.array(np.where(md.mesh.y < 1000000.1)),np.array(np.where(md.mesh.y > 999999.9))) +posB = np.intersect1d(np.array(np.where(md.mesh.y < 0.1)),np.array(np.where(md.mesh.y > -0.1))) +pos = np.unique(np.concatenate((posA,posB))) +md.stressbalance.spcvy[pos] = 0. +pos2 = np.intersect1d(np.array(np.where(md.mesh.x < 0.1)), np.array(np.where(md.mesh.x > -0.1))) +md.stressbalance.spcvx[pos2] = 0. +md.stressbalance.spcvy[pos2] = 0. + +md.materials.rheology_B = 1. / ((10**-25)**(1./3.)) * np.ones((md.mesh.numberofvertices,)) +md.materials.rheology_law = 'None' +md.friction.coefficient[:] = np.sqrt(1e7) * np.ones((md.mesh.numberofvertices,)) +md.friction.p = 3. * np.ones((md.mesh.numberofelements,)) +md.smb.mass_balance[:] = 1. +md.basalforcings.groundedice_melting_rate[:] = 0. +md.basalforcings.floatingice_melting_rate[:] = 30. +md.transient.isthermal = 0 +md.transient.isstressbalance = 1 +md.transient.isgroundingline = 1 +md.transient.ismasstransport = 1 +md.transient.issmb = 1 +md.transient.requested_outputs = ['default','BasalforcingsFloatingiceMeltingRate'] +md.groundingline.migration = 'SubelementMigration2' +md.timestepping.final_time = 30 +md.timestepping.time_step = 10 + +md.cluster = generic('name',gethostname(),'np',3) +md = solve(md,'Transient') + +#Fields and tolerances to track changes +field_names = [ + 'Bed1','Surface1','Thickness1','Floatingice1','Vx1','Vy1','Vz1','Pressure1','FloatingiceMeltingrate1', + 'Bed2','Surface2','Thickness2','Floatingice2','Vx2','Vy2','Vz2','Pressure2','FloatingiceMeltingrate2', + 'Bed3','Surface3','Thickness3','Floatingice3','Vx3','Vy3','Vz3','Pressure3','FloatingiceMeltingrate3'] +field_tolerances = [ + 2e-11,5e-12,2e-11,1e-11,5e-10,3e-08,6e-10,1e-13,1e-13, + 3e-11,3e-11,9e-10,7e-11,6e-09,1e-07,1e-09,1e-10,1e-13, + 1e-8,2e-08,7e-09,2e-7 ,1e-03,8e-04,2e-09,1e-10,1e-13] +field_values = [ + md.results.TransientSolution[0].Base, + md.results.TransientSolution[0].Surface, + md.results.TransientSolution[0].Thickness, + md.results.TransientSolution[0].MaskGroundediceLevelset, + md.results.TransientSolution[0].Vx, + md.results.TransientSolution[0].Vy, + md.results.TransientSolution[0].Vz, + md.results.TransientSolution[0].Pressure, + md.results.TransientSolution[0].BasalforcingsFloatingiceMeltingRate, + md.results.TransientSolution[1].Base, + md.results.TransientSolution[1].Surface, + md.results.TransientSolution[1].Thickness, + md.results.TransientSolution[1].MaskGroundediceLevelset, + md.results.TransientSolution[1].Vx, + md.results.TransientSolution[1].Vy, + md.results.TransientSolution[1].Vz, + md.results.TransientSolution[1].Pressure, + md.results.TransientSolution[1].BasalforcingsFloatingiceMeltingRate, + md.results.TransientSolution[2].Base, + md.results.TransientSolution[2].Surface, + md.results.TransientSolution[2].Thickness, + md.results.TransientSolution[2].MaskGroundediceLevelset, + md.results.TransientSolution[2].Vx, + md.results.TransientSolution[2].Vy, + md.results.TransientSolution[2].Vz, + md.results.TransientSolution[2].Pressure, + md.results.TransientSolution[2].BasalforcingsFloatingiceMeltingRate, + ] Index: /issm/trunk-jpl/test/NightlyRun/test460.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test460.py (revision 22267) +++ /issm/trunk-jpl/test/NightlyRun/test460.py (revision 22267) @@ -0,0 +1,38 @@ +#Test Name: SquareSheetShelfStressFSEstar +import numpy as np +from model import * +from socket import gethostname +from triangle import * +from setmask import * +from parameterize import * +from setflowequation import * +from solve import * +from matestar import * + +md=triangle(model(),'../Exp/Square.exp',180000.) +md=setmask(md,'../Exp/SquareShelf.exp','') +md=parameterize(md,'../Par/SquareSheetShelf.py') +md = md.extrude(3,1.) +md.materials = matestar() +md.materials.rheology_B = 3.15e8 * np.ones((md.mesh.numberofvertices,)) +md.materials.rheology_Ec = np.ones((md.mesh.numberofvertices,)) +md.materials.rheology_Es = 3 * np.ones((md.mesh.numberofvertices,)) +md.cluster = generic('name',gethostname(),'np',3) + +#Go solve +field_names=[] +field_tolerances=[] +field_values=[] +#md.initialization.pressure=md.constants.g*md.materials.rho_ice*(md.geometry.surface-md.mesh.y); +for i in ['SSA','HO','FS']: + md = setflowequation(md,i,'all') + md = solve(md,'Stressbalance') + field_names = field_names + ['Vx'+i,'Vy'+i,'Vz'+i,'Vel'+i,'Pressure'+i] + field_tolerances = field_tolerances + [6e-07,6e-07,2e-06,1e-06,5e-07] + field_values = field_values + [ + md.results.StressbalanceSolution.Vx, + md.results.StressbalanceSolution.Vy, + md.results.StressbalanceSolution.Vz, + md.results.StressbalanceSolution.Vel, + md.results.StressbalanceSolution.Pressure, + ] Index: /issm/trunk-jpl/test/NightlyRun/test461.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test461.py (revision 22267) +++ /issm/trunk-jpl/test/NightlyRun/test461.py (revision 22267) @@ -0,0 +1,53 @@ +#Test Name: SquareSheetShelfThermalFSEstar +import numpy as np +from model import * +from socket import gethostname +from triangle import * +from setmask import * +from parameterize import * +from setflowequation import * +from solve import * +from matestar import * + +md = triangle(model(),'../Exp/Square.exp',180000.) +md = setmask(md,'../Exp/SquareShelf.exp','') +md = parameterize(md,'../Par/SquareSheetShelf.py') +md = md.extrude(3,1.) +md.materials = matestar() +md.materials.rheology_B = 3.15e8 * np.ones((md.mesh.numberofvertices,)) +md.materials.rheology_Ec = np.ones((md.mesh.numberofvertices,)) +md.materials.rheology_Es = 3. * np.ones((md.mesh.numberofvertices,)) + +md = setflowequation(md,'FS','all') +md.initialization.waterfraction = np.zeros((md.mesh.numberofvertices,)) +md.initialization.watercolumn = np.zeros((md.mesh.numberofvertices,)) +md.transient.isstressbalance = 0 +md.transient.ismasstransport = 0 +md.transient.issmb = 1 +md.transient.isthermal = 1 +md.transient.isgroundingline = 0 +md.thermal.isenthalpy = 1 +md.thermal.isdynamicbasalspc = 1 +md.cluster = generic('name',gethostname(),'np',3) +md = solve(md,'Transient') + +#Fields and tolerances to track changes +field_names = [ + 'Enthalpy1','Waterfraction1','Temperature1', + 'Enthalpy2','Waterfraction2','Temperature2', + 'Enthalpy3','Waterfraction3','Temperature3'] +field_tolerances = [ + 1e-12,1e-11,1e-12, + 1e-12,1e-10,1e-12, + 1e-12,1e-9,1e-12] +field_values = [ + md.results.TransientSolution[0].Enthalpy, + md.results.TransientSolution[0].Waterfraction, + md.results.TransientSolution[0].Temperature, + md.results.TransientSolution[1].Enthalpy, + md.results.TransientSolution[1].Waterfraction, + md.results.TransientSolution[1].Temperature, + md.results.TransientSolution[2].Enthalpy, + md.results.TransientSolution[2].Waterfraction, + md.results.TransientSolution[2].Temperature, + ] Index: /issm/trunk-jpl/test/NightlyRun/test462.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test462.py (revision 22267) +++ /issm/trunk-jpl/test/NightlyRun/test462.py (revision 22267) @@ -0,0 +1,49 @@ +#Test Name: SquareSheetShelfAmrBamgField +import numpy as np +from model import * +from socket import gethostname +from triangle import * +from setmask import * +from parameterize import * +from setflowequation import * +from solve import * + +md = triangle(model(),'../Exp/Square.exp',150000.) +md = setmask(md,'../Exp/SquareShelf.exp','') +md = parameterize(md,'../Par/SquareSheetShelf.py') +md = setflowequation(md,'SSA','all') +md.cluster = generic('name',gethostname(),'np',3) +md.transient.isstressbalance = 1 +md.transient.ismasstransport = 1 +md.transient.issmb = 0 +md.transient.isthermal = 0 +md.transient.isgroundingline = 0 +#amr bamg settings, just field +md.amr.hmin = 10000 +md.amr.hmax = 100000 +md.amr.fieldname = 'Vel' +md.amr.keepmetric = 1 +md.amr.gradation = 1.2 +md.amr.groundingline_resolution = 2000 +md.amr.groundingline_distance = 0 +md.amr.icefront_resolution = 1000 +md.amr.icefront_distance = 0 +md.amr.thicknesserror_resolution = 1000 +md.amr.thicknesserror_threshold = 0 +md.amr.deviatoricerror_resolution = 1000 +md.amr.deviatoricerror_threshold = 0 +md.transient.amr_frequency = 1 +md.timestepping.start_time = 0 +md.timestepping.final_time = 3 +md.timestepping.time_step = 1 +md = solve(md,'Transient') + +#Fields and tolerances to track changes +field_names = ['Vx','Vy','Vel','Pressure'] +field_tolerances = [1e-13,1e-13,1e-13,1e-13] +field_values = [ + md.results.TransientSolution[2].Vx, + md.results.TransientSolution[2].Vy, + md.results.TransientSolution[2].Vel, + md.results.TransientSolution[2].Pressure, + ] Index: /issm/trunk-jpl/test/NightlyRun/test463.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test463.py (revision 22267) +++ /issm/trunk-jpl/test/NightlyRun/test463.py (revision 22267) @@ -0,0 +1,49 @@ +#Test Name: SquareSheetShelfAmrBamgGroundingline +import numpy as np +from model import * +from socket import gethostname +from triangle import * +from setmask import * +from parameterize import * +from setflowequation import * +from solve import * + +md = triangle(model(),'../Exp/Square.exp',150000.) +md = setmask(md,'../Exp/SquareShelf.exp','') +md = parameterize(md,'../Par/SquareSheetShelf.py') +md = setflowequation(md,'SSA','all') +md.cluster = generic('name',gethostname(),'np',3) +md.transient.isstressbalance = 1 +md.transient.ismasstransport = 1 +md.transient.issmb = 0 +md.transient.isthermal = 0 +md.transient.isgroundingline = 0 +#amr bamg settings, just grounding line +md.amr.hmin = 10000 +md.amr.hmax = 100000 +md.amr.fieldname = 'None' +md.amr.keepmetric = 0 +md.amr.gradation = 1.2 +md.amr.groundingline_resolution = 12000 +md.amr.groundingline_distance = 100000 +md.amr.icefront_resolution = 1000 +md.amr.icefront_distance = 0 +md.amr.thicknesserror_resolution = 1000 +md.amr.thicknesserror_threshold = 0 +md.amr.deviatoricerror_resolution = 1000 +md.amr.deviatoricerror_threshold = 0 +md.transient.amr_frequency = 1 +md.timestepping.start_time = 0 +md.timestepping.final_time = 3 +md.timestepping.time_step = 1 +md = solve(md,'Transient') + +#Fields and tolerances to track changes +field_names = ['Vx','Vy','Vel','Pressure'] +field_tolerances = [1e-8,1e-8,1e-8,1e-8] +field_values = [ + md.results.TransientSolution[2].Vx, + md.results.TransientSolution[2].Vy, + md.results.TransientSolution[2].Vel, + md.results.TransientSolution[2].Pressure, + ] Index: /issm/trunk-jpl/test/NightlyRun/test464.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test464.py (revision 22267) +++ /issm/trunk-jpl/test/NightlyRun/test464.py (revision 22267) @@ -0,0 +1,49 @@ +#Test Name: SquareSheetShelfAmrBamgIceFront +import numpy as np +from model import * +from socket import gethostname +from triangle import * +from setmask import * +from parameterize import * +from setflowequation import * +from solve import * + +md = triangle(model(),'../Exp/Square.exp',150000.) +md = setmask(md,'../Exp/SquareShelf.exp','') +md = parameterize(md,'../Par/SquareSheetShelf.py') +md = setflowequation(md,'SSA','all') +md.cluster = generic('name',gethostname(),'np',3) +md.transient.isstressbalance = 1 +md.transient.ismasstransport = 1 +md.transient.issmb = 0 +md.transient.isthermal = 0 +md.transient.isgroundingline = 0 +#amr bamg settings, just ice front +md.amr.hmin = 10000 +md.amr.hmax = 100000 +md.amr.fieldname = 'None' +md.amr.keepmetric = 0 +md.amr.gradation = 1.2 +md.amr.groundingline_resolution = 12000 +md.amr.groundingline_distance = 0 +md.amr.icefront_resolution = 12000 +md.amr.icefront_distance = 100000 +md.amr.thicknesserror_resolution = 1000 +md.amr.thicknesserror_threshold = 0 +md.amr.deviatoricerror_resolution = 1000 +md.amr.deviatoricerror_threshold = 0 +md.transient.amr_frequency = 1 +md.timestepping.start_time = 0 +md.timestepping.final_time = 3 +md.timestepping.time_step = 1 +md = solve(md,'Transient') + +#Fields and tolerances to track changes +field_names = ['Vx','Vy','Vel','Pressure'] +field_tolerances = [1e-13,1e-13,1e-13,1e-13] +field_values = [ + md.results.TransientSolution[2].Vx, + md.results.TransientSolution[2].Vy, + md.results.TransientSolution[2].Vel, + md.results.TransientSolution[2].Pressure, + ] Index: /issm/trunk-jpl/test/NightlyRun/test465.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test465.py (revision 22267) +++ /issm/trunk-jpl/test/NightlyRun/test465.py (revision 22267) @@ -0,0 +1,49 @@ +#Test Name: SquareSheetShelfAmrBamgAll +import numpy as np +from model import * +from socket import gethostname +from triangle import * +from setmask import * +from parameterize import * +from setflowequation import * +from solve import * + +md = triangle(model(),'../Exp/Square.exp',150000.) +md = setmask(md,'../Exp/SquareShelf.exp','') +md = parameterize(md,'../Par/SquareSheetShelf.py') +md = setflowequation(md,'SSA','all') +md.cluster = generic('name',gethostname(),'np',3) +md.transient.isstressbalance = 1 +md.transient.ismasstransport = 1 +md.transient.issmb = 0 +md.transient.isthermal = 0 +md.transient.isgroundingline = 0 +#amr bamg settings, field, grounding line and ice front +md.amr.hmin = 20000 +md.amr.hmax = 100000 +md.amr.fieldname = 'Vel' +md.amr.keepmetric = 0 +md.amr.gradation = 1.2 +md.amr.groundingline_resolution = 10000 +md.amr.groundingline_distance = 100000 +md.amr.icefront_resolution = 15000 +md.amr.icefront_distance = 100000 +md.amr.thicknesserror_resolution = 1000 +md.amr.thicknesserror_threshold = 0 +md.amr.deviatoricerror_resolution = 1000 +md.amr.deviatoricerror_threshold = 0 +md.transient.amr_frequency = 1 +md.timestepping.start_time = 0 +md.timestepping.final_time = 3 +md.timestepping.time_step = 1 +md = solve(md,'Transient') + +#Fields and tolerances to track changes +field_names = ['Vx','Vy','Vel','Pressure'] +field_tolerances = [1e-8,1e-8,1e-8,1e-8] +field_values = [ + md.results.TransientSolution[2].Vx, + md.results.TransientSolution[2].Vy, + md.results.TransientSolution[2].Vel, + md.results.TransientSolution[2].Pressure, + ] Index: /issm/trunk-jpl/test/NightlyRun/test540.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test540.py (revision 22267) +++ /issm/trunk-jpl/test/NightlyRun/test540.py (revision 22267) @@ -0,0 +1,66 @@ +#Test Name: PigTranCalvingDevSSA2d +from model import * +from socket import gethostname +from triangle import * +from setmask import * +from parameterize import * +from setflowequation import * +from solve import * +from calvingvonmises import * + +md = triangle(model(),'../Exp/Pig.exp',10000.) +md = setmask(md,'../Exp/PigShelves.exp','../Exp/PigIslands.exp') +md = parameterize(md,'../Par/Pig.py') +md = setflowequation(md,'SSA','all') +md.timestepping.time_step = 2 +md.timestepping.final_time = 50 + +#calving parameters +md.mask.ice_levelset = 1e4 * (md.mask.ice_levelset + 0.5) +md.calving = calvingvonmises() +md.calving.meltingrate = np.zeros((md.mesh.numberofvertices,)) +md.transient.ismovingfront = 1 +md.levelset.spclevelset = float('NaN') * np.ones((md.mesh.numberofvertices,)) +pos = np.where(md.mesh.vertexonboundary) +md.levelset.spclevelset[pos] = md.mask.ice_levelset[pos] + +md.cluster = generic('name',gethostname(),'np',2) +md = solve(md,'Transient') + +#Fields and tolerances to track changes +field_names = [ + 'Vx1' ,'Vy1' ,'Vel1' ,'Pressure1' ,'Bed1' ,'Surface1' ,'Thickness1' ,'MaskIceLevelset1' , + 'Vx2' ,'Vy2' ,'Vel2' ,'Pressure2' ,'Bed2' ,'Surface2' ,'Thickness2' ,'MaskIceLevelset2' , + 'Vx10','Vy10','Vel10','Pressure10','Bed10','Surface10','Thickness10','MaskIceLevelset10', + ] +field_tolerances = [ + 1e-12,2e-12,2e-12,1e-13,1e-13,1e-13,1e-13,1e-13, + 1e-12,1e-12,1e-12,1e-13,1e-13,1e-13,1e-13,1e-12, + 1e-11,1e-11,1e-11,1e-11,1e-11,1e-11,1e-11,1e-9, + ] +field_values = [ + md.results.TransientSolution[0].Vx, + md.results.TransientSolution[0].Vy, + md.results.TransientSolution[0].Vel, + md.results.TransientSolution[0].Pressure, + md.results.TransientSolution[0].Base, + md.results.TransientSolution[0].Surface, + md.results.TransientSolution[0].Thickness, + md.results.TransientSolution[0].MaskIceLevelset, + md.results.TransientSolution[1].Vx, + md.results.TransientSolution[1].Vy, + md.results.TransientSolution[1].Vel, + md.results.TransientSolution[1].Pressure, + md.results.TransientSolution[1].Base, + md.results.TransientSolution[1].Surface, + md.results.TransientSolution[1].Thickness, + md.results.TransientSolution[1].MaskIceLevelset, + md.results.TransientSolution[9].Vx, + md.results.TransientSolution[9].Vy, + md.results.TransientSolution[9].Vel, + md.results.TransientSolution[9].Pressure, + md.results.TransientSolution[9].Base, + md.results.TransientSolution[9].Surface, + md.results.TransientSolution[9].Thickness, + md.results.TransientSolution[9].MaskIceLevelset, + ] Index: /issm/trunk-jpl/test/NightlyRun/test701.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test701.py (revision 22267) +++ /issm/trunk-jpl/test/NightlyRun/test701.py (revision 22267) @@ -0,0 +1,75 @@ +#Test Name: FlowbandFSshelf +import numpy as np +from scipy.interpolate import interp1d +from model import * +from solve import * +from setflowequation import * +from bamgflowband import * +from paterson import * + +x = np.arange(1,3001,100).T +h = np.linspace(1000,300,np.size(x)).T +b = -917./1023.*h + +md = bamgflowband(model(),x,b+h,b,'hmax',80.) +print md.mesh.numberofvertices + +#Geometry #interp1d returns a function to be called on md.mesh.x +md.geometry.surface = interp1d(x,b+h)(md.mesh.x) +md.geometry.base = interp1d(x,b)(md.mesh.x) +md.geometry.thickness = md.geometry.surface-md.geometry.base + +#mask +md.mask.ice_levelset = -np.ones((md.mesh.numberofvertices,)) +md.mask.ice_levelset[np.where(md.mesh.vertexflags(2))] = 0. +md.mask.groundedice_levelset = np.zeros((md.mesh.numberofvertices,)) - 0.5 + +#materials +md.initialization.temperature = (273.-20.)*np.ones((md.mesh.numberofvertices,)) +md.materials.rheology_B = paterson(md.initialization.temperature) +md.materials.rheology_n = 3.*np.ones((md.mesh.numberofelements,)) + +#friction +md.friction.coefficient = np.zeros((md.mesh.numberofvertices,)) +md.friction.coefficient[np.where(md.mesh.vertexflags(1))] = 20. +md.friction.p = np.ones((md.mesh.numberofelements,)) +md.friction.q = np.ones((md.mesh.numberofelements,)) + +#Boundary conditions +md.stressbalance.referential = float('NaN')*np.ones((md.mesh.numberofvertices,6)) +md.stressbalance.loadingforce = 0. * np.ones((md.mesh.numberofvertices,3)) +md.stressbalance.spcvx = float('NaN') * np.ones((md.mesh.numberofvertices,)) +md.stressbalance.spcvy = float('NaN') * np.ones((md.mesh.numberofvertices,)) +md.stressbalance.spcvz = float('NaN') * np.ones((md.mesh.numberofvertices,)) +md.stressbalance.spcvx[np.where(md.mesh.vertexflags(4))] = 0. +md.stressbalance.spcvy[np.where(md.mesh.vertexflags(4))] = 0. + +#Misc +md = setflowequation(md,'FS','all') +md.stressbalance.abstol = float('NaN') +#md.stressbalance.reltol = 10**-16 +md.stressbalance.FSreconditioning = 1. +md.stressbalance.maxiter = 20 +md.flowequation.augmented_lagrangian_r = 10000. +md.miscellaneous.name = 'test701' +md.verbose = verbose('convergence',True) +md.cluster = generic('np',2) + +#Go solve +field_names = [] +field_tolerances = [] +field_values = [] +#md.initialization.pressure = md.constants.g*md.materials.rho_ice*(md.geometry.surface-md.mesh.y) +for i in ['MINI','MINIcondensed','TaylorHood','LATaylorHood','CrouzeixRaviart','LACrouzeixRaviart']: + print ' ' + print '======Testing ' +i+ ' Full-Stokes Finite element=====' + md.flowequation.fe_FS = i + md = solve(md,'Stressbalance') + field_names = field_names + [['Vx'+i],['Vy'+i],['Vel'+i],['Pressure'+i]] + field_tolerances = field_tolerances + [9e-5,9e-5,9e-5,1e-10] + field_values = field_values + [ + md.results.StressbalanceSolution.Vx, + md.results.StressbalanceSolution.Vy, + md.results.StressbalanceSolution.Vel, + md.results.StressbalanceSolution.Pressure + ] Index: /issm/trunk-jpl/test/NightlyRun/test703.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test703.py (revision 22267) +++ /issm/trunk-jpl/test/NightlyRun/test703.py (revision 22267) @@ -0,0 +1,160 @@ +#Test Name: FlowbandFSsheetshelfTrans +import numpy as np +from scipy.interpolate import interp1d +import copy +from model import * +from socket import gethostname +from triangle import * +from meshconvert import * +from setmask import * +from parameterize import * +from setflowequation import * +from solve import * +from NowickiProfile import * +from bamgflowband import * +from paterson import * + +#mesh parameters +x = np.arange(-5,5.5,.5).T +[b,h,sea] = NowickiProfile(x) +x = x * 10**3 +h = h * 10**3 +b = (b-sea) * 10**3 + +#mesh domain +md = bamgflowband(model(),x,b+h,b,'hmax',150) +print md.mesh.numberofvertices + +#geometry +md.geometry.surface = interp1d(x,b+h)(md.mesh.x) +md.geometry.base = interp1d(x,b)(md.mesh.x) +md.geometry.thickness = md.geometry.surface - md.geometry.base + +#mask +md.mask.ice_levelset = -np.ones((md.mesh.numberofvertices,)) +md.mask.ice_levelset[np.where(md.mesh.vertexflags(2))] = 0 +md.mask.groundedice_levelset = np.zeros((md.mesh.numberofvertices,)) - 0.5 + +#materials +md.initialization.temperature = (273. - 20.) * np.ones((md.mesh.numberofvertices,)) +md.materials.rheology_B = paterson(md.initialization.temperature) +md.materials.rheology_n = 3 * np.ones((md.mesh.numberofelements,)) + +#friction +md.friction.coefficient = np.zeros((md.mesh.numberofvertices,)) +md.friction.coefficient[np.where(md.mesh.vertexflags(1))] = 20 +md.friction.p = np.ones((md.mesh.numberofelements,)) +md.friction.q = np.ones((md.mesh.numberofelements,)) + +#boundary conditions +md.stressbalance.spcvx = float('NaN') * np.ones((md.mesh.numberofvertices,)) +md.stressbalance.spcvy = float('NaN') * np.ones((md.mesh.numberofvertices,)) +md.stressbalance.spcvz = float('NaN') * np.ones((md.mesh.numberofvertices,)) +md.stressbalance.referential = float('NaN') * np.ones((md.mesh.numberofvertices,6)) +md.stressbalance.loadingforce = 0 * np.ones((md.mesh.numberofvertices,3)) +md.stressbalance.spcvx[np.where(md.mesh.vertexflags(4))] = 800. +md.stressbalance.spcvy[np.where(md.mesh.vertexflags(4))] = 0. + +#print md.stressbalance.referential +#print md.stressbalance.loadingforce +#print md.stressbalance.spcvx +#print md.stressbalance.spcvy +#print md.stressbalance.spcvz +#print np.shape(md.stressbalance.referential) +#print np.shape(md.stressbalance.loadingforce) +#print np.shape(md.stressbalance.spcvx) +#print np.shape(md.stressbalance.spcvy) +#print np.shape(md.stressbalance.spcvz) + +#Misc +md = setflowequation(md,'FS','all') +md.flowequation.fe_FS = 'TaylorHood' +md.stressbalance.abstol = float('NaN') +md.miscellaneous.name = 'test703' + +#Transient settings +md.timestepping.time_step = 0.000001 +md.timestepping.final_time = 0.000005 +md.stressbalance.shelf_dampening = 1. +md.smb.mass_balance = np.zeros((md.mesh.numberofvertices,)) +md.basalforcings.groundedice_melting_rate = np.zeros((md.mesh.numberofvertices,)) +md.basalforcings.floatingice_melting_rate = np.zeros((md.mesh.numberofvertices,)) +md.basalforcings.geothermalflux = np.zeros((md.mesh.numberofvertices,)) +posb = np.intersect1d(np.where(md.mesh.x > 0.), np.where(md.mesh.vertexonbase)) +md.basalforcings.groundedice_melting_rate[posb] = 18. +md.basalforcings.floatingice_melting_rate[posb] = 18. +md.initialization.vx = np.zeros((md.mesh.numberofvertices,)) +md.initialization.vy = np.zeros((md.mesh.numberofvertices,)) +md.initialization.pressure = np.zeros((md.mesh.numberofvertices,)) +md.masstransport.spcthickness = float('NaN') * np.ones((md.mesh.numberofvertices,)) +md.thermal.spctemperature = float('NaN') * np.ones((md.mesh.numberofvertices,)) +md.transient.isthermal = 0 +md.masstransport.isfreesurface = 1 + +#Go solve +md.cluster = generic('np',3) +md.stressbalance.shelf_dampening = 1 +md1 = solve(md,'Transient') + +md.stressbalance.shelf_dampening = 0 +md = solve(md,'Transient') + +#Fields and tolerances to track changes +field_names = [ + 'Vx1','Vy1','Vel1','Pressure1','Bed1','Surface1','Thickness1', + 'Vx2','Vy2','Vel2','Pressure2','Bed2','Surface2','Thickness2', + 'Vx3','Vy3','Vel3','Pressure3','Bed3','Surface3','Thickness3', + 'Vx1_damp','Vy1_damp','Vel1_damp','Pressure1_damp','Bed1_damp','Surface1_damp','Thickness1_damp', + 'Vx2_damp','Vy2_damp','Vel2_damp','Pressure2_damp','Bed2_damp','Surface2_damp','Thickness2_damp', + 'Vx3_damp','Vy3_damp','Vel3_damp','Pressure3_damp','Bed3_damp','Surface3_damp','Thickness3_damp'] +field_tolerances = [ + 2e-08,2e-08,2e-08,1e-08,1e-10,1e-10,1e-10, + 2e-08,2e-08,2e-08,1e-08,1e-10,1e-10,1e-10, + 2e-08,2e-08,2e-08,1e-08,1e-10,1e-10,1e-10, + 5e-08,5e-08,5e-08,1e-08,1e-10,1e-10,1e-10, + 5e-08,5e-08,5e-08,1e-08,1e-10,1e-10,1e-10, + 5e-08,5e-08,5e-08,1e-08,1e-10,1e-10,1e-10] +field_values = [ + md.results.TransientSolution[0].Vx, + md.results.TransientSolution[0].Vy, + md.results.TransientSolution[0].Vel, + md.results.TransientSolution[0].Pressure, + md.results.TransientSolution[0].Base, + md.results.TransientSolution[0].Surface, + md.results.TransientSolution[0].Thickness, + md.results.TransientSolution[1].Vx, + md.results.TransientSolution[1].Vy, + md.results.TransientSolution[1].Vel, + md.results.TransientSolution[1].Pressure, + md.results.TransientSolution[1].Base, + md.results.TransientSolution[1].Surface, + md.results.TransientSolution[1].Thickness, + md.results.TransientSolution[2].Vx, + md.results.TransientSolution[2].Vy, + md.results.TransientSolution[2].Vel, + md.results.TransientSolution[2].Pressure, + md.results.TransientSolution[2].Base, + md.results.TransientSolution[2].Surface, + md.results.TransientSolution[2].Thickness, + md1.results.TransientSolution[0].Vx, + md1.results.TransientSolution[0].Vy, + md1.results.TransientSolution[0].Vel, + md1.results.TransientSolution[0].Pressure, + md1.results.TransientSolution[0].Base, + md1.results.TransientSolution[0].Surface, + md1.results.TransientSolution[0].Thickness, + md1.results.TransientSolution[1].Vx, + md1.results.TransientSolution[1].Vy, + md1.results.TransientSolution[1].Vel, + md1.results.TransientSolution[1].Pressure, + md1.results.TransientSolution[1].Base, + md1.results.TransientSolution[1].Surface, + md1.results.TransientSolution[1].Thickness, + md1.results.TransientSolution[2].Vx, + md1.results.TransientSolution[2].Vy, + md1.results.TransientSolution[2].Vel, + md1.results.TransientSolution[2].Pressure, + md1.results.TransientSolution[2].Base, + md1.results.TransientSolution[2].Surface, + md1.results.TransientSolution[2].Thickness, + ] Index: /issm/trunk-jpl/test/NightlyRun/test808.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test808.py (revision 22267) +++ /issm/trunk-jpl/test/NightlyRun/test808.py (revision 22267) @@ -0,0 +1,76 @@ +#Test Name: SquareShelfLevelsetCalvingSSA2dMinThickness +from model import * +from socket import gethostname +from triangle import * +from setmask import * +from parameterize import * +from setflowequation import * +from solve import * +import numpy as np +from calvingminthickness import * + +md = triangle(model(),'../Exp/Square.exp',30000.) +md = setmask(md,'all','') +md = parameterize(md,'../Par/SquareShelf.py') +md = setflowequation(md,'SSA','all') +md.cluster = generic('name',gethostname(),'np',3) + +x = md.mesh.x +xmin = min(x) +xmax = max(x) +Lx = (xmax - xmin) +alpha = 2. / 3. +md.mask.ice_levelset = -1 + 2 * (md.mesh.y > 9e5) + +md.timestepping.time_step = 1 +md.timestepping.final_time = 30 + +#Transient +md.transient.isstressbalance = 1 +md.transient.ismasstransport = 1 +md.transient.issmb = 1 +md.transient.isthermal = 0 +md.transient.isgroundingline = 0 +md.transient.isgia = 0 +md.transient.ismovingfront = 1 + +md.calving = calvingminthickness() +md.calving.min_thickness = 400 +md.calving.meltingrate = np.zeros((md.mesh.numberofvertices,)) +md.levelset.spclevelset = float('NaN')* np.ones((md.mesh.numberofvertices,)) +md.levelset.reinit_frequency = 1 + +md = solve(md,'Transient') + +#Fields and tolerances to track changes +field_names = [ + 'Vx1','Vy1','Vel1','Pressure1','Thickness1','Surface1','MaskIceLevelset1' + 'Vx2','Vy2','Vel2','Pressure2','Thickness2','Surface2','MaskIceLevelset2' + 'Vx3','Vy3','Vel3','Pressure3','Thickness3','Surface3','MaskIceLevelset3'] +field_tolerances = [ + 1e-8,1e-8,1e-8,1e-9,1e-9,1e-9,3e-9, + 1e-8,1e-8,1e-8,1e-9,1e-9,1e-9,3e-9, + 1e-8,1e-8,1e-8,1e-9,1e-9,1e-9,3e-9] +field_values = [ + md.results.TransientSolution[0].Vx, + md.results.TransientSolution[0].Vy, + md.results.TransientSolution[0].Vel, + md.results.TransientSolution[0].Pressure, + md.results.TransientSolution[0].Thickness, + md.results.TransientSolution[0].Surface, + md.results.TransientSolution[0].MaskIceLevelset, + md.results.TransientSolution[1].Vx, + md.results.TransientSolution[1].Vy, + md.results.TransientSolution[1].Vel, + md.results.TransientSolution[1].Pressure, + md.results.TransientSolution[1].Thickness, + md.results.TransientSolution[1].Surface, + md.results.TransientSolution[1].MaskIceLevelset, + md.results.TransientSolution[2].Vx, + md.results.TransientSolution[2].Vy, + md.results.TransientSolution[2].Vel, + md.results.TransientSolution[2].Pressure, + md.results.TransientSolution[2].Thickness, + md.results.TransientSolution[2].Surface, + md.results.TransientSolution[2].MaskIceLevelset + ] Index: /issm/trunk-jpl/test/Par/ISMIPE.par =================================================================== --- /issm/trunk-jpl/test/Par/ISMIPE.par (revision 22266) +++ /issm/trunk-jpl/test/Par/ISMIPE.par (revision 22267) @@ -14,4 +14,5 @@ md.geometry.surface(i)=(1.-coeff)*data(point1,3)+coeff*data(point2,3); end + md.geometry.thickness=md.geometry.surface-md.geometry.base; md.geometry.thickness(find(~md.geometry.thickness))=0.01; Index: /issm/trunk-jpl/test/Par/ISMIPE.py =================================================================== --- /issm/trunk-jpl/test/Par/ISMIPE.py (revision 22266) +++ /issm/trunk-jpl/test/Par/ISMIPE.py (revision 22267) @@ -13,7 +13,8 @@ point1=numpy.floor(y/100.) point2=numpy.minimum(point1+1,50) - coeff=(y-(point1-1.)*100.)/100. + coeff=(y-(point1)*100.)/100. md.geometry.base[i]=(1.-coeff)*data[point1,1]+coeff*data[point2,1] md.geometry.surface[i]=(1.-coeff)*data[point1,2]+coeff*data[point2,2] + md.geometry.thickness=md.geometry.surface-md.geometry.base md.geometry.thickness[numpy.nonzero(numpy.logical_not(md.geometry.thickness))]=0.01 Index: /issm/trunk-jpl/test/Par/SquareThermal.py =================================================================== --- /issm/trunk-jpl/test/Par/SquareThermal.py (revision 22266) +++ /issm/trunk-jpl/test/Par/SquareThermal.py (revision 22267) @@ -26,4 +26,7 @@ print " creating temperatures" md.initialization.temperature=(273.-20.)*numpy.ones((md.mesh.numberofvertices)) +md.initialization.pressure=numpy.zeros((md.mesh.numberofvertices,)) +md.initialization.waterfraction=numpy.zeros((md.mesh.numberofvertices,)) +md.initialization.watercolumn=numpy.zeros((md.mesh.numberofvertices,)) print " creating flow law parameter" @@ -33,5 +36,7 @@ print " creating surface mass balance" md.smb.mass_balance=numpy.ones((md.mesh.numberofvertices))/md.constants.yts #1m/a -md.basalforcings.melting_rate=0.*numpy.ones((md.mesh.numberofvertices))/md.constants.yts #1m/a +#md.basalforcings.melting_rate=0.*numpy.ones((md.mesh.numberofvertices))/md.constants.yts #1m/a +md.basalforcings.groundedice_melting_rate=0.*numpy.ones((md.mesh.numberofvertices))/md.constants.yts #1m/a +md.basalforcings.floatingice_melting_rate=0.*numpy.ones((md.mesh.numberofvertices))/md.constants.yts #1m/a #Deal with boundary conditions: @@ -41,5 +46,5 @@ print " boundary conditions for thermal model" -md.thermal.spctemperature=md.initialization.temperature +md.thermal.spctemperature[:]=md.initialization.temperature md.basalforcings.geothermalflux=numpy.zeros((md.mesh.numberofvertices)) md.basalforcings.geothermalflux[numpy.nonzero(md.mask.groundedice_levelset>0.)[0]]=1.*10**-3 #1 mW/m^2