source: issm/oecreview/Archive/21724-22754/ISSM-22270-22271.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: 8.4 KB
RevLine 
[22755]1Index: ../trunk-jpl/test/NightlyRun/test703.py
2===================================================================
3--- ../trunk-jpl/test/NightlyRun/test703.py (revision 22270)
4+++ ../trunk-jpl/test/NightlyRun/test703.py (nonexistent)
5@@ -1,160 +0,0 @@
6-#Test Name: FlowbandFSsheetshelfTrans
7-import numpy as np
8-from scipy.interpolate import interp1d
9-import copy
10-from model import *
11-from socket import gethostname
12-from triangle import *
13-from meshconvert import *
14-from setmask import *
15-from parameterize import *
16-from setflowequation import *
17-from solve import *
18-from NowickiProfile import *
19-from bamgflowband import *
20-from paterson import *
21-
22-#mesh parameters
23-x = np.arange(-5,5.5,.5).T
24-[b,h,sea] = NowickiProfile(x)
25-x = x * 10**3
26-h = h * 10**3
27-b = (b-sea) * 10**3
28-
29-#mesh domain
30-md = bamgflowband(model(),x,b+h,b,'hmax',150)
31-print md.mesh.numberofvertices
32-
33-#geometry
34-md.geometry.surface = interp1d(x,b+h)(md.mesh.x)
35-md.geometry.base = interp1d(x,b)(md.mesh.x)
36-md.geometry.thickness = md.geometry.surface - md.geometry.base
37-
38-#mask
39-md.mask.ice_levelset = -np.ones((md.mesh.numberofvertices,))
40-md.mask.ice_levelset[np.where(md.mesh.vertexflags(2))] = 0
41-md.mask.groundedice_levelset = np.zeros((md.mesh.numberofvertices,)) - 0.5
42-
43-#materials
44-md.initialization.temperature = (273. - 20.) * np.ones((md.mesh.numberofvertices,))
45-md.materials.rheology_B = paterson(md.initialization.temperature)
46-md.materials.rheology_n = 3 * np.ones((md.mesh.numberofelements,))
47-
48-#friction
49-md.friction.coefficient = np.zeros((md.mesh.numberofvertices,))
50-md.friction.coefficient[np.where(md.mesh.vertexflags(1))] = 20
51-md.friction.p = np.ones((md.mesh.numberofelements,))
52-md.friction.q = np.ones((md.mesh.numberofelements,))
53-
54-#boundary conditions
55-md.stressbalance.spcvx = float('NaN') * np.ones((md.mesh.numberofvertices,))
56-md.stressbalance.spcvy = float('NaN') * np.ones((md.mesh.numberofvertices,))
57-md.stressbalance.spcvz = float('NaN') * np.ones((md.mesh.numberofvertices,))
58-md.stressbalance.referential = float('NaN') * np.ones((md.mesh.numberofvertices,6))
59-md.stressbalance.loadingforce = 0 * np.ones((md.mesh.numberofvertices,3))
60-md.stressbalance.spcvx[np.where(md.mesh.vertexflags(4))] = 800.
61-md.stressbalance.spcvy[np.where(md.mesh.vertexflags(4))] = 0.
62-
63-#print md.stressbalance.referential
64-#print md.stressbalance.loadingforce
65-#print md.stressbalance.spcvx
66-#print md.stressbalance.spcvy
67-#print md.stressbalance.spcvz
68-#print np.shape(md.stressbalance.referential)
69-#print np.shape(md.stressbalance.loadingforce)
70-#print np.shape(md.stressbalance.spcvx)
71-#print np.shape(md.stressbalance.spcvy)
72-#print np.shape(md.stressbalance.spcvz)
73-
74-#Misc
75-md = setflowequation(md,'FS','all')
76-md.flowequation.fe_FS = 'TaylorHood'
77-md.stressbalance.abstol = float('NaN')
78-md.miscellaneous.name = 'test703'
79-
80-#Transient settings
81-md.timestepping.time_step = 0.000001
82-md.timestepping.final_time = 0.000005
83-md.stressbalance.shelf_dampening = 1.
84-md.smb.mass_balance = np.zeros((md.mesh.numberofvertices,))
85-md.basalforcings.groundedice_melting_rate = np.zeros((md.mesh.numberofvertices,))
86-md.basalforcings.floatingice_melting_rate = np.zeros((md.mesh.numberofvertices,))
87-md.basalforcings.geothermalflux = np.zeros((md.mesh.numberofvertices,))
88-posb = np.intersect1d(np.where(md.mesh.x > 0.), np.where(md.mesh.vertexonbase))
89-md.basalforcings.groundedice_melting_rate[posb] = 18.
90-md.basalforcings.floatingice_melting_rate[posb] = 18.
91-md.initialization.vx = np.zeros((md.mesh.numberofvertices,))
92-md.initialization.vy = np.zeros((md.mesh.numberofvertices,))
93-md.initialization.pressure = np.zeros((md.mesh.numberofvertices,))
94-md.masstransport.spcthickness = float('NaN') * np.ones((md.mesh.numberofvertices,))
95-md.thermal.spctemperature = float('NaN') * np.ones((md.mesh.numberofvertices,))
96-md.transient.isthermal = 0
97-md.masstransport.isfreesurface = 1
98-
99-#Go solve
100-md.cluster = generic('np',3)
101-md.stressbalance.shelf_dampening = 1
102-md1 = solve(md,'Transient')
103-
104-md.stressbalance.shelf_dampening = 0
105-md = solve(md,'Transient')
106-
107-#Fields and tolerances to track changes
108-field_names = [
109- 'Vx1','Vy1','Vel1','Pressure1','Bed1','Surface1','Thickness1',
110- 'Vx2','Vy2','Vel2','Pressure2','Bed2','Surface2','Thickness2',
111- 'Vx3','Vy3','Vel3','Pressure3','Bed3','Surface3','Thickness3',
112- 'Vx1_damp','Vy1_damp','Vel1_damp','Pressure1_damp','Bed1_damp','Surface1_damp','Thickness1_damp',
113- 'Vx2_damp','Vy2_damp','Vel2_damp','Pressure2_damp','Bed2_damp','Surface2_damp','Thickness2_damp',
114- 'Vx3_damp','Vy3_damp','Vel3_damp','Pressure3_damp','Bed3_damp','Surface3_damp','Thickness3_damp']
115-field_tolerances = [
116- 2e-08,2e-08,2e-08,1e-08,1e-10,1e-10,1e-10,
117- 2e-08,2e-08,2e-08,1e-08,1e-10,1e-10,1e-10,
118- 2e-08,2e-08,2e-08,1e-08,1e-10,1e-10,1e-10,
119- 5e-08,5e-08,5e-08,1e-08,1e-10,1e-10,1e-10,
120- 5e-08,5e-08,5e-08,1e-08,1e-10,1e-10,1e-10,
121- 5e-08,5e-08,5e-08,1e-08,1e-10,1e-10,1e-10]
122-field_values = [
123- md.results.TransientSolution[0].Vx,
124- md.results.TransientSolution[0].Vy,
125- md.results.TransientSolution[0].Vel,
126- md.results.TransientSolution[0].Pressure,
127- md.results.TransientSolution[0].Base,
128- md.results.TransientSolution[0].Surface,
129- md.results.TransientSolution[0].Thickness,
130- md.results.TransientSolution[1].Vx,
131- md.results.TransientSolution[1].Vy,
132- md.results.TransientSolution[1].Vel,
133- md.results.TransientSolution[1].Pressure,
134- md.results.TransientSolution[1].Base,
135- md.results.TransientSolution[1].Surface,
136- md.results.TransientSolution[1].Thickness,
137- md.results.TransientSolution[2].Vx,
138- md.results.TransientSolution[2].Vy,
139- md.results.TransientSolution[2].Vel,
140- md.results.TransientSolution[2].Pressure,
141- md.results.TransientSolution[2].Base,
142- md.results.TransientSolution[2].Surface,
143- md.results.TransientSolution[2].Thickness,
144- md1.results.TransientSolution[0].Vx,
145- md1.results.TransientSolution[0].Vy,
146- md1.results.TransientSolution[0].Vel,
147- md1.results.TransientSolution[0].Pressure,
148- md1.results.TransientSolution[0].Base,
149- md1.results.TransientSolution[0].Surface,
150- md1.results.TransientSolution[0].Thickness,
151- md1.results.TransientSolution[1].Vx,
152- md1.results.TransientSolution[1].Vy,
153- md1.results.TransientSolution[1].Vel,
154- md1.results.TransientSolution[1].Pressure,
155- md1.results.TransientSolution[1].Base,
156- md1.results.TransientSolution[1].Surface,
157- md1.results.TransientSolution[1].Thickness,
158- md1.results.TransientSolution[2].Vx,
159- md1.results.TransientSolution[2].Vy,
160- md1.results.TransientSolution[2].Vel,
161- md1.results.TransientSolution[2].Pressure,
162- md1.results.TransientSolution[2].Base,
163- md1.results.TransientSolution[2].Surface,
164- md1.results.TransientSolution[2].Thickness,
165- ]
166Index: ../trunk-jpl/src/m/mesh/bamg.py
167===================================================================
168--- ../trunk-jpl/src/m/mesh/bamg.py (revision 22270)
169+++ ../trunk-jpl/src/m/mesh/bamg.py (revision 22271)
170@@ -1,6 +1,8 @@
171 import os.path
172 import numpy as np
173+from mesh2d import *
174 from mesh2dvertical import *
175+from mesh3dsurface import *
176 from collections import OrderedDict
177 from pairoptions import pairoptions
178 from bamggeom import bamggeom
179@@ -91,18 +93,18 @@
180 for i,domaini in enumerate(domain):
181
182 #Check that the domain is closed
183- if (domaini.x[0]!=domaini.x[-1] or domaini.y[0]!=domaini.y[-1]):
184+ if (domaini['x'][0]!=domaini['x'][-1] or domaini['y'][0]!=domaini['y'][-1]):
185 raise RuntimeError("bamg error message: all contours provided in ''domain'' should be closed")
186
187 #Checks that all holes are INSIDE the principle domain outline
188 if i:
189- flags=ContourToNodes(domaini.x,domaini.y,domainfile,0)[0]
190+ flags=ContourToNodes(domaini['x'],domaini['y'],domainfile,0)[0]
191 if np.any(np.logical_not(flags)):
192 raise RuntimeError("bamg error message: All holes should be strictly inside the principal domain")
193
194 #Add all points to bamg_geometry
195- nods=domaini.nods-1 #the domain are closed 0=end
196- bamg_geometry.Vertices=np.vstack((bamg_geometry.Vertices,np.vstack((domaini.x[0:nods],domaini.y[0:nods],np.ones((nods)))).T))
197+ nods=domaini['nods']-1 #the domain are closed 0=end
198+ bamg_geometry.Vertices=np.vstack((bamg_geometry.Vertices,np.vstack((domaini['x'][0:nods],domaini['y'][0:nods],np.ones((nods)))).T))
199 bamg_geometry.Edges =np.vstack((bamg_geometry.Edges,np.vstack((np.arange(count+1,count+nods+1),np.hstack((np.arange(count+2,count+nods+1),count+1)),1.*np.ones((nods)))).T))
200 if i:
201 bamg_geometry.SubDomains=np.vstack((bamg_geometry.SubDomains,[2,count+1,1,-subdomain_ref]))
Note: See TracBrowser for help on using the repository browser.