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

Last change on this file since 22755 was 22755, checked in by Mathieu Morlighem, 7 years ago

CHG: added 21724-22754
File size: 10.4 KB
RevLine 
[22755]1Index: ../trunk-jpl/src/m/classes/misfit.py
2===================================================================
3--- ../trunk-jpl/src/m/classes/misfit.py (nonexistent)
4+++ ../trunk-jpl/src/m/classes/misfit.py (revision 22104)
5@@ -0,0 +1,109 @@
6+import numpy as np
7+from project3d import project3d
8+from pairoptions import *
9+from collections import OrderedDict
10+from fielddisplay import fielddisplay
11+from checkfield import checkfield
12+from WriteData import WriteData
13+
14+
15+class misfit(object):
16+ """
17+ MISFIT class definition
18+
19+ Usage:
20+ misfit=misfit()
21+ misfit=misfit(name='SurfaceAltimetry',
22+ definitionstring='Outputdefinition1',
23+ model_string='Surface',
24+ observation_string='SurfaceObservations',
25+ observation=md.geometry.surface,
26+ timeinterpolation='nearestneighbor',
27+ local=1,
28+ weights=np.ones((md.mesh.numberofvertices,1)),
29+ weights_string='WeightsSurfaceObservations')
30+ """
31+
32+ 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):
33+ # {{{
34+ self.name = name if name is not None else ''
35+
36+ #string that identifies this output definition uniquely, from 'Outputdefinition[1-100]'
37+ self.definitionstring = definitionstring if definitionstring is not None else ''
38+
39+ #string for field that is modeled
40+ self.model_string = model_string if model_string is not None else ''
41+
42+ #observed field that we compare the model against
43+ self.observation = observation if observation is not None else float('NaN')
44+
45+ #string for observed field.
46+ self.observation_string = observation_string if observation_string is not None else ''
47+
48+ self.timeinterpolation = timeinterpolation if timeinterpolation is not None else 'nearestneighbor'
49+
50+ self.local = local if local is not None else 1
51+
52+ #weight coefficients for every vertex
53+ self.weights = weights if weights is not None else float('NaN')
54+
55+ #string to identify this particular set of weights
56+ self.weights_string = weights_string if weights_string is not None else ''
57+
58+ #do we cumulate misfit through time?
59+ self.cumulated = cumulated if cumulated is not None else float('NaN')
60+ #}}}
61+
62+ def __repr__(self): # {{{
63+ string=' Misfit:'
64+
65+ string="%s\n%s"%(string,fielddisplay(self,'name','identifier for this misfit response'))
66+ string="%s\n%s"%(string,fielddisplay(self,'definitionstring','string that identifies this output definition uniquely, from "Outputdefinition[1-10]"'))
67+ string="%s\n%s"%(string,fielddisplay(self,'model_string','string for field that is modeled'))
68+ string="%s\n%s"%(string,fielddisplay(self,'observation','observed field that we compare the model against'))
69+ string="%s\n%s"%(string,fielddisplay(self,'observation_string','observation string'))
70+ string="%s\n%s"%(string,fielddisplay(self,'local','is the response local to the elements, or global? (default is 1)'))
71+ string="%s\n%s"%(string,fielddisplay(self,'timeinterpolation','interpolation routine used to interpolate misfit between two time steps (default is "nearestneighbor"'))
72+ string="%s\n%s"%(string,fielddisplay(self,'weights','weights (at vertices) to apply to the misfit'))
73+ string="%s\n%s"%(string,fielddisplay(self,'weights_string','string for weights for identification purposes'))
74+ return string
75+ #}}}
76+
77+ def extrude(self,md): # {{{
78+ if not np.any(np.isnan(self.weights)):
79+ self.weights = project3d(md,'vector',self.weights,'type','node')
80+ if not np.any(np.isnan(self.observation)):
81+ self.observation = project3d(md,'vector',self.observation,'type','node')
82+ return self
83+ #}}}
84+
85+ def checkconsistency(self,md,solution,analyses): # {{{
86+ if type(self.name) != str:
87+ raise TypeError('misfit error message: "name" field should be a string!')
88+
89+ OutputdefinitionStringArray = []
90+ for i in range(100):
91+ OutputdefinitionStringArray.append('Outputdefinition' + str(i))
92+
93+ md = checkfield(md,'fieldname','self.definitionstring','field',self.definitionstring,'values',OutputdefinitionStringArray)
94+ if type(self.timeinterpolation) != str:
95+ raise TypeError('misfit error message: "timeinterpolation" field should be a string!')
96+
97+ md = checkfield(md,'fieldname','self.observation','field',self.observation,'timeseries',1,'NaN',1,'Inf',1)
98+ md = checkfield(md,'fieldname','self.timeinterpolation','field',self.timeinterpolation,'values',['nearestneighbor'])
99+ md = checkfield(md,'fieldname','self.weights','field',self.weights,'timeseries',1,'NaN',1,'Inf',1)
100+
101+ return md
102+ # }}}
103+
104+ def marshall(self,prefix,md,fid): # {{{
105+ WriteData(fid,prefix,'data',self.name,'name','md.misfit.name','format','String')
106+ WriteData(fid,prefix,'data',self.definitionstring,'name','md.misfit.definitionstring','format','String')
107+ WriteData(fid,prefix,'data',self.model_string,'name','md.misfit.model_string','format','String')
108+ WriteData(fid,prefix,'data',self.observation,'name','md.misfit.observation','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts)
109+ WriteData(fid,prefix,'data',self.observation_string,'name','md.misfit.observation_string','format','String')
110+ WriteData(fid,prefix,'data',self.local,'name','md.misfit.local','format','Integer')
111+ WriteData(fid,prefix,'data',self.timeinterpolation,'name','md.misfit.timeinterpolation','format','String')
112+ WriteData(fid,prefix,'data',self.weights,'name','md.misfit.weights','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts)
113+ WriteData(fid,prefix,'data',self.weights_string,'name','md.misfit.weights_string','format','String')
114+ # }}}
115Index: ../trunk-jpl/test/NightlyRun/test123.py
116===================================================================
117--- ../trunk-jpl/test/NightlyRun/test123.py (nonexistent)
118+++ ../trunk-jpl/test/NightlyRun/test123.py (revision 22104)
119@@ -0,0 +1,43 @@
120+#Test Name: SquareShelfConstrainedTranMisfitSurface
121+from model import *
122+from socket import gethostname
123+from triangle import *
124+from setmask import *
125+from parameterize import *
126+from setflowequation import *
127+from solve import *
128+from misfit import *
129+
130+
131+md=triangle(model(),'../Exp/Square.exp',180000)
132+md=setmask(md,'all','')
133+md=parameterize(md,'../Par/SquareShelfConstrained.py')
134+md=setflowequation(md,'SSA','all')
135+md.cluster=generic('name',gethostname(),'np',3)
136+
137+fake_surface = np.vstack((np.append(np.array(md.geometry.surface)+100,1.1),
138+ np.append(np.array(md.geometry.surface)+200,2.1),
139+ np.append(np.array(md.geometry.surface)+300,2.5))).T
140+
141+md.transient.requested_outputs = ['default','SurfaceMisfit']
142+md.outputdefinition.definitions = [misfit(
143+ name='SurfaceMisfit',
144+ definitionstring='Outputdefinition1',
145+ model_string='Surface',
146+ observation=fake_surface,
147+ observation_string='SurfaceObservation',
148+ timeinterpolation='nearestneighbor',
149+ weights=np.ones((md.mesh.numberofvertices,1)),
150+ weights_string='WeightsSurfaceObservation'
151+ )]
152+
153+md=solve(md,'Transient')
154+
155+#Fields and tolerances to track changes
156+field_names = ['SurfaceMisfitFirstStep','SurfaceMisfitSecondStep','SurfaceMisfitThirdStep']
157+field_tolerances = [1e-13,1e-13,1e-13]
158+field_values = [
159+ md.results.TransientSolution[0].SurfaceMisfit,
160+ md.results.TransientSolution[1].SurfaceMisfit,
161+ md.results.TransientSolution[2].SurfaceMisfit
162+ ]
163Index: ../trunk-jpl/test/NightlyRun/test124.py
164===================================================================
165--- ../trunk-jpl/test/NightlyRun/test124.py (nonexistent)
166+++ ../trunk-jpl/test/NightlyRun/test124.py (revision 22104)
167@@ -0,0 +1,56 @@
168+#Test Name: SquareShelfConstrainedTranFSFreeSurface
169+from model import *
170+from socket import gethostname
171+from triangle import *
172+from setmask import *
173+from parameterize import *
174+from setflowequation import *
175+from solve import *
176+
177+md=triangle(model(),'../Exp/Square.exp',150000)
178+md=setmask(md,'all','')
179+md=parameterize(md,'../Par/SquareShelfConstrained.py')
180+md.extrude(3,1)
181+md=setflowequation(md,'FS','all')
182+
183+#Free surface
184+md.masstransport.isfreesurface = 1
185+md.timestepping.time_step = 0.00001
186+md.timestepping.final_time = 0.00005
187+
188+#Go solve
189+md.cluster=generic('name',gethostname(),'np',3)
190+md=solve(md,'Transient')
191+
192+#Fields and tolerances to track changes
193+field_names = [
194+ 'Vx1','Vy1','Vel1','Pressure1','Bed1','Surface1','Thickness1',
195+ 'Vx2','Vy2','Vel2','Pressure2','Bed2','Surface2','Thickness2',
196+ 'Vx3','Vy3','Vel3','Pressure3','Bed3','Surface3','Thickness3']
197+field_tolerances = [
198+ 2e-09,3e-9,3e-9,3e-9,1e-13,1e-12,1e-12,
199+ 2e-09,3e-9,3e-9,3e-9,1e-10,1e-10,1e-10,
200+ 2e-09,3e-9,3e-9,3e-9,1e-10,1e-10,1e-10]
201+field_values = [
202+ md.results.TransientSolution[0].Vx,
203+ md.results.TransientSolution[0].Vy,
204+ md.results.TransientSolution[0].Vel,
205+ md.results.TransientSolution[0].Pressure,
206+ md.results.TransientSolution[0].Base,
207+ md.results.TransientSolution[0].Surface,
208+ md.results.TransientSolution[0].Thickness,
209+ md.results.TransientSolution[1].Vx,
210+ md.results.TransientSolution[1].Vy,
211+ md.results.TransientSolution[1].Vel,
212+ md.results.TransientSolution[1].Pressure,
213+ md.results.TransientSolution[1].Base,
214+ md.results.TransientSolution[1].Surface,
215+ md.results.TransientSolution[1].Thickness,
216+ md.results.TransientSolution[2].Vx,
217+ md.results.TransientSolution[2].Vy,
218+ md.results.TransientSolution[2].Vel,
219+ md.results.TransientSolution[2].Pressure,
220+ md.results.TransientSolution[2].Base,
221+ md.results.TransientSolution[2].Surface,
222+ md.results.TransientSolution[2].Thickness,
223+ ]
224Index: ../trunk-jpl/test/NightlyRun/runme.py
225===================================================================
226--- ../trunk-jpl/test/NightlyRun/runme.py (revision 22103)
227+++ ../trunk-jpl/test/NightlyRun/runme.py (revision 22104)
228@@ -179,7 +179,8 @@
229 #compare to archive
230 # Matlab uses base 1, so use base 1 in labels
231 archive=np.array(archread(archive_file,archive_name+'_field'+str(k+1)))
232- if archive == None:
233+ #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):
234+ if str(archive) == 'None':
235 raise NameError("Field name '"+archive_name+'_field'+str(k+1)+"' does not exist in archive file.")
236 error_diff=np.amax(np.abs(archive-field),axis=0)/(np.amax(np.abs(archive),axis=0)+float_info.epsilon)
237 if not np.isscalar(error_diff): error_diff=error_diff[0]
Note: See TracBrowser for help on using the repository browser.