/[escript]/trunk/downunder/test/python/run_comm4.py
ViewVC logotype

Contents of /trunk/downunder/test/python/run_comm4.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 5753 - (show annotations)
Fri Jul 17 01:06:22 2015 UTC (3 years, 9 months ago) by sshaw
File MIME type: text/x-python
File size: 28290 byte(s)
replacing xrange for py3 compatibility
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 from __future__ import print_function, division
17 """
18 Test script to run test model COMMEMI-4
19 """
20
21 try:
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 from matplotlib import pyplot
28 HAVE_MPL=True
29 except:
30 HAVE_MPL=False
31
32 import esys.escriptcore.utestselect as unittest
33 from esys.escriptcore.testing import *
34
35 import numpy
36 import datetime
37 import esys.downunder.magtel2d as mt2d
38 import esys.escript as escript
39 import esys.finley as finley
40 import esys.escript.pdetools as pdetools
41
42 # Matplotlib uses outdated code -- ignore the warnings until an update is available:
43 import warnings
44 warnings.filterwarnings("ignore") #, category=DeprecationWarning
45 import logging
46 # this is mainly to avoid warning messages
47 logging.basicConfig(format='%(name)s: %(message)s', level=logging.INFO)
48
49 try:
50 from esys.finley import Rectangle as fRect, Brick as fBrick
51 HAVE_FINLEY = True
52 except ImportError:
53 HAVE_FINLEY = False
54
55 # ==========================================================
56 # ==========================================================
57
58
59
60 def setupMesh(mode, coord, elem_sizes):
61 #---------------------------------------------------------------------------
62 # DESCRIPTION:
63 # -----------
64 # This is a utility function which setups the COMMEMI-4 mesh.
65 #
66 #
67 # ARGUMENTS:
68 # ----------
69 # mode :: TE or TM mode.
70 # coord :: dictionary with coordinate tuples.
71 # elem_sizes :: mesh element sizes, large, normal, small.
72 #
73 # RETURNS:
74 # --------
75 # <Nothing> A mesh file is written to the output folder.
76 #
77 #
78 # AUTHOR:
79 # -------
80 # Ralf Schaa,
81 # University of Queensland
82 #
83 #---------------------------------------------------------------------------
84
85
86
87 #---------------------------------------------------------------------------
88 # Imports.
89 #---------------------------------------------------------------------------
90
91 import esys.pycad as pycad # @UnresolvedImport
92 import esys.finley as finley # @UnresolvedImport
93 import esys.escript as escript # @UnresolvedImport
94 import esys.weipa as weipa # @UnresolvedImport
95 # <Note>: "@UnresolvedImport" ignores any warnings in Eclipse/PyDev (PyDev has trouble with external libraries).
96
97
98
99 model = "COMMEMI-4"
100
101 print("Preparing the mesh " + model + " ...")
102 print("")
103
104 # Warn about magnetotelluric TM mode:
105 if mode.lower() == 'tm':
106 print("TM mode not yet supported")
107 return
108
109
110
111 # Path to write the mesh:
112 outpath = "../out/commemi4"
113
114
115
116
117 # --------------------------------------------------------------------------
118 # Initialisations.
119 # --------------------------------------------------------------------------
120
121 # Get coordinates from dictionary as list of tuples
122 a0 = coord["air"]
123 l1 = coord["lyr1"]
124 s1 = coord["slab"]
125 b1 = coord["basin"]
126 l2 = coord["lyr2"]
127 l3 = coord["lyr3"]
128
129 # Mesh length from top-boundary.
130 x_extent = abs(a0[3][0]-a0[0][0])
131
132
133
134
135 # --------------------------------------------------------------------------
136 # Point definitions.
137 # --------------------------------------------------------------------------
138
139 #<Note>: define all points spanning the mesh, anomalies and layers;
140 # note also shared domain points must be defined only once.
141
142
143 # Mesh top boundary.
144 air = []
145 air.append( pycad.Point( *a0[0] ) ) # 0: left , top (@ boundary)
146 air.append( pycad.Point( *a0[3] ) ) # 3: right , top (@ boundary)
147
148
149 # First-layer.
150 ly1 = []
151 ly1.append( pycad.Point( *l1[0] ) ) # 0: left , top (@ air/earth interface)
152 ly1.append( pycad.Point( *l1[1] ) ) # 1: left , bottom (@ boundary)
153 ly1.append( pycad.Point( *l1[2] ) ) # 2: right , bottom (@ slab/basin)
154 ly1.append( pycad.Point( *l1[3] ) ) # 3: right , bottom (@ boundary)
155 ly1.append( pycad.Point( *l1[4] ) ) # 4: right , top (@ air/earth interface)
156
157
158 # Slab.
159 sl1 = []
160 sl1.append( ly1[1] ) # 0: left , top (@ boundary)
161 sl1.append( pycad.Point( *s1[1] ) ) # 1: left , bottom (@ boundary)
162 sl1.append( pycad.Point( *s1[2] ) ) # 2: right , bottom (@ slab/basin)
163 sl1.append( ly1[2] ) # 3: right , top (@ slab/basin)
164
165
166 # Basin.
167 bs1 = []
168 bs1.append( ly1[2] ) # 0: left , top (@ slab/basin)
169 bs1.append( sl1[2] ) # 1: left , centre (@ slab/basin)
170 bs1.append( pycad.Point( *b1[2] ) ) # 2: left , bottom (@ lyr1/basin)
171 bs1.append( pycad.Point( *b1[3] ) ) # 3: centre, bottom (@ lyr1/basin)
172 bs1.append( pycad.Point( *b1[4] ) ) # 4: edge , bottom (@ lyr1/basin)
173 bs1.append( pycad.Point( *b1[5] ) ) # 5: right , bottom (@ boundary)
174 bs1.append( ly1[3] ) # 6: right , top
175
176
177 # Second-Layer.
178 ly2 = []
179 ly2.append( sl1[1] ) # 0: left , top (@ lyr2/slab)
180 ly2.append( pycad.Point( *l2[1] ) ) # 1: left , bottom (@ boundary)
181 ly2.append( pycad.Point( *l2[2] ) ) # 2: right , bottom (@ boundary)
182 ly2.append( bs1[5] ) # 3: right , top (@ basin/boundary)
183 ly2.append( bs1[4] ) # 4: edge , top (@ lyr2/basin)
184 ly2.append( bs1[3] ) # 5: centre, top (@ lyr2/basin)
185 ly2.append( bs1[2] ) # 6: left , top (@ lyr2/basin)
186 ly2.append( sl1[2] ) # 7: left , centre (@ slab/basin)
187
188
189 # Basement layer.
190 ly3 = []
191 ly3.append( ly2[1] ) # 0: left , top (@ boundary)
192 ly3.append( pycad.Point( *l3[1] ) ) # 1: left , bottom (@ boundary)
193 ly3.append( pycad.Point( *l3[2] ) ) # 2: right , bottom (@ boundary)
194 ly3.append( ly2[2] ) # 3: right , top (@ boundary)
195 #___________________________________________________________________________
196
197
198
199
200 # --------------------------------------------------------------------------
201 # Line definitions.
202 # --------------------------------------------------------------------------
203
204 #<Note>: connects the points to define lines counterclockwise;
205 # shared lines are re-used to ensure that all domains
206 # are recognised as parts of the same mesh.
207
208 # Air.
209 ln0 = []
210 ln0.append( pycad.Line(air[0], ly1[0]) ) # 0 left-top to left-bottom.
211 ln0.append( pycad.Line(ly1[0], ly1[4]) ) # 1 left-bottom to right-bottom (air-earth interface).
212 ln0.append( pycad.Line(ly1[4], air[1]) ) # 2 right-bottom to right-top.
213 ln0.append( pycad.Line(air[1], air[0]) ) # 3 right-top to left-top.
214
215 # Top Layer.
216 ln1 = []
217 ln1.append( pycad.Line(ly1[0], ly1[1]) ) # 0 left-top to left-bottom.
218 ln1.append( pycad.Line(ly1[1], ly1[2]) ) # 1 left-bottom to start-slab/basin.
219 ln1.append( pycad.Line(ly1[2], ly1[3]) ) # 2 start-slab/basin to basin-boundary
220 ln1.append( pycad.Line(ly1[3], ly1[4]) ) # 3 basin-boundary to right-top.
221 ln1.append( -ln0[1] ) # 4 right-top to left-top.
222
223
224 # Slab.
225 ln2 = []
226 ln2.append( pycad.Line(sl1[0], sl1[1]) ) # 0 left-top to left-bottom.
227 ln2.append( pycad.Line(sl1[1], sl1[2]) ) # 1 left-bottom to right-bottom.
228 ln2.append( pycad.Line(sl1[2], sl1[3]) ) # 2 right-bottom to right-top.
229 ln2.append( -ln1[1] ) # 3 right-top to left-top
230
231
232 # Basin.
233 ln3 = []
234 ln3.append( -ln2[2] ) # 0 left-top to left-centre.
235 ln3.append( pycad.Line(bs1[1], bs1[2]) ) # 1 left-centre to left-bottom.
236 ln3.append( pycad.Line(bs1[2], bs1[3]) ) # 2 left-bottom to mid-bottom.
237 ln3.append( pycad.Line(bs1[3], bs1[4]) ) # 3 mid-bottom to right-mid-top.
238 ln3.append( pycad.Line(bs1[4], bs1[5]) ) # 4 right-mid-top to right-bottom.
239 ln3.append( pycad.Line(bs1[5], bs1[6]) ) # 5 right-bottom to right-top.
240 ln3.append( -ln1[2] ) # 6 right-top to right-slab/basin.
241
242
243 # Layer below.
244 ln4 = []
245 ln4.append( pycad.Line(ly2[0], ly2[1]) ) # 0 left-top to left-bottom.
246 ln4.append( pycad.Line(ly2[1], ly2[2]) ) # 1 left-bottom to right-bottom.
247 ln4.append( pycad.Line(ly2[2], ly2[3]) ) # 2 right-bottom to right-top.
248 ln4.append( -ln3[4] ) # 3 right-top to right-mid-top.
249 ln4.append( -ln3[3] ) # 4 right-mid-top to mid-bottom.
250 ln4.append( -ln3[2] ) # 5 mid-bottom to left-bottom.
251 ln4.append( -ln3[1] ) # 6 left-bottom to left-centre.
252 ln4.append( -ln2[1] ) # 7 left-centre to left-top.
253
254 # Basement layer.
255 ln5 = []
256 ln5.append( pycad.Line(ly3[0], ly3[1]) ) # 0 left-top to left-bottom.
257 ln5.append( pycad.Line(ly3[1], ly3[2]) ) # 1 left-bottom to right-bottom.
258 ln5.append( pycad.Line(ly3[2], ly3[3]) ) # 2 right-bottom to right-top.
259 ln5.append( -ln4[1] ) # 3 right-top to left-top.
260 #___________________________________________________________________________
261
262
263
264
265 # --------------------------------------------------------------------------
266 # Domain definitions.
267 # --------------------------------------------------------------------------
268
269
270 # First define all borders.
271 borders = []
272 borders.append( pycad.CurveLoop(*ln0) )
273 borders.append( pycad.CurveLoop(*ln1) )
274 borders.append( pycad.CurveLoop(*ln2) )
275 borders.append( pycad.CurveLoop(*ln3) )
276 borders.append( pycad.CurveLoop(*ln4) )
277 borders.append( pycad.CurveLoop(*ln5) )
278
279 # And next the domains.
280 domains = []
281 for i in range( len(borders) ):
282 domains.append( pycad.PlaneSurface(borders[i]) )
283 #___________________________________________________________________________
284
285
286
287
288 # --------------------------------------------------------------------------
289 # Set element sizes in domains.
290 # --------------------------------------------------------------------------
291
292 # Horizontal extents of segments along slab and basin:
293 x_extents = []
294 x_extents.append( l1[2][0] - l1[0][0] ) # 0
295 x_extents.append( l1[3][0] - l1[2][0] ) # 1
296
297 # Number of elements in the air-domain, first-layer as well as slab- and basin-domain.
298 domains[0].setElementDistribution( x_extent / elem_sizes["large"] )
299 domains[1].setElementDistribution( x_extent / (elem_sizes["small"]) )
300 domains[2].setElementDistribution( 0.4*x_extent / (elem_sizes["small"]) )
301 domains[3].setElementDistribution( 0.5*x_extent / (elem_sizes["small"]) )
302 #<Note> slab and basin multiplied by approximate ratio of their x_extent.
303 #___________________________________________________________________________
304
305
306
307
308 #---------------------------------------------------------------------------
309 # Now define the gmsh 'design' object.
310 #---------------------------------------------------------------------------
311
312 design2D = pycad.gmsh.Design(dim=2, element_size=elem_sizes['large'], keep_files=False)
313
314 # Also specify the domains for tagging with property values later on:
315 design2D.addItems(
316 pycad.PropertySet( "air" , domains[0]) ,
317 pycad.PropertySet( "lyr1" , domains[1]) ,
318 pycad.PropertySet( "slab" , domains[2]) ,
319 pycad.PropertySet( "basin" , domains[3]) ,
320 pycad.PropertySet( "lyr2" , domains[4]) ,
321 pycad.PropertySet( "lyr3" , domains[5]) )
322
323 # Now define the unstructured finley-mesh..
324 model2D = finley.MakeDomain(design2D)
325 #___________________________________________________________________________
326
327
328 return model2D
329
330 def generateCommemi4Mesh():
331 #---------------------------------------------------------------------------
332 # DESCRIPTION:
333 # ------------
334 # Script for preparing the COMMEMI-2 2D model.
335 #
336 # The COMMEMI-4 2D model consist of a 3-layered halfspace,
337 # hosting an anomalous horizontal slab and a basin-structure
338 # in the first layer.
339 #
340 # References:
341 # -----------
342 # See Franke A., p.89, 2003 (MSc. Thesis).
343 #
344 # Antje Franke, "Zweidimensionale Finite-Elemente-Modellierung
345 # niederfrequenter elektromagnetischer Felder in der Fernzone",
346 # Diplomarbeit (MSc.), 2003, Technische Universtitaet Freiberg.
347 #
348 # --------------------------------------------------------------------------
349
350
351 #---------------------------------------------------------------------------
352 # Geometric mesh parameters.
353 # --------------------------------------------------------------------------
354
355 # Horizontal extent and start point of mesh.
356 a_extent = 50000 # 50km - Vertical extent of air-layer in (m).
357 z_extent = 50000 # 50km - Vertical extent of subsurface in (m).
358 x_extent = 60000 # 60km - Horizontal extent of model in (m).
359
360 # Start point of mesh.
361 x_start = 0 #-x_extent/2.0
362
363 # Mesh elements sizes.
364 elem_sizes = {
365 'large' : 4.00 * x_extent/100.0, #
366 'normal': 2.00 * x_extent/100.0, #
367 'small' : 0.25 * x_extent/100.0 #
368 }
369 #____________________________________________________________________________
370
371
372
373
374
375 #---------------------------------------------------------------------------
376 # Coordinate definitions.
377 # --------------------------------------------------------------------------
378
379 # X-coordinates of all domain corners (in order of appearance, left to right).
380 x0 = x_start # left (@ boundary)
381 x1 = x_start + 24000 # centre (@ slab/basin)
382 x2 = x_start + 24000 + 8000 # edge-bottom (@ slab/lyr1)
383 x3 = x_start + 24000 + 8000 + 3000 # edge-top (@ slab/lyr1)
384 x4 = x_start + x_extent # right (@ boundary)
385
386 # Y-coordinates of all domain corners (in order of appearance, top to bottom).
387 y0 = a_extent # top
388 y1 = 0 # centre (@ air/earth)
389 y2 =-500 # lyr1-bottom (@ boundary-left)
390 y3 =-1000 # basin-bottom (@ boundary-right)
391 y4 =-2000 # slab-bottom (@ boundary-left)
392 y5 =-4000 # basin-bottom (@ centre)
393 y6 =-25000 # lyr1-bottom
394 y7 =-z_extent # bottom
395
396 # Save in dictionary as a list of tuples for each domain, from left-top corner, counterclockwise.
397 coord = {
398 'air' : ([x0, y0, 0], # 0: left , top
399 [x0, y1, 0], # 1: left , bottom (@ air/earth)
400 [x4, y1, 0], # 2: right , bottom (@ air/earth)
401 [x4, y0, 0]), # 3: right , top
402
403 'lyr1' : ([x0, y1, 0], # 0: left , top
404 [x0, y2, 0], # 1: left , bottom
405 [x1, y2, 0], # 2: right , bottom (@ slab/basin)
406 [x4, y2, 0], # 3: right , bottom (@ boundary)
407 [x4, y1, 0]), # 4: right , top
408
409 'slab' : ([x0, y2, 0], # 0: left , top
410 [x0, y4, 0], # 1: left , bottom
411 [x1, y4, 0], # 2: right , bottom (@ slab/basin)
412 [x1, y2, 0]), # 3: right , top (@ slab/basin)
413
414 'basin': ([x1, y2, 0], # 0: left , top (@ slab/basin)
415 [x1, y4, 0], # 1: left , centre (@ slab/basin)
416 [x1, y5, 0], # 2: left , bottom (@ lyr1/basin)
417 [x2, y5, 0], # 3: centre, bottom (@ lyr1/basin)
418 [x3, y3, 0], # 4: edge , bottom (@ lyr1/basin)
419 [x4, y3, 0], # 5: right , bottom (@ boundary)
420 [x4, y2, 0]), # 6: right , top
421
422 'lyr2' : ([x0, y4, 0], # 0: left , top
423 [x0, y6, 0], # 1: left , bottom
424 [x4, y6, 0], # 2: right , bottom
425 [x4, y3, 0], # 3: right , top (@ basin/boundary)
426 [x3, y3, 0], # 4: edge , top (@ lyr2/basin)
427 [x2, y5, 0], # 5: centre, top (@ lyr2/basin)
428 [x1, y5, 0], # 6: left , top (@ lyr2/basin)
429 [x1, y4, 0]), # 7: left , centre (@ slab/basin)
430
431 'lyr3' : ([x0, y6, 0], # 0: left , top
432 [x0, y7, 0], # 1: left , bottom
433 [x4, y7, 0], # 2: right , bottom
434 [x4, y6, 0]), # 3: right , top
435 }
436 #___________________________________________________________________________
437
438
439 #---------------------------------------------------------------------------
440 # Setup the COMMEMI-4 mesh.
441 #---------------------------------------------------------------------------
442
443 # This creates the mesh and saves it to the output folder.
444 return setupMesh("TE", coord, elem_sizes)
445 #___________________________________________________________________________
446
447
448
449 # ==========================================================
450 # ==========================================================
451
452 class Test_COMMEMI4(unittest.TestCase):
453 @unittest.skipIf(not HAVE_FINLEY, "Test requires finley to be available")
454 @unittest.skipIf(not escript.getEscriptParamInt("PASO_DIRECT"), "Missing direct solvers")
455 @unittest.skipIf(escript.getMPISizeWorld() > 1, "Direct solvers and MPI are currently incompatible")
456 def test_comm4(self):
457 # ---
458 # Initialisations
459 # ---
460
461 # Get timing:
462 startTime = datetime.datetime.now()
463
464 # Mode (TE includes air-layer, whereas TM does not):
465 mode = 'TE'
466
467 # Read the mesh file and define the 'finley' domain:
468 #mesh_file = "mesh/commemi-4/commemi4_tm.msh"
469 #domain = finley.ReadGmsh(mesh_file, numDim=2)
470 domain = generateCommemi4Mesh()
471
472 #mesh_file = "mesh/commemi-4/commemi4_tm.fly"
473 #domain = finley.ReadMesh(mesh_file)
474
475 # Sounding frequencies (in Hz):
476 freq_def = {"high":1.0e+0,"low":1.0e-0,"step":1}
477 # Frequencies will be mapped on a log-scale from
478 # 'high' to 'low' with 'step' points per decade.
479 # (also only one frequency must be passed via dict)
480
481 # Step sizes for sampling along vertical and horizontal axis (in m):
482 xstep=100
483 zstep=100
484
485
486
487 # ---
488 # Resistivity model
489 # ---
490
491 # Resistivity values assigned to tagged regions (in Ohm.m):
492 rho = [
493 1.0e+14, # 0: air 1.0e-30
494 25.0 , # 1: lyr1 0.04
495 10.0 , # 2: slab 0.1
496 2.5 , # 3: basin 0.4
497 1000.0 , # 4: lyr2 0.001
498 5.0 # 5: lyr3 0.2
499 ]
500
501 # Tags must match those in the file:
502 tags = ["air", "lyr1", "slab", "basin", "lyr2", "lyr3"]
503
504 # Optional user defined map of resistivity:
505 def f4(x,z,r): return escript.sqrt(escript.sqrt(x*x+z*z))/r
506 maps = [None, None, None, None, f4, None]
507
508
509
510 # ---
511 # Layer definitions for 1D response at boundaries.
512 # ---
513
514 # List with resistivity values for left and right boundary.
515 rho_1d_left = [ rho[0], rho[1], rho[2], rho[4], rho[5] ]
516 rho_1d_rght = [ rho[0], rho[1], rho[3], rho[4], rho[5] ]
517
518 # Associated interfaces for 1D response left and right (must match the mesh file).
519 ifc_1d_left = [ 50000, 0, -500, -2000, -25000, -50000]
520 ifc_1d_rght = [ 50000, 0, -500, -1000, -25000, -50000]
521
522 # Save in dictionary with layer interfaces and resistivities left and right:
523 ifc_1d = {"left":ifc_1d_left , "right":ifc_1d_rght}
524 rho_1d = {"left":rho_1d_left , "right":rho_1d_rght}
525
526
527
528 # ---
529 # Adjust parameters here for TM mode
530 # ---
531
532 # Simply delete first element from lists:
533 if mode.upper() == 'TM':
534 tags.pop(0)
535 rho.pop(0)
536 rho_1d['left'].pop(0)
537 rho_1d['right'].pop(0)
538 ifc_1d['left'].pop(0)
539 ifc_1d['right'].pop(0)
540 if maps is not None:
541 maps.pop(0)
542
543
544
545 # ---
546 # Run MT_2D
547 # ---
548
549 # Class options:
550 mt2d.MT_2D._solver = "DIRECT" #"ITERATIVE" #"CHOLEVSKY" #"CGLS " #"BICGSTAB" #"DIRECT" "ITERATIVE"
551 mt2d.MT_2D._debug = False
552
553 # Instantiate an MT_2D object with required & optional parameters:
554 obj_mt2d = mt2d.MT_2D(domain, mode, freq_def, tags, rho, rho_1d, ifc_1d,
555 xstep=xstep ,zstep=zstep, maps=None, plot=False)
556
557 # Solve for fields, apparent resistivity and phase:
558 mt2d_fields, arho_2d, aphi_2d = obj_mt2d.pdeSolve()
559
560
561 #import random
562
563 #mt2d_fields[0]['real']+=random.random()
564 #mt2d_fields[0]['imag']+=50*random.random()
565
566 #print(arho_2d[0][0])
567 #for i in range(len(aphi_2d[0])):
568 #aphi_2d[0][i]+=(50*random.random())
569
570 #for i in range(len(arho_2d[0])):
571 #arho_2d[0][i]-=7*(random.random())
572
573
574 # ---
575 # User defined plots
576 # ---
577
578 from scipy.interpolate import InterpolatedUnivariateSpline
579
580 # Setup abscissas/Ordinates for escript data:
581 x = numpy.array( obj_mt2d.loc.getX() )[:,0]
582 y0 = numpy.array( obj_mt2d.loc.getValue(arho_2d[0]) )
583 y1 = numpy.array( obj_mt2d.loc.getValue(aphi_2d[0]) )
584
585 # Zhdanov et al, 1997, -- Model 2D-1 Table B.33. Model2D-4 (T=1.0, z=0), see
586 # "Methods for modelling electromagnetic fields. Results from COMMEMI -- the
587 # international project on the comparison of modelling results for electromag-
588 # netic induction", Journal of Applied Geophysics, 133-271
589 rte = [12.70, 12.00, 8.80, 6.84, 6.67, 6.25] # TE rho_a (3 Canada)
590 rtm = [11.40, 11.50, 9.03, 6.78, 6.80, 5.71] # TM rho_a (3 Canada)
591 if mode.lower() == 'te':
592 ra = rte
593 else:
594 ra = rtm
595 # Associated stations shifted to match escript coordinates:
596 xs = numpy.array( [-10, -7, -6, -5, 2, 5] )*1000 + x.max()/2.0
597
598 # Setup interpolation to get values at specified stations (for comparison):
599 fi = InterpolatedUnivariateSpline(x, y0)
600 # Save esscript values at comparison points in text file:
601 # uncomment to investigate
602 #numpy.savetxt("mesh/commemi-4/commemi4_"+mode.lower()+".dat", numpy.column_stack((xs,fi(xs))), fmt='%g')
603
604 # X plot-limits:
605 x0lim = [2000,38000]
606 #y1lim = [0,120]
607 #y2lim = [40,85]
608
609 # Plot labels:
610 title = ' escript COMMEMI-4 MT-2D ' + '(' + mode.upper() + ')' + ' freq: ' + str(obj_mt2d.frequencies[0]) + ' Hz'
611 ylbl0 = r'resistivity $(\Omega m)$'
612 ylbl1 = r'phase $(\circ)$'
613 xlbl1 = 'X (m)'
614 if HAVE_MPL:
615 # Setup the plot window with app. res. on top and phase on bottom:
616 f, ax = pyplot.subplots(2, figsize=(3.33,3.33), dpi=1200,
617 facecolor='w', edgecolor='k', sharex=True) # Mind shared axis
618 f.subplots_adjust(hspace=0.1, top=0.95, left=0.135, bottom=0.125, right=0.975)
619 f.suptitle(title, y=0.99,fontsize=8) #
620
621 # Top: apparent resistivity and points from Weaver for comparison:
622 ax[0].plot(x, y0, color='red', label = 'escript')
623 ax[0].plot(xs,ra, linestyle='', markersize=3, marker='o',
624 color='blue', label = 'Weaver')
625 ax[0].grid(b=True, which='both', color='grey',linestyle=':')
626 ax[0].set_ylabel( ylbl0)
627 ax[0].yaxis.set_label_coords(-0.082, 0.5)
628 # Plot limits:
629 ax[0].set_xlim(x0lim)
630 #ax[0].set_ylim(y1lim)
631
632 # Bottom: phase on linear plot
633 ax[1].plot(x,y1, color='blue')
634 ax[1].grid(b=True, which='both', color='grey',linestyle=':')
635 ax[1].set_xlabel( xlbl1 )
636 ax[1].set_ylabel( ylbl1 )
637 # Plot limits:
638 ax[1].set_xlim(x0lim)
639 #ax[1].set_ylim(y2lim)
640
641 # ask matplotlib for the plotted objects and their labels
642 lna, la = ax[0].get_legend_handles_labels()
643 ax[0].legend(lna, la, bbox_to_anchor=(0.02, 0.325), loc=2,
644 borderaxespad=0.,prop={'size':8}, frameon=False)
645
646 pyplot.ticklabel_format(style='sci', axis='x', scilimits=(0,0),
647 useMathText=True)
648 ax[0].xaxis.major.formatter._useMathText = True
649 pyplot.rc('font', **{'size': 8,'family':'sans-serif'})
650 #uncomment to allow visual inspection
651 #f.savefig("mesh/commemi-4/commemi4_"+mode.lower()+".png", dpi=1200)
652
653 # Now let's see if the points match
654 # First, we need to find correspondance between xs and x
655 indices=[]
656 for i in range(len(xs)):
657 mindiff=40000
658 mindex=0
659 for j in range(len(x)):
660 if abs(xs[i]-x[j]) < mindiff:
661 mindiff=abs(xs[i]-x[j])
662 mindex=j
663 indices.append(mindex)
664
665 # The following are very simple checks based on the visual shape of the correct result
666 maxdiff=0
667 for i in range(len(indices)):
668 if abs(y0[indices[i]]-ra[i])>maxdiff:
669 maxdiff=abs(y0[indices[i]]-ra[i])
670
671 if maxdiff>5: #Threshold is pretty arbitrary
672 raise RuntimeError("Mismatch with reference data")
673
674 c=0
675 for y in y1:
676 if y<46:
677 c+=1
678
679 if not (62 < escript.Lsup(y1) < 64):
680 raise RuntimeError("Peak of bottom plot is off.")
681
682 if not (0.62 < c/len(y1) < 0.65):
683 print(c/len(y1))
684 raise RuntimeError("Bottom plot has too many high points")
685
686 #
687 print("Runtime:", datetime.datetime.now()-startTime)
688 print("Done!")
689
690
691 if __name__ == '__main__':
692 run_tests(__name__, exit_on_failure=True)

  ViewVC Help
Powered by ViewVC 1.1.26