source: issm/oecreview/Archive/21724-22754/ISSM-22103-22104.diff@ 22755

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

@@ +1 @@
+import numpy as np
+from project3d import project3d
+from pairoptions import *
+from collections import OrderedDict
+from fielddisplay import fielddisplay
+from checkfield import checkfield
+from WriteData import WriteData
+
+
+class misfit(object):
+        """
+        MISFIT class definition
+
+        Usage:
+                misfit=misfit()
+                misfit=misfit(name='SurfaceAltimetry',
+                        definitionstring='Outputdefinition1',
+                        model_string='Surface',
+                        observation_string='SurfaceObservations',
+                        observation=md.geometry.surface,
+                        timeinterpolation='nearestneighbor',
+                        local=1,
+                        weights=np.ones((md.mesh.numberofvertices,1)),
+                        weights_string='WeightsSurfaceObservations')
+        """
+
+        def __init__(self, name = None, definitionstring = None, model_string = None, observation = None, observation_string = None, timeinterpolation = None, local = None, weights = None, weights_string = None, cumulated = None):
+                # {{{
+                self.name = name if name is not None else ''
+
+                #string that identifies this output definition uniquely, from 'Outputdefinition[1-100]'
+                self.definitionstring = definitionstring if definitionstring is not None else ''
+
+                #string for field that is modeled
+                self.model_string = model_string if model_string is not None else ''
+
+                #observed field that we compare the model against
+                self.observation = observation if observation is not None else float('NaN')
+
+                #string for observed field.
+                self.observation_string = observation_string if observation_string is not None else ''
+
+                self.timeinterpolation = timeinterpolation if timeinterpolation is not None else 'nearestneighbor'
+
+                self.local = local if local is not None else 1
+
+                #weight coefficients for every vertex
+                self.weights = weights if weights is not None else float('NaN')
+
+                #string to identify this particular set of weights
+                self.weights_string = weights_string if weights_string is not None else ''
+
+                #do we cumulate misfit through time?
+                self.cumulated = cumulated if cumulated is not None else float('NaN')           
+                #}}}
+
+        def __repr__(self): # {{{
+                string='   Misfit:'
+
+                string="%s\n%s"%(string,fielddisplay(self,'name','identifier for this misfit response'))
+                string="%s\n%s"%(string,fielddisplay(self,'definitionstring','string that identifies this output definition uniquely, from "Outputdefinition[1-10]"'))
+                string="%s\n%s"%(string,fielddisplay(self,'model_string','string for field that is modeled'))
+                string="%s\n%s"%(string,fielddisplay(self,'observation','observed field that we compare the model against'))
+                string="%s\n%s"%(string,fielddisplay(self,'observation_string','observation string'))
+                string="%s\n%s"%(string,fielddisplay(self,'local','is the response local to the elements, or global? (default is 1)'))
+                string="%s\n%s"%(string,fielddisplay(self,'timeinterpolation','interpolation routine used to interpolate misfit between two time steps (default is "nearestneighbor"'))
+                string="%s\n%s"%(string,fielddisplay(self,'weights','weights (at vertices) to apply to the misfit'))
+                string="%s\n%s"%(string,fielddisplay(self,'weights_string','string for weights for identification purposes'))
+                return string
+                #}}}
+
+        def extrude(self,md): # {{{
+                if not np.any(np.isnan(self.weights)):
+                        self.weights = project3d(md,'vector',self.weights,'type','node')
+                if not np.any(np.isnan(self.observation)):
+                        self.observation = project3d(md,'vector',self.observation,'type','node')
+                return self
+        #}}}
+
+        def checkconsistency(self,md,solution,analyses):    # {{{
+                if type(self.name) != str:
+                        raise TypeError('misfit error message: "name" field should be a string!')
+
+                OutputdefinitionStringArray = []
+                for i in range(100):
+                        OutputdefinitionStringArray.append('Outputdefinition' + str(i))
+
+                md = checkfield(md,'fieldname','self.definitionstring','field',self.definitionstring,'values',OutputdefinitionStringArray)
+                if type(self.timeinterpolation) != str:
+                        raise TypeError('misfit error message: "timeinterpolation" field should be a string!')
+
+                md = checkfield(md,'fieldname','self.observation','field',self.observation,'timeseries',1,'NaN',1,'Inf',1)
+                md = checkfield(md,'fieldname','self.timeinterpolation','field',self.timeinterpolation,'values',['nearestneighbor'])
+                md = checkfield(md,'fieldname','self.weights','field',self.weights,'timeseries',1,'NaN',1,'Inf',1)
+
+                return md
+        # }}}
+
+        def marshall(self,prefix,md,fid):    #  {{{
+                WriteData(fid,prefix,'data',self.name,'name','md.misfit.name','format','String')
+                WriteData(fid,prefix,'data',self.definitionstring,'name','md.misfit.definitionstring','format','String')
+                WriteData(fid,prefix,'data',self.model_string,'name','md.misfit.model_string','format','String')
+                WriteData(fid,prefix,'data',self.observation,'name','md.misfit.observation','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts)
+                WriteData(fid,prefix,'data',self.observation_string,'name','md.misfit.observation_string','format','String')
+                WriteData(fid,prefix,'data',self.local,'name','md.misfit.local','format','Integer')
+                WriteData(fid,prefix,'data',self.timeinterpolation,'name','md.misfit.timeinterpolation','format','String')
+                WriteData(fid,prefix,'data',self.weights,'name','md.misfit.weights','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts)
+                WriteData(fid,prefix,'data',self.weights_string,'name','md.misfit.weights_string','format','String')
+        # }}}

../trunk-jpl/test/NightlyRun/test123.py

@@ +1 @@
+#Test Name: SquareShelfConstrainedTranMisfitSurface
+from model import *
+from socket import gethostname
+from triangle import *
+from setmask import *
+from parameterize import *
+from setflowequation import *
+from solve import *
+from misfit import *
+
+
+md=triangle(model(),'../Exp/Square.exp',180000)
+md=setmask(md,'all','')
+md=parameterize(md,'../Par/SquareShelfConstrained.py')
+md=setflowequation(md,'SSA','all')
+md.cluster=generic('name',gethostname(),'np',3)
+
+fake_surface = np.vstack((np.append(np.array(md.geometry.surface)+100,1.1),
+                np.append(np.array(md.geometry.surface)+200,2.1),
+                np.append(np.array(md.geometry.surface)+300,2.5))).T
+
+md.transient.requested_outputs = ['default','SurfaceMisfit']
+md.outputdefinition.definitions = [misfit(
+                name='SurfaceMisfit',
+                definitionstring='Outputdefinition1',
+                model_string='Surface',
+                observation=fake_surface,
+                observation_string='SurfaceObservation',
+                timeinterpolation='nearestneighbor',
+                weights=np.ones((md.mesh.numberofvertices,1)),
+                weights_string='WeightsSurfaceObservation'
+                )]
+
+md=solve(md,'Transient')
+
+#Fields and tolerances to track changes
+field_names = ['SurfaceMisfitFirstStep','SurfaceMisfitSecondStep','SurfaceMisfitThirdStep']
+field_tolerances = [1e-13,1e-13,1e-13]
+field_values = [
+        md.results.TransientSolution[0].SurfaceMisfit,
+        md.results.TransientSolution[1].SurfaceMisfit,
+        md.results.TransientSolution[2].SurfaceMisfit
+        ]

../trunk-jpl/test/NightlyRun/test124.py

@@ +1 @@
+#Test Name: SquareShelfConstrainedTranFSFreeSurface
+from model import *
+from socket import gethostname
+from triangle import *
+from setmask import *
+from parameterize import *
+from setflowequation import *
+from solve import *
+
+md=triangle(model(),'../Exp/Square.exp',150000)
+md=setmask(md,'all','')
+md=parameterize(md,'../Par/SquareShelfConstrained.py')
+md.extrude(3,1)
+md=setflowequation(md,'FS','all')
+
+#Free surface
+md.masstransport.isfreesurface = 1
+md.timestepping.time_step = 0.00001
+md.timestepping.final_time = 0.00005
+
+#Go solve
+md.cluster=generic('name',gethostname(),'np',3)
+md=solve(md,'Transient')
+
+#Fields and tolerances to track changes
+field_names      = [
+        'Vx1','Vy1','Vel1','Pressure1','Bed1','Surface1','Thickness1',
+        'Vx2','Vy2','Vel2','Pressure2','Bed2','Surface2','Thickness2',
+        'Vx3','Vy3','Vel3','Pressure3','Bed3','Surface3','Thickness3']
+field_tolerances = [
+        2e-09,3e-9,3e-9,3e-9,1e-13,1e-12,1e-12,
+        2e-09,3e-9,3e-9,3e-9,1e-10,1e-10,1e-10,
+        2e-09,3e-9,3e-9,3e-9,1e-10,1e-10,1e-10]
+field_values = [
+        md.results.TransientSolution[0].Vx,
+        md.results.TransientSolution[0].Vy,
+        md.results.TransientSolution[0].Vel,
+        md.results.TransientSolution[0].Pressure,
+        md.results.TransientSolution[0].Base,
+        md.results.TransientSolution[0].Surface,
+        md.results.TransientSolution[0].Thickness,
+        md.results.TransientSolution[1].Vx,
+        md.results.TransientSolution[1].Vy,
+        md.results.TransientSolution[1].Vel,
+        md.results.TransientSolution[1].Pressure,
+        md.results.TransientSolution[1].Base,
+        md.results.TransientSolution[1].Surface,
+        md.results.TransientSolution[1].Thickness,
+        md.results.TransientSolution[2].Vx,
+        md.results.TransientSolution[2].Vy,
+        md.results.TransientSolution[2].Vel,
+        md.results.TransientSolution[2].Pressure,
+        md.results.TransientSolution[2].Base,
+        md.results.TransientSolution[2].Surface,
+        md.results.TransientSolution[2].Thickness,
+        ]

../trunk-jpl/test/NightlyRun/runme.py

@@ -179 +179 @@
                                                 #compare to archive
                                                 # Matlab uses base 1, so use base 1 in labels
                                                 archive=np.array(archread(archive_file,archive_name+'_field'+str(k+1)))
-                                                if archive == None:
+                                                #Because np.array is weird (str(np.array(None)) becomes 'None' but np.array(None) is never equal to None, it basically becomes a type of string in an array):
+                                                if str(archive) == 'None':
                                                         raise NameError("Field name '"+archive_name+'_field'+str(k+1)+"' does not exist in archive file.")
                                                 error_diff=np.amax(np.abs(archive-field),axis=0)/(np.amax(np.abs(archive),axis=0)+float_info.epsilon)
                                                 if not np.isscalar(error_diff): error_diff=error_diff[0]