Index: /issm/trunk-jpl/src/m/classes/misfit.py =================================================================== --- /issm/trunk-jpl/src/m/classes/misfit.py (revision 22104) +++ /issm/trunk-jpl/src/m/classes/misfit.py (revision 22104) @@ -0,0 +1,109 @@ +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') + # }}} Index: /issm/trunk-jpl/test/NightlyRun/runme.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/runme.py (revision 22103) +++ /issm/trunk-jpl/test/NightlyRun/runme.py (revision 22104) @@ -180,5 +180,6 @@ # 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) Index: /issm/trunk-jpl/test/NightlyRun/test123.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test123.py (revision 22104) +++ /issm/trunk-jpl/test/NightlyRun/test123.py (revision 22104) @@ -0,0 +1,43 @@ +#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 + ] Index: /issm/trunk-jpl/test/NightlyRun/test124.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test124.py (revision 22104) +++ /issm/trunk-jpl/test/NightlyRun/test124.py (revision 22104) @@ -0,0 +1,56 @@ +#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, + ]