Changeset 26013

Timestamp:

02/21/21 15:55:44 (4 years ago)

Author:

bulthuis

Message:

NEW: add sampling capability for python

Location:

issm/trunk-jpl/src/m

Files:

• classes/initialization.py (modified) (4 diffs)
• classes/model.py (modified) (5 diffs)
• consistency/ismodelselfconsistent.py (modified) (1 diff)
• sampling.py (modified) (4 diffs)
• solve/solve.py (modified) (2 diffs)

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

@@ -29 +29 @@
         self.hydraulic_potential = np.nan
         self.channelarea = np.nan
+        self.sample = np.nan
 
         #set defaults
@@ -50 +51 @@
         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 +126 @@
                 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 +146 @@
         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):

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

@@ -76 +76 @@
 from contourenvelope import contourenvelope
 from DepthAverage import DepthAverage
+from sampling import sampling
 #}}}
 
@@ -116 +117 @@
         self.love = None
         self.esa = None
+        self.sampling = None
         self.autodiff = None
         self.inversion = None
@@ -173 +175 @@
                 'love',
                 'esa',
+                'sampling',
                 'autodiff',
                 'inversion',
@@ -224 +227 @@
         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 +275 @@
         self.love = fourierlove()
         self.esa = esa()
+        self.sampling = sampling()
         self.autodiff = autodiff()
         self.inversion = inversion()

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

@@ -79 +79 @@
     elif 'DamageEvolutionSolution':
         analyses = ['DamageEvolutionAnalysis']
+    elseif strcmp(solutiontype,'SamplingSolution')
+                analyses=['SamplingAnalysis']  
     else:
         raise TypeError('solution type: {} not supported yet!'.format(solutiontype))

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

@@ -55 +55 @@
         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 +68 @@
             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 +81 @@
 
     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 +98 @@
 
         # 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')
     # }}}

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

@@ -35 +35 @@
     - 'Sealevelchange'     or 'slc'
     - 'Love'               or 'lv'
+    - 'Sampling'           or 'smp'
 
     Extra options:
@@ -80 +81 @@
     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)