source: issm/oecreview/Archive/25834-26739/ISSM-26012-26013.diff

../trunk-jpl/src/m/classes/initialization.py

@@ -28 +28 @@
         self.epl_thickness = np.nan
         self.hydraulic_potential = np.nan
         self.channelarea = np.nan
+        self.sample = np.nan
 
         #set defaults
         self.setdefaultparameters()
@@ -49 +50 @@
         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 +125 @@
                 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 +145 @@
         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

../trunk-jpl/src/m/classes/model.py

@@ -75 +75 @@
 from ElementConnectivity import ElementConnectivity
 from contourenvelope import contourenvelope
 from DepthAverage import DepthAverage
+from sampling import sampling
 #}}}
 
 
@@ -115 +116 @@
         self.gia = None
         self.love = None
         self.esa = None
+        self.sampling = None
         self.autodiff = None
         self.inversion = None
         self.qmu = None
@@ -172 +174 @@
                 'gia',
                 'love',
                 'esa',
+                'sampling',
                 'autodiff',
                 'inversion',
                 'qmu',
@@ -223 +226 @@
         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 +274 @@
         self.gia = giamme()
         self.love = fourierlove()
         self.esa = esa()
+        self.sampling = sampling()
         self.autodiff = autodiff()
         self.inversion = inversion()
         self.qmu = qmu()

../trunk-jpl/src/m/consistency/ismodelselfconsistent.py

@@ -78 +78 @@
         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))
 

../trunk-jpl/src/m/sampling.py

@@ -54 +54 @@
         # 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 +67 @@
         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 +97 @@
         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

../trunk-jpl/src/m/solve/solve.py

@@ -34 +34 @@
     - 'Esa'                or 'esa'
     - 'Sealevelchange'     or 'slc'
     - 'Love'               or 'lv'
+    - 'Sampling'           or 'smp'
 
     Extra options:
     - loadonly         : does not solve. only load results
@@ -79 +80 @@
         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)