/[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 6651 - (show annotations)
Wed Feb 7 02:12:08 2018 UTC (21 months, 1 week ago) by jfenwick
File MIME type: text/x-python
File size: 16565 byte(s)
Make everyone sad by touching all the files

Copyright dates update

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

  ViewVC Help
Powered by ViewVC 1.1.26