/[escript]/branches/subworld2/inversionSurveyDCResJ.py
ViewVC logotype

Contents of /branches/subworld2/inversionSurveyDCResJ.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 5530 - (show annotations)
Wed Mar 11 23:39:58 2015 UTC (3 years, 11 months ago) by jfenwick
File MIME type: text/x-python
File size: 10040 byte(s)
Finally making the modded script available.   Some other tweaks
1
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 """
18 this a very simple test for inversion of Dc resistivity Data
19 """
20
21 #For this hacked version, you will need to create a (non-version controlled)
22 #subdir called pklMsh with the mesh in it
23
24
25 __copyright__="""Copyright (c) 2003-2014 by University of Queensland
26 http://www.uq.edu.au
27 Primary Business: Queensland, Australia"""
28 __license__="""Licensed under the Open Software License version 3.0
29 http://www.opensource.org/licenses/osl-3.0.php"""
30 __url__="https://launchpad.net/escript-finley"
31
32 import os
33 import sys
34 from esys.escript import *
35 from esys.downunder import *
36 import esys.ripley as ripley
37 from esys.escript import unitsSI as U
38 from esys.weipa import saveSilo
39 import numpy as np
40 import esys.escript.pdetools as pdetools
41 import pickle
42
43 import logging
44 logger=logging.getLogger('inv')
45 #logger.setLevel(logging.INFO)
46 logger.setLevel(logging.DEBUG)
47
48 """
49 Inversion script for 3d dc resistivity data
50 It is highly recommended to use the same .msh file that was used for generating the
51 sythetic if you are using synthetic measurements.
52 """
53
54 def getApparentResistivity(delPhi, n, a, current):
55 resistivity=(-delPhi/current)*(n*(n+1))*pi*a
56 return resistivity
57
58 current=0.5 #500 milli amps
59 #load data to be inverted, the pkl file is to be generated using the forward
60 #modeling code adjust path to suit
61
62
63 pklFileName = "dcResResistiveLayer20-2-24-10.pkl"
64 pklFilePtr = open(pklFileName, "r")
65 pklFile = pickle.load(pklFilePtr)
66 #sourceInfo descibes the current electode setup, it containts pairs of tags as a list of tupples.
67 sourceInfo = pklFile["srcInfo"]
68 #VPSQ is a list of potential differences across different measurement pairs s and source pairs p and q
69 vpqs = pklFile["delPotArr"]
70 sampleInfo = pklFile["samples"]
71 #run name contains the name of the msh file to be used
72 runName = pklFile["runName"]
73 interval_a= pklFile["interval_a"]
74 numElectrodes=pklFile["numElectrodes"]
75 tags=[]
76 points=[]
77 electrodeDict = pklFile["electrodeDict"]
78 for i in range(numElectrodes):
79 tags.append("e%d"%i)
80 points.append(electrodeDict["e%d"%i][:-1]) #first three elements of electrode dict provide point info
81
82
83 sw=SplitWorld(getMPISizeWorld()) # Depending on how big the meshes are, we may want to be more efficient about number of worlds
84 buildDomains(sw,ReadGmsh, runName+".msh", 3 , diracTags=tags, diracPoints=points)
85
86
87 addVariable(sw,"regularization", makeLocalOnly)
88 addVariable(sw,"mappings", makeLocalOnly)
89 addVariable(sw,"fwdmodels", makeLocalOnly)
90
91 # model count is how many models this Job should create, modelnumstart is where in the global numbering should
92 # these jobs be?
93 # This imposes a constraint on the arity of functions accepted for this call
94 #Note that self, is needed because this will be executed in the context of a job and we want to be able to access it
95 def makeParamBlob(self, **kwargs):
96 dom=self.domain
97 DIM=dom.getDim()
98 z=ContinuousFunction(dom).getX()[1]
99 try:
100 modelcount=kwargs["rangestart"] #Where does our part of the range start
101 modelnumstart=kwargs["rangelen"] #How many things in our part of the range
102 nLS=kwargs["numLevelSets"]
103 except KeyError as e:
104 raise ValueError("The following expected value was not supplied to the function: "+e.message)
105
106 #====================== THIS IS WHERE THE INVERSION STARTS:
107
108 # create regularization with two level set functions:
109 z=Solution(dom).getX()[DIM-1]
110 x=Solution(dom).getX()[0]
111 # ask luts about this
112 reg_mask = whereZero(z-inf(z)) # fix the top resistivity only.
113 #====================== SETUP INITIAL VALUES
114 primaryRes=100.0
115 Sigma0 = Scalar(0,Function(dom))
116 SigmaPrimary = Scalar(0,Function(dom))
117 Sigma0.setTaggedValue("volume-1",1/primaryRes)
118 Sigma0.setTaggedValue("volume-2",1/primaryRes)
119 Sigma0.setTaggedValue("volume-3",1/primaryRes)
120 Sigma0.setTaggedValue("volume-4",1/primaryRes)
121 SigmaPrimary.setTaggedValue("volume-1",1/primaryRes)
122 SigmaPrimary.setTaggedValue("volume-2",1/primaryRes)
123 SigmaPrimary.setTaggedValue("volume-3",1/primaryRes)
124 SigmaPrimary.setTaggedValue("volume-4",1/primaryRes)
125 Sigma0.expand()
126 SigmaPrimary.expand()
127 #===================== SETUP SAMPLES
128 samples=[]
129 for i in range(len(sampleInfo)):
130 samples.append([electrodeDict[sampleInfo[i][0]],electrodeDict[sampleInfo[i][1]]])
131
132 #===================== SETUP PDE FOR PRIMARY POTENTIAL
133 pde=LinearPDE(dom, numEquations=1)
134 kro=kronecker(dom)
135 tol=1e-8
136 DIM=dom.getDim()
137 x = dom.getX()
138 q=whereZero(x[DIM-1]-inf(x[DIM-1]))
139 for i in xrange(DIM-1):
140 xi=x[i]
141 q+=whereZero(xi-inf(xi))+whereZero(xi-sup(xi))
142 A=pde.createCoefficient('A')
143 A=SigmaPrimary * kro
144 dcm=mappings.DcResMapping(2*Sigma0)
145
146
147 # iterate over source info
148 n=1
149 j=0
150
151 # Need to make sure that the numLevelSets here is the same as the one given as a parameter
152 # to the overall InversionCostFunction but don't want to pull values out of the subworld
153 regularization=Regularization(dom, numLevelSets=nLS,
154 w1=np.ones(DIM), # consider gradient terms
155 location_of_set_m=reg_mask)
156
157
158 # Now create a number of forward models to be processed on this subworld
159 forwardmodels=[] # this will hold the list of forward models
160 for k in range(0, modelcount):
161 i=modelnumstart+k
162 nCount = (numElectrodes-3) - (2*(n-1))
163 loc = pdetools.Locator(ContinuousFunction(dom), samples[i])
164 y_dirac=Scalar(0,DiracDeltaFunctions(dom))
165 y_dirac.setTaggedValue(sourceInfo[i][0],current)
166 y_dirac.setTaggedValue(sourceInfo[i][1],-current)
167 pde.getSolverOptions().setTolerance(tol)
168 pde.setSymmetryOn()
169 pde.setValue(A=A,y_dirac=y_dirac,q=q)
170 phiPrimary=pde.getSolution()
171 if (vpqs[n-1][j]!=0):
172 DCRES=DcRes(dom, loc, (vpqs[n-1][j],), [sampleInfo[i]], phiPrimary, SigmaPrimary, w=(1./vpqs[n-1][j])**2)
173 else:
174 DCRES=DcRes(dom, loc, (vpqs[n-1][j],), [sampleInfo[i]], phiPrimary, SigmaPrimary, w=(1.))
175 j+=1
176 if j == nCount:
177 n+=1
178 j=0
179 forwardmodels.append((DCRES, 0))
180 dcm=(mappings.DcResMapping(2*Sigma0),)
181
182 dcm, forwardmodels=InversionCostFunction.processMapsAndModels(dcm, forwardmodels, nLS)
183
184 print "dcm is",dcm
185
186
187 trafo = regularization.getCoordinateTransformation()
188 for m in forwardmodels:
189 print "--------------------------------"
190 print type(m[0])
191 print "--------------------------------"
192 if not m[0].getCoordinateTransformation() == trafo:
193 raise ValueError("Coordinate transformation for regularization and model don't match.")
194
195
196 muVal=5000
197 mu_c=None
198 #cf.setTradeOffFactorsRegularization(mu=muVal, mu_c=None)
199 regularization.setTradeOffFactorsForVariation(muVal)
200 regularization.setTradeOffFactorsForCrossGradient(mu_c)
201
202
203 # Need to export forward models as well
204 self.exportValue("fwdmodels", forwardmodels)
205 # Need to export regularization
206 self.exportValue("regularization",regularization)
207
208 # Need to export dcm
209
210 self.exportValue("mappings", dcm)
211
212
213
214 # finally we can set up the cost function:
215 # Need to pass in subworld and world init function here
216 # Also need to modify so we don't take in regularization or mappings .... or forward models
217 #cf=InversionCostFunction(regularization, (dcm,), tuple(forwardmodels) , splitw=sw, worldsinit_fn=makeParamBlob)
218 cf=InversionCostFunction(None, None, None , splitw=sw, worldsinit_fn=makeParamBlob, numLevelSets=1,
219 numModels=len(sourceInfo), numMappings=1)
220
221 #============== set regularization
222 #muVal=5000
223 #cf.setTradeOffFactorsRegularization(mu=muVal, mu_c=None)
224
225 # this code can be used to save silo at every step
226 #def minimizerCallBack(k, x, Jx, g_Jxx):
227 # saveSilo("/scratch/jplessis/mtinversion/silos/%d-%d.silo"%(siloCount,k),est=x,sigma=cf.getProperties(x))
228
229 #==================================== setup and run solver
230 # and then start running the show:
231 solver=MinimizerLBFGS()
232 solver.setCostFunction(cf)
233 # solver.setCallback(minimizerCallBack)
234 solver.setTolerance(1e-6)
235 solver.setMaxIterations(180)
236
237 try:
238 solver.run(Scalar(0, Solution(dom)))
239 except (MinimizerMaxIterReached ,RuntimeError) as exceptionVal:
240 print "Caught exception saving silo and rethrowing."
241 print "silo save to,","/scratch/uqjdupl1/dcres/silos/sigmaResistiveLayerfailed%d-%d.silo"%(numElectrodes, muVal)
242 m=solver.getResult()
243 sigma = cf.getProperties(m)
244 saveSilo("/scratch/uqjdupl1/dcres/silos/sigmaResistiveLayerfailed%d-%d-%d.silo"%(interval_a, numElectrodes ,muVal), sigma=sigma, con=1./sigma)
245 raise exceptionVal
246
247 m=solver.getResult()
248 sigma = cf.getProperties(m)
249 saveSilo("/scratch/uqjdupl1/dcres/silos/sigmaResistiveLayer%d-%d-%d.silo"%(interval_a, numElectrodes ,muVal), sigma=sigma, con=1./sigma)
250 print ("output silo in:","silos/sigmaResistiveLayer%d-%d-%d.silo"%(interval_a, numElectrodes ,muVal))
251 resIn=[]
252 resFinal=[]
253 n=1
254 j=0
255 for i in forwardmodels:
256 nCount = (numElectrodes-3) - (2*(n-1))
257 phi, loc_phi=i[0].getArguments(sigma)
258 delPot=loc_phi[1]-loc_phi[0]
259 delPhiIn = (i[0].delphiIn)[0]
260 resIn.append(primaryRes + getApparentResistivity(delPhiIn, n, interval_a, current))
261 resFinal.append(primaryRes + getApparentResistivity(delPot, n, interval_a, current))
262 j+=1
263 if j == nCount:
264 n+=1
265 j=0
266
267 pklF=open(pklFileName[:-4]+"-%d-inverted.pkl"%muVal,"w")
268 print ("output in:",pklFileName[:-4]+"-%d-inverted.pkl"%muVal)
269 pklDict={"resIn":resIn,"resFinal":resFinal,"regVal":muVal}
270 pickle.dump(pklDict,pklF)
271 print "m=",m
272 print "sigma=", sigma

  ViewVC Help
Powered by ViewVC 1.1.26