/[escript]/trunk/downunder/py_src/magtel2d.py
ViewVC logotype

Contents of /trunk/downunder/py_src/magtel2d.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 5748 - (show annotations)
Wed Jul 15 06:49:05 2015 UTC (3 years, 9 months ago) by sshaw
File MIME type: text/x-python
File size: 32364 byte(s)
python3 prints for 2D MT source
1 # -*- coding: utf-8 -*-
2 ##############################################################################
3 #
4 # Copyright (c) 2003-2014 by University of Queensland
5 # http://www.uq.edu.au
6 #
7 # Primary Business: Queensland, Australia
8 # Licensed under the Open Software License version 3.0
9 # http://www.opensource.org/licenses/osl-3.0.php
10 #
11 # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 # Development 2012-2013 by School of Earth Sciences
13 # Development from 2014 by Centre for Geoscience Computing (GeoComp)
14 #
15 ##############################################################################
16
17 __copyright__="""Copyright (c) 2003-2014 by University of Queensland
18 http://www.uq.edu.au
19 Primary Business: Queensland, Australia"""
20 __license__="""Licensed under the Open Software License version 3.0
21 http://www.opensource.org/licenses/osl-3.0.php"""
22 __url__="https://launchpad.net/escript-finley"
23
24 """
25 2D Magnetotelluric modelling for TE and TM mode.
26
27 :var __author__: name of author
28 :var __copyright__: copyrights
29 :var __license__: licence agreement
30 :var __url__: url entry point on documentation
31 :var __version__: version
32 :var __date__: date of the version
33 """
34
35 __author__="Ralf Schaa, r.schaa@uq.edu.au"
36
37 import os
38 import numpy
39 import math
40 import cmath
41 import types
42 from . import magtel1d as mt1d
43 import matplotlib.pyplot as plt
44 import esys.weipa as weipa
45 import esys.escript as escript
46 import esys.finley as finley
47 import esys.escript.pdetools as pdetools
48 import esys.escript.linearPDEs as pde
49
50
51 class MT_2D(object):
52
53 # class options:
54 _debug = False #
55 _solver = "DEFAULT" #
56
57 # 'private' field:
58 __version = 0.1 #
59
60 """
61 DESCRIPTION:
62 ------------
63 solves the scalar 2-D electromagnetic diffusion equation,
64 (where 'u' is the electric field E or magnetic field H).
65
66 [1] -div( k*grad(u) ) + q*u = 0 (+ Boundary Conditions)
67
68 In 2D the equation is solved for the transverse electric
69 field (TE mode) or transverse magnetic field (TM mode).
70 These fields are parallel to the 2D strike direction.
71 Based on the actual mode, the coefficients are given by:
72
73 TE: k = 1/mu , q = i*w*sigma
74 TM: k = 1/sigma, q = i*w*mu
75
76 'mu' is the vacuum permeability,
77 'i' is the imaginary unit
78 'w' is the angular frequency
79 'sigma' is the conductivity
80
81 The EM diffusion equation is complex and is solved as
82 a coupled PDE for the real and imaginary parts. The
83 coupled PDE is given by the following equations, with
84 Er, Ei and Hr, Hi are the real and imaginary components
85 of the electric and magnetic field, respectively:
86
87 TE:
88 [2] div( grad(Er) ) + w*mu*sigma*Ei = 0
89 [3] div( grad(Ei) ) - w*mu*sigma*Er = 0
90
91 the complementary magnetic fields
92 are calculated via Faraday's Law:
93
94 [4] Hr =-d/dz(Ei) / (w*mu)
95 [5] Hi = d/dz(Er) / (w*mu)
96
97
98 TM:
99 [6] div( rho*grad(Hr) ) + w*mu*Hi = 0
100 [7] div( rho*grad(Hi) ) - w*mu*Hr = 0
101
102 (resistivity 'rho' is 1/sigma)
103 the complementary magnetic fields
104 are calculated via Ampere's Law:
105
106 [8] Er = d/dz(Hr) * rho
107 [9] Ei = d/dz(Hi) * rho
108
109
110 Based on the ratio of electric to magnetic field
111 apparent resistivity and phase is calculated, viz:
112
113 rho_a = (1/w*mu) * [ (Er)^2 + (Ei)^2 ] / [ (Hr)^2 + (Hi)^2 ]
114 phase = arctan( [Ei*Hr - Er*Hi] / [Er*Hr + Ei*Hi] )
115
116
117 Boundary conditions:
118 --------------------
119 the source term on the right-hand-side of equation [1] is zero,
120 i.e. no artificial source is employed but instead the 'source'
121 is provided via the boundary conditions of the PDE which are
122 given as Dirichlet conditions at all boundaries. To calculate
123 the Dirichlet values, a 1D response is calculated at the left
124 and right boundary (based on the 1D recursion formula for MT).
125 Interpolation from the left to the right sides then provides
126 the values at the top and bottom boundary. See module 'mt1d'
127 for details of the computation of the 1D response. Once the
128 values on the boundaries have been calculated, the values
129 inside the domain are solved in this class.
130 """
131
132 def __init__(self, domain, mode, freq_def, tags, rho, rho_1d, ifc_1d, xstep=100, zstep=100, maps=None, plot=False):
133 """
134 DESCRIPTION:
135 -----------
136 Constructor which initialises the 2D magnetotelluric class:
137 (*) check for argument type
138 (*) check for valid argument values
139 (*) initialises required values
140
141 ARGUMENTS:
142 ----------
143 param domain :: the 2d mesh domain
144 type domain :: ``escript data object``
145
146 param mode :: TE or TM mode
147 type mode :: ``string``
148
149 param freq_def :: highest/lowest frequency & points per decade
150 type freq_def :: ``dictionary``
151
152 param tags :: the tag names of the regions defined in the mesh
153 type tags :: ``list``
154
155 param rho :: the resistivity values of the regions in the mesh
156 type rho :: ``list``
157
158 param rho_1d :: the resistivity values at the left & right boundary
159 type rho_1d :: ``dictionary``
160
161 param ifc_1d :: the layer interface depths of the left & right boundary
162 type ifc_1d :: ``dictionary``
163
164 param xstep :: user-defined step size for horizontal sample list
165 type xstep :: ``number`` (optional)
166
167 param zstep :: user-defined step size for vertical sample list
168 type zstep :: ``number`` (optional)
169
170 param maps :: list with user-defined functions which map the resistivity for each region
171 type maps :: ``list`` (optional)
172
173 param plot :: user-defined flag to show a plot of apparent resistivity and phase at each frequency
174 type plot :: ``boolean`` (optional)
175
176
177
178 DATA ATTRIBUTES:
179 ---------------
180 self.domain :: escript data object of mesh
181 self.X :: escript data object with all mesh coordinates
182 self.mode :: string with TE or TM mode
183 self.xmin :: float with x-coordinate minimum
184 self.xmax :: float with x-coordinate maximum
185 self.zmin :: float with z-coordinate minimum
186 self.zmax :: float with z-coordinate maximum
187 self.zstep :: number with sample step in vertical direction
188 self.xstep :: number with sample step in horizontal direction
189 self.rho :: list with resistivity values of all regions
190 self.rho_1d :: dictionary with resistivity values at boundaries left/right
191 self.ifc_1d :: dictionary with interface depths at boundaries left/right
192 self.plot :: boolean flag to show plots of apparent resistivity and phase
193 self.sigma :: escript data object with the conductivity model (based on 'rho' and 'maps')
194 self.frequencies :: list of sounding frequencies
195 self.boundary_mask :: Dirichlet mask at boundaries
196 """
197
198
199 # ---
200 # Checks
201 # ---
202
203 # Types:
204 if not isinstance(domain, finley.finleycpp.MeshAdapter ):
205 raise ValueError("Input parameter DOMAIN must be an Escript mesh")
206 if not isinstance(mode, str):
207 raise ValueError("Input parameter MODE must be a string")
208 if not isinstance(freq_def, dict) or len(freq_def) != 3:
209 raise ValueError("Input parameter FREQ_DEF must be a dictionary with length 3")
210 if not isinstance(tags, list) or not all(isinstance(t,str) for t in tags):
211 raise ValueError("Input parameter TAGS must be a list of strings")
212 if not isinstance(rho, list) or not all(isinstance(d,(int,long,float)) for d in rho):
213 raise ValueError("Input parameter RHO must be a list of numbers")
214 if not isinstance(rho_1d, dict) or len(rho_1d) != 2:
215 raise ValueError("Input parameter RHO_1D must be a dictionary with length 2")
216 if not isinstance(ifc_1d, dict) or len(ifc_1d) != 2:
217 raise ValueError("Input parameter IFC_1D must be a dictionary with length 2")
218 if not isinstance(xstep, (int,long,float)):
219 raise ValueError("Optional input parameter XSTEP must be a number")
220 if not isinstance(zstep, (int,long,float)):
221 raise ValueError("Optional input parameter ZSTEP must be a number")
222 if maps is not None:
223 if not isinstance(maps, list) or not all(isinstance(m,(types.FunctionType, types.NoneType)) for m in maps):
224 raise ValueError("Optional input parameter MAPS must be a list of Functions or Nones")
225 if plot is not None:
226 if not isinstance(plot, bool):
227 raise ValueError("Optional input parameter PLOT must be True or False")
228
229 # Values:
230 if mode.upper() != "TE" and mode.upper() != "TM": # Check mode:
231 raise ValueError("Input parameter mode must be either 'TE' or 'TM'")
232 if not 'high' in freq_def and not 'low' in freq_def and not 'step' in freq_def:
233 raise ValueError("Input dictionary FREQ_DEF must have keys 'high', 'low' and 'step' defined" )
234 if freq_def['high'] < freq_def['low']:
235 raise ValueError("High frequency value is smaller than low frequency value in input parameter FREQ_DEF")
236 if freq_def['step'] < 1:
237 raise ValueError("Step frequency value is smaller than 1 in input parameter FREQ_DEF")
238 if not all(r>0 for r in rho): # Check resistivity values:
239 raise ValueError("Input parameter RHO must be all positive")
240 if len(rho) != len(tags): # Check resistivity list-length:
241 raise ValueError("Input parameter RHO must have the same length as input parameter TAGS")
242 if not 'left' in rho_1d and not 'right' in rho_1d:
243 raise ValueError("Input dictionary RHO_1D must have keys 'left' and 'right' defined" )
244 if not 'left' in ifc_1d and not 'right' in ifc_1d:
245 raise ValueError("Input dictionary IFC_1D must have keys 'left' and 'right' defined" )
246 if len(ifc_1d['left'])-1 != len(rho_1d['left']) and len(ifc_1d['right'])-1 != len(rho_1d['right']):
247 raise ValueError("Lists with values in input dictionary RHO_1D must have length equal to IFC_1D" )
248 if xstep < 0.5: # Step size should be non-zero but should not be smaller than 0.5m:
249 raise ValueError("Input parameter XSTEP must be at least 0.5" )
250 if zstep < 0.5: # Step size should be non-zero but should not be smaller than 0.5m:
251 raise ValueError("Input parameter ZSTEP must be at least 0.5" )
252
253
254
255 # ---
256 # Domain coordinates & tags:
257 # ---
258
259 # Extract the model coordinates..
260 X = escript.Solution(domain).getX()
261
262 # Get the Min/Max coordinates:
263 xmin = escript.inf(X[0])
264 xmax = escript.sup(X[0])
265 zmin = escript.inf(X[1])
266 zmax = escript.sup(X[1])
267
268 # Get the tag names from the mesh file
269 mesh_tags = escript.getTagNames(domain)
270
271 if xmin >= xmax or zmin >= zmax: raise ValueError("The mesh limits are not valid (min >= max)" )
272 if zmin >= 0 : raise ValueError("The mesh must be defined with a negative vertical axis" )
273 if not set(mesh_tags) == set(tags) :
274 print("user-tags:", tags)
275 print("mesh-tags:", mesh_tags)
276 raise ValueError("Input parameter TAGS does not match the tags defined in the mesh")
277
278
279
280 # ---
281 # Define the boundary mask:
282 # ---
283
284 boundary_mask = self.__setBoundaryMask(X)
285
286
287 # ---
288 # Calculate list of sounding frequencies:
289 # ---
290
291 frequencies = self.__getSoundingFrequencies(freq_def)
292
293
294
295 # ---
296 # Tag the domain with conductivity maps:
297 # ---
298
299 sigma_model = self.__tagDomain(domain, X, tags, rho, maps)
300
301 # Check for valid values
302 if escript.inf(sigma_model) < 0 or escript.sup(sigma_model) < 0:
303 raise ValueError("Negative conductivity encountered" )
304 if cmath.isnan( escript.inf(sigma_model) ) or \
305 cmath.isnan( escript.sup(sigma_model) ) or \
306 cmath.isinf( escript.sup(sigma_model) ):
307 raise ValueError("The conductivity model contains NaNs or INFs" )
308
309
310
311 # ---
312 # Print information:
313 # ---
314
315 print("")
316 print("="*72)
317 print("Escript MT2D, version", self.__version)
318 print("="*72)
319 print("Mesh XMin/XMax : ", xmin, xmax)
320 print("Mesh ZMin/ZMax : ", zmin, zmax)
321 print("Number of Tags : ", len( tags ))
322 print("Mapping : ", {True: 'Yes', False: 'No'}[maps is not None])
323 print("Conductivity Model : ", sigma_model)
324 print("Nr of Frequencies : ", len( frequencies ))
325 print("Start/End/Step (Hz) : ", freq_def["high"], freq_def["low"], freq_def["step"])
326 print("Mode : ", mode.upper())
327 print("Solver : ", MT_2D._solver)
328 print("Show plots : ", plot)
329 print("="*72)
330 print("")
331
332 if self._debug:
333 print("Mesh-Info : ")
334 print(domain.print_mesh_info(full=False))
335
336
337
338 # ---
339 # Set all required variables as data attributes
340 # ---
341
342 self.domain = domain
343 self.X = X
344 self.mode = mode
345 self.xmin = xmin
346 self.xmax = xmax
347 self.zmin = zmin
348 self.zmax = zmax
349 self.zstep = zstep
350 self.xstep = xstep
351 self.rho = rho
352 self.rho_1d = rho_1d
353 self.ifc_1d = ifc_1d
354 self.plot = plot
355 self.sigma = sigma_model
356 self.frequencies = frequencies
357 self.boundary_mask = boundary_mask
358
359 #__________________________________________________________________________________________________
360
361
362 def __interpolLinear(self,dx,x0,x1,y0,y1):
363 """
364 DESCRIPTION:
365 -----------
366 Function for simple 1D interpolation using the line-equation.
367
368 ARGUMENTS:
369 ----------
370 dx :: interpolation step.
371 x0 :: first coordinate point of known value y0.
372 x1 :: last coordinate point of known value y1.
373 y0 :: known value at first coordinate.
374 y1 :: known value at last coordinate.
375
376 RETURNS:
377 --------
378 y :: list with interpolated values
379 """
380 # Initialise return lists.
381 y = []
382
383 # Test for long enough interval.
384 if abs(x1-x0) <= dx: return y
385 # Test for correct abscissae.
386 if x0 >= x1: return y
387
388 x = x0
389 while x <= x1:
390 y.append( y0 + (y1-y0)*(x-x0)/(x1-x0) )
391 x = x + dx
392
393 return y
394 #__________________________________________________________________________________________________
395
396
397 def __getSamplePoints(self, min,max,step,constant=None):
398 """
399 DESCRIPTION:
400 -----------
401 Function to setup a list with sample points. If a
402 constant value was passed a 2D list is returned
403 where the second column is set to the constant.
404
405 ARGUMENTS:
406 ----------
407 min :: minimum value.
408 max :: maximum value.
409 step :: step value.
410 constant :: optional constant value for 2nd column.
411
412 RETURNS:
413 --------
414 sample :: list with samples.
415 """
416
417 # Initialise return list.
418 sample = []
419
420 # Cycle with step-size and fill sample list.
421 dp = min
422 while dp <= max:
423 if constant is not None:
424 sample.append([dp,constant])
425 else:
426 sample.append(dp)
427 # Increment the step.
428 dp = dp + step
429
430 # Return the list:
431 return sample
432 #______________________________________________________________________________________________
433
434
435 def __getSoundingFrequencies(self, frequencies):
436 """
437 DESCRIPTION:
438 -----------
439 Defines the sounding frequencies in Hz.
440
441 ARGUMENTS:
442 ----------
443 frequencies :: dictionary with frequency start/stop/step
444
445 RETURNS:
446 --------
447 sounding_frequencies :: list with frequency values
448 """
449 # Output list with frequencies in Hertz:
450 sounding_frequencies = []
451
452 # Period definition (from freq to time):
453 tme_1 = 1.0/frequencies["high"]
454 tme_n = 1.0/frequencies["low"]
455
456 # Number of points per decade:
457 tme_p = frequencies["step"]
458
459 # Number of periods in range:
460 nt = int(math.log10(tme_n/tme_1) * tme_p) + 1
461 # Fill list with times:
462 for n in xrange(nt):
463 # Sounding period in seconds:
464 period = tme_1*10**( (n)/float(tme_p))
465 # And store as frequency in Hertz:
466 sounding_frequencies.append( 1.0/period )
467
468 return sounding_frequencies
469 #__________________________________________________________________________________________________
470
471
472 def __getGradField(self, proj, mt2d_field, wm):
473 """
474 DESCRIPTION:
475 -----------
476 Calculates the complementary fields via Faraday's Law (TE-mode)
477 or via Ampere's Law (TM-mode). Partial derivative w.r.t. the
478 vertical coordinate are taken at the surface for which an Escript
479 'Projector' object is used to calculate the gradient.
480
481 ARGUMENTS:
482 ----------
483 proj :: escript Projection object
484 mt2d_field :: calculated magnetotelluric field
485 wm :: number with actual angular sounding frequency * mu
486
487 RETURNS:
488 --------
489 mt2d_grad :: dictionary with computed gradient fields
490 """
491
492 # Define the derived gradient fields:
493 if self.mode.upper() == 'TE':
494 # H = -(dE/dz) / iwm
495 grad_real =-proj.getValue( escript.grad(mt2d_field["imag"])/wm )
496 grad_imag = proj.getValue( escript.grad(mt2d_field["real"])/wm )
497 #<Note the coupled dependency on real/imaginary part>:
498 else:
499 # E = (dH/dz) / sigma
500 grad_real = proj.getValue( escript.grad(mt2d_field["real"])/self.sigma )
501 grad_imag = proj.getValue( escript.grad(mt2d_field["imag"])/self.sigma )
502 #<'sigma' is an Escript data-object and as such the division
503 # will use the tagged sigma values of the associated domains>
504
505
506 # And return as dictionary for real and imaginary parts:
507 mt2d_grad = {"real": grad_real[1], "imag":grad_imag[1] }
508 #<Note>: the derivative w.r.t. 'z' is used (i.e. '[1]').
509
510 return mt2d_grad
511 #__________________________________________________________________________________________________
512
513
514 def __tagDomain(self, domain, X, tags, rho, maps):
515 """
516 DESCRIPTION:
517 -----------
518 Defines the conductivity model. Conductivities of tagged regions can be mapped
519 via user-defined functions passed in 'maps'. If no user-defined functions are
520 passed, a constant value is applied as provided in list 'rho' for each region.
521 User-defined functions have 3 arguments: x-coordinate, z-coordinate, resistivity.
522
523 ARGUMENTS:
524 ----------
525 domain :: escript object of mesh
526 X :: escript object with all coordinates
527 tags :: list with domain tags
528 rho :: list with domain resistivities
529 maps :: list with user-defined resistivity mappings
530
531 RETURNS:
532 --------
533 sigma :: escript object of conductivity model
534
535 """
536 # Setup the conductivity structure (acts on elements and can be discontinuous).
537 sigma = escript.Scalar(0, escript.Function(domain))
538
539 # Setup conductivity domains.
540 for i in xrange( len(tags) ):
541
542 # Default: assign conductivity which is the inverse of resistivity:
543 m = 1.0/rho[i]
544
545 # Map a user-defined conductivity distribution if given:
546 if maps is not None:
547 # Guard against undefined elements:
548 if maps[i] is not None:
549 # Map the conductivity according to the defined functions:
550 m = maps[i]( X[0], X[1], rho[i] )
551
552 # Tag the mesh with the conductivity distributions at each iteration:
553 sigma += m * escript.insertTaggedValues(escript.Scalar(0,escript.Function(domain)),**{ tags[i] : 1})
554
555
556 if self._debug == True:
557 sigma.expand()
558 mydir = os.getcwd()
559 dbgfl = mydir + os.sep + "mt2d_sigma_dbg.silo"
560 print("")
561 print("DEBUG: writing SILO debug output of conductivity model:")
562 print(dbgfl)
563 print("")
564 weipa.saveSilo(dbgfl, sigma = sigma)
565
566
567 # All done:
568 return sigma
569 #__________________________________________________________________________________________________
570
571
572 def __setBoundaryMask(self, X):
573 """
574 DESCRIPTION:
575 -----------
576 Define Dirichlet model boundaries conditions.
577
578 ARGUMENTS:
579 ----------
580 X :: escript object with all coordinates
581
582 RETURNS:
583 --------
584 boundary_mask :: escript object with mask values at boundaries
585
586 """
587 # Boundaries are defined as masks (1 or 0) for all mesh coordinates;
588 # values at the boundary are '1', whereas all other values are '0'.
589 mask_l = escript.whereZero( X[0] - escript.inf(X[0]) )
590 mask_r = escript.whereZero( X[0] - escript.sup(X[0]) )
591 mask_t = escript.whereZero( X[1] - escript.inf(X[1]) )
592 mask_b = escript.whereZero( X[1] - escript.sup(X[1]) )
593
594 # Combine the mask for all boundaries:
595 boundary_mask = mask_t + mask_b + mask_l + mask_r
596
597 return boundary_mask
598 #<Note>: this boundary mask is used later on as PDE coefficient 'q'.
599 #__________________________________________________________________________________________________
600
601
602 def __getBoundaryValues(self, mode, X, rho_1d, ifc_1d, xstep, zstep, frequency):
603 """
604 DESCRIPTION:
605 -----------
606 Returns a list with boundary values along each Dirichlet boundary.
607 Values at the left and right side of the domain are evaluated at
608 sample points and interpolated across the domain. The subroutine
609 expects that layers at the right- and left-hand-side are defined.
610
611 ARGUMENTS:
612 ----------
613 mode :: string with TE or TM mode
614 X :: escript object with all coordinates
615 rho_1d :: dictionary with resistivities at the left/right boundary
616 ifc_1d :: dictionary with layer interfaces at the left/right boundary
617 xstep :: number with step size for horizontal sample list
618 zstep :: number with step size for vertical sample list
619 frequency :: number with actual sounding frequency
620
621 RETURNS:
622 --------
623 bondary_value :: dictionary with lists of boundary values at sample points
624 """
625
626 # ---
627 # Sample lists at vertical and horizontal boundaries.
628 # ---
629
630 # Horizontal extents:
631 xmin = escript.inf(X[0])
632 xmax = escript.sup(X[0])
633
634 # Vertical extents:
635 zmin = escript.inf(X[1])
636 zmax = escript.sup(X[1])
637
638 # Setup a list with sample points along the vertical mesh extent, bottom to top:
639 zsamples = self.__getSamplePoints(-zmax,-zmin,zstep)
640
641
642 # ---
643 # Calculate the 1D response at the left- and right-hand-side boundaries
644 # ---
645
646 # Instantiate an 'mt1d' object for the left- and right-hand-sides:
647 mt1d_left = mt1d.MT_1D( frequency, ifc_1d['left'] , rho_1d['left'] , zsamples )
648 mt1d_rght = mt1d.MT_1D( frequency, ifc_1d['right'], rho_1d['right'], zsamples )
649
650 # Compute the 1D field values at the sample nodes:
651 te1d_left, tm1d_left = mt1d_left.mt1d( )
652 te1d_rght, tm1d_rght = mt1d_rght.mt1d( )
653
654 # Distinguish TE and TM mode and save 1D values in dictionary:
655 if mode.upper() == "TE":
656 mt_1d = {"left":te1d_left, "right":te1d_rght}
657 else:
658 mt_1d = {"left":tm1d_left, "right":tm1d_rght}
659
660
661 # ---
662 # Interpolation across mesh.
663 # ---
664
665 # Now setup a 2D-table from left to right at each sampled depth for mesh-interpolation.
666 table2d_real = []
667 table2d_imag = []
668
669 # 1D-interpolation of values from left to right at different depths 'i':
670 for i in xrange( len(zsamples)):
671 table2d_real.append( self.__interpolLinear(xstep, xmin, xmax, mt_1d["left"].real[i], mt_1d["right"].real[i]) )
672 table2d_imag.append( self.__interpolLinear(xstep, xmin, xmax, mt_1d["left"].imag[i], mt_1d["right"].imag[i]) )
673
674 # 2D-interpolation to map the values on the mesh coordinates:
675 bondary_value_real = escript.interpolateTable( table2d_real, X, (xmin,zmin), (xstep,zstep) )
676 bondary_value_imag = escript.interpolateTable( table2d_imag, X, (xmin,zmin), (xstep,zstep) )
677
678 # Return the real and imaginary values as a dictionary:
679 boundary_value = {"real":bondary_value_real, "imag":bondary_value_imag}
680
681
682 return boundary_value
683 #__________________________________________________________________________________________________
684
685
686 def __getAppResPhase(self, mt2d_field, mt2d_grad, wm):
687 """
688 DESCRIPTION:
689 -----------
690 Calculates the apparent resistivity and phase.
691
692 ARGUMENTS:
693 ----------
694 mt2d_field :: dictionary with real/imaginary field values
695 mt2d_grad :: dictionary with real/imaginary gradient values
696
697 RETURNS:
698 --------
699 apparent resistivity and phase
700 """
701
702 # Define the associated modelled fields in readable variables:
703 if self.mode.upper() == 'TE':
704 # Transverse electric field:
705 Er = mt2d_field["real"]
706 Ei = mt2d_field["imag"]
707 Hr = mt2d_grad["real"]
708 Hi = mt2d_grad["imag"]
709 else:
710 # Transverse magnetic field :
711 Hr = mt2d_field["real"]
712 Hi = mt2d_field["imag"]
713 Er = mt2d_grad["real"]
714 Ei = mt2d_grad["imag"]
715
716
717 # Return apparent Resistivity and Phase:
718 arho_2d = ( (Er**2 + Ei**2)/(Hr**2 + Hi**2) ) / wm
719 aphi_2d = escript.atan( (Ei*Hr - Er*Hi)/(Er*Hr + Ei*Hi) ) * 180.0/cmath.pi
720
721 return arho_2d, aphi_2d
722 #__________________________________________________________________________________________________
723
724
725 def __showPlot(self, loc, rho_2d, phi_2d, f):
726 """
727 DESCRIPTION:
728 -----------
729 Plot of apparent resistivity and phase.
730
731 ARGUMENTS:
732 ----------
733 loc :: escript Locator object
734 rho_2d :: list with computed apparent resistivities
735 phi_2d :: list with computed phase values
736 f :: sounding frequency
737
738 RETURNS:
739 --------
740 Plot in window.
741
742 """
743
744 # Abscissas/Ordinates:
745 x = numpy.array( loc.getX() )[:,0]
746 y0 = numpy.array( loc.getValue(rho_2d) )
747 y1 = numpy.array( loc.getValue(phi_2d) )
748
749 # Plot labels:
750 title = 'Escript MT-2D ' + '(' + self.mode.upper() + ')' + ' freq: ' + str(f) + ' Hz'
751 ylbl0 = r'Apparent Resistivity $(\Omega\cdot\,m)$'
752 ylbl1 = r'Phase $(^{\circ})$'
753 xlbl1 = 'Easting (m)'
754
755
756 # Setup the plot window with app. res. on top and phase on bottom:
757 f, ax = plt.subplots(2, figsize=(8,8), sharex=True) # Mind shared axis
758 f.subplots_adjust(top=0.9) # Little extra space for 'suptitle'
759 f.suptitle(title) # This is actually the plot-title
760
761 # Top: apparent resistivity on semi-log plot
762 ax[0].semilogy(x,y0, color='red')
763 ax[0].grid(b=True, which='both', color='grey',linestyle=':')
764 ax[0].set_title( ylbl0 )
765
766 # Bottom: phase on linear plot
767 ax[1].plot(x,y1, color='blue')
768 ax[1].grid(b=True, which='both', color='grey',linestyle=':')
769 ax[1].set_xlabel( xlbl1 )
770 ax[1].set_title( ylbl1 )
771
772 plt.show()
773
774 #__________________________________________________________________________________________________
775
776
777 def __setSolver(self, mode, domain, sigma, boundary_mask, boundary_value, f):
778 """
779 DESCRIPTION:
780 -----------
781 Setups the coupled PDE for real and complex part.
782
783 ARGUMENTS:
784 ----------
785 mode :: string with TE or TM mode
786 domain :: escript object with mesh domain
787 sigma :: escript object with conductivity model
788 boundary_mask :: escript object with boundary mask
789 boundary_value :: dictionary with real/imag boundary values
790 f :: sounding frequency
791
792 RETURNS:
793 --------
794 mt2d_fields :: dictionary with solved PDE, magnetotelluric fields real/imag
795 """
796
797 # Constants:
798 pi = cmath.pi # Ratio circle circumference to diameter.
799 mu0 = 4*pi*1e-7 # Free space permeability in V.s/(A.m).
800 wm = 2*pi*f*mu0 # Angular frequency times mu0.
801
802
803 # ---
804 # Setup the coupled PDE for real/imaginary parts:
805 # ---
806
807 # Initialise the PDE object for two coupled equations (real/imaginary).
808 mtpde = pde.LinearPDE(domain, numEquations=2)
809
810 # If set, solve the 2D case using the direct solver:
811 if MT_2D._solver.upper() == "DIRECT":
812 mtpde.getSolverOptions().setSolverMethod(pde.SolverOptions().DIRECT)
813 else:
814 mtpde.getSolverOptions().setSolverMethod(pde.SolverOptions().DEFAULT)
815
816 # Now initialise the PDE coefficients 'A' and 'D',
817 # as well as the Dirichlet variables 'q' and 'r':
818 A = mtpde.createCoefficient("A")
819 D = mtpde.createCoefficient("D")
820 q = mtpde.createCoefficient("q")
821 r = mtpde.createCoefficient("r")
822
823 # Set the appropriate values for the coefficients depending on the mode:
824 if mode.upper() == "TE":
825 a_val = 1.0
826 d_val = wm*sigma
827 elif mode.upper() == "TM":
828 a_val = 1.0/sigma
829 d_val = wm
830
831
832 # ---
833 # Define the PDE parameters, mind boundary conditions.
834 # ---
835
836 # Now define the rank-4 coefficient A:
837 for i in range(domain.getDim()):
838 A[0,i,0,i] = a_val
839 A[1,i,1,i] = a_val
840
841 # And define the elements of 'D' which are decomposed into real/imaginary values:
842 D[0,0] = 0 ; D[1,0] = d_val
843 D[0,1] =-d_val ; D[1,1] = 0
844
845 # Set Dirichlet boundaries and values:
846 q[0] = boundary_mask ; r[0] = boundary_value['real']
847 q[1] = boundary_mask ; r[1] = boundary_value['imag']
848
849 # ---
850 # Solve the PDE
851 # ---
852
853 mtpde.setValue(A=A, D=D, q=q, r=r )
854 pde_solution = mtpde.getSolution()
855
856 # And return the real and imaginary parts individually:
857 mt2d_fields = {"real":pde_solution[0], "imag":pde_solution[1] }
858 #<Note>: the electric field is returned for TE-mode.
859 # the magnetic field is returned for TM-mode.
860
861 return mt2d_fields
862 #__________________________________________________________________________________________________
863
864
865 def pdeSolve(self):
866 """
867 DESCRIPTION:
868 -----------
869 Solves the PDE for either the TE or the TM mode.
870 (TE mode is the transverse Electric field).
871 (TM mode is the transverse Magnetic field).
872
873 ARGUMENTS:
874 ----------
875 (uses `self`)
876
877 RETURNS:
878 --------
879 mt2d :: list with real/imag fields for each sounding frequency
880 arho :: list with apparent resistivities for each sounding frequency
881 aphi :: list with phase values for each sounding frequency
882
883 """
884
885
886 # ---
887 # Constants.
888 # ---
889
890 # Pi & vacuum permeability:
891 pi = cmath.pi
892 mu = 4*pi*1e-7
893
894 # Number of frequencies:
895 nfreq = len(self.frequencies)
896
897
898 # ---
899 # Projector and Locator objects.
900 # ---
901
902 print("Setting up Escript Locator and Projector objects...")
903
904 # Setup a list with sample points along the vertical mesh extent, bottom to top:
905 xsample = self.__getSamplePoints(escript.inf(self.X[0]),escript.sup(self.X[0]),self.xstep, constant=0.0)
906
907 # Get the locations of mesh points at the surface via the Locator object
908 # operating on the continuous function-space (i.e. nodes) of the domain.
909 loc = pdetools.Locator(escript.ContinuousFunction(self.domain),xsample )
910
911 # Instantiate the Projector class with smoothing on (fast=False);
912 # the Projector is used to calculate the gradient correctly.
913 proj = pdetools.Projector(self.domain, reduce=False, fast=False)
914
915
916 # ---
917 # Solve the PDE for all frequencies.
918 # ---
919
920 # Prepare lists to store the values at each frequency:
921 arho = []
922 aphi = []
923 mt2d = []
924
925 # Cycle over all frequencies:
926 print("Solving for frequency: ...")
927 for n in xrange( nfreq ):
928
929 f = self.frequencies[n] # actual frequency (Hz)
930 wm = (2*pi*f)*mu # angular frequency (rad/s)
931 T = 1.0 / f # sounding period (s)
932
933 print(n+1,":", f, "(Hz)")
934
935 # Calculate 1D Dirichlet boundary values:
936 boundary_value = self.__getBoundaryValues(self.mode.upper(), self.X, self.rho_1d, self.ifc_1d, self.xstep, self.zstep, f)
937
938 # Solve the 2D-MT PDE:
939 fld_2d = self.__setSolver(self.mode.upper(),self.domain, self.sigma, self.boundary_mask, boundary_value, f)
940
941 # Calculate the field gradients:
942 grd_2d = self.__getGradField(proj, fld_2d, wm)
943
944 # Calculate the apparent resistivity and phase:
945 rho_2d, phi_2d = self.__getAppResPhase(fld_2d, grd_2d, wm)
946
947 # Save in lists for each frequency:
948 mt2d.append( fld_2d )
949 arho.append( rho_2d )
950 aphi.append( phi_2d )
951
952 # Optionally plot the apparent resistivity and phase:
953 if self.plot:
954 self.__showPlot(loc, rho_2d, phi_2d, f)
955
956
957 # ---
958 # All done
959 # ---
960
961 print("field calculations finished.")
962 return mt2d, arho, aphi
963 #__________________________________________________________________________________________________
964
965
966
967
968
969
970

  ViewVC Help
Powered by ViewVC 1.1.26