Index: ../trunk-jpl/src/m/classes/initialization.py =================================================================== --- ../trunk-jpl/src/m/classes/initialization.py (revision 26012) +++ ../trunk-jpl/src/m/classes/initialization.py (revision 26013) @@ -28,6 +28,7 @@ self.epl_thickness = np.nan self.hydraulic_potential = np.nan self.channelarea = np.nan + self.sample = np.nan #set defaults self.setdefaultparameters() @@ -49,6 +50,7 @@ s += '{}\n'.format(fielddisplay(self, 'epl_thickness', 'thickness of the epl [m]')) 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 #}}} @@ -123,6 +125,9 @@ md = checkfield(md, 'fieldname', 'initialization.watercolumn', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices]) 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 # }}} @@ -140,6 +145,7 @@ 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) if md.thermal.isenthalpy: if (np.size(self.enthalpy) <= 1): tpmp = md.materials.meltingpoint - md.materials.beta * md.initialization.pressure Index: ../trunk-jpl/src/m/classes/model.py =================================================================== --- ../trunk-jpl/src/m/classes/model.py (revision 26012) +++ ../trunk-jpl/src/m/classes/model.py (revision 26013) @@ -75,6 +75,7 @@ from ElementConnectivity import ElementConnectivity from contourenvelope import contourenvelope from DepthAverage import DepthAverage +from sampling import sampling #}}} @@ -115,6 +116,7 @@ self.gia = None self.love = None self.esa = None + self.sampling = None self.autodiff = None self.inversion = None self.qmu = None @@ -172,6 +174,7 @@ 'gia', 'love', 'esa', + 'sampling', 'autodiff', 'inversion', 'qmu', @@ -223,6 +226,7 @@ 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")) s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("autodiff", "[%s %s]" % ("1x1", obj.autodiff.__class__.__name__), "automatic differentiation parameters")) s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("inversion", "[%s %s]" % ("1x1", obj.inversion.__class__.__name__), "parameters for inverse methods")) @@ -270,6 +274,7 @@ self.gia = giamme() self.love = fourierlove() self.esa = esa() + self.sampling = sampling() self.autodiff = autodiff() self.inversion = inversion() self.qmu = qmu() Index: ../trunk-jpl/src/m/consistency/ismodelselfconsistent.py =================================================================== --- ../trunk-jpl/src/m/consistency/ismodelselfconsistent.py (revision 26012) +++ ../trunk-jpl/src/m/consistency/ismodelselfconsistent.py (revision 26013) @@ -78,6 +78,8 @@ analyses = ['L2ProjectionBaseAnalysis', 'HydrologyShreveAnalysis', 'HydrologyDCInefficientAnalysis', 'HydrologyDCEfficientAnalysis'] elif 'DamageEvolutionSolution': analyses = ['DamageEvolutionAnalysis'] + elseif strcmp(solutiontype,'SamplingSolution') + analyses=['SamplingAnalysis'] else: raise TypeError('solution type: {} not supported yet!'.format(solutiontype)) Index: ../trunk-jpl/src/m/sampling.py =================================================================== --- ../trunk-jpl/src/m/sampling.py (revision 26012) +++ ../trunk-jpl/src/m/sampling.py (revision 26013) @@ -54,9 +54,9 @@ # Temporal correlation factor 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'] return self @@ -67,26 +67,26 @@ if (SamplingAnalysis' not in analyses): 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 # }}} 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 outputs = self.requested_outputs @@ -97,10 +97,10 @@ WriteData(fid, prefix, 'data', outputs, 'name', 'md.masstransport.requested_outputs', 'format', 'StringArray') # 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') # }}} \ No newline at end of file Index: ../trunk-jpl/src/m/solve/solve.py =================================================================== --- ../trunk-jpl/src/m/solve/solve.py (revision 26012) +++ ../trunk-jpl/src/m/solve/solve.py (revision 26013) @@ -34,6 +34,7 @@ - 'Esa' or 'esa' - 'Sealevelchange' or 'slc' - 'Love' or 'lv' + - 'Sampling' or 'smp' Extra options: - loadonly : does not solve. only load results @@ -79,6 +80,8 @@ solutionstring = 'EsaSolution' 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) options = pairoptions('solutionstring', solutionstring, *args)