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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 5707 - (show annotations)
Mon Jun 29 03:59:06 2015 UTC (3 years, 8 months ago) by sshaw
File MIME type: text/x-python
File size: 4197 byte(s)
adding copyright headers to files without copyright info, moved header to top of file in some cases where it wasn't
1
2 ##############################################################################
3 #
4 # Copyright (c) 2003-2015 by The 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 # -------------------------------------------------------------------------------------------------
18 # DESCRIPTION:
19 # ------------
20 # Script to calculate the DC Resistivity response for a half-sphere
21 # embedded in a half-space (see example from Rucker et al (2006), Figure 4).
22 #
23 # REFERENCES:
24 # -----------
25 # "Three-dimensional modelling and inversion of DC resistivity data incorporating
26 # topography -- I. Modelling", Geophysical Journal International (2006) 166, 495-505
27 # Carsten Rucker, Thomas Gunther and Klaus Spitzer
28 # -------------------------------------------------------------------------------------------------
29
30 from __future__ import division, print_function
31
32 import esys.finley as finley
33 import esys.escript as escript
34 from esys.downunder import *
35 import math
36
37
38 # -------------------------------------------------------------------------------------------------
39 # Constants
40 # -------------------------------------------------------------------------------------------------
41
42 pi = math.pi
43
44 # -------------------------------------------------------------------------------------------------
45 # Input
46 # -------------------------------------------------------------------------------------------------
47
48 mesh_file = "data/HalfSphere_v1.4.msh"
49
50 # Tag volume names and conductivity values (Sm/m) for primary and secondary potential:
51 tag_p = {"domain" : 1/10.0, "sphere" : 1/10.0} # Primary (homogeneous).
52 tag_s = {"domain" : 1/10.0, "sphere" : 1/1.0 } # Secondary.
53
54 xe_0 = -5.0 # start X-coordinate
55 numEle = 21 # number of electrodes
56 a = 0.5 # step size
57 n=9 #
58 midPoint = [xe_0 + (((numEle-1)*a)/2.), 0, 0]
59 current = 1.0 # (Ampere)
60 domain = finley.ReadGmsh(mesh_file, 3)
61 mesh_tags = escript.getTagNames(domain)
62 directionVector = [1,0]
63
64 # -------------------------------------------------------------------------------------------------
65 # Define the primary and secondary resistivity model.
66 # -------------------------------------------------------------------------------------------------
67
68 #<Note>: the mesh was created with domains, which were all 'tagged' with unique id-s;
69 # based on the tagged domain id-s, the conductivity values are assigned to the domains.
70 #<Note>: the primary conductivity is the conductivity in the vicinity of the electrodes;
71 # the secondary conductivity is defined as the difference between primary and the
72 # actual input conductivity.
73
74 # Setup primary and secondary conductivity objects:
75 sig_p = escript.Scalar(0,escript.ContinuousFunction(domain))
76 sig_s = escript.Scalar(0,escript.ContinuousFunction(domain))
77
78 # Cycle tag_values dictionary and assign conductivity values;
79 # this is done for both, primary and secondary, potentials as
80 # both dictionaries are setup with identical dictionary-keys:
81 for tag in tag_p:
82 # All initially defined tags must be in the tag list.
83 # Print an error if it doesn't and exit the program.
84 if tag in mesh_tags:
85 # Assign value:
86 sig_p.setTaggedValue( tag, tag_p[tag] )
87 sig_s.setTaggedValue( tag, tag_s[tag] )
88 else:
89 print("Error: the defined tag is not defined in the mesh: " & tag)
90 sys.exit()
91
92 # Expand the data objects for output.
93 sig_p.expand()
94 sig_s.expand()
95
96
97 schs=SchlumbergerSurvey(domain, sig_p, sig_s, current, a, n, midPoint, directionVector, numEle)
98 pot=schs.getPotentialAnalytic()
99 primaryApparentRes=schs.getApparentResistivity(pot[0])
100 SecondaryApparentRes=schs.getApparentResistivity(pot[1])
101 totalApparentRes=schs.getApparentResistivity(pot[2])
102
103
104 n=1
105 print ("Total:\n")
106 for i in totalApparentRes:
107 print ("n = %d:"%n)
108 print (i,"\n")
109 n=n+1

  ViewVC Help
Powered by ViewVC 1.1.26