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

Contents of /release/4.1/doc/examples/inversion/test_commemi4.py

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.26