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

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)