Index: /issm/trunk-jpl/src/m/classes/initialization.py =================================================================== --- /issm/trunk-jpl/src/m/classes/initialization.py (revision 26012) +++ /issm/trunk-jpl/src/m/classes/initialization.py (revision 26013) @@ -29,4 +29,5 @@ self.hydraulic_potential = np.nan self.channelarea = np.nan + self.sample = np.nan #set defaults @@ -50,4 +51,5 @@ s += '{}\n'.format(fielddisplay(self, 'hydraulic_potential', 'Hydraulic potential (for GlaDS) [Pa]')) s += '{}\n'.format(fielddisplay(self, 'channelarea', 'subglaciale water channel area (for GlaDS) [m2]')) + s += '{}\n'.format(fielddisplay(self,'sample','Realization of a Gaussian random field')) return s #}}} @@ -124,4 +126,7 @@ md = checkfield(md, 'fieldname', 'initialization.hydraulic_potential', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices]) md = checkfield(md, 'fieldname', 'initialization.channelarea', 'NaN', 1, 'Inf', 1, '>=', 0, 'size', [md.mesh.numberofelements]) + if 'SamplingAnalysis' in analyses and not solution == 'TransientSolution' and not md.transient.issampling: + if np.any(np.isnan(md.initialization.sample)): + md = checkfield(md,'fieldname','initialization.sample','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]) return md # }}} @@ -141,4 +146,5 @@ 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) if md.thermal.isenthalpy: if (np.size(self.enthalpy) <= 1): Index: /issm/trunk-jpl/src/m/classes/model.py =================================================================== --- /issm/trunk-jpl/src/m/classes/model.py (revision 26012) +++ /issm/trunk-jpl/src/m/classes/model.py (revision 26013) @@ -76,4 +76,5 @@ from contourenvelope import contourenvelope from DepthAverage import DepthAverage +from sampling import sampling #}}} @@ -116,4 +117,5 @@ self.love = None self.esa = None + self.sampling = None self.autodiff = None self.inversion = None @@ -173,4 +175,5 @@ 'love', 'esa', + 'sampling', 'autodiff', 'inversion', @@ -224,4 +227,5 @@ 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")) s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("autodiff", "[%s %s]" % ("1x1", obj.autodiff.__class__.__name__), "automatic differentiation parameters")) @@ -271,4 +275,5 @@ self.love = fourierlove() self.esa = esa() + self.sampling = sampling() self.autodiff = autodiff() self.inversion = inversion() Index: /issm/trunk-jpl/src/m/consistency/ismodelselfconsistent.py =================================================================== --- /issm/trunk-jpl/src/m/consistency/ismodelselfconsistent.py (revision 26012) +++ /issm/trunk-jpl/src/m/consistency/ismodelselfconsistent.py (revision 26013) @@ -79,4 +79,6 @@ elif 'DamageEvolutionSolution': analyses = ['DamageEvolutionAnalysis'] + elseif strcmp(solutiontype,'SamplingSolution') + analyses=['SamplingAnalysis'] else: raise TypeError('solution type: {} not supported yet!'.format(solutiontype)) Index: /issm/trunk-jpl/src/m/sampling.py =================================================================== --- /issm/trunk-jpl/src/m/sampling.py (revision 26012) +++ /issm/trunk-jpl/src/m/sampling.py (revision 26013) @@ -55,7 +55,7 @@ self.phi = 0 # Exponent in fraction SPDE (default=2, biLaplacian covariance operator) - self.alpha=2; + self.alpha = 2 # Seed for pseudorandom number generator (default -1 for random seed) - self.alpha=-1; + self.alpha = -1 # Default output self.requested_outputs = ['default'] @@ -68,12 +68,12 @@ return md - md = checkfield(md,'fieldname','sampling.kappa','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1],'>',0); - md = checkfield(md,'fieldname','sampling.tau','NaN',1,'Inf',1,'numel',1,'>',0); - md = checkfield(md,'fieldname','sampling.beta','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1],'>',0); - md = checkfield(md,'fieldname','sampling.phi','NaN',1,'Inf',1,'numel',1,'>=',0); - md = checkfield(md,'fieldname','sampling.alpha','NaN',1,'Inf',1,'numel',1,'>',0); - md = checkfield(md,'fieldname','sampling.robin','numel',1,'values',[0 1]); - md = checkfield(md,'fieldname','sampling.seed','NaN',1,'Inf',1,'numel',1); - md = checkfield(md,'fieldname','sampling.requested_outputs','stringrow',1); + md = checkfield(md,'fieldname','sampling.kappa','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1],'>',0) + md = checkfield(md,'fieldname','sampling.tau','NaN',1,'Inf',1,'numel',1,'>',0) + md = checkfield(md,'fieldname','sampling.beta','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1],'>',0) + md = checkfield(md,'fieldname','sampling.phi','NaN',1,'Inf',1,'numel',1,'>=',0) + md = checkfield(md,'fieldname','sampling.alpha','NaN',1,'Inf',1,'numel',1,'>',0) + md = checkfield(md,'fieldname','sampling.robin','numel',1,'values',[0 1]) + md = checkfield(md,'fieldname','sampling.seed','NaN',1,'Inf',1,'numel',1) + md = checkfield(md,'fieldname','sampling.requested_outputs','stringrow',1) return md @@ -81,11 +81,11 @@ def marshall(self, prefix, md, fid): # {{{ - WriteData(fid,prefix,'object',self,'fieldname','kappa','format','DoubleMat','mattype',1); - WriteData(fid,prefix,'object',self,'fieldname','tau','format','Double'); - WriteData(fid,prefix,'object',self,'fieldname','beta','format','DoubleMat','mattype',1); - WriteData(fid,prefix,'object',self,'fieldname','phi','format','Double'); - WriteData(fid,prefix,'object',self,'fieldname','alpha','format','Integer'); - WriteData(fid,prefix,'object',self,'fieldname','robin','format','Boolean'); - WriteData(fid,prefix,'object',self,'fieldname','seed','format','Integer'); + WriteData(fid,prefix,'object',self,'fieldname','kappa','format','DoubleMat','mattype',1) + WriteData(fid,prefix,'object',self,'fieldname','tau','format','Double') + WriteData(fid,prefix,'object',self,'fieldname','beta','format','DoubleMat','mattype',1) + WriteData(fid,prefix,'object',self,'fieldname','phi','format','Double') + WriteData(fid,prefix,'object',self,'fieldname','alpha','format','Integer') + WriteData(fid,prefix,'object',self,'fieldname','robin','format','Boolean') + WriteData(fid,prefix,'object',self,'fieldname','seed','format','Integer') # Process requested outputs @@ -98,9 +98,9 @@ # Process requested outputs - outputs = self.requested_outputs; + outputs = self.requested_outputs indices = [i for i, x in enumerate(outputs) if x == 'default'] if len(indices) > 0: outputscopy = outputs[0:max(0, indices[0] - 1)] + self.defaultoutputs(md) + outputs[indices[0] + 1:] outputs = outputscopy - WriteData(fid,prefix,'data',outputs,'name','md.sampling.requested_outputs','format','StringArray'); + WriteData(fid,prefix,'data',outputs,'name','md.sampling.requested_outputs','format','StringArray') # }}} Index: /issm/trunk-jpl/src/m/solve/solve.py =================================================================== --- /issm/trunk-jpl/src/m/solve/solve.py (revision 26012) +++ /issm/trunk-jpl/src/m/solve/solve.py (revision 26013) @@ -35,4 +35,5 @@ - 'Sealevelchange' or 'slc' - 'Love' or 'lv' + - 'Sampling' or 'smp' Extra options: @@ -80,4 +81,6 @@ elif solutionstring.lower() == 'slc' or solutionstring.lower() == 'sealevelchange': solutionstring = 'SealevelchangeSolution' + elif solutionstring.lower() == 'smp' or solutionstring.lower() == 'sampling': + solutionstring = 'SamplingSolution' else: raise ValueError("solutionstring '%s' not supported!" % solutionstring)