Index: ../trunk-jpl/src/m/boundaryconditions/SetIceSheetBC.m =================================================================== --- ../trunk-jpl/src/m/boundaryconditions/SetIceSheetBC.m (revision 26058) +++ ../trunk-jpl/src/m/boundaryconditions/SetIceSheetBC.m (revision 26059) @@ -30,10 +30,10 @@ %Initialize surface and basal forcings md.smb = initialize(md.smb,md); -md.basalforcings = initialize(md.basalforcings,md); +md.basalforcings = initialize(md.basalforcings,md); %Initialize ocean forcings and sealevel -md.dsl= initialize(md.dsl,md); +md.dsl = initialize(md.dsl,md); %Deal with other boundary conditions if isnan(md.balancethickness.thickening_rate), Index: ../trunk-jpl/src/m/boundaryconditions/SetIceSheetBC.py =================================================================== --- ../trunk-jpl/src/m/boundaryconditions/SetIceSheetBC.py (revision 26058) +++ ../trunk-jpl/src/m/boundaryconditions/SetIceSheetBC.py (revision 26059) @@ -2,13 +2,12 @@ def SetIceSheetBC(md): - """ - SETICESHEETBC - Create the boundary conditions for stressbalance and thermal models for an IceSheet with no Ice Front + """SETICESHEETBC - Create the boundary conditions for stressbalance and thermal models for an IceSheet with no Ice Front - Usage: - md = SetIceSheetBC(md) + Usage: + md = SetIceSheetBC(md) - See also: SETICESHELFBC, SETMARINEICESHEETBC + See also: SETICESHELFBC, SETMARINEICESHEETBC """ #node on Dirichlet @@ -32,10 +31,13 @@ #No ice front -> do nothing - #Create zeros basalforcings and smb + #Initialize surface and basal forcings md.smb.initialize(md) md.basalforcings.initialize(md) + #Initialize ocean forcings and sealevel + md.dsl.initialize(md) + #Deal with other boundary conditions if np.all(np.isnan(md.balancethickness.thickening_rate)): md.balancethickness.thickening_rate = np.zeros((md.mesh.numberofvertices)) Index: ../trunk-jpl/src/m/classes/hydrologydc.py =================================================================== --- ../trunk-jpl/src/m/classes/hydrologydc.py (revision 26058) +++ ../trunk-jpl/src/m/classes/hydrologydc.py (revision 26059) @@ -6,11 +6,10 @@ class hydrologydc(object): - """ - Hydrologydc class definition + """HYDROLOGYDC class definition Usage: - hydrologydc = hydrologydc() + hydrologydc = hydrologydc() """ def __init__(self): # {{{ @@ -48,11 +47,14 @@ self.epl_conductivity = 0 self.eplflip_lock = 0 - #set defaults self.setdefaultparameters() # }}} def __repr__(self): # {{{ + # TODO: + # - Convert all formatting to calls to .format (see any + # already converted .__repr__ method for examples) + # string = ' hydrology Dual Porous Continuum Equivalent parameters:' string = ' - general parameters' string = "%s\n%s" % (string, fielddisplay(self, 'water_compressibility', 'compressibility of water [Pa^ - 1]')) Index: ../trunk-jpl/src/m/classes/hydrologytws.m =================================================================== --- ../trunk-jpl/src/m/classes/hydrologytws.m (revision 26058) +++ ../trunk-jpl/src/m/classes/hydrologytws.m (revision 26059) @@ -23,7 +23,7 @@ end % }}} function list = defaultoutputs(self,md) % {{{ list = {''}; - end % }}} + end % }}} function self = setdefaultparameters(self) % {{{ self.requested_outputs = {'default'}; @@ -45,16 +45,16 @@ function marshall(self,prefix,md,fid) % {{{ WriteData(fid,prefix,'name','md.hydrology.model','data',2,'format','Integer'); WriteData(fid,prefix,'object',self,'fieldname','spcwatercolumn','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); - outputs = self.requested_outputs; - pos = find(ismember(outputs,'default')); - if ~isempty(pos), - outputs(pos) = []; %remove 'default' from outputs - outputs = [outputs defaultoutputs(self,md)]; %add defaults - end - WriteData(fid,prefix,'data',outputs,'name','md.hydrology.requested_outputs','format','StringArray'); + outputs = self.requested_outputs; + pos = find(ismember(outputs,'default')); + if ~isempty(pos), + outputs(pos) = []; %remove 'default' from outputs + outputs = [outputs defaultoutputs(self,md)]; %add defaults + end + WriteData(fid,prefix,'data',outputs,'name','md.hydrology.requested_outputs','format','StringArray'); end % }}} function savemodeljs(self,fid,modelname) % {{{ - + writejs1Darray(fid,[modelname '.hydrology.spcwatercolumn'],self.spcwatercolumn); end % }}} Index: ../trunk-jpl/src/m/classes/hydrologytws.py =================================================================== --- ../trunk-jpl/src/m/classes/hydrologytws.py (nonexistent) +++ ../trunk-jpl/src/m/classes/hydrologytws.py (revision 26059) @@ -0,0 +1,64 @@ +import numpy as np + +from structtoobj import structtoobj + +class hydrologytws(object): + """HYDROLOGYTWS class definition + + Usage: + hydrologytws = hydrologytws() + """ + + def __init__(self): # {{{ + self.spcwatercolumn = np.nan + self.requested_outputs = np.nan + + nargs = len(args) + if nargs == 0: + self.setdefaultparameters() + elif nargs == 1: + self = structtoobj(self, args[0]) + else: + raise RuntimeError('constructor not supported') + #}}} + + def __repr__(self): # {{{ + s = ' hydrologytws solution parameters:\n' + s += '{}\n'.format(fielddisplay(self, 'spcwatercolumn', 'water thickness constraints (NaN means no constraint) [m]')) + s += '{}\n'.format(fielddisplay(self, 'requested_outputs', 'additional outputs requested')) + return s + #}}} + + def defaultoutputs(self, md): # {{{ + return [''] + #}}} + + def setdefaultparameters(self): # {{{ + self.requested_outputs = ['defualt'] + return self + #}}} + + def extrude(self, md): # {{{ + return self + #}}} + + def checkconsistency(self, md, solution, analyses): # {{{ + # Early return + if 'HydrologyTwsAnalysis' not in analyses: + return + md = checkfield(md, 'fieldname', 'hydrology.spcwatercolumn', 'Inf', 1, 'timeseries', 1) + #}}} + + def marshall(self, prefix, md, fid): # {{{ + WriteData(fid, prefix, 'name', 'md.hydrology.model', 'data', 2, 'format', 'Integer') + WriteData(fid, prefix, 'object', self, 'fieldname', 'spcwatercolumn', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts) + outputs = self.requested_outputs + pos = find(ismember(outputs,'default')) + if not len(pos), + outputs[pos] = []; # remove 'default' from outputs + outputs.extend(defaultoutputs(self, md)) # add defaults + end + WriteData(fid, prefix, 'data', outputs, 'name', 'md.hydrology.requested_outputs', 'format', 'StringArray') + # }}} + + Index: ../trunk-jpl/src/m/classes/matenhancedice.m =================================================================== --- ../trunk-jpl/src/m/classes/matenhancedice.m (revision 26058) +++ ../trunk-jpl/src/m/classes/matenhancedice.m (revision 26059) @@ -5,26 +5,26 @@ classdef matenhancedice properties (SetAccess=public) - rho_ice = 0.; - rho_water = 0.; - rho_freshwater = 0.; - mu_water = 0.; - heatcapacity = 0.; - latentheat = 0.; - thermalconductivity = 0.; - temperateiceconductivity = 0.; + rho_ice = 0.; + rho_water = 0.; + rho_freshwater = 0.; + mu_water = 0.; + heatcapacity = 0.; + latentheat = 0.; + thermalconductivity = 0.; + temperateiceconductivity = 0.; effectiveconductivity_averaging = 0.; - meltingpoint = 0.; - beta = 0.; - mixed_layer_capacity = 0.; - thermal_exchange_velocity = 0.; - rheology_E = NaN; - rheology_B = NaN; - rheology_n = NaN; - rheology_law = ''; + meltingpoint = 0.; + beta = 0.; + mixed_layer_capacity = 0.; + thermal_exchange_velocity = 0.; + rheology_E = NaN; + rheology_B = NaN; + rheology_n = NaN; + rheology_law = ''; - %slr - earth_density = 0; + %SLC + earth_density = 0; end methods Index: ../trunk-jpl/src/m/classes/matestar.m =================================================================== --- ../trunk-jpl/src/m/classes/matestar.m (revision 26058) +++ ../trunk-jpl/src/m/classes/matestar.m (revision 26059) @@ -5,26 +5,26 @@ classdef matestar properties (SetAccess=public) - rho_ice = 0.; - rho_water = 0.; - rho_freshwater = 0.; - mu_water = 0.; - heatcapacity = 0.; - latentheat = 0.; - thermalconductivity = 0.; - temperateiceconductivity = 0.; + rho_ice = 0.; + rho_water = 0.; + rho_freshwater = 0.; + mu_water = 0.; + heatcapacity = 0.; + latentheat = 0.; + thermalconductivity = 0.; + temperateiceconductivity = 0.; effectiveconductivity_averaging = 0; - meltingpoint = 0.; - beta = 0.; - mixed_layer_capacity = 0.; - thermal_exchange_velocity = 0.; - rheology_B = NaN; - rheology_Ec = NaN; - rheology_Es = NaN; - rheology_law = ''; + meltingpoint = 0.; + beta = 0.; + mixed_layer_capacity = 0.; + thermal_exchange_velocity = 0.; + rheology_B = NaN; + rheology_Ec = NaN; + rheology_Es = NaN; + rheology_law = ''; %slc - earth_density = 0; + earth_density = 0; end methods Index: ../trunk-jpl/src/m/classes/matice.py =================================================================== --- ../trunk-jpl/src/m/classes/matice.py (revision 26058) +++ ../trunk-jpl/src/m/classes/matice.py (revision 26059) @@ -7,12 +7,11 @@ class matice(object): - ''' - MATICE class definition + """MATICE class definition Usage: matice = matice() - ''' + """ def __init__(self): #{{{ self.rho_ice = 0. @@ -32,13 +31,7 @@ self.rheology_n = np.nan self.rheology_law = '' - #giaivins - self.lithosphere_shear_modulus = 0. - self.lithosphere_density = 0. - self.mantle_shear_modulus = 0. - self.mantle_density = 0. - - #slr + #slc self.earth_density = 0 self.setdefaultparameters() #}}} @@ -61,10 +54,6 @@ s = "%s\n%s" % (s, fielddisplay(self, "rheology_B", "flow law parameter [Pa s^(1/n)]")) s = "%s\n%s" % (s, fielddisplay(self, "rheology_n", "Glen's flow law exponent")) s = "%s\n%s" % (s, fielddisplay(self, "rheology_law", "law for the temperature dependance of the rheology: 'None', 'BuddJacka', 'Cuffey', 'CuffeyTemperate', 'Paterson', 'Arrhenius', 'LliboutryDuval', 'NyeCO2', or 'NyeH2O'")) - s = "%s\n%s" % (s, fielddisplay(self, "lithosphere_shear_modulus", "Lithosphere shear modulus [Pa]")) - s = "%s\n%s" % (s, fielddisplay(self, "lithosphere_density", "Lithosphere density [g/cm^-3]")) - s = "%s\n%s" % (s, fielddisplay(self, "mantle_shear_modulus", "Mantle shear modulus [Pa]")) - s = "%s\n%s" % (s, fielddisplay(self, "mantle_density", "Mantle density [g/cm^-3]")) s = "%s\n%s" % (s, fielddisplay(self, "earth_density", "Mantle density [kg/m^-3]")) return s @@ -107,11 +96,9 @@ #available: none, paterson and arrhenius self.rheology_law = 'Paterson' - # GIA: - self.lithosphere_shear_modulus = 6.7e10 # (Pa) - self.lithosphere_density = 3.32 # (g/cm^-3) - self.mantle_shear_modulus = 1.45e11 # (Pa) - self.mantle_density = 3.34 # (g/cm^-3) + # Rheology for ice + self.rheology_B = 2.1 * 1e8 + self.rheology_n = 3 # SLR self.earth_density = 5512 # average density of the Earth, (kg/m^3) @@ -118,25 +105,18 @@ #}}} def checkconsistency(self, md, solution, analyses): #{{{ - if solution != 'SealevelriseSolution': + if solution == 'TransientSolution' and md.transient.isslc: + md = checkfield(md, 'fieldname', 'materials.earth_density', '>', 0, 'numel', [1]) + else: 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, 'universal', 1, 'NaN', 1, 'Inf', 1) - md = checkfield(md, 'fieldname', 'materials.rheology_n', '>', 0, 'size', [md.mesh.numberofelements]) + md = checkfield(md, 'fieldname', 'materials.rheology_n', '>', 0, 'universal',1, 'NaN', 1, 'Inf', 1) md = checkfield(md, 'fieldname', 'materials.rheology_law', 'values', ['None', 'BuddJacka', 'Cuffey', 'CuffeyTemperate', 'Paterson', 'Arrhenius', 'LliboutryDuval', 'NyeCO2', 'NyeH2O']) - md = checkfield(md, 'fieldname', 'materials.effectiveconductivity_averaging', 'numel', [1], 'values', [0, 1, 2]) + md = checkfield(md,'fieldname','materials.effectiveconductivity_averaging', 'numel', [1], 'values', [0, 1, 2]) - 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 #}}} @@ -158,9 +138,5 @@ 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: ../trunk-jpl/src/m/classes/surfaceload.m =================================================================== --- ../trunk-jpl/src/m/classes/surfaceload.m (revision 26058) +++ ../trunk-jpl/src/m/classes/surfaceload.m (revision 26059) @@ -19,13 +19,13 @@ end end % }}} function self = setdefaultparameters(self) % {{{ - + icethicknesschange=[]; waterheightchange=[]; otherchange=[]; - + end % }}} - function md = checkconsistency(self,md,solution,analyses) % {{{ + function md = checkconsistency(self,md,solution,analyses) % {{{ if ~ismember('SealevelchangeAnalysis',analyses) | (strcmp(solution,'TransientSolution') & md.transient.isslc==0), return; Index: ../trunk-jpl/src/m/classes/transient.m =================================================================== --- ../trunk-jpl/src/m/classes/transient.m (revision 26058) +++ ../trunk-jpl/src/m/classes/transient.m (revision 26059) @@ -7,7 +7,7 @@ properties (SetAccess=public) issmb = 0; ismasstransport = 0; - isoceantransport = 0; + isoceantransport = 0; isstressbalance = 0; isthermal = 0; isgroundingline = 0; @@ -15,7 +15,7 @@ isdamageevolution = 0; ismovingfront = 0; ishydrology = 0; - issampling = 0; + issampling = 0; isslc = 0; amr_frequency = 0; isoceancoupling = 0; @@ -43,7 +43,7 @@ self.isdamageevolution = 0; self.ismovingfront =0; self.ishydrology = 0; - self.issampling = 0; + self.issampling = 0; self.isslc = 0; self.isoceancoupling = 0; self.amr_frequency = 0; @@ -64,7 +64,7 @@ self.isdamageevolution = 0; self.ismovingfront = 0; self.ishydrology = 0; - self.issampling = 0; + self.issampling = 0; self.isslc = 0; self.isoceancoupling = 0; self.amr_frequency = 0; @@ -97,7 +97,7 @@ md = checkfield(md,'fieldname','transient.requested_outputs','stringrow',1); md = checkfield(md,'fieldname','transient.isslc','numel',[1],'values',[0 1]); md = checkfield(md,'fieldname','transient.isoceancoupling','numel',[1],'values',[0 1]); - md = checkfield(md,'fieldname','transient.issampling','numel',[1],'values',[0 1]); + md = checkfield(md,'fieldname','transient.issampling','numel',[1],'values',[0 1]); md = checkfield(md,'fieldname','transient.amr_frequency','numel',[1],'>=',0,'NaN',1,'Inf',1); if (~strcmp(solution,'TransientSolution') & md.transient.iscoupling==1), @@ -120,7 +120,7 @@ fielddisplay(self,'isdamageevolution','indicates whether damage evolution is used in the transient'); fielddisplay(self,'ismovingfront','indicates whether a moving front capability is used in the transient'); fielddisplay(self,'ishydrology','indicates whether an hydrology model is used'); - fielddisplay(self,'issampling','indicates whether sampling is used in the transient') + fielddisplay(self,'issampling','indicates whether sampling is used in the transient') fielddisplay(self,'isslc','indicates whether a sea-level change solution is used in the transient'); fielddisplay(self,'isoceancoupling','indicates whether a coupling with an ocean model is used in the transient'); fielddisplay(self,'amr_frequency','frequency at which mesh is refined in simulations with multiple time_steps'); @@ -138,7 +138,7 @@ WriteData(fid,prefix,'object',self,'fieldname','isdamageevolution','format','Boolean'); WriteData(fid,prefix,'object',self,'fieldname','ishydrology','format','Boolean'); WriteData(fid,prefix,'object',self,'fieldname','ismovingfront','format','Boolean'); - WriteData(fid,prefix,'object',self,'fieldname','issampling','format','Boolean'); + WriteData(fid,prefix,'object',self,'fieldname','issampling','format','Boolean'); WriteData(fid,prefix,'object',self,'fieldname','isslc','format','Boolean'); WriteData(fid,prefix,'object',self,'fieldname','isoceancoupling','format','Boolean'); WriteData(fid,prefix,'object',self,'fieldname','amr_frequency','format','Integer'); @@ -164,7 +164,7 @@ writejsdouble(fid,[modelname '.trans.isdamageevolution'],self.isdamageevolution); writejsdouble(fid,[modelname '.trans.ismovingfront'],self.ismovingfront); writejsdouble(fid,[modelname '.trans.ishydrology'],self.ishydrology); - writejsdouble(fid,[modelname '.trans.issampling'],self.issampling); + writejsdouble(fid,[modelname '.trans.issampling'],self.issampling); writejsdouble(fid,[modelname '.trans.isslc'],self.isslc); writejsdouble(fid,[modelname '.trans.isoceancoupling'],self.isoceancoupling); writejsdouble(fid,[modelname '.trans.amr_frequency'],self.amr_frequency); Index: ../trunk-jpl/src/m/consistency/ismodelselfconsistent.m =================================================================== --- ../trunk-jpl/src/m/consistency/ismodelselfconsistent.m (revision 26058) +++ ../trunk-jpl/src/m/consistency/ismodelselfconsistent.m (revision 26059) @@ -2,7 +2,7 @@ %ISMODELSELFCONSISTENT - check that model forms a closed form solvable problem. % % Usage: -% ismodelselfconsistent(md), +% ismodelselfconsistent(md); %initialize consistency as true md.private.isconsistent=true; Index: ../trunk-jpl/src/m/plot/radarpower.py =================================================================== --- ../trunk-jpl/src/m/plot/radarpower.py (nonexistent) +++ ../trunk-jpl/src/m/plot/radarpower.py (revision 26059) @@ -0,0 +1,16 @@ +def solveslm(slm, solutionstringi, *args): + """RADARPOWER - overlay a power radar image on an existing mesh + + This routine will overlay a power radar image on an existing mesh. + The power amplitude will be output to vel for now. + In the future, think about a field to hold this value. + + Usage: + md=radarpower(md,options); + md=radarpower(md) + + TODO: + - Translate from MATLAB API as we bring Python plotting capabilities online + """ + + return md Index: ../trunk-jpl/src/m/solve/solveslm.m =================================================================== --- ../trunk-jpl/src/m/solve/solveslm.m (revision 26058) +++ ../trunk-jpl/src/m/solve/solveslm.m (revision 26059) @@ -1,5 +1,5 @@ function slm=solveslm(slm,solutionstringi,varargin) -%SOLVESLR - apply solution sequence for this sealevel model +%SOLVESLM - apply solution sequence for this sealevel model % % Usage: % slm=solve(slm,solutionstring,varargin) @@ -42,7 +42,7 @@ slm.private.solution=solutionstring; cluster=slm.cluster; batch=0; -%now, go through icecaps, glacies and earth, and upload all the data independently: +%now, go through icecaps, glaciers and earth, and upload all the data independently: disp('solving ice caps first'); for i=1:length(slm.icecaps), slm.icecaps{i}=solve(slm.icecaps{i},solutionstringi,'batch','yes'); @@ -50,7 +50,7 @@ disp('solving earth now'); slm.earth=solve(slm.earth,solutionstringi,'batch','yes'); -%Firs, build a runtime name that is unique +%First, build a runtime name that is unique c=clock; slm.private.runtimename=sprintf('%s-%02i-%02i-%04i-%02i-%02i-%02i-%i',slm.miscellaneous.name,c(2),c(3),c(1),c(4),c(5),floor(c(6)),feature('GetPid')); Index: ../trunk-jpl/src/m/solve/solveslm.py =================================================================== --- ../trunk-jpl/src/m/solve/solveslm.py (nonexistent) +++ ../trunk-jpl/src/m/solve/solveslm.py (revision 26059) @@ -0,0 +1,95 @@ +from datetime import datetime +import os + +import numpy as np + +from loadresultsfromcluster import loadresultsfromcluster +from pairoptions import pairoptions +from waitonlock import waitonlock + + +def solveslm(slm, solutionstringi, *args): + """SOLVESLM - apply solution sequence for this sealevel model + + Usage: + slm=solve(slm,solutionstring,varargin) + where varargin is a lit of paired arguments of string OR enums + + solution types available comprise: + - 'Transient' + + extra options: + + Examples: + slm=solve(slm,'Transient'); + """ + + # Recover and process solve options + if solutionstringi.lower() == 'tr' or solutionstringi.lower() == 'transient': + solutionstring = 'TransientSolution' + else: + raise RuntimeError('solutionstring {} not supported!'.format(solutionstringi)) + + # Default settings for debugging + valgrind = 0 + #slm.cluster.interactive = 0 + #valgrind = 1 + + # Check consistency + slm.checkconsistency(solutionstring) + + # Process options + options = pairoptions('solutionstring', solutionstring, *args) + + # Make sure we request sum of cluster processors + totalnp = 0 + for i in range(len(slm.icecaps)): + totalnp = totalnp + slm.icecaps[i].cluster.np + totalnp = totalnp + slm.earth.cluster.np + if totalnp != slm.cluster.np: + raise RuntimeError('sum of all icecaps and earch cluster processors requestes should be equal to slm.cluster.np') + + # Recover some fields + slm.private.solution = solutionstring + cluster = slm.cluster + batch = 0 + # Now, go through icecaps, glaciers and earth, and upload all the data independently + print('solving ice caps first') + for i in range(len(slm.icecaps)): + slm.icecaps[i] = solve(slm.icecaps[i], solutionastringi,'batch','yes') + print('solving earth now') + slm.earth = solve(slm.earth, solutionstringi, 'batch', 'yes') + + # First, build a runtime name that is unique + c = datetime.now() + md.private.runtimename = "%s-%02i-%02i-%04i-%02i-%02i-%02i-%i" % (md.miscellaneous.name, c.month, c.day, c.year, c.hour, c.minute, c.second, os.getpid()) + + # Write all input files + privateruntimenames = [] + miscellaneousnames = [] + nps = [] + for i in range(len(slm.icecaps)): + privateruntimenames.append(slm.icecaps[i],private.runtimename) + miscellaneousnames.append(slm.earth.miscellaneous.name) + nps.append(slm.earth.cluster.np) + + BuildQueueScriptMultipleModels(cluster, slm.private.runtimename, slm.miscellaneous.name, slm.private.solution, valgrind, privateruntimenames, miscellaneousnames, nps) + + # Upload all required files, given that each individual solution for icecaps and earth model already did + filelist = [slm.miscellaneous.name + '.queue'] + UploadQueueJob(cluster, slm.miscellaneous.name, slm.private.runtimename, filelist) + + # Launch queue job + LaunchQueueJob(cluster, slm.miscellaneous.name, slm.private.runtimename, filelist, '', batch) + + # Wait on lock + if slm.settings.waitonlock > 0: + islock = waitonlock(slm) + if islock == 0: # no results to be loaded + print('The results must be loaded manually with md = loadresultsfromcluster(md).') + else: # load results + if slm.verbose.solution: + print('loading results from cluster') + slm = loadresultsfromcluster(slm) + + return slm Index: ../trunk-jpl/src/m/classes/dsl.m =================================================================== --- ../trunk-jpl/src/m/classes/dsl.m (revision 26058) +++ ../trunk-jpl/src/m/classes/dsl.m (revision 26059) @@ -6,9 +6,9 @@ classdef dsl properties (SetAccess=public) - global_average_thermosteric_sea_level; %corresponds to zostoga field in CMIP5 archives. Specified as a temporally variable quantity (in m) - sea_surface_height_above_geoid; %corresponds to zos field in CMIP5 archives. Spatial average is 0. Specified as a spatio-temporally variable (in m) - sea_water_pressure_at_sea_floor; %corresponds to bpo field in CMIP5 archives. Specified as a spatio-temporally variable quantity (in m equivalent, not in Pa!) + global_average_thermosteric_sea_level; %Corresponds to zostoga field in CMIP5 archives. Specified as a temporally variable quantity (in m). + sea_surface_height_above_geoid; %Corresponds to zos field in CMIP5 archives. Spatial average is 0. Specified as a spatio-temporally variable quantity (in m). + sea_water_pressure_at_sea_floor; %Corresponds to bpo field in CMIP5 archives. Specified as a spatio-temporally variable quantity (in m equivalent, not in Pa!). end methods @@ -50,9 +50,9 @@ function disp(self) % {{{ disp(sprintf(' dsl parameters:')); - fielddisplay(self,'global_average_thermosteric_sea_level','corresponds to zostoga field in CMIP5 archives. Specified as a temporally variable quantity (in m)'); - fielddisplay(self,'sea_surface_height_above_geoid','corresponds to zos field in CMIP5 archives. Spatial average is 0. Specified as a spatio-temporally quantity (in m)'); - fielddisplay(self,'sea_water_pressure_at_sea_floor','corresponds to bpo field in CMIP5 archives. Specified as a spatio-temporally variable quantity (in m equivalent)'); + fielddisplay(self,'global_average_thermosteric_sea_level','Corresponds to zostoga field in CMIP5 archives. Specified as a temporally variable quantity (in m).'); + fielddisplay(self,'sea_surface_height_above_geoid','Corresponds to zos field in CMIP5 archives. Spatial average is 0. Specified as a spatio-temporally variable quantity (in m).'); + fielddisplay(self,'sea_water_pressure_at_sea_floor','Corresponds to bpo field in CMIP5 archives. Specified as a spatio-temporally variable quantity (in m equivalent, not in Pa!).'); end % }}} function marshall(self,prefix,md,fid) % {{{ Index: ../trunk-jpl/src/m/classes/dsl.py =================================================================== --- ../trunk-jpl/src/m/classes/dsl.py (revision 26058) +++ ../trunk-jpl/src/m/classes/dsl.py (revision 26059) @@ -14,24 +14,22 @@ """ def __init__(self): #{{{ - self.global_average_thermosteric_sea_level_change = 0 #corresponds to zostoga field in CMIP5 archives. Specified as a temporally variable global rate (mm/yr) - self.sea_surface_height_change_above_geoid = float('NaN') #corresponds to zos field in CMIP5 archives. Spatial average is 0. Specified as a spatio-temporally variable rate (mm/yr) - self.sea_water_pressure_change_at_sea_floor = float('NaN') #corresponds to bpo field in CMIP5 archives. Specified as a spatio-temporally variable rate (in Pa/yr) - self.compute_fingerprints = 0; #do we use the sea water pressure change to compute fingerprints and correct sea_surface_height_change_above_geoid + self.global_average_thermosteric_sea_level = np.nan # Corresponds to zostoga field in CMIP5 archives. Specified as a temporally variable quantity (in m). + self.sea_surface_height_above_geoid = np.nan # Corresponds to zos field in CMIP5 archives. Spatial average is 0. Specified as a spatio-temporally variable quantity (in m). + self.sea_water_pressure_at_sea_floor = np.nan # Corresponds to bpo field in CMIP5 archives. Specified as a spatio-temporally variable quantity (in m equivalent, not in Pa!). #}}} def __repr__(self): #{{{ s = ' dsl parameters:\n' - s += '{}\n'.format(fielddisplay(self, 'global_average_thermosteric_sea_level_change', 'corresponds to zostoga field in CMIP5 archives. Specified as a temporally variable global rate (mm/yr)')) - s += '{}\n'.format(fielddisplay(self, 'sea_surface_height_change_above_geoid', 'corresponds to zos field in CMIP5 archives. Spatial average is 0. Specified as a spatio-temporally variable rate (mm/yr)')) - s += '{}\n'.format(fielddisplay(self, 'sea_water_pressure_change_at_sea_floor', 'corresponds to bpo field in CMIP5 archives. Specified as a spatio-temporally variable rate (in Pa/yr)')) - s += '{}\n'.format(fielddisplay(self, 'compute_fingerprints', 'do we use the sea water pressure change to compute fingerprints and correct sea_surface_height_change_above_geoid')) + s += '{}\n'.format(fielddisplay(self, 'global_average_thermosteric_sea_level', 'Corresponds to zostoga field in CMIP5 archives. Specified as a temporally variable quantity (in m).')) + s += '{}\n'.format(fielddisplay(self, 'sea_surface_height_above_geoid', 'Corresponds to zos field in CMIP5 archives. Spatial average is 0. Specified as a spatio-temporally variable quantity (in m).')) + s += '{}\n'.format(fielddisplay(self, 'sea_water_pressure_at_sea_floor', 'Corresponds to bpo field in CMIP5 archives. Specified as a spatio-temporally variable quantity (in m equivalent, not in Pa!).')) return s #}}} def extrude(self, md): #{{{ - self.sea_surface_height_change_above_geoid = project3d(md, 'vector', self.sea_surface_height_change_above_geoid, 'type', 'node') - self.sea_water_pressure_change_at_sea_floor = project3d(md, 'vector', self.sea_water_pressure_change_at_sea_floor, 'type', 'node') + self.sea_surface_height_above_geoid = project3d(md, 'vector', self.sea_surface_height_above_geoid, 'type', 'node') + self.sea_water_pressure_at_sea_floor = project3d(md, 'vector', self.sea_water_pressure_at_sea_floor, 'type', 'node') return self #}}} @@ -41,18 +39,16 @@ def checkconsistency(self, md, solution, analyses): #{{{ # Early return - if 'SealevelriseAnalysis' not in analyses: + if ('SealevelriseAnalysis' not in analyses) or (solution == 'TransientSolution' and not md.transient.isslc): return md - if solution == 'TransientSolution' and md.transient.isslc == 0: - return md - md = checkfield(md, 'fieldname', 'dsl.global_average_thermosteric_sea_level_change', 'NaN', 1, 'Inf', 1) - md = checkfield(md, 'fieldname', 'dsl.sea_surface_height_change_above_geoid', 'NaN', 1, 'Inf', 1, 'timeseries', 1) - md = checkfield(md, 'fieldname', 'dsl.sea_water_pressure_change_at_sea_floor', 'NaN', 1, 'Inf', 1, 'timeseries', 1) - md = checkfield(md, 'fieldname', 'dsl.compute_fingerprints', 'NaN', 1, 'Inf', 1, 'values', [0, 1]) - if self.compute_fingerprints: - # Check if geodetic flag of slr is on - if not md.slr.geodetic: - raise RuntimeError('DSL checkconsistency error message: if bottom pressure fingerprints computations are requested, slr class should have geodetic flag on') + + md = checkfield(md, 'fieldname', 'dsl.global_average_thermosteric_sea_level', 'NaN', 1, 'Inf', 1) + md = checkfield(md, 'fieldname', 'dsl.sea_surface_height_above_geoid', 'NaN', 1, 'Inf', 1, 'timeseries', 1) + md = checkfield(md, 'fieldname', 'dsl.sea_water_pressure_at_sea_floor', 'NaN', 1, 'Inf', 1, 'timeseries', 1) + + if md.solidearth.settings.compute_bp_grd: + md = checkfield(md, 'fieldname', dsl.sea_water_pressure_at_sea_floor, 'empty', 1) + return md # }}} @@ -59,8 +55,7 @@ def marshall(self, prefix, md, fid): #{{{ yts = md.constants.yts WriteData(fid, prefix, 'name', 'md.dsl.model', 'data', 1, 'format', 'Integer') - WriteData(fid, prefix, 'object', self, 'class', 'dsl', 'fieldname', 'compute_fingerprints', 'format', 'Integer') - WriteData(fid, prefix, 'object', self, 'class', 'dsl', 'fieldname', 'global_average_thermosteric_sea_level_change', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', 1+1, 'yts', md.constants.yts, 'scale', 1e-3/md.constants.yts) - WriteData(fid, prefix, 'object', self, 'class', 'dsl', 'fieldname', 'sea_water_pressure_change_at_sea_floor', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices+1, 'yts', md.constants.yts, 'scale', 1e-3/md.constants.yts) - WriteData(fid, prefix, 'object', self, 'class', 'dsl', 'fieldname', 'sea_surface_height_change_above_geoid', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices+1, 'yts', md.constants.yts) + WriteData(fid, prefix, 'object', self, 'class', 'dsl', 'fieldname', 'global_average_thermosteric_sea_level', 'format', 'DoubleMat', 'mattype', 2, 'timeserieslength', 1+1, 'yts', md.constants.yts, 'scale', 1e-3/md.constants.yts) # mattype 2, because we are sending a GMSL value identical everywhere on each element. + WriteData(fid, prefix, 'object', self, 'class', 'dsl', 'fieldname', 'sea_water_pressure_at_sea_floor', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices+1, 'yts', md.constants.yts, 'scale', 1e-3/md.constants.yts) # mattype 1 because we specify DSL at vertex locations. + WriteData(fid, prefix, 'object', self, 'class', 'dsl', 'fieldname', 'sea_surface_height_above_geoid', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices+1, 'yts', md.constants.yts) # mattype 1 because we specify bottom pressure at vertex locations. # }}} Index: ../trunk-jpl/src/m/classes/dslmme.m =================================================================== --- ../trunk-jpl/src/m/classes/dslmme.m (revision 26058) +++ ../trunk-jpl/src/m/classes/dslmme.m (revision 26059) @@ -7,9 +7,9 @@ properties (SetAccess=public) modelid; %index into the multi-model ensemble, determine which field will be used. - global_average_thermosteric_sea_level; %corresponds to zostoga fields in CMIP5 archives. Specified as a temporally quantity (in m) for each ensemble. - sea_surface_height_above_geoid; %corresponds to zos fields in CMIP5 archives. Spatial average is 0. Specified as a spatio-temporally quantity (in m) for each ensemble - sea_water_pressure_at_sea_floor; %corresponds to bpo field in CMIP5 archives. Specified as a spatio-temporally variable quantity (in m equivalent, not in Pa!) for each ensemble + global_average_thermosteric_sea_level; %Corresponds to zostoga field in CMIP5 archives. Specified as a temporally variable quantity (in m) for each ensemble. + sea_surface_height_above_geoid; %Corresponds to zos field in CMIP5 archives. Spatial average is 0. Specified as a spatio-temporally variable quantity (in m) for each ensemble. + sea_water_pressure_at_sea_floor; %Corresponds to bpo field in CMIP5 archives. Specified as a spatio-temporally variable quantity (in m equivalent, not in Pa!) for each ensemble. end methods @@ -52,9 +52,9 @@ disp(sprintf(' dsl mme parameters:')); fielddisplay(self,'modelid','index into the multi-model ensemble, determine which field will be used.'); - fielddisplay(self,'global_average_thermosteric_sea_level','corresponds to zostoga fields in CMIP5 archives. Specified as a temporally variable quantity (in m) for each ensemble.'); - fielddisplay(self,'sea_surface_height_above_geoid','corresponds to zos fields in CMIP5 archives. Spatial average is 0. Specified as a spatio-temporally quantity (in m) for each ensemble.'); - fielddisplay(self,'sea_water_pressure_at_sea_floor','corresponds to bpo fields in CMIP5 archives. Specified as a spatio-temporally quantity (in m) for each ensemble.'); + fielddisplay(self,'global_average_thermosteric_sea_level','Corresponds to zostoga field in CMIP5 archives. Specified as a temporally variable quantity (in m) for each ensemble.'); + fielddisplay(self,'sea_surface_height_above_geoid','Corresponds to zos field in CMIP5 archives. Spatial average is 0. Specified as a spatio-temporally variable quantity (in m) for each ensemble.'); + fielddisplay(self,'sea_water_pressure_at_sea_floor','Corresponds to bpo field in CMIP5 archives. Specified as a spatio-temporally variable quantity (in m equivalent, not in Pa!) for each ensemble.'); end % }}} function marshall(self,prefix,md,fid) % {{{ Index: ../trunk-jpl/src/m/classes/dslmme.py =================================================================== --- ../trunk-jpl/src/m/classes/dslmme.py (revision 26058) +++ ../trunk-jpl/src/m/classes/dslmme.py (revision 26059) @@ -14,15 +14,14 @@ """ def __init__(self, *args): #{{{ - self.modelid = 0 #index into the multi-model ensemble - self.global_average_thermosteric_sea_level_change = [] #corresponds to zostoga fields in CMIP5 archives. Specified as a temporally variable global rate (mm/yr) for each ensemble. - self.sea_surface_height_change_above_geoid = [] #corresponds to zos fields in CMIP5 archives. Spatial average is 0. Specified as a spatio-temporally variable rate (mm/yr) for each ensemble - self.sea_water_pressure_change_at_sea_floor = [] #corresponds to bpo fields in CMIP5 archives. Specified as a spatio-temporally variable rate (in mm/yr equivalent, not in Pa/yr!) for each ensemble - self.compute_fingerprints = 0 #corresponds to zos fields in CMIP5 archives. Spatial average is 0. Specified as a spatio-temporally variable rate (mm/yr) for each ensemble - - nargin = len(args) + self.modelid = 0 # Index into the multi-model ensemble + self.global_average_thermosteric_sea_level = [] # Corresponds to zostoga field in CMIP5 archives. Specified as a temporally variable quantity (in m) for each ensemble. + self.sea_surface_height_above_geoid = [] # Corresponds to zos field in CMIP5 archives. Spatial average is 0. Specified as a spatio-temporally variable quantity (in m) for each ensemble. + self.sea_water_pressure_at_sea_floor = [] #Corresponds to bpo field in CMIP5 archives. Specified as a spatio-temporally variable quantity (in m equivalent, not in Pa!) for each ensemble. - if nargin == 0: + nargs = len(args) + + if nargs == 0: self.setdefaultparameters() else: raise Exception('constructor not supported') @@ -31,10 +30,9 @@ def __repr__(self): # {{{ s = ' dsl mme parameters:\n' s += '{}\n'.format(fielddisplay(self, 'modelid', 'index into the multi-model ensemble, determines which field will be used.')) - s += '{}\n'.format(fielddisplay(self, 'global_average_thermosteric_sea_level_change', 'corresponds to zostoga fields in CMIP5 archives. Specified as a temporally variable global rate (mm/yr) for each ensemble.')) - s += '{}\n'.format(fielddisplay(self, 'sea_surface_height_change_above_geoid', 'corresponds to zos fields in CMIP5 archives. Spatial average is 0. Specified as a spatio-temporally variable rate (mm/yr) for each ensemble.')) - s += '{}\n'.format(fielddisplay(self, 'sea_water_pressure_change_at_sea_floor', 'corresponds to bpo fields in CMIP5 archives. Specified as a spatio-temporally variable rate (in mm/yr) for each ensemble.')) - s += '{}\n'.format(fielddisplay(self, 'compute_fingerprints', 'do we use the sea water pressure change to compute fingerprints and correct sea_surface_height_change_above_geoid')) + s += '{}\n'.format(fielddisplay(self, 'global_average_thermosteric_sea_level', 'Corresponds to zostoga field in CMIP5 archives. Specified as a temporally variable quantity (in m) for each ensemble.')) + s += '{}\n'.format(fielddisplay(self, 'sea_surface_height_above_geoid', 'Corresponds to zos field in CMIP5 archives. Spatial average is 0. Specified as a spatio-temporally variable quantity (in m) for each ensemble.')) + s += '{}\n'.format(fielddisplay(self, 'sea_water_pressure_at_sea_floor', 'Corresponds to bpo field in CMIP5 archives. Specified as a spatio-temporally variable quantity (in m equivalent, not in Pa!) for each ensemble.')) return s #}}} @@ -43,34 +41,35 @@ #}}} def checkconsistency(self, md, solution, analyses): # {{{ - if ('SealevelriseAnalysis' not in analyses) or (solution == 'TransientSolution' and not md.transient.isslr): + # Early return + if ('SealevelriseAnalysis' not in analyses) or (solution == 'TransientSolution' and not md.transient.isslc): return md - for i in range(len(self.global_average_thermosteric_sea_level_change)): - md = checkfield(md, 'field', self.global_average_thermosteric_sea_level_change[i], 'NaN', 1, 'Inf', 1) - md = checkfield(md, 'field', self.sea_surface_height_change_above_geoid[i], 'NaN', 1, 'Inf', 1, 'timeseries', 1) - md = checkfield(md, 'field', self.sea_water_pressure_change_at_sea_floor[i], 'NaN', 1, 'Inf', 1, 'timeseries', 1) - md = checkfield(md, 'field', self.modelid, 'NaN', 1, 'Inf', 1, '>=', 1, '<=',len(self.global_average_thermosteric_sea_level_change)) - if self.compute_fingerprints: - #check geodetic flag of slr is on - if not md.solidearth.settings.computesealevelchange: - raise Exception('DSL checkconsistency error message: if bottom pressure fingerprints computations are requested, slr class should have geodetic flag on') + + for i in range(len(self.global_average_thermosteric_sea_level)): + md = checkfield(md, 'field', self.global_average_thermosteric_sea_level[i], 'NaN', 1, 'Inf', 1) + md = checkfield(md, 'field', self.sea_surface_height_above_geoid[i], 'NaN', 1, 'Inf', 1, 'timeseries', 1) + md = checkfield(md, 'field', self.sea_water_pressure_at_sea_floor[i], 'NaN', 1, 'Inf', 1, 'timeseries', 1) + md = checkfield(md, 'field', self.modelid, 'NaN', 1, 'Inf', 1, '>=', 1, '<=',len(self.global_average_thermosteric_sea_level)) + + if self.solidearth.settings.compute_bp_grd: + md = checkfield(md, 'field', self.sea_water_pressure_at_sea_floor, 'empty', 1) + return md #}}} def marshall(self, prefix, md, fid): #{{{ WriteData(fid, prefix, 'name', 'md.dsl.model', 'data', 2, 'format', 'Integer') - WriteData(fid, prefix, 'object', self, 'fieldname', 'compute_fingerprints', 'format', 'Integer') WriteData(fid, prefix, 'object', self, 'fieldname', 'modelid', 'format', 'Double') - WriteData(fid, prefix, 'name', 'md.dsl.nummodels', 'data', len(self.global_average_thermosteric_sea_level_change), 'format', 'Integer') - WriteData(fid, prefix, 'object', self, 'fieldname', 'global_average_thermosteric_sea_level_change', 'format', 'MatArray', 'timeseries', 1, 'timeserieslength', 2, 'yts', md.constants.yts, 'scale', 1e-3 / md.constants.yts) - WriteData(fid, prefix, 'object', self, 'fieldname', 'sea_water_pressure_change_at_sea_floor', 'format', 'MatArray', 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts, 'scale', 1e-3) - WriteData(fid, prefix, 'object', self, 'fieldname', 'sea_surface_height_change_above_geoid', 'format', 'MatArray', 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts, 'scale', 1e-3 / md.constants.yts) + WriteData(fid, prefix, 'name', 'md.dsl.nummodels', 'data', len(self.global_average_thermosteric_sea_level), 'format', 'Integer') + WriteData(fid, prefix, 'object', self, 'fieldname', 'global_average_thermosteric_sea_level', 'format', 'MatArray', 'timeseries', 1, 'timeserieslength', 2, 'yts', md.constants.yts, 'scale', 1e-3 / md.constants.yts) + WriteData(fid, prefix, 'object', self, 'fieldname', 'sea_water_pressure_at_sea_floor', 'format', 'MatArray', 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts, 'scale', 1e-3) + WriteData(fid, prefix, 'object', self, 'fieldname', 'sea_surface_height_above_geoid', 'format', 'MatArray', 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts, 'scale', 1e-3 / md.constants.yts) #}}} def extrude(self, md): #{{{ - for i in range(len(self.global_average_thermosteric_sea_level_change)): - self.sea_surface_height_change_above_geoid[i] = project3d(md, 'vector', self.self.sea_surface_height_change_above_geoid[i], 'type', 'node', 'layer', 1) - self.sea_water_pressure_change_at_sea_floor[i] = project3d(md, 'vector', self.sea_water_pressure_change_at_sea_floor[i], 'type', 'node', 'layer', 1) + for i in range(len(self.global_average_thermosteric_sea_level)): + self.sea_surface_height_above_geoid[i] = project3d(md, 'vector', self.self.sea_surface_height_above_geoid[i], 'type', 'node', 'layer', 1) + self.sea_water_pressure_at_sea_floor[i] = project3d(md, 'vector', self.sea_water_pressure_at_sea_floor[i], 'type', 'node', 'layer', 1) return self #}}} Index: ../trunk-jpl/src/m/classes/fourierlove.m =================================================================== --- ../trunk-jpl/src/m/classes/fourierlove.m (revision 26058) +++ ../trunk-jpl/src/m/classes/fourierlove.m (revision 26059) @@ -79,12 +79,8 @@ end %need 'litho' material: - if ~isa(md.materials,'materials') + if ~isa(md.materials,'materials') | ~sum(strcmpi(md.materials.nature,'litho')) error('Need a ''litho'' material to run a Fourier Love number analysis'); - else - if ~sum(strcmpi(md.materials.nature,'litho')), - error('Need a ''litho'' material to run a Fourier Love number analysis'); - end end end % }}} Index: ../trunk-jpl/src/m/classes/fourierlove.py =================================================================== --- ../trunk-jpl/src/m/classes/fourierlove.py (revision 26058) +++ ../trunk-jpl/src/m/classes/fourierlove.py (revision 26059) @@ -4,11 +4,10 @@ class fourierlove(object): - """ - Fourier Love Number class definition + """FOURIERLOVE - Fourier Love Number class definition - Usage: - fourierlove = fourierlove() + Usage: + fourierlove = fourierlove() """ def __init__(self): # {{{ self.nfreq = 0 @@ -22,11 +21,14 @@ self.love_kernels = 0 self.forcing_type = 0 - #set defaults self.setdefaultparameters() #}}} def __repr__(self): # {{{ + # TODO: + # - Convert all formatting to calls to .format (see any + # already converted .__repr__ method for examples) + # string = ' Fourier Love class:' string = "%s\n%s" % (string, fielddisplay(self, 'nfreq', 'number of frequencies sampled (default 1, elastic) [Hz]')) string = "%s\n%s" % (string, fielddisplay(self, 'frequencies', 'frequencies sampled (convention defaults to 0 for the elastic case) [Hz]')) @@ -75,6 +77,8 @@ #}}} def checkconsistency(self, md, solution, analyses): # {{{ + if 'LoveAnalysis' not in analyses: + return md md = checkfield(md, 'fieldname', 'love.nfreq', 'NaN', 1, 'Inf', 1, 'numel', [1], '>', 0) md = checkfield(md, 'fieldname', 'love.frequencies', 'NaN', 1, 'Inf', 1, 'numel', [md.love.nfreq]) md = checkfield(md, 'fieldname', 'love.sh_nmax', 'NaN', 1, 'Inf', 1, 'numel', [1], '>', 0) @@ -88,6 +92,9 @@ if md.love.sh_nmin <= 1 and md.love.forcing_type == 9: raise RuntimeError("Degree 1 not supported for Volumetric Potential forcing. Use sh_min >= 2 for this kind of calculation.") + # need 'litho' material + if not hasattr(md.materials, 'materials') or 'litho' not in md.materials.nature: + raise RuntimeError('Need a \'litho\' material to run a Fourier Love number analysis') return md # }}} Index: ../trunk-jpl/src/m/classes/geometry.m =================================================================== --- ../trunk-jpl/src/m/classes/geometry.m (revision 26058) +++ ../trunk-jpl/src/m/classes/geometry.m (revision 26059) @@ -84,10 +84,10 @@ end % }}} function marshall(self,prefix,md,fid) % {{{ - - if (size(self.thickness==md.mesh.numberofvertices) | (self.thickness==md.mesh.numberofvertices+1)), + + if (size(self.thickness)==md.mesh.numberofvertices) | (size(self.thickness)==md.mesh.numberofvertices+1)), WriteData(fid,prefix,'object',self,'fieldname','thickness','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); - elseif (size(self.thickness==md.mesh.numberofelements) | (self.thickness==md.mesh.numberofelements+1)), + elseif (size(self.thickness)==md.mesh.numberofelements) | (size(self.thickness)==md.mesh.numberofelements+1)), WriteData(fid,prefix,'object',self,'fieldname','thickness','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts); else error('geometry thickness time series should be a vertex or element time series'); Index: ../trunk-jpl/src/m/classes/geometry.py =================================================================== --- ../trunk-jpl/src/m/classes/geometry.py (revision 26058) +++ ../trunk-jpl/src/m/classes/geometry.py (revision 26059) @@ -50,12 +50,8 @@ #}}} def checkconsistency(self, md, solution, analyses): #{{{ - if (solution == 'TransientSolution' and md.transient.isgia) or (solution == 'GiaSolution'): - md = checkfield(md, 'fieldname', 'geometry.thickness', 'timeseries', 1, 'NaN', 1, 'Inf', 1) - elif solution == 'SealevelriseSolution': - md = checkfield(md, 'fieldname', 'geometry.bed', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices]) - elif solution == 'LoveSolution': - return + if solution == 'LoveSolution': + return md else: md = checkfield(md, 'fieldname', 'geometry.surface', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices]) md = checkfield(md, 'fieldname', 'geometry.base', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices]) @@ -70,13 +66,18 @@ if np.any(np.abs(self.bed[pos] - self.base[pos]) > 10**-9): md.checkmessage('equality base = bed on grounded ice violated') md = checkfield(md, 'fieldname', 'geometry.bed', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices]) - return md # }}} def marshall(self, prefix, md, fid): #{{{ + if (len(self.thickness) == md.mesh.numberofvertices) or (len(self.thickness) == md.mesh.numberofvertices + 1): + WriteData(fid, prefix, 'object', self, 'fieldname', 'thickness', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts) + elif (len(self.thickness) == md.mesh.numberofelements) or (len(self.thickness) == md.mesh.numberofelements + 1): + WriteData(fid, prefix, 'object', self, 'fieldname', 'thickness', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts) + else: + raise RuntimeError('geometry thickness time series should be a vertex or element time series') + WriteData(fid, prefix, 'object', self, 'fieldname', 'surface', 'format', 'DoubleMat', 'mattype', 1) - WriteData(fid, prefix, 'object', self, 'fieldname', 'thickness', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts) WriteData(fid, prefix, 'object', self, 'fieldname', 'base', 'format', 'DoubleMat', 'mattype', 1) WriteData(fid, prefix, 'object', self, 'fieldname', 'bed', 'format', 'DoubleMat', 'mattype', 1) WriteData(fid, prefix, 'object', self, 'fieldname', 'hydrostatic_ratio', 'format', 'DoubleMat', 'mattype', 1) Index: ../trunk-jpl/src/m/classes/hydrologyshreve.m =================================================================== --- ../trunk-jpl/src/m/classes/hydrologyshreve.m (revision 26058) +++ ../trunk-jpl/src/m/classes/hydrologyshreve.m (revision 26059) @@ -30,7 +30,7 @@ %Type of stabilization to use 0:nothing 1:artificial_diffusivity self.stabilization = 1; - self.requested_outputs = {'default'}; + self.requested_outputs = {'default'}; end % }}} function md = checkconsistency(self,md,solution,analyses) % {{{ Index: ../trunk-jpl/src/m/classes/hydrologyshreve.py =================================================================== --- ../trunk-jpl/src/m/classes/hydrologyshreve.py (revision 26058) +++ ../trunk-jpl/src/m/classes/hydrologyshreve.py (revision 26059) @@ -4,11 +4,10 @@ class hydrologyshreve(object): - """ - HYDROLOGYSHREVE class definition + """HYDROLOGYSHREVE class definition - Usage: - hydrologyshreve = hydrologyshreve() + Usage: + hydrologyshreve = hydrologyshreve() """ def __init__(self): # {{{ @@ -15,13 +14,15 @@ self.spcwatercolumn = float('NaN') self.stabilization = 0 self.requested_outputs = [] - #set defaults + self.setdefaultparameters() - #}}} def __repr__(self): # {{{ - + # TODO: + # - Convert all formatting to calls to .format (see any + # already converted .__repr__ method for examples) + # string = ' hydrologyshreve solution parameters:' string = "%s\n%s" % (string, fielddisplay(self, 'spcwatercolumn', 'water thickness constraints (NaN means no constraint) [m]')) string = "%s\n%s" % (string, fielddisplay(self, 'stabilization', 'artificial diffusivity (default is 1). can be more than 1 to increase diffusivity.')) @@ -47,7 +48,7 @@ def checkconsistency(self, md, solution, analyses): # {{{ #Early return - if 'HydrologyShreveAnalysis' not in analyses: + if 'HydrologyShreveAnalysis' not in analyses or (solution == 'TransientSolution' and not md.transient.ishydrology): return md md = checkfield(md, 'fieldname', 'hydrology.spcwatercolumn', 'Inf', 1, 'timeseries', 1) @@ -60,7 +61,7 @@ WriteData(fid, prefix, 'name', 'md.hydrology.model', 'data', 2, 'format', 'Integer') WriteData(fid, prefix, 'object', self, 'fieldname', 'spcwatercolumn', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts) WriteData(fid, prefix, 'object', self, 'fieldname', 'stabilization', 'format', 'Double') - #process requested outputs + #process requested outputs outputs = self.requested_outputs indices = [i for i, x in enumerate(outputs) if x == 'default'] if len(indices) > 0: Index: ../trunk-jpl/src/m/classes/initialization.m =================================================================== --- ../trunk-jpl/src/m/classes/initialization.m (revision 26058) +++ ../trunk-jpl/src/m/classes/initialization.m (revision 26059) @@ -21,7 +21,7 @@ channelarea = NaN; sealevel = NaN; bottompressure = NaN; - sample = NaN; + sample = NaN; end methods function self = extrude(self,md) % {{{ @@ -51,7 +51,6 @@ end end % }}} function self = setdefaultparameters(self) % {{{ - end % }}} function md = checkconsistency(self,md,solution,analyses) % {{{ if ismember('StressbalanceAnalysis',analyses) & ~(strcmp(solution,'TransientSolution') & md.transient.isstressbalance == 0), @@ -97,12 +96,9 @@ end if ismember('HydrologyShreveAnalysis',analyses), if isa(md.hydrology,'hydrologyshreve'), - if strcmp(solution,'TransientSolution') & md.transient.ishydrology, + if (strcmp(solution,'TransientSolution') & md.transient.ishydrology) | strcmp(solution,'HydrologySolution'), md = checkfield(md,'fieldname','initialization.watercolumn','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]); end - if strcmp(solution,'HydrologySolution'), - md = checkfield(md,'fieldname','initialization.watercolumn','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]); - end end end if ismember('HydrologyTwsAnalysis',analyses), @@ -110,7 +106,7 @@ md = checkfield(md,'fieldname','initialization.watercolumn','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]); end end - if ismember('SealevelchangeAnalysis',analyses) & ~(strcmp(solution,'TransientSolution') & md.transient.isslc == 0), + if ismember('SealevelchangeAnalysis',analyses), if strcmp(solution,'TransientSolution') & md.transient.isslc, md = checkfield(md,'fieldname','initialization.sealevel','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]); end @@ -133,13 +129,13 @@ md = checkfield(md,'fieldname','initialization.epl_head','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]); md = checkfield(md,'fieldname','initialization.epl_thickness','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]); end - end - end - if ismember('SamplingAnalysis',analyses) & ~(strcmp(solution,'TransientSolution') & md.transient.issampling == 0), - if ~isnan(md.initialization.sample) - md = checkfield(md,'fieldname','initialization.sample','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]); - end + end end + if ismember('SamplingAnalysis',analyses) & ~(strcmp(solution,'TransientSolution') & md.transient.issampling == 0), + if ~isnan(md.initialization.sample) + md = checkfield(md,'fieldname','initialization.sample','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]); + end + end end % }}} function disp(self) % {{{ disp(sprintf(' initial field values:')); @@ -158,7 +154,7 @@ fielddisplay(self,'watercolumn','subglacial water sheet thickness (for Shreve and GlaDS) [m]'); fielddisplay(self,'hydraulic_potential','Hydraulic potential (for GlaDS) [Pa]'); fielddisplay(self,'channelarea','subglacial water channel area (for GlaDS) [m2]'); - fielddisplay(self,'sample','Realization of a Gaussian random field'); + fielddisplay(self,'sample','Realization of a Gaussian random field'); end % }}} function marshall(self,prefix,md,fid) % {{{ @@ -178,8 +174,8 @@ WriteData(fid,prefix,'object',self,'fieldname','watercolumn','format','DoubleMat','mattype',1); WriteData(fid,prefix,'object',self,'fieldname','channelarea','format','DoubleMat','mattype',1); WriteData(fid,prefix,'object',self,'fieldname','hydraulic_potential','format','DoubleMat','mattype',1); - WriteData(fid,prefix,'object',self,'fieldname','sample','format','DoubleMat','mattype',1); - + WriteData(fid,prefix,'object',self,'fieldname','sample','format','DoubleMat','mattype',1); + if md.thermal.isenthalpy, if numel(self.enthalpy) <= 1, %reconstruct enthalpy @@ -194,7 +190,7 @@ end end % }}} function savemodeljs(self,fid,modelname) % {{{ - + writejs1Darray(fid,[modelname '.initialization.vx'],self.vx); writejs1Darray(fid,[modelname '.initialization.vy'],self.vy); writejs1Darray(fid,[modelname '.initialization.vz'],self.vz); @@ -209,8 +205,8 @@ writejs1Darray(fid,[modelname '.initialization.watercolumn'],self.watercolumn); writejs1Darray(fid,[modelname '.initialization.hydraulic_potential'],self.hydraulic_potential); writejs1Darray(fid,[modelname '.initialization.channel'],self.channelarea); - writejs1Darray(fid,[modelname '.initialization.sample'],self.sample); - + writejs1Darray(fid,[modelname '.initialization.sample'],self.sample); + end % }}} end end Index: ../trunk-jpl/src/m/classes/initialization.py =================================================================== --- ../trunk-jpl/src/m/classes/initialization.py (revision 26058) +++ ../trunk-jpl/src/m/classes/initialization.py (revision 26059) @@ -28,6 +28,8 @@ self.epl_thickness = np.nan self.hydraulic_potential = np.nan self.channelarea = np.nan + self.sealevel = np.nan + self.bottompressure = np.nan self.sample = np.nan #set defaults @@ -66,6 +68,8 @@ self.sediment_head = project3d(md, 'vector', self.sediment_head, 'type', 'node', 'layer', 1) self.epl_head = project3d(md, 'vector', self.epl_head, 'type', 'node', 'layer', 1) self.epl_thickness = project3d(md, 'vector', self.epl_thickness, 'type', 'node', 'layer', 1) + self.sealevel = project3d(md, 'vector', self.sealevel, 'type', 'node', 'layer', 1) + self.bottompressure = project3d(md, 'vector', self.bottompressure, 'type', 'node', 'layer', 1) #Lithostatic pressure by default if np.ndim(md.geometry.surface) == 2: @@ -89,6 +93,9 @@ if 'MasstransportAnalysis' in analyses and not solution == 'TransientSolution' and not md.transient.ismasstransport: md = checkfield(md, 'fieldname', 'initialization.vx', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices]) md = checkfield(md, 'fieldname', 'initialization.vy', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices]) + if 'OceantransportAnalysis' in analyses: + if solution == 'TransientSolution' and md.transient.isslc and md.transient.isoceantransport: + md = checkfield(md, 'fieldname', 'initialization.bottompressure', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices]) if 'BalancethicknessAnalysis' in analyses and solution == 'BalancethicknessSolution': md = checkfield(md, 'fieldname', 'initialization.vx', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices]) md = checkfield(md, 'fieldname', 'initialization.vy', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices]) @@ -111,15 +118,14 @@ md = checkfield(md, 'fieldname', 'delta Tpmp', 'field', np.absolute(md.initialization.temperature[pos] - (md.materials.meltingpoint - md.materials.beta * md.initialization.pressure[pos])), '<', 1e-11, 'message', 'set temperature to pressure melting point at locations with waterfraction > 0') if 'HydrologyShreveAnalysis' in analyses: if hasattr(md.hydrology, 'hydrologyshreve'): + if (solution == 'TransientSolution' and md.transient.ishydrology) or solution == 'HydrologySolution': + md = checkfield(md, 'fieldname', 'initialization.watercolumn', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices]) + if 'HydrologyTwsAnalysis' in analyses: + if hasattr(md.hydrology, 'hydrologytws'): md = checkfield(md, 'fieldname', 'initialization.watercolumn', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices]) - if 'HydrologyDCInefficientAnalysis' in analyses: - if hasattr(md.hydrology, 'hydrologydc'): - md = checkfield(md, 'fieldname', 'initialization.sediment_head', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices]) - if 'HydrologyDCEfficientAnalysis' in analyses: - if hasattr(md.hydrology, 'hydrologydc'): - if md.hydrology.isefficientlayer == 1: - md = checkfield(md, 'fieldname', 'initialization.epl_head', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices]) - md = checkfield(md, 'fieldname', 'initialization.epl_thickness', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices]) + if 'SealevelchangeAnalysis' in analyses: + if solution == 'TransientSolution' and md.transient.isslc: + md = checkfield(md, 'fieldname', 'initialization.sealevel', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices]) if 'HydrologyGlaDSAnalysis' in analyses: if hasattr(md.hydrology, 'hydrologyglads'): md = checkfield(md, 'fieldname', 'initialization.watercolumn', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices]) @@ -137,6 +143,8 @@ WriteData(fid, prefix, 'object', self, 'fieldname', 'vy', 'format', 'DoubleMat', 'mattype', 1, 'scale', 1. / yts) WriteData(fid, prefix, 'object', self, 'fieldname', 'vz', 'format', 'DoubleMat', 'mattype', 1, 'scale', 1. / yts) WriteData(fid, prefix, 'object', self, 'fieldname', 'pressure', 'format', 'DoubleMat', 'mattype', 1) + WriteData(fid, prefix, 'object', self, 'fieldname', 'sealevel', 'format', 'DoubleMat', 'mattype', 1) + WriteData(fid, prefix, 'object', self, 'fieldname', 'bottompressure', 'format', 'DoubleMat', 'mattype', 1) WriteData(fid, prefix, 'object', self, 'fieldname', 'temperature', 'format', 'DoubleMat', 'mattype', 1) WriteData(fid, prefix, 'object', self, 'fieldname', 'waterfraction', 'format', 'DoubleMat', 'mattype', 1) WriteData(fid, prefix, 'object', self, 'fieldname', 'sediment_head', 'format', 'DoubleMat', 'mattype', 1) Index: ../trunk-jpl/src/m/consistency/ismodelselfconsistent.py =================================================================== --- ../trunk-jpl/src/m/consistency/ismodelselfconsistent.py (revision 26058) +++ ../trunk-jpl/src/m/consistency/ismodelselfconsistent.py (revision 26059) @@ -1,10 +1,9 @@ def ismodelselfconsistent(md): #{{{ - ''' - ISMODELSELFCONSISTENT - check that model forms a closed form solvable problem. + """ISMODELSELFCONSISTENT - check that model forms a closed form solvable problem. - Usage: - ismodelselfconsistent(md), - ''' + Usage: + ismodelselfconsistent(md) + """ #initialize consistency as true md.private.isconsistent = True @@ -52,6 +51,8 @@ analyses = ['EnthalpyAnalysis', 'ThermalAnalysis', 'MeltingAnalysis'] elif solutiontype == 'MasstransportSolution': analyses = ['MasstransportAnalysis'] + elif solutiontype == 'OceantransportSolution': + analyses = ['OceantransportAnalysis'] elif solutiontype == 'BalancethicknessSolution': analyses = ['BalancethicknessAnalysis'] elif solutiontype == 'Balancethickness2Solution': @@ -71,11 +72,11 @@ elif solutiontype == 'EsaSolution': analyses = ['EsaAnalysis'] elif solutiontype == 'TransientSolution': - analyses = ['StressbalanceAnalysis', 'StressbalanceVerticalAnalysis', 'StressbalanceSIAAnalysis', 'L2ProjectionBaseAnalysis', 'ThermalAnalysis', 'MeltingAnalysis', 'EnthalpyAnalysis', 'MasstransportAnalysis', 'HydrologyShaktiAnalysis', 'HydrologyGladsAnalysis', 'HydrologyDCInefficientAnalysis', 'HydrologyDCEfficientAnalysis', 'SealevelriseAnalysis'] + analyses = ['StressbalanceAnalysis', 'StressbalanceVerticalAnalysis', 'StressbalanceSIAAnalysis', 'L2ProjectionBaseAnalysis', 'ThermalAnalysis', 'MeltingAnalysis', 'EnthalpyAnalysis', 'MasstransportAnalysis', 'OceantransportAnalysis', 'HydrologyShaktiAnalysis', 'HydrologyGladsAnalysis', 'HydrologyShreveAnalysis', 'HydrologyTwsAnalysis', 'HydrologyDCInefficientAnalysis', 'HydrologyDCEfficientAnalysis', 'SealevelriseAnalysis'] elif solutiontype == 'SealevelriseSolution': analyses = ['SealevelriseAnalysis'] elif solutiontype == 'HydrologySolution': - analyses = ['L2ProjectionBaseAnalysis', 'HydrologyShreveAnalysis', 'HydrologyDCInefficientAnalysis', 'HydrologyDCEfficientAnalysis'] + analyses = ['L2ProjectionBaseAnalysis', 'HydrologyShreveAnalysis', 'HydrologyDCInefficientAnalysis', 'HydrologyDCEfficientAnalysis', 'HydrologyGladsAnalysis', 'HydrologyShaktiAnalysis', 'HydrologyTwsAnalysis'] elif 'DamageEvolutionSolution': analyses = ['DamageEvolutionAnalysis'] elif 'SamplingSolution': Index: ../trunk-jpl/src/m/solve/solve.m =================================================================== --- ../trunk-jpl/src/m/solve/solve.m (revision 26058) +++ ../trunk-jpl/src/m/solve/solve.m (revision 26059) @@ -9,7 +9,7 @@ % Solution types available comprise: % - 'Stressbalance' or 'sb' % - 'Masstransport' or 'mt' -% - 'Oceantransport' or 'oceant' +% - 'Oceantransport' or 'oceant' % - 'Thermal' or 'th' % - 'Steadystate' or 'ss' % - 'Transient' or 'tr' Index: ../trunk-jpl/test/NightlyRun/test2001.m =================================================================== --- ../trunk-jpl/test/NightlyRun/test2001.m (revision 26058) +++ ../trunk-jpl/test/NightlyRun/test2001.m (revision 26059) @@ -21,8 +21,8 @@ %Loading history md.timestepping.start_time=-2400000; %4,800 kyr :: EVALUATION TIME md.timestepping.time_step= 2400000; %2,400 kyr :: EVALUATION TIME -% to get rid of default final_time: make sure final_time>start_time -md.timestepping.final_time=2400000; %2,500 kyr +% to get rid of default final_time, make sure final_time > start_time +md.timestepping.final_time=2400000; %2,400 kyr md.masstransport.spcthickness=[... [md.geometry.thickness; 0],... [md.geometry.thickness; 2400000]... @@ -45,7 +45,7 @@ %md.cluster=generic('name',oshostname(),'np',3); md.verbose=verbose('11111111111'); md.verbose.solver=0; -md=solve(md,'tr'); +md=solve(md,'Transient'); %Fields and tolerances to track changes field_names ={'UGrd'}; Index: ../trunk-jpl/src/m/classes/matdamageice.m =================================================================== --- ../trunk-jpl/src/m/classes/matdamageice.m (revision 26058) +++ ../trunk-jpl/src/m/classes/matdamageice.m (revision 26059) @@ -5,25 +5,25 @@ classdef matdamageice properties (SetAccess=public) - rho_ice = 0.; - rho_water = 0.; - rho_freshwater = 0.; - mu_water = 0.; - heatcapacity = 0.; - latentheat = 0.; - thermalconductivity = 0.; - temperateiceconductivity = 0.; + rho_ice = 0.; + rho_water = 0.; + rho_freshwater = 0.; + mu_water = 0.; + heatcapacity = 0.; + latentheat = 0.; + thermalconductivity = 0.; + temperateiceconductivity = 0.; effectiveconductivity_averaging = 0.; - meltingpoint = 0.; - beta = 0.; - mixed_layer_capacity = 0.; - thermal_exchange_velocity = 0.; - rheology_B = NaN; - rheology_n = NaN; - rheology_law = ''; + meltingpoint = 0.; + beta = 0.; + mixed_layer_capacity = 0.; + thermal_exchange_velocity = 0.; + rheology_B = NaN; + rheology_n = NaN; + rheology_law = ''; %slc - earth_density = 0; + earth_density = 0; end methods Index: ../trunk-jpl/src/m/classes/matdamageice.py =================================================================== --- ../trunk-jpl/src/m/classes/matdamageice.py (revision 26058) +++ ../trunk-jpl/src/m/classes/matdamageice.py (revision 26059) @@ -30,18 +30,17 @@ 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. + #SLC + self.earth_density = 5512 # average density of the Earth, (kg / m^3) - #SLR - self.earth_density = 5512 # average density of the Earth, (kg / m^3) self.setdefaultparameters() #}}} def __repr__(self): # {{{ + # TODO: + # - Convert all formatting to calls to .format (see any + # already converted .__repr__ method for examples) + # 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]")) @@ -59,10 +58,6 @@ 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 #}}} @@ -104,13 +99,7 @@ #available: none, paterson and arrhenius self.rheology_law = 'Paterson' - # GIA: - self.lithosphere_shear_modulus = 6.7e10 # (Pa) - self.lithosphere_density = 3.32 # (g / cm^ - 3) - self.mantle_shear_modulus = 1.45e11 # (Pa) - self.mantle_density = 3.34 # (g / cm^ - 3) - - #SLR + #SLC self.earth_density = 5512 #average density of the Earth, (kg / m^3) return self #}}} @@ -130,12 +119,6 @@ md = checkfield(md, 'fieldname', 'materials.mantle_density', '>', 0, 'numel', [1]) md = checkfield(md, 'fieldname', 'materials.earth_density', '>', 0, 'numel', [1]) - 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) @@ -160,10 +143,6 @@ WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'rheology_B', 'format', 'DoubleMat', 'mattype', 1) 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: ../trunk-jpl/src/m/classes/matenhancedice.py =================================================================== --- ../trunk-jpl/src/m/classes/matenhancedice.py (revision 26058) +++ ../trunk-jpl/src/m/classes/matenhancedice.py (revision 26059) @@ -1,15 +1,14 @@ +from checkfield import checkfield from fielddisplay import fielddisplay from project3d import project3d -from checkfield import checkfield from WriteData import WriteData class matenhancedice(object): - """ - MATICE class definition + """MATICE class definition - Usage: - matenhancedice = matenhancedice() + Usage: + matenhancedice = matenhancedice() """ def __init__(self): #{{{ @@ -31,18 +30,17 @@ self.rheology_n = float('NaN') self.rheology_law = '' - #GIA - self.lithosphere_shear_modulus = 0. - self.lithosphere_density = 0. - self.mantle_shear_modulus = 0. - self.mantle_density = 0. + #SLC + self.earth_density = 0 # average density of the Earth, (kg/m^3) - #SLR - self.earth_density = 0 # average density of the Earth, (kg/m^3) self.setdefaultparameters() #}}} def __repr__(self): #{{{ + # TODO: + # - Convert all formatting to calls to .format (see any + # already converted .__repr__ method for examples) + # s = " Materials:" s = "%s\n%s" % (s, fielddisplay(self, "rho_ice", "ice density [kg/m^3]")) s = "%s\n%s" % (s, fielddisplay(self, "rho_water", "water density [kg/m^3]")) @@ -61,10 +59,6 @@ s = "%s\n%s" % (s, fielddisplay(self, "rheology_B", "flow law parameter [Pa s^(1/n)]")) s = "%s\n%s" % (s, fielddisplay(self, "rheology_n", "Glen's flow law exponent")) s = "%s\n%s" % (s, fielddisplay(self, "rheology_law", "law for the temperature dependance of the rheology: 'None', 'BuddJacka', 'Cuffey', 'CuffeyTemperate', 'Paterson', 'Arrhenius' or 'LliboutryDuval'")) - s = "%s\n%s" % (s, fielddisplay(self, "lithosphere_shear_modulus", "Lithosphere shear modulus [Pa]")) - s = "%s\n%s" % (s, fielddisplay(self, "lithosphere_density", "Lithosphere density [g/cm^-3]")) - s = "%s\n%s" % (s, fielddisplay(self, "mantle_shear_modulus", "Mantle shear modulus [Pa]")) - s = "%s\n%s" % (s, fielddisplay(self, "mantle_density", "Mantle density [g/cm^-3]")) s = "%s\n%s" % (s, fielddisplay(self, "earth_density", "Mantle density [kg/m^-3]")) return s @@ -108,13 +102,13 @@ #available: none, paterson and arrhenius self.rheology_law = 'Paterson' - # GIA: + #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 + #SLC self.earth_density = 5512 #average density of the Earth, (kg / m^3) return self @@ -131,12 +125,6 @@ md = checkfield(md, 'fieldname', 'materials.rheology_law', 'values', ['None', 'BuddJacka', 'Cuffey', 'CuffeyTemperate', 'Paterson', 'Arrhenius', 'LliboutryDuval']) md = checkfield(md, 'fieldname', 'materials.effectiveconductivity_averaging', 'numel', [1], 'values', [0, 1, 2]) - 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 @@ -161,9 +149,5 @@ 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: ../trunk-jpl/src/m/classes/materials.m =================================================================== --- ../trunk-jpl/src/m/classes/materials.m (revision 26058) +++ ../trunk-jpl/src/m/classes/materials.m (revision 26059) @@ -122,7 +122,6 @@ self.rheology_B = 1*1e8; self.rheology_n = 3; - case 'litho' %we default to a configuration that enables running GIA solutions using giacaron and/or giaivins. self.numlayers=2; Index: ../trunk-jpl/src/m/classes/materials.py =================================================================== --- ../trunk-jpl/src/m/classes/materials.py (revision 26058) +++ ../trunk-jpl/src/m/classes/materials.py (revision 26059) @@ -7,12 +7,11 @@ class materials(object): - ''' - MATERIALS class definition + """MATERIALS class definition - Usage: - materials = materials() - ''' + Usage: + materials = materials() + """ def __init__(self, *args): #{{{ self.nature = [] @@ -37,6 +36,7 @@ setattr(self, 'latentheat', 0) setattr(self, 'thermalconductivity', 0) setattr(self, 'temperateiceconductivity', 0) + setattr(self, 'effectiveconductivity_averaging', 0) setattr(self, 'meltingpoint', 0) setattr(self, 'beta', 0) setattr(self, 'mixed_layer_capacity', 0) @@ -59,9 +59,9 @@ setattr(self, 'rho_ice', 0) setattr(self, 'rho_water', 0) setattr(self, 'rho_freshwater', 0) - setattr(self, 'earth_density', 0) else: raise RuntimeError("materials constructor error message: nature of the material not supported yet! ('ice' or 'litho' or 'hydro')") + setattr(self, 'earth_density', 0) #set default parameters: self.setdefaultparameters() @@ -73,37 +73,58 @@ if nat == 'ice': #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.34e5 + #ice thermal conductivity (W/m/K) self.thermalconductivity = 2.4 + #wet ice thermal conductivity (W/m/K) self.temperateiceconductivity = 0.24 + + #computation of effective conductivity + self.effectiveconductivity_averaging = 1 + #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.8e-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.00e-4 + #Rheology law: what is the temperature dependence of B with T #available: none, paterson and arrhenius self.rheology_law = 'Paterson' + + #Rheology fields default + self.rheology_B = 1e8 + self.rheology_n = 3 elif nat == 'litho': - #we default to a configuration that enables running GIA solutions using giacaron and / or giaivins. + #we default to a configuration that enables running GIA solutions using giacaron and/or giaivins. self.numlayers = 2 + #center of the earth (approximation, must not be 0), then the lab (lithosphere / asthenosphere boundary) then the surface #(with 1d3 to avoid numerical singularities) self.radius = [1e3, 6278e3, 6378e3] + self.viscosity = [1e21, 1e40] #mantle and lithosphere viscosity (respectively) [Pa.s] self.lame_mu = [1.45e11, 6.7e10] # (Pa) #lithosphere and mantle shear modulus (respectively) [Pa] self.lame_lambda = self.lame_mu # (Pa) #mantle and lithosphere lamba parameter (respectively) [Pa] @@ -115,14 +136,17 @@ elif nat == 'hydro': #ice density (kg/m^3) self.rho_ice = 917. + #ocean water density (kg/m^3) self.rho_water = 1023. - #average density of the Earth (kg/m^3) - self.earth_density = 5512 + #fresh water density (kg/m^3) self.rho_freshwater = 1000. else: raise RuntimeError("materials setdefaultparameters error message: nature of the material not supported yet! ('ice' or 'litho' or 'hydro')") + + #average density of the Earth (kg/m^3) + self.earth_density = 5512 #}}} def __repr__(self): #{{{ @@ -222,7 +246,6 @@ #1: MatdamageiceEnum 2: MatestarEnum 3: MaticeEnum 4: MatenhancediceEnum 5: MaterialsEnum WriteData(fid, prefix, 'name', 'md.materials.nature', 'data', naturetointeger(self.nature), 'format', 'IntMat', 'mattype', 3) WriteData(fid, prefix, 'name', 'md.materials.type', 'data', 5, 'format', 'Integer') #DANGER, this can evolve if you have classes - for i in range(len(self.nature)): nat = self.nature[i] if nat == 'ice': @@ -235,6 +258,7 @@ 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', 'effectiveconductivity_averaging', 'format', 'Integer') 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') @@ -253,13 +277,19 @@ WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'isburgers', 'format', 'DoubleMat', 'mattype', 3) WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'burgers_viscosity', 'format', 'DoubleMat', 'mattype', 3) WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'burgers_mu', 'format', 'DoubleMat', 'mattype', 3) + # Compute earth density compatible with our layer density distribution + earth_density = 0 + for i in range(len(self.numlayers)): + earth_density = earth_density + (self.radius[i + 1] ** 3 - self.radius[i] ** 3) * self.density[i] + earth_density = earth_density / self.radius[self.numlayers + 1] ** 3 + self.earth_density = earth_density elif nat == 'hydro': 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', 'earth_density', 'format', 'Double') WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'rho_freshwater', 'format', 'Double') else: raise RuntimeError("materials constructor error message: nature of the material not supported yet! ('ice' or 'litho' or 'hydro')") + WriteData(fid, prefix, 'data', self.earth_density, 'name', 'md.materials.earth_density', 'format', 'Double') #}}} def extrude(self, md): #{{{ Index: ../trunk-jpl/src/m/classes/matestar.py =================================================================== --- ../trunk-jpl/src/m/classes/matestar.py (revision 26058) +++ ../trunk-jpl/src/m/classes/matestar.py (revision 26059) @@ -33,13 +33,7 @@ self.rheology_Es = np.nan self.rheology_law = '' - #GIA - self.lithosphere_shear_modulus = 0. - self.lithosphere_density = 0. - self.mantle_shear_modulus = 0. - self.mantle_density = 0. - - #SLR + #slc self.earth_density = 0 #set default parameters @@ -65,10 +59,6 @@ s = "%s\n%s" % (s, fielddisplay(self, 'rheology_Ec', 'compressive enhancement factor')) s = "%s\n%s" % (s, fielddisplay(self, 'rheology_Es', 'shear enhancement factor')) s = "%s\n%s" % (s, fielddisplay(self, 'rheology_law', ['law for the temperature dependance of the rheology: ''None'', ''BuddJacka'', ''Cuffey'', ''CuffeyTemperate'', ''Paterson'', ''Arrhenius'' or ''LliboutryDuval'''])) - s = "%s\n%s" % (s, fielddisplay(self, 'lithosphere_shear_modulus', 'Lithosphere shear modulus [Pa]')) - s = "%s\n%s" % (s, fielddisplay(self, 'lithosphere_density', 'Lithosphere density [g/cm^-3]')) - s = "%s\n%s" % (s, fielddisplay(self, 'mantle_shear_modulus', 'Mantle shear modulus [Pa]')) - s = "%s\n%s" % (s, fielddisplay(self, 'mantle_density', 'Mantle density [g/cm^-3]')) s = "%s\n%s" % (s, fielddisplay(self, 'earth_density', 'Mantle density [kg/m^-3]')) return s @@ -112,12 +102,7 @@ #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.7e10 # (Pa) - self.lithosphere_density = 3.32 # (g / cm^ - 3) - self.mantle_shear_modulus = 1.45e11 # (Pa) - self.mantle_density = 3.34 # (g / cm^ - 3) - #SLR + #slc self.earth_density = 5512 # average density of the Earth, (kg / m^3) return self @@ -134,12 +119,6 @@ md = checkfield(md, 'fieldname', 'materials.rheology_law', 'values', ['None', 'BuddJacka', 'Cuffey', 'CuffeyTemperate', 'Paterson', 'Arrhenius', 'LliboutryDuval']) md = checkfield(md, 'fieldname', 'materials.effectiveconductivity_averaging', 'numel', [1], 'values', [0, 1, 2]) - 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) @@ -165,9 +144,5 @@ 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', 1.0e3) - 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', 1.0e3) WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'earth_density', 'format', 'Double') # }}} Index: ../trunk-jpl/src/m/classes/matice.m =================================================================== --- ../trunk-jpl/src/m/classes/matice.m (revision 26058) +++ ../trunk-jpl/src/m/classes/matice.m (revision 26059) @@ -103,20 +103,18 @@ end % }}} function md = checkconsistency(self,md,solution,analyses) % {{{ - + if strcmpi(solution,'TransientSolution') & md.transient.isslc, md = checkfield(md,'fieldname','materials.earth_density','>',0,'numel',1); - return; - end - - 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,'universal',1,'NaN',1,'Inf',1); - md = checkfield(md,'fieldname','materials.rheology_n','>',0,'universal',1,'NaN',1,'Inf',1); - md = checkfield(md,'fieldname','materials.rheology_law','values',{'None' 'BuddJacka' 'Cuffey' 'CuffeyTemperate' 'Paterson' 'Arrhenius' 'LliboutryDuval' 'NyeCO2' 'NyeH2O'}); - md = checkfield(md,'fieldname','materials.effectiveconductivity_averaging','numel',[1],'values',[0 1 2]); + else + 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,'universal',1,'NaN',1,'Inf',1); + md = checkfield(md,'fieldname','materials.rheology_n','>',0,'universal',1,'NaN',1,'Inf',1); + md = checkfield(md,'fieldname','materials.rheology_law','values',{'None' 'BuddJacka' 'Cuffey' 'CuffeyTemperate' 'Paterson' 'Arrhenius' 'LliboutryDuval' 'NyeCO2' 'NyeH2O'}); + md = checkfield(md,'fieldname','materials.effectiveconductivity_averaging','numel',[1],'values',[0 1 2]); end % }}} function disp(self) % {{{ Index: ../trunk-jpl/src/m/classes/model.m =================================================================== --- ../trunk-jpl/src/m/classes/model.m (revision 26058) +++ ../trunk-jpl/src/m/classes/model.m (revision 26059) @@ -43,7 +43,7 @@ frontalforcings = 0; love = 0; esa = 0; - sampling = 0; + sampling = 0; autodiff = 0; inversion = 0; @@ -155,7 +155,7 @@ %2019 Jan.. if isa(md.frontalforcings,'double'); if(isprop('meltingrate',md.calving) & ~isnan(md.calving.meltingrate)) - disp('Warning: md.calving.meltingrate is now in md.frontalforcings'); + gia disp('Warning: md.calving.meltingrate is now in md.frontalforcings'); end md.frontalforcings=frontalforcings(md.calving); end @@ -186,8 +186,8 @@ md.smb.isconstrainsurfaceT = 0; end end - end - %2021 February 17 + end + %2021 February 17 if isa(md.sampling,'double'); md.sampling=sampling(); end end% }}} end @@ -241,7 +241,7 @@ md.frontalforcings = frontalforcings(); md.love = fourierlove(); md.esa = esa(); - md.sampling = sampling(); + md.sampling = sampling(); md.autodiff = autodiff(); md.inversion = inversion(); md.qmu = qmu(); @@ -614,8 +614,8 @@ %size = number of elements * n elseif fieldsize(1)==numberofelements1 md2.(model_fields{i}).(object_fields{j})=field(pos_elem,:); - elseif (fieldsize(1)==numberofelements1+1) - md2.(model_fields{i}).(object_fields{j})=[field(pos_elem,:); field(end,:)]; + elseif (fieldsize(1)==numberofelements1+1) + md2.(model_fields{i}).(object_fields{j})=[field(pos_elem,:); field(end,:)]; end end else @@ -627,7 +627,7 @@ %size = number of elements * n elseif fieldsize(1)==numberofelements1 md2.(model_fields{i})=field(pos_elem,:); - elseif (fieldsize(1)==numberofelements1+1) + elseif (fieldsize(1)==numberofelements1+1) md2.(model_fields{i})=[field(pos_elem,:); field(end,:)]; end end @@ -774,18 +774,18 @@ %get subfields % loop over time steps for p=1:length(md1.results.(solutionfields{i})) - current = md1.results.(solutionfields{i})(p); - solutionsubfields=fields(current); - for j=1:length(solutionsubfields), + current = md1.results.(solutionfields{i})(p); + solutionsubfields=fields(current); + for j=1:length(solutionsubfields), field=md1.results.(solutionfields{i})(p).(solutionsubfields{j}); if length(field)==numberofvertices1, - md2.results.(solutionfields{i})(p).(solutionsubfields{j})=field(pos_node); + md2.results.(solutionfields{i})(p).(solutionsubfields{j})=field(pos_node); elseif length(field)==numberofelements1, - md2.results.(solutionfields{i})(p).(solutionsubfields{j})=field(pos_elem); + md2.results.(solutionfields{i})(p).(solutionsubfields{j})=field(pos_elem); else - md2.results.(solutionfields{i})(p).(solutionsubfields{j})=field; + md2.results.(solutionfields{i})(p).(solutionsubfields{j})=field; end - end + end end else field=md1.results.(solutionfields{i}); @@ -1500,7 +1500,7 @@ disp(sprintf('%19s: %-22s -- %s','frontalforcings' ,['[1x1 ' class(self.frontalforcings) ']'],'parameters for frontalforcings')); disp(sprintf('%19s: %-22s -- %s','esa' ,['[1x1 ' class(self.esa) ']'],'parameters for elastic adjustment solution')); disp(sprintf('%19s: %-22s -- %s','love' ,['[1x1 ' class(self.love) ']'],'parameters for love solution')); - disp(sprintf('%19s: %-22s -- %s','sampling' ,['[1x1 ' class(self.sampling) ']'],'parameters for stochastic sampler')); + disp(sprintf('%19s: %-22s -- %s','sampling' ,['[1x1 ' class(self.sampling) ']'],'parameters for stochastic sampler')); disp(sprintf('%19s: %-22s -- %s','autodiff' ,['[1x1 ' class(self.autodiff) ']'],'automatic differentiation parameters')); disp(sprintf('%19s: %-22s -- %s','inversion' ,['[1x1 ' class(self.inversion) ']'],'parameters for inverse methods')); disp(sprintf('%19s: %-22s -- %s','qmu' ,['[1x1 ' class(self.qmu) ']'],'Dakota properties')); Index: ../trunk-jpl/src/m/classes/model.py =================================================================== --- ../trunk-jpl/src/m/classes/model.py (revision 26058) +++ ../trunk-jpl/src/m/classes/model.py (revision 26059) @@ -54,8 +54,6 @@ from thermal import thermal from steadystate import steadystate from transient import transient -from giaivins import giaivins -from giamme import giamme from esa import esa from autodiff import autodiff from inversion import inversion @@ -113,7 +111,6 @@ self.levelset = None self.calving = None self.frontalforcings = None - self.gia = None self.love = None self.esa = None self.sampling = None @@ -171,7 +168,6 @@ 'levelset', 'calving', 'frontalforcings', - 'gia', 'love', 'esa', 'sampling', @@ -224,7 +220,6 @@ s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("levelset", "[%s %s]" % ("1x1", obj.levelset.__class__.__name__), "parameters for moving boundaries (level - set method)")) s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("calving", "[%s %s]" % ("1x1", obj.calving.__class__.__name__), "parameters for calving")) s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("frontalforcings", "[%s %s]" % ("1x1", obj.frontalforcings.__class__.__name__), "parameters for frontalforcings")) - s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("gia", "[%s %s]" % ("1x1", obj.gia.__class__.__name__), "parameters for gia solution")) s = "%s\n%s" % (s, '%19s: %-22s -- %s' % ("esa", "[%s %s]" % ("1x1", obj.esa.__class__.__name__), "parameters for elastic adjustment solution")) s = "%s\n%s" % (s, '%19s: %-22s -- %s' % ("sampling", "[%s %s]" % ("1x1", obj.sampling.__class__.__name__), "parameters for stochastic sampler")) s = "%s\n%s" % (s, '%19s: %-22s -- %s' % ("love", "[%s %s]" % ("1x1", obj.love.__class__.__name__), "parameters for love solution")) @@ -271,7 +266,6 @@ self.levelset = levelset() self.calving = calving() self.frontalforcings = frontalforcings() - self.gia = giamme() self.love = fourierlove() self.esa = esa() self.sampling = sampling() @@ -851,13 +845,6 @@ if not np.isnan(md.initialization.watercolumn).all(): md.initialization.watercolumn = project2d(md, md.initialization.watercolumn, 1) - # giaivins - if md.gia.__class__.__name__ == 'giaivins': - if not np.isnan(md.gia.mantle_viscosity).all(): - md.gia.mantle_viscosity = project2d(md, md.gia.mantle_viscosity, 1) - if not np.isnan(md.gia.lithosphere_thickness).all(): - md.gia.lithosphere_thickness = project2d(md, md.gia.lithosphere_thickness, 1) - # elementstype if not np.isnan(md.flowequation.element_equation).all(): md.flowequation.element_equation = project2d(md, md.flowequation.element_equation, 1) Index: ../trunk-jpl/src/m/classes/sealevelmodel.m =================================================================== --- ../trunk-jpl/src/m/classes/sealevelmodel.m (revision 26058) +++ ../trunk-jpl/src/m/classes/sealevelmodel.m (revision 26059) @@ -90,15 +90,15 @@ for i=1:length(slm.icecaps), md= slm.icecaps{i}; if ~isempty(md.solidearth.external), - error(sprintf('cannot run external forcings on an ice sheet when runing a coupling earth/ice sheet model'); + error('cannot run external forcings on an ice sheet when running a coupling earth/ice sheet model'); end end - %make sure that we have the rigth grd model for computing ouf sealevel patterns: + %make sure that we have the right grd model for computing out sealevel patterns: for i=1:length(slm.icecaps), md= slm.icecaps{i}; if md.solidearth.settings.grdmodel~=0 - error(sprintf('sealevelmodel checkconsistenty error message: ice sheets do not run GRD module, specify solidearth.settings.grdmodel=0 on ice cap %i',i)); + error(sprintf('sealevelmodel checkconsistency error message: ice sheets do not run GRD module, specify solidearth.settings.grdmodel=0 on ice cap %i',i)); end end Index: ../trunk-jpl/src/m/classes/sealevelmodel.py =================================================================== --- ../trunk-jpl/src/m/classes/sealevelmodel.py (revision 26058) +++ ../trunk-jpl/src/m/classes/sealevelmodel.py (revision 26059) @@ -110,6 +110,22 @@ md = slm.icecaps[i] if np.nonzero(md.slr.steric_rate - slm.earth.slr.steric_rate[slm.earth.slr.transitions[i]]) != []: raise Exception('steric rate on ice cap {} is not the same as for the earth'.format(md.miscellaneous.name)) + + # Make sure grd is the same everywhere + for i in range(len(slm.icecaps)): + md = slm.icecaps[i] + if md.solidearthsettings.isgrd != slm.earth.solidearthsettings.isgrd: + raise RuntimeError('isgrd on ice cap {} is not the same as for the earth\n'.format(md.miscellaneous.name)) + # Make sure that there is no solid earth external forcing on the basins + for i in range(len(slm.icecaps)): + md = slm.icecaps[i] + if md.solidearth.external: + raise RuntimeError('cannot run external forcings on an ice sheet when running a coupling earth/ice sheet model') + # Make sure that we have the right grd model for computing our sealevel patterns + for i in range(len(slm.icecaps)): + md = slm.icecaps[i] + if md.solidearth.settings.grdmodel != 0: + raise RuntimeError('sealevelmodel checkconsistency error message: ice sheets do not run GRD module, specify solidearth.settings.grdmodel=0 on ice cap {}'.format(i)) #}}} def mergeresults(self): # {{{ Index: ../trunk-jpl/src/m/classes/transient.py =================================================================== --- ../trunk-jpl/src/m/classes/transient.py (revision 26058) +++ ../trunk-jpl/src/m/classes/transient.py (revision 26059) @@ -13,10 +13,10 @@ def __init__(self): # {{{ self.issmb = 0 self.ismasstransport = 0 + self.isoceantransport = 0 self.isstressbalance = 0 self.isthermal = 0 self.isgroundingline = 0 - self.isgia = 0 self.isesa = 0 self.isdamageevolution = 0 self.ismovingfront = 0 @@ -23,7 +23,6 @@ self.ishydrology = 0 self.issampling = 0 self.isslc = 0 - self.iscoupler = 0 self.amr_frequency = 0 self.isoceancoupling = 0 self.requested_outputs = [] @@ -35,10 +34,10 @@ s = ' transient solution parameters:\n' s += '{}\n'.format(fielddisplay(self, 'issmb', 'indicates if a surface mass balance solution is used in the transient')) s += '{}\n'.format(fielddisplay(self, 'ismasstransport', 'indicates if a masstransport solution is used in the transient')) + s += '{}\n'.format(fielddisplay(self, 'isoceantransport', 'indicates whether an ocean masstransport solution is used in the transient')) s += '{}\n'.format(fielddisplay(self, 'isstressbalance', 'indicates if a stressbalance solution is used in the transient')) s += '{}\n'.format(fielddisplay(self, 'isthermal', 'indicates if a thermal solution is used in the transient')) s += '{}\n'.format(fielddisplay(self, 'isgroundingline', 'indicates if a groundingline migration is used in the transient')) - s += '{}\n'.format(fielddisplay(self, 'isgia', 'indicates if a postglacial rebound is used in the transient')) s += '{}\n'.format(fielddisplay(self, 'isesa', 'indicates whether an elastic adjustment model is used in the transient')) s += '{}\n'.format(fielddisplay(self, 'isdamageevolution', 'indicates whether damage evolution is used in the transient')) s += '{}\n'.format(fielddisplay(self, 'ismovingfront', 'indicates whether a moving front capability is used in the transient')) @@ -46,7 +45,6 @@ s += '{}\n'.format(fielddisplay(self, 'issampling', 'indicates whether sampling is used in the transient')) s += '{}\n'.format(fielddisplay(self, 'isslc', 'indicates if a sea level change solution is used in the transient')) s += '{}\n'.format(fielddisplay(self, 'isoceancoupling', 'indicates whether coupling with an ocean model is used in the transient')) - s += '{}\n'.format(fielddisplay(self, 'iscoupler', 'indicates whether different models are being run with need for coupling')) s += '{}\n'.format(fielddisplay(self, 'amr_frequency', 'frequency at which mesh is refined in simulations with multiple time_steps')) s += '{}\n'.format(fielddisplay(self, 'requested_outputs', 'list of additional outputs requested')) return s @@ -62,10 +60,10 @@ def deactivateall(self): #{{{ self.issmb = 0 self.ismasstransport = 0 + self.isoceantransport = 0 self.isstressbalance = 0 self.isthermal = 0 self.isgroundingline = 0 - self.isgia = 0 self.isesa = 0 self.isdamageevolution = 0 self.ismovingfront = 0 @@ -73,7 +71,6 @@ self.issampling = 0 self.isslc = 0 self.isoceancoupling = 0 - self.iscoupler = 0 self.amr_frequency = 0 # Default output @@ -85,10 +82,10 @@ #full analysis: Stressbalance, Masstransport and Thermal but no groundingline migration for now self.issmb = 1 self.ismasstransport = 1 + self.isoceantransport = 0 self.isstressbalance = 1 self.isthermal = 1 self.isgroundingline = 0 - self.isgia = 0 self.isesa = 0 self.isdamageevolution = 0 self.ismovingfront = 0 @@ -96,7 +93,6 @@ self.issampling = 0 self.isslc = 0 self.isoceancoupling = 0 - self.iscoupler = 0 self.amr_frequency = 0 # Default output @@ -111,10 +107,10 @@ md = checkfield(md, 'fieldname', 'transient.issmb', 'numel', [1], 'values', [0, 1]) md = checkfield(md, 'fieldname', 'transient.ismasstransport', 'numel', [1], 'values', [0, 1]) + md = checkfield(md, 'fieldname', 'transient.isocealtransport', 'numel', [1], 'values', [0, 1]) md = checkfield(md, 'fieldname', 'transient.isstressbalance', 'numel', [1], 'values', [0, 1]) md = checkfield(md, 'fieldname', 'transient.isthermal', 'numel', [1], 'values', [0, 1]) md = checkfield(md, 'fieldname', 'transient.isgroundingline', 'numel', [1], 'values', [0, 1]) - md = checkfield(md, 'fieldname', 'transient.isgia', 'numel', [1], 'values', [0, 1]) md = checkfield(md, 'fieldname', 'transient.isesa', 'numel', [1], 'values', [0, 1]) md = checkfield(md, 'fieldname', 'transient.isdamageevolution', 'numel', [1], 'values', [0, 1]) md = checkfield(md, 'fieldname', 'transient.ishydrology', 'numel', [1], 'values', [0, 1]) @@ -122,7 +118,6 @@ md = checkfield(md, 'fieldname', 'transient.ismovingfront', 'numel', [1], 'values', [0, 1]) md = checkfield(md, 'fieldname', 'transient.isslc', 'numel', [1], 'values', [0, 1]) md = checkfield(md, 'fieldname', 'transient.isoceancoupling', 'numel', [1], 'values', [0, 1]) - md = checkfield(md, 'fieldname', 'transient.iscoupler', 'numel', [1], 'values', [0, 1]) md = checkfield(md, 'fieldname', 'transient.amr_frequency', 'numel', [1], '>=', 0, 'NaN', 1, 'Inf', 1) if solution != 'TransientSolution' and md.transient.iscoupling: @@ -135,10 +130,10 @@ def marshall(self, prefix, md, fid): # {{{ WriteData(fid, prefix, 'object', self, 'fieldname', 'issmb', 'format', 'Boolean') WriteData(fid, prefix, 'object', self, 'fieldname', 'ismasstransport', 'format', 'Boolean') + WriteData(fid, prefix, 'object', self, 'fieldname', 'isoceantransport', 'format', 'Boolean') WriteData(fid, prefix, 'object', self, 'fieldname', 'isstressbalance', 'format', 'Boolean') WriteData(fid, prefix, 'object', self, 'fieldname', 'isthermal', 'format', 'Boolean') WriteData(fid, prefix, 'object', self, 'fieldname', 'isgroundingline', 'format', 'Boolean') - WriteData(fid, prefix, 'object', self, 'fieldname', 'isgia', 'format', 'Boolean') WriteData(fid, prefix, 'object', self, 'fieldname', 'isesa', 'format', 'Boolean') WriteData(fid, prefix, 'object', self, 'fieldname', 'isdamageevolution', 'format', 'Boolean') WriteData(fid, prefix, 'object', self, 'fieldname', 'ishydrology', 'format', 'Boolean') @@ -146,7 +141,6 @@ WriteData(fid, prefix, 'object', self, 'fieldname', 'issampling', 'format', 'Boolean') WriteData(fid, prefix, 'object', self, 'fieldname', 'isslc', 'format', 'Boolean') WriteData(fid, prefix, 'object', self, 'fieldname', 'isoceancoupling', 'format', 'Boolean') - WriteData(fid, prefix, 'object', self, 'fieldname', 'iscoupler', 'format', 'Boolean') WriteData(fid, prefix, 'object', self, 'fieldname', 'amr_frequency', 'format', 'Integer') # Process requested outputs Index: ../trunk-jpl/src/m/solve/solve.py =================================================================== --- ../trunk-jpl/src/m/solve/solve.py (revision 26058) +++ ../trunk-jpl/src/m/solve/solve.py (revision 26059) @@ -1,6 +1,5 @@ from datetime import datetime import os -import shutil from ismodelselfconsistent import ismodelselfconsistent from loadresultsfromcluster import loadresultsfromcluster @@ -21,6 +20,7 @@ Solution types available comprise: - 'Stressbalance' or 'sb' - 'Masstransport' or 'mt' + - 'Oceantransport' or 'oceant' - 'Thermal' or 'th' - 'Steadystate' or 'ss' - 'Transient' or 'tr' @@ -31,9 +31,8 @@ - 'Hydrology' or 'hy' - 'DamageEvolution' or 'da' - 'Gia' or 'gia' + - 'Love' or 'lv' - 'Esa' or 'esa' - - 'Sealevelchange' or 'slc' - - 'Love' or 'lv' - 'Sampling' or 'smp' Extra options: @@ -54,6 +53,8 @@ solutionstring = 'StressbalanceSolution' elif solutionstring.lower() == 'mt' or solutionstring.lower() == 'masstransport': solutionstring = 'MasstransportSolution' + elif solutionstring.lower() == 'oceant' or solutionstring.lower() == 'oceantransport': + solutionstring = 'OceantransportSolution' elif solutionstring.lower() == 'th' or solutionstring.lower() == 'thermal': solutionstring = 'ThermalSolution' elif solutionstring.lower() == 'st' or solutionstring.lower() == 'steadystate': @@ -78,8 +79,6 @@ solutionstring = 'LoveSolution' elif solutionstring.lower() == 'esa': solutionstring = 'EsaSolution' - elif solutionstring.lower() == 'slc' or solutionstring.lower() == 'sealevelchange': - solutionstring = 'SealevelchangeSolution' elif solutionstring.lower() == 'smp' or solutionstring.lower() == 'sampling': solutionstring = 'SamplingSolution' else: @@ -150,9 +149,9 @@ # Wait on lock if md.settings.waitonlock > 0: islock = waitonlock(md) - if islock == 0: # no results to be loaded + if islock == 0: # no results to be loaded print('The results must be loaded manually with md = loadresultsfromcluster(md).') - else: # load results + else: # load results if md.verbose.solution: print('loading results from cluster') md = loadresultsfromcluster(md) Index: ../trunk-jpl/test/NightlyRun/test2001.py =================================================================== --- ../trunk-jpl/test/NightlyRun/test2001.py (revision 26058) +++ ../trunk-jpl/test/NightlyRun/test2001.py (revision 26059) @@ -3,7 +3,7 @@ import numpy as np -from giaivins import giaivins +from materials import * from model import * from parameterize import * from setflowequation import * @@ -16,35 +16,50 @@ md = setmask(md, '', '') md = parameterize(md, '../Par/SquareSheetConstrained.py') -#GIA -md.gia = giaivins() -md.gia.lithosphere_thickness = 100. * np.ones(md.mesh.numberofvertices) # in km -md.gia.mantle_viscosity = 1.0e21 * np.ones(md.mesh.numberofvertices) # in Pa.s -md.materials.lithosphere_shear_modulus = 6.7e10 # in Pa -md.materials.lithosphere_density = 3.32 # in g/cm^3 -md.materials.mantle_shear_modulus = 1.45e11 # in Pa -md.materials.mantle_density = 3.34 # in g/cm^3 +# GIA Ivins, 2 layer model +md.solidearth.settings.grdmodel = 2 -#Indicate what you want to compute -md.gia.cross_section_shape = 1 # for square-edged x-section +md.materials = materials('litho','ice') +md.materials.numlayers = 2; +md.materials.radius = [10, 6271e3, 6371e3] +md.materials.density = [3.34e3, 3.32e3] +md.materials.lame_mu = [1.45e11, 6.7e10] +md.materials.viscosity = [1e21, 0] +md.initialization.sealevel = np.zeros(md.mesh.numberofvertices) +md.solidearth.settings.cross_section_shape = 1 # for square-edged x-section +md.solidearth.settings.computesealevelchange = 0 # do not compute sea level, only deformation +md.solidearth.requested_outputs = ['Sealevel', 'SealevelUGrd'] -#Define loading history (see test2001.m for the description) -md.timestepping.start_time = 2400000 # 2, 400 kyr -md.timestepping.final_time = 2500000 # 2, 500 kyr -md.geometry.thickness = np.array([ - np.append(md.geometry.thickness * 0.0, 0.0), - np.append(md.geometry.thickness / 2.0, 0.1), - np.append(md.geometry.thickness, 0.2), - np.append(md.geometry.thickness, 1.0), - np.append(md.geometry.thickness, md.timestepping.start_time) +# Loading history +md.timestepping.start_time = -2400000 # 4,800 kyr :: EVALUATION TIME +md.timestepping.time_step = 2400000 # 2,400 kyr :: EVALUATION TIME +# To get rid of default final_time, make sure final_time > start_time +md.timestepping.final_time = 2400000 # 2,400 kyr +md.masstransport.spcthickness np.array([ + np.append(md.geometry.thickness, 0), + np.append(md.geometry.thickness, 2400000) ]).T +# Geometry at 0 initially +md.geometry.thickness = np.zeros(md.mesh.numberofvertices) +md.geometry.surface = np.zeros(md.mesh.numberofvertices) +md.geometry.base = np.zeros(md.mesh.numberofvertices) + +# Physics +md.transient.issmb = 0 +md.transient.isstressbalance = 0 +md.transient.isthermal = 0 +md.transient.ismasstransport = 0 +md.transient.isslc = 0 + #Solve for GIA deflection +md.cluster = generic('name', gethostname(), 'np', 1) md.cluster = generic('name', gethostname(), 'np', 3) -md.verbose = verbose('1111111') -md = solve(md, 'Gia') +md.verbose = verbose('11111111111') +md.verbose.solver = 0 +md = solve(md, 'Transient') #Fields and tolerances to track changes -field_names = ['UGia', 'UGiaRate'] -field_tolerances = [1e-13, 1e-13] -field_values = [md.results.GiaSolution.UGia, md.results.GiaSolution.UGiaRate] +field_names = ['UGrd'] +field_tolerances = [1e-13] +field_values = [md.results.TransientSolution[0].SealevelUGrd]