/[escript]/trunk/doc/examples/inversion/test_commemi1.py
ViewVC logotype

Contents of /trunk/doc/examples/inversion/test_commemi1.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 5775 - (show annotations)
Thu Jul 30 08:01:06 2015 UTC (3 years, 10 months ago) by sshaw
File MIME type: text/x-python
File size: 15971 byte(s)
pushing release to trunk
1 ##############################################################################
2 #
3 # Copyright (c) 2015 by The University of Queensland
4 # http://www.uq.edu.au
5 #
6 # Primary Business: Queensland, Australia
7 # Licensed under the Open Software License version 3.0
8 # http://www.opensource.org/licenses/osl-3.0.php
9 #
10 # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
11 # Development 2012-2013 by School of Earth Sciences
12 # Development from 2014 by Centre for Geoscience Computing (GeoComp)
13 #
14 ##############################################################################
15
16
17 """
18 Test script to run test model COMMEMI-1
19 """
20
21 from __future__ import print_function, division
22
23 import matplotlib
24 # The following line is here to allow automated testing. Remove or comment if
25 # you would like to display the final plot in a window instead.
26 matplotlib.use('agg')
27
28
29 import esys.downunder.magtel2d as mt2d
30 import numpy
31 import datetime
32 import esys.escript as escript
33 import esys.finley as finley
34 import esys.escript.pdetools as pdetools
35
36 #-------------------------------------------------------------
37 # The following functions create the mesh used by this example
38 #-------------------------------------------------------------
39
40
41 def makeLayerCake(x_start,x_extent,z_layers):
42 # --------------------------------------------------------------------------
43 # DESCRIPTION:
44 # -----------
45 # This is a utility function which sets up a 2D model with N layers.
46 #
47 # ARGUMENTS:
48 # ----------
49 # x_start :: start coordinate of mesh.
50 # x_extent :: horizontal extent of mesh.
51 # z_layers :: list with interface coordinates.
52 #
53 # RETURNS:
54 # --------
55 # borders :: borders of layers.
56 # air_earth_interface :: line at the air/earth interface.
57 #
58 # AUTHOR:
59 # -------
60 # Ralf Schaa,
61 # University of Queensland
62 #
63 #
64 # HISTORY:
65 # --------
66 #
67 # --------------------------------------------------------------------------
68
69 import esys.pycad as pycad # @UnresolvedImport
70 import esys.weipa as weipa # @UnresolvedImport
71 import esys.finley as finley # @UnresolvedImport
72 import esys.escript as escript # @UnresolvedImport
73
74
75 # --------------------------------------------------------------------------
76 # Point definitions.
77 # --------------------------------------------------------------------------
78
79 # Loop through all layers and define the vertices at all interfaces.
80 scale = 1.0
81 points = []
82 for i in range(0,len(z_layers)):
83 # Adjust scale at corners of air/earth interface:
84 if z_layers[i] == 0:
85 scale = 0.15
86 else:
87 scale = 1.0
88 points.append( pycad.Point(x_start , z_layers[i], 0.0, local_scale = scale) ) # Left-Corner.
89 points.append( pycad.Point(x_start + x_extent, z_layers[i], 0.0, local_scale = scale) ) # Right-Corner.
90
91
92 # --------------------------------------------------------------------------
93 # Line definitions.
94 # --------------------------------------------------------------------------
95
96 # Now connect the points to define the horizontal lines for all interfaces:
97 hlines = []
98 for i in range(0,len(points),2):
99 if i <= len(points)-1:
100 hlines.append( pycad.Line(points[i],points[i+1]) )
101
102 # Now connect the points to define the vertical lines for all interfaces:
103 vlines_left = []
104 for i in range(0,len(points),2):
105 if i <= len(points)-3:
106 vlines_left.append( pycad.Line(points[i],points[i+2]) )
107
108 vlines_right = []
109 for i in range(0,len(points),2):
110 if i <= len(points)-4:
111 vlines_right.append( pycad.Line(points[i+1],points[i+3]) )
112
113
114
115 # --------------------------------------------------------------------------
116 # Curveloop and Area definitions.
117 # --------------------------------------------------------------------------
118
119 # Join line segments for each layer.
120 borders = []
121 for i in range(0,len(z_layers)-1):
122 border = [ hlines[i],vlines_right[i],-hlines[i+1],-vlines_left[i] ]
123 borders.append( pycad.CurveLoop( border) )
124
125
126
127 # --------------------------------------------------------------------------
128 # Return values.
129 # --------------------------------------------------------------------------
130
131 # Explicitly specify the air-earth-boundary:
132 air_earth_interface = hlines[1]
133
134 return borders, air_earth_interface
135
136 #_______________________________________________________________________________
137
138
139
140
141 def setupMesh(mode, x_start, x_extent, a_extent, z_layers, anomaly_coord, elem_sizes):
142 # --------------------------------------------------------------------------
143 # DESCRIPTION:
144 # -----------
145 # This is a utility function which sets up the COMMEMI-1 mesh.
146 #
147 #
148 # ARGUMENTS:
149 # ----------
150 # mode :: TE or TM mode.
151 # x_start :: horizontal start-point mesh.
152 # x_extent :: horizontal extent of mesh.
153 # a_extent :: vertical extent of air-layer.
154 # z_layers :: list with coordinates of top-interfaces in Z-direction, incl. basement.
155 # anomaly_coord :: dictionary with coordinate tuples of anomalies, counterclockwise.
156 # elem_sizes :: mesh element sizes, large, normal, small.
157 #
158 # RETURNS:
159 # --------
160 # <Nothing> A mesh file is written to the output folder.
161 #
162 #
163 # AUTHOR:
164 # -------
165 # Ralf Schaa,
166 # The University of Queensland
167 #
168 #
169 # HISTORY:
170 # --------
171 #
172 # --------------------------------------------------------------------------
173
174
175
176 # --------------------------------------------------------------------------
177 # Imports.
178 # --------------------------------------------------------------------------
179
180 # System imports.
181 import math
182
183 # Escript modules.
184 import esys.pycad as pycad # @UnresolvedImport
185 import esys.finley as finley # @UnresolvedImport
186 import esys.escript as escript # @UnresolvedImport
187 import esys.weipa as weipa # @UnresolvedImport
188 # <Note>: "@UnresolvedImport" ignores any warnings in Eclipse/PyDev (PyDev has trouble with external libraries).
189
190 # Warn about magnetotelluric TM mode:
191 if mode.lower() == 'tm':
192 print("TM mode not yet supported")
193 return None
194
195 # --------------------------------------------------------------------------
196 # Anomaly border.
197 # --------------------------------------------------------------------------
198
199 #<Note>: define the anomaly which must be 'cut out' in the main mesh.
200
201
202 # Prepare list to store the anomaly borders:
203 border_anomaly = []
204
205 # Cycle anomaly dictionary and define the border for each.
206 for anomaly in anomaly_coord:
207
208 # Extract the coordinates for current key:
209 coord = anomaly_coord[anomaly]
210
211 # Points defining the anomaly from left-top.
212 points0 = []
213 for i in range( 0, len(coord) ):
214 points0.append(pycad.Point(coord[i][0], coord[i][1], 0.0))
215
216 # Define the line segments connecting the points.
217 lines0 = []
218 for i in range( 0, len(points0)-1 ):
219 lines0.append(pycad.Line(points0[i],points0[i+1]))
220 # Connect the last segment from end to start:
221 lines0.append(pycad.Line(points0[-1], points0[0]))
222
223 # And define the border of the anomalous area.
224 border_anomaly.append( pycad.CurveLoop(*lines0) )
225 #___________________________________________________________________________
226
227
228
229
230 # --------------------------------------------------------------------------
231 # Get the borders for each layer (air & host).
232 # --------------------------------------------------------------------------
233
234 # Borders around layers and the air/earth interface.
235 borders, air_earth_interface = makeLayerCake(x_start,x_extent,z_layers)
236
237
238
239
240
241 # --------------------------------------------------------------------------
242 # Specification of number of elements in domains.
243 # --------------------------------------------------------------------------
244
245 #<Note>: specifying the number of mesh elements is somewhat heuristic
246 # and is dependent on the mesh size and the anomaly sizes.
247
248 coord = anomaly_coord["anomaly_1"]
249
250 # First get the max-length of the anomaly to specify the number of elements.
251 length = max(( abs(coord[2][0]-coord[0][0]) ), # X-length
252 ( abs(coord[2][1]-coord[0][1]) )) # Y-length
253
254 # Specify number of elements in air, anomaly and on air/earth interface:
255 nr_elements_air = 1 * x_extent / elem_sizes["large"]
256 nr_elements_anomaly = 2 * length / elem_sizes["small"]
257 nr_elements_interface = 4 * x_extent / elem_sizes["small"]
258 #___________________________________________________________________________
259
260
261
262
263 # --------------------------------------------------------------------------
264 # Domain definitions.
265 # --------------------------------------------------------------------------
266
267 # Define the air & layer areas; note the 'holes' specifiers.
268 domain_air = pycad.PlaneSurface( borders[0] )
269 domain_host = pycad.PlaneSurface( borders[1] , holes = [ border_anomaly[0] ] )
270 domain_anomaly = pycad.PlaneSurface( border_anomaly[0] )
271
272 # Specify the element sizes in the domains and along the interface.
273 #<Note>: Sizes must be assigned in the order as they appear below:
274 domain_air.setElementDistribution( nr_elements_air )
275 domain_anomaly.setElementDistribution( nr_elements_anomaly )
276 air_earth_interface.setElementDistribution( nr_elements_interface )
277
278 # Ready to define the mesh-design..
279 design2D = pycad.gmsh.Design(dim=2, element_size=elem_sizes["normal"] , keep_files=False)
280 # ..and also specify the domains for tagging with property values later on:
281 design2D.addItems( pycad.PropertySet("domain_air" , domain_air),
282 pycad.PropertySet("domain_host" , domain_host),
283 pycad.PropertySet("domain_anomaly", domain_anomaly) )
284
285 # Now define the unstructured finley-mesh..
286 model2D = finley.MakeDomain(design2D)
287 #___________________________________________________________________________
288
289
290 return model2D
291 #___________________________________________________________________________
292
293 def generateCommemi1Mesh():
294 # --------------------------------------------------------------------------
295 # Geometric mesh parameters.
296 # --------------------------------------------------------------------------
297
298 # Mesh extents.
299 a_extent = 20000 # 20km - Vertical extent of air-layer in (m).
300 z_extent = 20000 # 20km - Vertical extent of subsurface in (m).
301 x_extent = 40000 # 40km - Horizontal extent of mesh in (m).
302
303 # Start point of mesh.
304 x_start = 0 #-x_extent/2.0
305
306 # Define interface locations in z-direction: top, air/earth, basement.
307 z_layers = [ a_extent, 0, -z_extent]
308
309 # Mesh elements sizes.
310 elem_sizes = {
311 'large' : 10.00 * x_extent/100.0, # 5.00% of x_extent.
312 'normal': 05.00 * x_extent/100.0, # 2.50% of x_extent.
313 'small' : 00.50 * x_extent/100.0 # 0.25% of x_extent.
314 }
315 #___________________________________________________________________________
316
317
318
319
320 # --------------------------------------------------------------------------
321 # Geometric anomaly parameters.
322 # --------------------------------------------------------------------------
323
324 # Extents of the rectangular 2D anomaly.
325 x_anomaly = 1000 # 1km - Horizontal extent of anomaly in (m).
326 z_anomaly = 2000 # 2km - Vertical extent of anomaly in (m).
327
328 # Coordinates of the rectangular 2D anomaly.
329 ya1 = -250 # Top
330 ya2 = -z_anomaly + ya1 # Bottom
331 xa1 = x_start + x_extent/2.0 - x_anomaly/2.0 # Left
332 xa2 = x_start + x_extent/2.0 + x_anomaly/2.0 # Right
333
334 # Save in dictionary as a list of tuples from left-top corner, counterclockwise.
335 anomaly_coord = {
336 'anomaly_1': ([xa1,ya1],[xa1,ya2],[xa2,ya2],[xa2,ya1])
337 }
338 #___________________________________________________________________________
339
340
341
342
343 # --------------------------------------------------------------------------
344 # Setup the COMMEMI-1 mesh.
345 # --------------------------------------------------------------------------
346
347 # This creates the mesh and saves it to the output folder.
348 return setupMesh("TE", x_start, x_extent, a_extent, z_layers, anomaly_coord, elem_sizes)
349 #___________________________________________________________________________
350
351
352
353 #-------------------------------------------------------------
354 # End of mesh set up functions
355 #-------------------------------------------------------------
356
357
358
359 # ---
360 # Initialisations
361 # ---
362
363 # Get timing:
364 startTime = datetime.datetime.now()
365
366 # Mode (TE includes air-layer, whereas TM does not):
367 mode = 'TE'
368
369 # Read the mesh file and define the 'finley' domain:
370 #mesh_file = "data/commemi1_te.fly"
371 #domain = finley.ReadMesh(mesh_file, numDim=2)
372 if escript.getEscriptParamInt('GMSH_SUPPORT'):
373 domain = generateCommemi1Mesh()
374
375 # Sounding frequencies (in Hz):
376 freq_def = {"high":1.0e+1,"low":1.0e+1,"step":1}
377 # Frequencies will be mapped on a log-scale from
378 # 'high' to 'low' with 'step' points per decade.
379 # (also only one frequency must be passed via dict)
380
381 # Step sizes for sampling along vertical and horizontal axis (in m):
382 xstep=400
383 zstep=200
384
385
386
387 # ---
388 # Resistivity model
389 # ---
390
391 # Resistivity values assigned to tagged regions (in Ohm.m):
392 rho = [
393 1.0e+14, # 0: air
394 100.0 , # 1: host
395 0.5 # 2: anomaly
396 ]
397
398 # Tags must match those in the file:
399 tags = ["domain_air", "domain_host", "domain_anomaly"]
400
401
402 # ---
403 # Layer definitions for 1D response at boundaries.
404 # ---
405
406 # List with resistivity values for left and right boundary.
407 rho_1d_left = [ rho[0], rho[1] ]
408 rho_1d_rght = [ rho[0], rho[1] ]
409
410 # Associated interfaces for 1D response left and right (must match the mesh file).
411 ifc_1d_left = [ 20000, 0, -20000]
412 ifc_1d_rght = [ 20000, 0, -20000]
413
414 # Save in dictionary with layer interfaces and resistivities left and right:
415 ifc_1d = {"left":ifc_1d_left , "right":ifc_1d_rght}
416 rho_1d = {"left":rho_1d_left , "right":rho_1d_rght}
417
418
419
420 # ---
421 # Run MT_2D
422 # ---
423
424 # Class options:
425 mt2d.MT_2D._solver = "DIRECT"
426 mt2d.MT_2D._debug = False
427
428 if mt2d.MT_2D._solver == "DIRECT" and escript.getMPISizeWorld() > 1:
429 print("Direct solvers and multiple MPI processes are not currently supported")
430 elif mt2d.MT_2D._solver == "DIRECT" and not escript.getEscriptParamInt('PASO_DIRECT'):
431 print("escript was not built with support for direct solvers, aborting")
432 elif not escript.getEscriptParamInt('GMSH_SUPPORT'):
433 print("This example requires gmsh")
434 else:
435 # Instantiate an MT_2D object with required & optional parameters:
436 obj_mt2d = mt2d.MT_2D(domain, mode, freq_def, tags, rho, rho_1d, ifc_1d,
437 xstep=xstep ,zstep=zstep, maps=None, plot=True)
438
439 # Solve for fields, apparent resistivity and phase:
440 mt2d_fields, arho_2d, aphi_2d = obj_mt2d.pdeSolve()
441
442
443 #
444 print(datetime.datetime.now()-startTime)
445
446
447 print("Done!")
448
449
450

  ViewVC Help
Powered by ViewVC 1.1.26