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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6651 - (show annotations)
Wed Feb 7 02:12:08 2018 UTC (2 years, 3 months ago) by jfenwick
File MIME type: text/x-python
File size: 22167 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 Test script to run test model COMMEMI-4
18 """
19
20 from __future__ import print_function, division
21
22 import matplotlib
23 # The following line is here to allow automated testing. Remove or comment if
24 # you would like to display the final plot in a window instead.
25 matplotlib.use('agg')
26
27 import datetime
28 import numpy
29
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 def setupMesh(mode, coord, elem_sizes):
43 #---------------------------------------------------------------------------
44 # DESCRIPTION:
45 # -----------
46 # This is a utility function which setups the COMMEMI-4 mesh.
47 #
48 #
49 # ARGUMENTS:
50 # ----------
51 # mode :: TE or TM mode.
52 # coord :: dictionary with coordinate tuples.
53 # elem_sizes :: mesh element sizes, large, normal, small.
54 #
55 # RETURNS:
56 # --------
57 # <Nothing> A mesh file is written to the output folder.
58 #
59 #
60 # AUTHOR:
61 # -------
62 # Ralf Schaa,
63 # University of Queensland
64 #
65 #---------------------------------------------------------------------------
66
67
68
69 #---------------------------------------------------------------------------
70 # Imports.
71 #---------------------------------------------------------------------------
72
73 import esys.pycad as pycad # @UnresolvedImport
74 import esys.finley as finley # @UnresolvedImport
75 import esys.escript as escript # @UnresolvedImport
76 import esys.weipa as weipa # @UnresolvedImport
77 # <Note>: "@UnresolvedImport" ignores any warnings in Eclipse/PyDev (PyDev has trouble with external libraries).
78
79
80
81 model = "COMMEMI-4"
82
83 print("Preparing the mesh " + model + " ...")
84 print("")
85
86 # Warn about magnetotelluric TM mode:
87 if mode.lower() == 'tm':
88 print("TM mode not yet supported")
89 return
90
91
92
93 # Path to write the mesh:
94 outpath = "../out/commemi4"
95
96
97
98
99 # --------------------------------------------------------------------------
100 # Initialisations.
101 # --------------------------------------------------------------------------
102
103 # Get coordinates from dictionary as list of tuples
104 a0 = coord["air"]
105 l1 = coord["lyr1"]
106 s1 = coord["slab"]
107 b1 = coord["basin"]
108 l2 = coord["lyr2"]
109 l3 = coord["lyr3"]
110
111 # Mesh length from top-boundary.
112 x_extent = abs(a0[3][0]-a0[0][0])
113
114
115
116
117 # --------------------------------------------------------------------------
118 # Point definitions.
119 # --------------------------------------------------------------------------
120
121 #<Note>: define all points spanning the mesh, anomalies and layers;
122 # note also shared domain points must be defined only once.
123
124
125 # Mesh top boundary.
126 air = []
127 air.append( pycad.Point( *a0[0] ) ) # 0: left , top (@ boundary)
128 air.append( pycad.Point( *a0[3] ) ) # 3: right , top (@ boundary)
129
130
131 # First-layer.
132 ly1 = []
133 ly1.append( pycad.Point( *l1[0] ) ) # 0: left , top (@ air/earth interface)
134 ly1.append( pycad.Point( *l1[1] ) ) # 1: left , bottom (@ boundary)
135 ly1.append( pycad.Point( *l1[2] ) ) # 2: right , bottom (@ slab/basin)
136 ly1.append( pycad.Point( *l1[3] ) ) # 3: right , bottom (@ boundary)
137 ly1.append( pycad.Point( *l1[4] ) ) # 4: right , top (@ air/earth interface)
138
139
140 # Slab.
141 sl1 = []
142 sl1.append( ly1[1] ) # 0: left , top (@ boundary)
143 sl1.append( pycad.Point( *s1[1] ) ) # 1: left , bottom (@ boundary)
144 sl1.append( pycad.Point( *s1[2] ) ) # 2: right , bottom (@ slab/basin)
145 sl1.append( ly1[2] ) # 3: right , top (@ slab/basin)
146
147
148 # Basin.
149 bs1 = []
150 bs1.append( ly1[2] ) # 0: left , top (@ slab/basin)
151 bs1.append( sl1[2] ) # 1: left , centre (@ slab/basin)
152 bs1.append( pycad.Point( *b1[2] ) ) # 2: left , bottom (@ lyr1/basin)
153 bs1.append( pycad.Point( *b1[3] ) ) # 3: centre, bottom (@ lyr1/basin)
154 bs1.append( pycad.Point( *b1[4] ) ) # 4: edge , bottom (@ lyr1/basin)
155 bs1.append( pycad.Point( *b1[5] ) ) # 5: right , bottom (@ boundary)
156 bs1.append( ly1[3] ) # 6: right , top
157
158
159 # Second-Layer.
160 ly2 = []
161 ly2.append( sl1[1] ) # 0: left , top (@ lyr2/slab)
162 ly2.append( pycad.Point( *l2[1] ) ) # 1: left , bottom (@ boundary)
163 ly2.append( pycad.Point( *l2[2] ) ) # 2: right , bottom (@ boundary)
164 ly2.append( bs1[5] ) # 3: right , top (@ basin/boundary)
165 ly2.append( bs1[4] ) # 4: edge , top (@ lyr2/basin)
166 ly2.append( bs1[3] ) # 5: centre, top (@ lyr2/basin)
167 ly2.append( bs1[2] ) # 6: left , top (@ lyr2/basin)
168 ly2.append( sl1[2] ) # 7: left , centre (@ slab/basin)
169
170
171 # Basement layer.
172 ly3 = []
173 ly3.append( ly2[1] ) # 0: left , top (@ boundary)
174 ly3.append( pycad.Point( *l3[1] ) ) # 1: left , bottom (@ boundary)
175 ly3.append( pycad.Point( *l3[2] ) ) # 2: right , bottom (@ boundary)
176 ly3.append( ly2[2] ) # 3: right , top (@ boundary)
177 #___________________________________________________________________________
178
179
180
181
182 # --------------------------------------------------------------------------
183 # Line definitions.
184 # --------------------------------------------------------------------------
185
186 #<Note>: connects the points to define lines counterclockwise;
187 # shared lines are re-used to ensure that all domains
188 # are recognised as parts of the same mesh.
189
190 # Air.
191 ln0 = []
192 ln0.append( pycad.Line(air[0], ly1[0]) ) # 0 left-top to left-bottom.
193 ln0.append( pycad.Line(ly1[0], ly1[4]) ) # 1 left-bottom to right-bottom (air-earth interface).
194 ln0.append( pycad.Line(ly1[4], air[1]) ) # 2 right-bottom to right-top.
195 ln0.append( pycad.Line(air[1], air[0]) ) # 3 right-top to left-top.
196
197 # Top Layer.
198 ln1 = []
199 ln1.append( pycad.Line(ly1[0], ly1[1]) ) # 0 left-top to left-bottom.
200 ln1.append( pycad.Line(ly1[1], ly1[2]) ) # 1 left-bottom to start-slab/basin.
201 ln1.append( pycad.Line(ly1[2], ly1[3]) ) # 2 start-slab/basin to basin-boundary
202 ln1.append( pycad.Line(ly1[3], ly1[4]) ) # 3 basin-boundary to right-top.
203 ln1.append( -ln0[1] ) # 4 right-top to left-top.
204
205
206 # Slab.
207 ln2 = []
208 ln2.append( pycad.Line(sl1[0], sl1[1]) ) # 0 left-top to left-bottom.
209 ln2.append( pycad.Line(sl1[1], sl1[2]) ) # 1 left-bottom to right-bottom.
210 ln2.append( pycad.Line(sl1[2], sl1[3]) ) # 2 right-bottom to right-top.
211 ln2.append( -ln1[1] ) # 3 right-top to left-top
212
213
214 # Basin.
215 ln3 = []
216 ln3.append( -ln2[2] ) # 0 left-top to left-centre.
217 ln3.append( pycad.Line(bs1[1], bs1[2]) ) # 1 left-centre to left-bottom.
218 ln3.append( pycad.Line(bs1[2], bs1[3]) ) # 2 left-bottom to mid-bottom.
219 ln3.append( pycad.Line(bs1[3], bs1[4]) ) # 3 mid-bottom to right-mid-top.
220 ln3.append( pycad.Line(bs1[4], bs1[5]) ) # 4 right-mid-top to right-bottom.
221 ln3.append( pycad.Line(bs1[5], bs1[6]) ) # 5 right-bottom to right-top.
222 ln3.append( -ln1[2] ) # 6 right-top to right-slab/basin.
223
224
225 # Layer below.
226 ln4 = []
227 ln4.append( pycad.Line(ly2[0], ly2[1]) ) # 0 left-top to left-bottom.
228 ln4.append( pycad.Line(ly2[1], ly2[2]) ) # 1 left-bottom to right-bottom.
229 ln4.append( pycad.Line(ly2[2], ly2[3]) ) # 2 right-bottom to right-top.
230 ln4.append( -ln3[4] ) # 3 right-top to right-mid-top.
231 ln4.append( -ln3[3] ) # 4 right-mid-top to mid-bottom.
232 ln4.append( -ln3[2] ) # 5 mid-bottom to left-bottom.
233 ln4.append( -ln3[1] ) # 6 left-bottom to left-centre.
234 ln4.append( -ln2[1] ) # 7 left-centre to left-top.
235
236 # Basement layer.
237 ln5 = []
238 ln5.append( pycad.Line(ly3[0], ly3[1]) ) # 0 left-top to left-bottom.
239 ln5.append( pycad.Line(ly3[1], ly3[2]) ) # 1 left-bottom to right-bottom.
240 ln5.append( pycad.Line(ly3[2], ly3[3]) ) # 2 right-bottom to right-top.
241 ln5.append( -ln4[1] ) # 3 right-top to left-top.
242 #___________________________________________________________________________
243
244
245
246
247 # --------------------------------------------------------------------------
248 # Domain definitions.
249 # --------------------------------------------------------------------------
250
251
252 # First define all borders.
253 borders = []
254 borders.append( pycad.CurveLoop(*ln0) )
255 borders.append( pycad.CurveLoop(*ln1) )
256 borders.append( pycad.CurveLoop(*ln2) )
257 borders.append( pycad.CurveLoop(*ln3) )
258 borders.append( pycad.CurveLoop(*ln4) )
259 borders.append( pycad.CurveLoop(*ln5) )
260
261 # And next the domains.
262 domains = []
263 for i in range( len(borders) ):
264 domains.append( pycad.PlaneSurface(borders[i]) )
265 #___________________________________________________________________________
266
267
268
269
270 # --------------------------------------------------------------------------
271 # Set element sizes in domains.
272 # --------------------------------------------------------------------------
273
274 # Horizontal extents of segments along slab and basin:
275 x_extents = []
276 x_extents.append( l1[2][0] - l1[0][0] ) # 0
277 x_extents.append( l1[3][0] - l1[2][0] ) # 1
278
279 # Number of elements in the air-domain, first-layer as well as slab- and basin-domain.
280 domains[0].setElementDistribution( x_extent / elem_sizes["large"] )
281 domains[1].setElementDistribution( x_extent / (elem_sizes["small"]) )
282 domains[2].setElementDistribution( 0.4*x_extent / (elem_sizes["small"]) )
283 domains[3].setElementDistribution( 0.5*x_extent / (elem_sizes["small"]) )
284 #<Note> slab and basin multiplied by approximate ratio of their x_extent.
285 #___________________________________________________________________________
286
287
288
289
290 #---------------------------------------------------------------------------
291 # Now define the gmsh 'design' object.
292 #---------------------------------------------------------------------------
293
294 design2D = pycad.gmsh.Design(dim=2, element_size=elem_sizes['large'], keep_files=False)
295
296 # Also specify the domains for tagging with property values later on:
297 design2D.addItems(
298 pycad.PropertySet( "air" , domains[0]) ,
299 pycad.PropertySet( "lyr1" , domains[1]) ,
300 pycad.PropertySet( "slab" , domains[2]) ,
301 pycad.PropertySet( "basin" , domains[3]) ,
302 pycad.PropertySet( "lyr2" , domains[4]) ,
303 pycad.PropertySet( "lyr3" , domains[5]) )
304
305 # Now define the unstructured finley-mesh..
306 model2D = finley.MakeDomain(design2D)
307 #___________________________________________________________________________
308
309
310 return model2D
311
312 def generateCommemi4Mesh():
313 #---------------------------------------------------------------------------
314 # DESCRIPTION:
315 # ------------
316 # Script for preparing the COMMEMI-2 2D model.
317 #
318 # The COMMEMI-4 2D model consist of a 3-layered halfspace,
319 # hosting an anomalous horizontal slab and a basin-structure
320 # in the first layer.
321 #
322 # References:
323 # -----------
324 # See Franke A., p.89, 2003 (MSc. Thesis).
325 #
326 # Antje Franke, "Zweidimensionale Finite-Elemente-Modellierung
327 # niederfrequenter elektromagnetischer Felder in der Fernzone",
328 # Diplomarbeit (MSc.), 2003, Technische Universtitaet Freiberg.
329 #
330 # --------------------------------------------------------------------------
331
332
333 #---------------------------------------------------------------------------
334 # Geometric mesh parameters.
335 # --------------------------------------------------------------------------
336
337 # Horizontal extent and start point of mesh.
338 a_extent = 50000 # 50km - Vertical extent of air-layer in (m).
339 z_extent = 50000 # 50km - Vertical extent of subsurface in (m).
340 x_extent = 60000 # 60km - Horizontal extent of model in (m).
341
342 # Start point of mesh.
343 x_start = 0 #-x_extent/2.0
344
345 # Mesh elements sizes.
346 elem_sizes = {
347 'large' : 4.00 * x_extent/100.0, #
348 'normal': 2.00 * x_extent/100.0, #
349 'small' : 0.25 * x_extent/100.0 #
350 }
351 #____________________________________________________________________________
352
353
354
355
356
357 #---------------------------------------------------------------------------
358 # Coordinate definitions.
359 # --------------------------------------------------------------------------
360
361 # X-coordinates of all domain corners (in order of appearance, left to right).
362 x0 = x_start # left (@ boundary)
363 x1 = x_start + 24000 # centre (@ slab/basin)
364 x2 = x_start + 24000 + 8000 # edge-bottom (@ slab/lyr1)
365 x3 = x_start + 24000 + 8000 + 3000 # edge-top (@ slab/lyr1)
366 x4 = x_start + x_extent # right (@ boundary)
367
368 # Y-coordinates of all domain corners (in order of appearance, top to bottom).
369 y0 = a_extent # top
370 y1 = 0 # centre (@ air/earth)
371 y2 =-500 # lyr1-bottom (@ boundary-left)
372 y3 =-1000 # basin-bottom (@ boundary-right)
373 y4 =-2000 # slab-bottom (@ boundary-left)
374 y5 =-4000 # basin-bottom (@ centre)
375 y6 =-25000 # lyr1-bottom
376 y7 =-z_extent # bottom
377
378 # Save in dictionary as a list of tuples for each domain, from left-top corner, counterclockwise.
379 coord = {
380 'air' : ([x0, y0, 0], # 0: left , top
381 [x0, y1, 0], # 1: left , bottom (@ air/earth)
382 [x4, y1, 0], # 2: right , bottom (@ air/earth)
383 [x4, y0, 0]), # 3: right , top
384
385 'lyr1' : ([x0, y1, 0], # 0: left , top
386 [x0, y2, 0], # 1: left , bottom
387 [x1, y2, 0], # 2: right , bottom (@ slab/basin)
388 [x4, y2, 0], # 3: right , bottom (@ boundary)
389 [x4, y1, 0]), # 4: right , top
390
391 'slab' : ([x0, y2, 0], # 0: left , top
392 [x0, y4, 0], # 1: left , bottom
393 [x1, y4, 0], # 2: right , bottom (@ slab/basin)
394 [x1, y2, 0]), # 3: right , top (@ slab/basin)
395
396 'basin': ([x1, y2, 0], # 0: left , top (@ slab/basin)
397 [x1, y4, 0], # 1: left , centre (@ slab/basin)
398 [x1, y5, 0], # 2: left , bottom (@ lyr1/basin)
399 [x2, y5, 0], # 3: centre, bottom (@ lyr1/basin)
400 [x3, y3, 0], # 4: edge , bottom (@ lyr1/basin)
401 [x4, y3, 0], # 5: right , bottom (@ boundary)
402 [x4, y2, 0]), # 6: right , top
403
404 'lyr2' : ([x0, y4, 0], # 0: left , top
405 [x0, y6, 0], # 1: left , bottom
406 [x4, y6, 0], # 2: right , bottom
407 [x4, y3, 0], # 3: right , top (@ basin/boundary)
408 [x3, y3, 0], # 4: edge , top (@ lyr2/basin)
409 [x2, y5, 0], # 5: centre, top (@ lyr2/basin)
410 [x1, y5, 0], # 6: left , top (@ lyr2/basin)
411 [x1, y4, 0]), # 7: left , centre (@ slab/basin)
412
413 'lyr3' : ([x0, y6, 0], # 0: left , top
414 [x0, y7, 0], # 1: left , bottom
415 [x4, y7, 0], # 2: right , bottom
416 [x4, y6, 0]), # 3: right , top
417 }
418 #___________________________________________________________________________
419
420
421
422
423
424
425
426
427
428 #---------------------------------------------------------------------------
429 # Setup the COMMEMI-4 mesh.
430 #---------------------------------------------------------------------------
431
432 # This creates the mesh and saves it to the output folder.
433 return setupMesh("TE", coord, elem_sizes)
434 #___________________________________________________________________________
435
436
437
438 if HAVE_FINLEY:
439 # ---
440 # Initialisations
441 # ---
442
443 # Get timing:
444 startTime = datetime.datetime.now()
445
446 # Mode (TE includes air-layer, whereas TM does not):
447 mode = 'TE'
448
449 # Read the mesh file and define the 'finley' domain:
450 #mesh_file = "commemi4_tm.fly"
451 #domain = finley.ReadMesh(mesh_file, numDim=2)
452 if escript.hasFeature('gmsh'):
453 domain=generateCommemi4Mesh()
454
455 # Sounding frequencies (in Hz):
456 freq_def = {"high":1.0e+0,"low":1.0e-0,"step":1}
457 # Frequencies will be mapped on a log-scale from
458 # 'high' to 'low' with 'step' points per decade.
459 # (also only one frequency must be passed via dict)
460
461 # Step sizes for sampling along vertical and horizontal axis (in m):
462 xstep=300
463 zstep=250
464
465
466 # ---
467 # Resistivity model
468 # ---
469
470 # Resistivity values assigned to tagged regions (in Ohm.m):
471 rho = [
472 1.0e+14, # 0: air 1.0e-30
473 25.0 , # 1: lyr1 0.04
474 10.0 , # 2: slab 0.1
475 2.5 , # 3: basin 0.4
476 1000.0 , # 4: lyr2 0.001
477 5.0 # 5: lyr3 0.2
478 ]
479
480 # Tags must match those in the file:
481 tags = ["air", "lyr1", "slab", "basin", "lyr2", "lyr3"]
482
483 # Optional user defined map of resistivity:
484 def f4(x,z,r): return escript.sqrt(escript.sqrt(x*x+z*z))/r
485 maps = [None, None, None, None, f4, None]
486
487
488 # ---
489 # Layer definitions for 1D response at boundaries.
490 # ---
491
492 # List with resistivity values for left and right boundary.
493 rho_1d_left = [ rho[0], rho[1], rho[2], rho[4], rho[5] ]
494 rho_1d_rght = [ rho[0], rho[1], rho[3], rho[4], rho[5] ]
495
496 # Associated interfaces for 1D response left and right (must match the mesh file).
497 ifc_1d_left = [ 50000, 0, -500, -2000, -25000, -50000]
498 ifc_1d_rght = [ 50000, 0, -500, -1000, -25000, -50000]
499
500 # Save in dictionary with layer interfaces and resistivities left and right:
501 ifc_1d = {"left":ifc_1d_left , "right":ifc_1d_rght}
502 rho_1d = {"left":rho_1d_left , "right":rho_1d_rght}
503
504
505 # ---
506 # Adjust parameters here for TM mode
507 # ---
508
509 # Simply delete first element from lists:
510 if mode.upper() == 'TM':
511 tags.pop(0)
512 rho.pop(0)
513 rho_1d['left'].pop(0)
514 rho_1d['right'].pop(0)
515 ifc_1d['left'].pop(0)
516 ifc_1d['right'].pop(0)
517 if maps is not None:
518 maps.pop(0)
519
520
521 # ---
522 # Run MT_2D
523 # ---
524
525 # Class options:
526 mt2d.MT_2D._solver = "DIRECT"
527 mt2d.MT_2D._debug = False
528
529 if mt2d.MT_2D._solver == "DIRECT" and not escript.hasFeature('paso'):
530 print("Trilinos direct solvers cannot currently handle PDE systems. Please compile with Paso.")
531 elif mt2d.MT_2D._solver == "DIRECT" and not HAVE_DIRECT:
532 if escript.getMPISizeWorld() > 1:
533 print("Direct solvers and multiple MPI processes are not currently supported.")
534 else:
535 print("escript was not built with support for direct solvers, aborting.")
536 elif not escript.hasFeature('gmsh'):
537 print("This example requires gmsh, aborting.")
538 else:
539 # Instantiate an MT_2D object with required & optional parameters:
540 obj_mt2d = mt2d.MT_2D(domain, mode, freq_def, tags, rho, rho_1d, ifc_1d,
541 xstep=xstep ,zstep=zstep, maps=None, plot=True)
542
543 # Solve for fields, apparent resistivity and phase:
544 mt2d_fields, arho_2d, aphi_2d = obj_mt2d.pdeSolve()
545
546 #
547 print("Runtime:", datetime.datetime.now()-startTime)
548 print("Done!")
549
550 else: # no finley
551 print("Finley module not available.")
552

  ViewVC Help
Powered by ViewVC 1.1.26