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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6651 - (show annotations)
Wed Feb 7 02:12:08 2018 UTC (14 months, 2 weeks ago) by jfenwick
File MIME type: text/x-python
File size: 27149 byte(s)
Make everyone sad by touching all the files

Copyright dates update

1 ##############################################################################
2 #
3 # Copyright (c) 2003-2018 by The University of Queensland
4 # http://www.uq.edu.au
5 #
6 # Primary Business: Queensland, Australia
7 # Licensed under the Apache License, version 2.0
8 # http://www.apache.org/licenses/LICENSE-2.0
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 __copyright__="""Copyright (c) 2003-2018 by The University of Queensland
19 http://www.uq.edu.au
20 Primary Business: Queensland, Australia"""
21 __license__="""Licensed under the Apache License, version 2.0
22 http://www.apache.org/licenses/LICENSE-2.0"""
23 __url__="https://launchpad.net/escript-finley"
24
25 import logging
26 import esys.escriptcore.utestselect as unittest
27 from esys.escriptcore.testing import *
28 import numpy as np
29 import os
30 import sys
31 import cmath
32 from esys.downunder import *
33 from esys.escript import unitsSI as U
34 from esys.escript import *
35 from esys.weipa import saveSilo
36 from esys.escript.linearPDEs import LinearSinglePDE, LinearPDE
37 from esys.escript.pdetools import Locator
38
39 try:
40 from esys.ripley import Rectangle as ripRectangle, Brick as ripBrick
41 HAVE_RIPLEY = True
42 except ImportError as e:
43 HAVE_RIPLEY = False
44
45 try:
46 from esys.finley import Rectangle as finRectangle
47 HAVE_FINLEY = True
48 except ImportError as e:
49 HAVE_FINLEY = False
50
51 mpisize = getMPISizeWorld()
52 # this is mainly to avoid warning messages
53 logging.basicConfig(format='%(name)s: %(message)s', level=logging.INFO)
54
55 HAVE_DIRECT = hasFeature("PASO_DIRECT") or hasFeature('trilinos')
56
57 @unittest.skipUnless(HAVE_RIPLEY, "Ripley module not available")
58 @unittest.skipUnless(HAVE_DIRECT, "more than 1 MPI rank or missing direct solver")
59 class TestAcousticInversion(unittest.TestCase):
60 def test_API(self):
61 NE=20
62 for x in [int(sqrt(mpisize)),2,3,5,7,1]:
63 NX=x
64 NY=mpisize//x
65 if NX*NY == mpisize:
66 break
67 domain=ripRectangle(n0=NE*NX-1, n1=NE*NY-1, l0=1., l1=1., d0=NX, d1=NY, diracPoints=[(0.5,1.)], diracTags=['sss'])
68
69 omega=2.
70
71 data=Data([1,2], FunctionOnBoundary(domain))
72 F=Data([2,3], Function(domain))
73 w=1.
74 self.assertRaises(ValueError, AcousticWaveForm, domain, omega, w, data, 1.) # F is a scalar
75 self.assertRaises(ValueError, AcousticWaveForm, domain, omega, w, [1,2], F) # data is not Data
76 self.assertRaises(ValueError, AcousticWaveForm, domain, omega, w, Data([1,2], Function(domain)), F) # data is not on boundary
77 self.assertRaises(ValueError, AcousticWaveForm, domain, omega, w, Scalar(1, Function(domain)), F) # data is not of shape (2,)
78 self.assertRaises(ValueError, AcousticWaveForm, domain, omega, [1,2], data, F) # w is not a scalar
79 self.assertRaises(ValueError, AcousticWaveForm, domain, omega, Scalar(1, Function(domain)), data, F) # w is not a scalar
80
81 # now we do a real one
82 acw=AcousticWaveForm(domain, omega, w, data, F)
83 self.assertEqual(acw.getDomain(), domain)
84 pde=acw.setUpPDE()
85 self.assertIsInstance(pde, LinearPDE)
86 self.assertEqual(pde.getNumEquations(), 2)
87 self.assertEqual(pde.getNumSolutions(), 2)
88 self.assertEqual(pde.getDomain(), domain)
89
90 def test_numeric2DscaleF(self):
91 NE=40
92 for x in [int(sqrt(mpisize)),2,3,5,7,1]:
93 NX=x
94 NY=mpisize//x
95 if NX*NY == mpisize:
96 break
97 domain=ripRectangle(n0=NE*NX-1, n1=NE*NY-1, l0=1., l1=1., d0=NX, d1=NY, diracPoints=[(0.5,1.)], diracTags=['sss'])
98 #domain=ripRectangle(100,100, diracPoints=[(0.5,1.)], diracTags=['sss'])
99 omega=2.
100
101 # test solution is u = a * z where a is complex
102 a=complex(3.45, 0.56)
103 sigma=complex(1e-3, 0.056)
104
105 data=Data([a.real, a.imag], FunctionOnBoundary(domain))
106 mydata=data.copy()
107
108 z=FunctionOnBoundary(domain).getX()[1]
109 w=whereZero(z-1.)
110 # source:
111 F=Data( [1,0],Function(domain))
112
113 acw=AcousticWaveForm(domain, omega, w, data, F, coordinates=None, fixAtBottom=False, tol=1e-8, saveMemory=True, scaleF=True)
114 # check rescaled data
115 surv=acw.getSurvey()
116 self.assertAlmostEqual( integrate(length(surv[0])**2 * surv[1]), 1.)
117
118 mydata_scale=sqrt( integrate(w*length(mydata)**2) )
119 self.assertAlmostEqual( acw.getSourceScaling(z*[1, 0.]) , a/mydata_scale )
120 self.assertAlmostEqual( acw.getSourceScaling(mydata) , 1./mydata_scale )
121
122 # this should be zero:
123 sigma_comps=[sigma.real, sigma.imag]
124 args=acw.getArguments(sigma_comps)
125 d=acw.getDefect(sigma_comps, *args)
126 self.assertTrue(isinstance(d, float))
127 self.assertLess(abs(d), 1e-10)
128
129 dg=acw.getGradient(sigma_comps, *args)
130 self.assertTrue(isinstance(dg, Data))
131 self.assertTrue(dg.getShape()==(2,))
132 self.assertTrue(dg.getFunctionSpace()==Solution(domain))
133 self.assertLess(Lsup(dg), 1e-10)
134
135 # this shuld be zero' too
136 sigma_comps=[2*sigma.real, sigma.imag/2.]
137 args=acw.getArguments(sigma_comps)
138 d=acw.getDefect(sigma_comps, *args)
139 self.assertTrue(isinstance(d, float))
140 self.assertLess(abs(d), 1e-10)
141
142 dg=acw.getGradient(sigma_comps, *args)
143 self.assertTrue(isinstance(dg, Data))
144 self.assertTrue(dg.getShape()==(2,))
145 self.assertTrue(dg.getFunctionSpace()==Solution(domain))
146 self.assertLess(Lsup(dg), 1e-10)
147
148 # this shouldn't be zero:
149 sigma0=[2*sigma.real, 10*a.imag]*(27*Function(domain).getX()[0]-Function(domain).getX()[1])
150 args=acw.getArguments(sigma0)
151 d0=acw.getDefect(sigma0, *args)
152 self.assertTrue(isinstance(d0, float))
153 self.assertTrue(d0 >= 0)
154 self.assertTrue(d0 > 1e-10)
155
156 dg0=acw.getGradient(sigma0, *args)
157 self.assertTrue(isinstance(dg0, Data))
158 self.assertTrue(dg0.getShape()==(2,))
159 self.assertTrue(dg0.getFunctionSpace()==Solution(domain))
160 self.assertTrue(Lsup(dg0) > 1e-10)
161
162 # test the gradient numerrically:
163 h=0.002
164 X=Function(domain).getX()
165 # .. increment:
166 p=h*exp(-(length(X-[0.6,0.6])/10)**2)*Lsup(length(sigma0))
167
168 sigma1=sigma0+p*[1,0]
169 args=acw.getArguments(sigma1)
170 d1=acw.getDefect(sigma1, *args)
171 self.assertLess( abs( d1-d0-integrate(dg0[0]*p) ), 1e-2*abs(d1-d0) )
172
173 sigma2=sigma0+p*[0,1]
174 args=acw.getArguments(sigma2)
175 d2=acw.getDefect(sigma2, *args)
176 self.assertLess( abs(d2-d0-integrate(dg0[1]*p)), 1e-2*abs(d2-d0) )
177
178 def test_numeric2DnoscaleF(self):
179 domain=ripRectangle(n0=10, n1=20, l0=1., l1=1., diracPoints=[(0.5,1.)], diracTags=['sss'])
180 omega=1.5
181
182 # test solution is u = a * z where a is complex
183 a=complex(3.45, 0.56)
184 sigma=complex(1e-3, 0.056)
185
186 data=Data([a.real, a.imag], FunctionOnBoundary(domain))
187 z=FunctionOnBoundary(domain).getX()[1]
188 w=whereZero(z-1.)
189 # F = - a*omega* sigma
190 F=Data( [-(a*omega**2*sigma).real, -(a*omega**2*sigma).imag ],Function(domain))
191
192 acw=AcousticWaveForm(domain, omega, w, data, F, coordinates=None, fixAtBottom=False, tol=1e-8, saveMemory=True, scaleF=False)
193 # this should be zero:
194 sigma_comps=[sigma.real, sigma.imag]
195 args=acw.getArguments(sigma_comps)
196 d=acw.getDefect(sigma_comps, *args)
197 self.assertTrue(isinstance(d, float))
198 self.assertLess(Lsup(d), 1e-10)
199 #self.assertTrue(d >= 0)
200
201 dg=acw.getGradient(sigma_comps, *args)
202
203 self.assertTrue(isinstance(dg, Data))
204 self.assertTrue(dg.getShape()==(2,))
205 self.assertTrue(dg.getFunctionSpace()==Solution(domain))
206 self.assertLess(Lsup(dg), 5e-10)
207 # this shouldn't be zero:
208 sigma0=Data([2*sigma.real, sigma.imag/2], Function(domain) )
209 args=acw.getArguments(sigma0)
210 d0=acw.getDefect(sigma0, *args)
211 self.assertTrue(isinstance(d0, float))
212 self.assertGreaterEqual(d0, 0)
213 self.assertGreater(d0, 1e-10)
214
215 dg0=acw.getGradient(sigma0, *args)
216 self.assertTrue(isinstance(dg0, Data))
217 self.assertTrue(dg0.getShape()==(2,))
218 self.assertTrue(dg0.getFunctionSpace()==Solution(domain))
219 self.assertTrue(Lsup(dg0) > 1e-10)
220 # test the gradient numerically:
221 h=0.001
222 X=Function(domain).getX()
223 p=h*sin(length(X)*np.pi)*Lsup(length(sigma0))
224
225 sigma1=sigma0+p*[1,0]
226 args=acw.getArguments(sigma1)
227 d1=acw.getDefect(sigma1, *args)
228
229 self.assertLess( abs( d1-d0-integrate(dg0[0]*p) ), 1e-2*abs(d1-d0) )
230
231 sigma2=sigma0+p*[0,1]
232 args=acw.getArguments(sigma2)
233 d2=acw.getDefect(sigma2, *args)
234 self.assertLess( abs(d2-d0-integrate(dg0[1]*p)), 1e-2*abs(d2-d0) )
235
236
237 @unittest.skipIf(not HAVE_RIPLEY, "Ripley module not available")
238 class TestSubsidence(unittest.TestCase):
239 def test_PDE(self):
240 lam=2.
241 mu=1.
242
243 domain=ripBrick(20,20,max(19,2*mpisize-1), d2=mpisize)
244
245 xb=FunctionOnBoundary(domain).getX()
246 m=whereZero(xb[2]-1)
247 w=m*[0,0,1]
248 d=m*2.5
249 acw=Subsidence(domain, w,d, lam, mu )
250
251 P0=10.
252 args0=acw.getArguments(P0)
253 u=args0[0]
254 self.assertLess(Lsup(u[0]), 1.e-8)
255 self.assertLess(Lsup(u[1]), 1.e-8)
256 self.assertLess(Lsup(u[2]-2.5*domain.getX()[2]), 1.e-8)
257
258 dd=acw.getDefect(P0, *args0)
259
260 self.assertTrue( dd >= 0.)
261 self.assertTrue( dd <= 1e-7 * 2.5 )
262
263 def test_Differential(self):
264 lam=2.
265 mu=1.
266 INC=0.01
267 domain=ripBrick(20,20,min(99,20*mpisize-1) , d2=mpisize)
268
269 xb=FunctionOnBoundary(domain).getX()
270 m=whereZero(xb[2]-1)
271 w=m*[0,0,1]
272 d=m*2.5
273 acw=Subsidence(domain, w,d, lam, mu )
274
275 x=Function(domain).getX()
276 P0=x[0]*x[1]
277 args0=acw.getArguments(P0)
278 d0=acw.getDefect(P0, *args0)
279 grad_d=acw.getGradient(P0, *args0)
280
281 dP=exp(-(length(x-[0.5,0.5,0.5])/0.06)**2)
282 P1=P0+INC*dP
283 args1=acw.getArguments(P1)
284 d1=acw.getDefect(P1, *args1)
285 ref=abs((d1-d0)/INC)
286 self.assertLess(abs((d1-d0)/INC-integrate(grad_d* dP)), ref * 1.e-5)
287
288 dP=exp(-(length(x-[0.3,0.3,0.5])/0.06)**2)
289 P2=P0-INC*dP
290 args2=acw.getArguments(P2)
291 d2=acw.getDefect(P2, *args2)
292 ref=abs((d2-d0)/INC)
293 self.assertLess(abs((d2-d0)/INC+integrate(grad_d* dP)), ref * 1.e-5)
294
295 @unittest.skipUnless(HAVE_FINLEY, "Finley module not available")
296 class TestDCResistivity(unittest.TestCase):
297
298 def test_PDE2D(self):
299 dx_tests=0.1
300 sigma0=1.
301 electrodes=[(0.5-2*dx_tests,1.), (0.5-dx_tests,1.), (0.5+dx_tests,1.), (0.5+2*dx_tests,1.)]
302 domain=finRectangle(20,20, d1=mpisize, diracPoints=electrodes, diracTags=["sl0", "sl1", "sr0", "sr1"] )
303 loc=Locator(domain,electrodes[2:])
304
305 # this creates some reference Data:
306 x=domain.getX()
307 q=whereZero(x[0]-inf(x[0]))+whereZero(x[0]-sup(x[0]))+whereZero(x[1]-inf(x[1]))
308 ppde=LinearPDE(domain, numEquations=1)
309 s=Scalar(0.,DiracDeltaFunctions(domain))
310 s.setTaggedValue("sl0" ,1.)
311 s.setTaggedValue("sl1",-1.)
312 ppde.setValue(A=kronecker(2)*sigma0, q=q, y_dirac=s)
313 pp=ppde.getSolution()
314 uu=loc(pp)
315
316 # arguments for DcRes
317 current = 10.
318 sourceInfo = [ "sl0", "sl1" ]
319 sampleTags = [ ("sr0", "sr1") ]
320
321 sigmaPrimary=7.
322 phiPrimary=pp*current*sigma0/sigmaPrimary
323
324 uuscale=1-current*sigma0/sigmaPrimary
325 delphi_in = [ (uu[1]-uu[0]) * uuscale]
326
327 acw=DcRes(domain, loc, delphi_in, sampleTags, phiPrimary, sigmaPrimary)
328
329 self.assertLess(Lsup(phiPrimary-acw.getPrimaryPotential()), 1.e-10 * Lsup(acw.getPrimaryPotential()))
330
331 SIGMA=10. # matches current
332 args0=acw.getArguments(SIGMA)
333 p=args0[0]
334 u=args0[1]
335
336 # true secondary potential
337 pps=pp-phiPrimary
338 self.assertLess(Lsup(p-pps), 1.e-6*Lsup(pps))
339
340 # test return values at electrodes:
341 self.assertLess(abs(u[0]-uu[0]*uuscale), 1.e-6 * abs(uu[0]*uuscale))
342 self.assertLess(abs(u[1]-uu[1]*uuscale), 1.e-6 * abs(uu[1]*uuscale))
343
344 # this sould be zero
345 dd=acw.getDefect(SIGMA, *args0)
346 self.assertTrue( dd >= 0.)
347 self.assertTrue( dd <= 1e-7 )
348
349 def test_Differential2D(self):
350 INC=0.001
351 sigma0=1.
352 dx_tests=0.1
353 electrodes=[(0.5-2*dx_tests,1.), (0.5-dx_tests,1.), (0.5+dx_tests,1.), (0.5+2*dx_tests,1.)]
354 domain=finRectangle(20,20, d1=mpisize, diracPoints=electrodes, diracTags=["sl0", "sl1", "sr0", "sr1"] )
355 loc=Locator(domain,electrodes[2:])
356
357 # arguments for DcRes
358 #current = 10.
359 sampleTags = [ ("sr0", "sr1") ]
360
361 delphi_in = [ 0.05 ]
362
363 sigmaPrimary=1
364 x=domain.getX()
365 phiPrimary=(x[0]-inf(x[0]))*(x[1]-inf(x[1]))*(x[0]-sup(x[0]))
366
367 acw=DcRes(domain, loc, delphi_in, sampleTags, phiPrimary, sigmaPrimary)
368
369 #=====================================================================
370 x=Function(domain).getX()
371 SIGMA0=x[0]*x[1]+1
372 args0=acw.getArguments(SIGMA0)
373 d0=acw.getDefect(SIGMA0, *args0)
374 grad_d=acw.getGradient(SIGMA0, *args0)
375
376 dS=exp(-(length(x-[0.5,0.5])/0.2)**2)
377 SIGMA1=SIGMA0+INC*dS
378 args1=acw.getArguments(SIGMA1)
379 d1=acw.getDefect(SIGMA1, *args1)
380 ref=abs((d1-d0)/INC)
381 self.assertLess(abs((d1-d0)/INC-integrate(grad_d* dS)), ref * 1.e-3)
382
383 dS=-exp(-(length(x-[0.5,0.5])/0.2)**2)
384 SIGMA2=SIGMA0+INC*dS
385 args2=acw.getArguments(SIGMA2)
386 d2=acw.getDefect(SIGMA2, *args2)
387 ref=abs((d2-d0)/INC)
388 self.assertLess(abs((d2-d0)/INC-integrate(grad_d* dS)), ref * 1.e-3)
389
390 dS=-1
391 SIGMA3=SIGMA0+INC*dS
392 args3=acw.getArguments(SIGMA3)
393 d3=acw.getDefect(SIGMA3, *args3)
394 ref=abs((d3-d0)/INC)
395 self.assertLess(abs((d3-d0)/INC-integrate(grad_d* dS)), ref * 1.e-3)
396
397 dS=1
398 SIGMA4=SIGMA0+INC*dS
399 args4=acw.getArguments(SIGMA4)
400 d4=acw.getDefect(SIGMA4, *args4)
401 ref=abs((d4-d0)/INC)
402 self.assertLess(abs((d4-d0)/INC-integrate(grad_d* dS)), ref * 1.e-3)
403
404 class TestIsostaticPressure(unittest.TestCase):
405 @unittest.skipIf(not HAVE_RIPLEY, "Ripley module not available")
406 def test_all(self):
407 domain=ripBrick(50,50,20*mpisize-1, d2=mpisize)
408
409 ps=IsostaticPressure(domain, level0=1., coordinates=None)
410 g=Vector(0., Function(domain))
411 rho=Scalar(100, Function(domain))
412 p0=ps.getPressure(g, rho)
413 p_ref=-(1.-domain.getX()[2])*981.
414 self.assertLess(Lsup(p0-p_ref), 1e-6 * Lsup(p_ref))
415
416 g=Vector([0,0,-10], Function(domain))
417 rho=Scalar(0, Function(domain))
418 p0=ps.getPressure(g, rho)
419 p_ref=-(1.-domain.getX()[2])*26700
420 self.assertLess(Lsup(p0-p_ref), 1e-6 * Lsup(p_ref))
421
422 g=Vector([0,0,-10], Function(domain))
423 rho=Scalar(100, Function(domain))
424 p0=ps.getPressure(g, rho)
425 p_ref=-(1.-domain.getX()[2])*(981.+26700+1000)
426 self.assertLess(Lsup(p0-p_ref), 1e-6 * Lsup(p_ref))
427
428 @unittest.skipUnless(HAVE_RIPLEY, "Ripley module not available")
429 @unittest.skipUnless(HAVE_DIRECT, "more than 1 MPI rank or missing direct solver")
430 class TestMT2DModelTEMode(unittest.TestCase):
431 def test_API(self):
432 domain=ripRectangle(25, 25, d1=mpisize)
433 omega=2.
434 SIGMA=15.
435
436 x=[ [0.2,0.5], [0.3,0.5] ]
437 Z_XY=[ complex(1.2,1.5), complex(1.3,2.5) ]
438
439 eta=1.
440 w0=1.
441 # now we do a real one
442 model=MT2DModelTEMode(domain, omega, x, Z_XY, eta, sigma0=SIGMA, w0=w0)
443 self.assertEqual(model.getDomain(), domain)
444 pde=model.setUpPDE()
445 self.assertIsInstance(pde, LinearPDE)
446 self.assertEqual(pde.getNumEquations(), 2)
447 self.assertEqual(pde.getNumSolutions(), 2)
448 self.assertEqual(pde.getDomain(), domain)
449
450 # other things that should work
451 model=MT2DModelTEMode(domain, omega, x, Z_XY, eta=None, w0=[2.,3.] , sigma0=SIGMA, airLayerLevel=1.)
452
453 # these shouldn't work
454 self.assertRaises(ValueError, MT2DModelTEMode, domain, omega, x, [3.], eta=[1.,1.], w0=[2.,3.], sigma0=SIGMA)
455 self.assertRaises(ValueError, MT2DModelTEMode, domain, omega, x, Z_XY, eta=[1.], w0=[2.,3.] , sigma0=SIGMA)
456 self.assertRaises(ValueError, MT2DModelTEMode, domain, omega, [(6.7,5)], Z_XY, eta=[1.,1.], w0=[2.,3.] , sigma0=SIGMA)
457
458 def test_PDE(self):
459 omega=1.
460 mu0=0.123
461 SIGMA=15.
462 k=cmath.sqrt(1j*omega*mu0*SIGMA) # Ex=exp(k*z)
463 NE=101
464 domain=ripRectangle(max(NE,30*mpisize-1),max(NE,30*mpisize-1), d1=mpisize)
465
466 Z0=0.5
467 H=1./NE
468 X1=(int(0.3/H)+0.5)*H
469 X2=(int(0.6/H)+0.5)*H
470
471
472 IMP=-(1j*omega*mu0)/k*cmath.exp(k*Z0)/cmath.exp(k*Z0)
473 Z_XY=[ IMP, IMP ]
474
475 x=[ [X1,Z0], [X2,Z0] ]
476 eta=None
477
478 z=domain.getX()[1]
479 Ex0_ex=cos(k.imag*(z-1))*exp(k.real*(z-1))
480 Ex1_ex=sin(k.imag*(z-1))*exp(k.real*(z-1))
481 Ex0_ex_z=-sin(k.imag*(z-1))*k.imag*exp(k.real*(z-1))+cos(k.imag*(z-1))*exp(k.real*(z-1))*k.real
482 Ex1_ex_z=cos(k.imag*(z-1))*k.imag*exp(k.real*(z-1))+sin(k.imag*(z-1))*exp(k.real*(z-1))*k.real
483
484 model=MT2DModelTEMode(domain, omega, x, Z_XY, eta, mu=mu0, tol=1e-9, directSolver=True, sigma0=SIGMA)
485
486 args=model.getArguments(SIGMA)
487 Ex=args[0]
488 Exz=args[1]
489 self.assertLess(Lsup(Ex[0]-Ex0_ex), 1e-4 * Lsup(Ex0_ex))
490 self.assertLess(Lsup(Ex[1]-Ex1_ex), 1e-4 * Lsup(Ex1_ex))
491 self.assertLess(Lsup(Exz[0]-Ex0_ex_z), 1e-2 * Lsup(Ex0_ex_z))
492 self.assertLess(Lsup(Exz[1]-Ex1_ex_z), 1e-2 * Lsup(Ex1_ex_z))
493
494 argsr=model.getArguments(0.)
495 ref=model.getDefect(0., *argsr)
496
497 # this should be almost zero:
498 args=model.getArguments(SIGMA)
499 d=model.getDefect(SIGMA, *args)
500 self.assertTrue( d > 0.)
501 self.assertTrue( ref > 0.)
502 self.assertLess( d, 3e-3 * ref ) # d should be zero (some sort of)
503
504 z=ReducedFunction(domain).getX()[1]
505 Ex0_ex=cos(k.imag*(z-1))*exp(k.real*(z-1))
506 Ex1_ex=sin(k.imag*(z-1))*exp(k.real*(z-1))
507 Ex0_ex_z=-sin(k.imag*(z-1))*k.imag*exp(k.real*(z-1))+cos(k.imag*(z-1))*exp(k.real*(z-1))*k.real
508 Ex1_ex_z=cos(k.imag*(z-1))*k.imag*exp(k.real*(z-1))+sin(k.imag*(z-1))*exp(k.real*(z-1))*k.real
509 # and this should be zero
510 d0=model.getDefect(SIGMA, Ex0_ex*[1.,0]+ Ex1_ex*[0,1.], Ex0_ex_z*[1.,0]+ Ex1_ex_z*[0,1.])
511 self.assertTrue( d0 <= 1e-8 * ref ) # d should be zero (some sort of)
512
513 # and this too
514 dg=model.getGradient(SIGMA, Ex0_ex*[1.,0]+ Ex1_ex*[0,1.], Ex0_ex_z*[1.,0]+ Ex1_ex_z*[0,1.])
515 self.assertTrue(isinstance(dg, Data))
516 self.assertTrue(dg.getShape()==())
517 self.assertLess(Lsup(dg), 1e-10)
518
519 def test_Differential(self):
520 INC=0.001
521 omega=5.
522 mu0=0.123
523 SIGMA=15.
524 k=cmath.sqrt(1j*omega*mu0*SIGMA) # Ex=exp(k*z)
525
526 NE=101
527 domain=ripRectangle(max(NE,50*mpisize-1), max(NE,50*mpisize-1), d1=mpisize)
528
529 Z0=0.5
530 IMP=-(1j*omega*mu0)/k*(cmath.exp(k*Z0)-cmath.exp(-k*Z0))/(cmath.exp(k*Z0)+cmath.exp(-k*Z0))
531 Z_XY=[ IMP, IMP ]
532 H=1./NE
533 X1=(int(0.3/H)+0.5)*H
534 X2=(int(0.6/H)+0.5)*H
535 x=[ [X1,Z0], [X2,Z0] ]
536 eta=None
537
538 z=domain.getX()[1]
539 Ex0_ex=cos(k.imag*z)*(exp(k.real*z)-exp(-k.real*z))
540 Ex0_ex_z=-sin(k.imag*z)*k.imag*(exp(k.real*z)-exp(-k.real*z))+cos(k.imag*z)*(exp(k.real*z)+exp(-k.real*z))*k.real
541 Ex1_ex=sin(k.imag*z)*(exp(k.real*z)+exp(-k.real*z))
542 Ex1_ex_z=cos(k.imag*z)*k.imag*(exp(k.real*z)+exp(-k.real*z))+sin(k.imag*z)*(exp(k.real*z)-exp(-k.real*z))*k.real
543
544 model=MT2DModelTEMode(domain, omega, x, Z_XY, eta, mu=mu0, sigma0=SIGMA, tol=1e-9, directSolver=True, airLayerLevel=1.)
545
546 # this is the base line:
547 xx=domain.getX()[0]
548 SIGMA0=3.*(xx+0.3)
549 args0=model.getArguments(SIGMA0)
550 d0=model.getDefect(SIGMA0, *args0)
551 dg0=model.getGradient(SIGMA0, *args0)
552 self.assertTrue(isinstance(dg0, Data))
553 self.assertTrue(dg0.getShape()==())
554
555 X=Function(domain).getX()
556
557 # test 1
558 p=INC
559 SIGMA1=SIGMA0+p
560 args1=model.getArguments(SIGMA1)
561 d1=model.getDefect(SIGMA1, *args1)
562 self.assertLess( abs( d1-d0-integrate(dg0*p) ), 1e-2*abs(d1-d0) )
563
564 # test 2
565 p=exp(-length(X-(0.2,0.2))**2/10)*INC
566 SIGMA1=SIGMA0+p
567 args1=model.getArguments(SIGMA1)
568 d1=model.getDefect(SIGMA1, *args1)
569 self.assertLess( abs( d1-d0-integrate(dg0*p) ), 1e-2*abs(d1-d0) )
570
571 # test 3
572 p=sin(length(X)*3*3.14)*INC
573 SIGMA1=SIGMA0+p
574 args1=model.getArguments(SIGMA1)
575 d1=model.getDefect(SIGMA1, *args1)
576 self.assertLess( abs( d1-d0-integrate(dg0*p) ), 1e-2*abs(d1-d0) )
577
578
579 @unittest.skipUnless(HAVE_RIPLEY, "Ripley module not available")
580 @unittest.skipUnless(HAVE_DIRECT, "more than 1 MPI rank or missing direct solver")
581 class TestMT2DModelTMMode(unittest.TestCase):
582 def test_API(self):
583 domain=ripRectangle(25, 25, d0=mpisize)
584 omega=2.
585 SIGMA=15.
586 x=[ [0.2,0.5], [0.3,0.5] ]
587 Z_XY=[ complex(1.2,1.5), complex(1.3,2.5) ]
588 eta=1.
589 w0=1.
590 Hx0=1.
591 # now we do a real one
592 model=MT2DModelTMMode(domain, omega, x, Z_XY, eta, w0=w0, sigma0=SIGMA)
593 self.assertEqual(model.getDomain(), domain)
594 pde=model.setUpPDE()
595 self.assertIsInstance(pde, LinearPDE)
596 self.assertEqual(pde.getNumEquations(), 2)
597 self.assertEqual(pde.getNumSolutions(), 2)
598 self.assertEqual(pde.getDomain(), domain)
599
600 # other things that should work
601 model=MT2DModelTMMode(domain, omega, x, Z_XY, eta=None, w0=[2.,3.], sigma0=SIGMA, airLayerLevel=1.)
602
603 # these shouldn't work
604 self.assertRaises(ValueError, MT2DModelTMMode, domain, omega, x, [3.], eta=[1.,1.], w0=[2.,3.], sigma0=SIGMA)
605 self.assertRaises(ValueError, MT2DModelTMMode, domain, omega, x, Z_XY, eta=[1.], w0=[2.,3.], sigma0=SIGMA)
606 self.assertRaises(ValueError, MT2DModelTMMode, domain, omega, [(6.7,5)], Z_XY, eta=[1.,1.], w0=[2.,3.], sigma0=SIGMA)
607
608 def test_PDE(self):
609 omega=10.
610 mu0=0.123
611 RHO=0.15
612 SIGMA=1/RHO
613 k=cmath.sqrt(1j*omega*mu0/RHO) # Hx=exp(k*z)
614 NE=151
615 L=1
616 domain=ripRectangle(NE,NE, d0=mpisize)
617
618 Z0=0.5
619 H=1./NE
620 X1=(int(0.3/H)+0.5)*H
621 X2=(int(0.6/H)+0.5)*H
622
623 IMP=RHO*k*cmath.exp(k*(Z0-L))/cmath.exp(k*(Z0-L))
624 Z_XY=[ IMP, IMP ]
625
626 x=[ [X1,Z0], [X2,Z0] ]
627 eta=None
628
629 z=domain.getX()[1]
630 Hx0_ex=cos(k.imag*(z-L))*exp(k.real*(z-L))
631 Hx1_ex=sin(k.imag*(z-L))*exp(k.real*(z-L))
632 Hx0_ex_z=-sin(k.imag*(z-L))*k.imag*exp(k.real*(z-L))+cos(k.imag*(z-L))*exp(k.real*(z-L))*k.real
633 Hx1_ex_z=cos(k.imag*(z-L))*k.imag*exp(k.real*(z-L))+sin(k.imag*(z-L))*exp(k.real*(z-L))*k.real
634
635 model=MT2DModelTMMode(domain, omega, x, Z_XY, eta, mu=mu0, sigma0=SIGMA, tol=1e-9, directSolver=True, airLayerLevel=1.)
636
637 args=model.getArguments(RHO)
638 Hx=args[0]
639 g_Hx=args[1]
640 Hxz=g_Hx[:,1]
641 self.assertLess(Lsup(Hx[0]-Hx0_ex), 1e-4 * Lsup(Hx0_ex))
642 self.assertLess(Lsup(Hx[1]-Hx1_ex), 1e-4 * Lsup(Hx1_ex))
643 self.assertLess(Lsup(Hxz[0]-Hx0_ex_z), 1e-2 * Lsup(Hx0_ex_z))
644 self.assertLess(Lsup(Hxz[1]-Hx1_ex_z), 1e-2 * Lsup(Hx1_ex_z))
645
646 argsr=model.getArguments(1.)
647 ref=model.getDefect(1., *argsr)
648
649 # this should be almost zero:
650 args=model.getArguments(RHO)
651 d=model.getDefect(RHO, *args)
652 self.assertTrue( d > 0.)
653 self.assertTrue( ref > 0.)
654 self.assertTrue( d <= 3e-3 * ref ) # d should be zero (some sort of)
655
656 z=ReducedFunction(domain).getX()[1]
657 Hx0_ex=cos(k.imag*(z-L))*exp(k.real*(z-L))
658 Hx1_ex=sin(k.imag*(z-L))*exp(k.real*(z-L))
659 Hx0_ex_z=-sin(k.imag*(z-L))*k.imag*exp(k.real*(z-L))+cos(k.imag*(z-L))*exp(k.real*(z-L))*k.real
660 Hx1_ex_z=cos(k.imag*(z-L))*k.imag*exp(k.real*(z-L))+sin(k.imag*(z-L))*exp(k.real*(z-L))*k.real
661 g_Hx = Data(0, (2,2), Hx0_ex_z.getFunctionSpace())
662 g_Hx[0,1] = Hx0_ex_z
663 g_Hx[1,1] = Hx1_ex_z
664 # and this should be zero
665 d0=model.getDefect(RHO, Hx0_ex*[1.,0]+ Hx1_ex*[0,1.], g_Hx)
666 self.assertLess( d0, 1e-8 * ref ) # d should be zero (some sort of)
667
668 # and this too
669 dg=model.getGradient(RHO, Hx0_ex*[1.,0]+Hx1_ex*[0,1.], g_Hx)
670 self.assertTrue(isinstance(dg, Data))
671 self.assertTrue(dg.getShape()==())
672 self.assertLess(Lsup(dg), 1e-10)
673
674 def test_Differential(self):
675 INC=0.001
676 omega=5.
677 mu0=0.123
678 RHO=0.15
679 SIGMA=1/RHO
680 k=cmath.sqrt(1j*omega*mu0*1/RHO) # Hx=exp(k*z)
681
682 L=1
683 NE=101
684 domain=ripRectangle(max(NE,50*mpisize-1), max(NE,50*mpisize-1), d1=mpisize)
685
686 Z0=0.5
687 IMP=RHO*k*(cmath.exp(k*(Z0-L))-cmath.exp(-k*(Z0-L)))/(cmath.exp(k*(Z0-L))+cmath.exp(-k*(Z0-L)))
688 Z_XY=[ IMP, IMP ]
689 H=1./NE
690 X1=(int(0.3/H)+0.5)*H
691 X2=(int(0.6/H)+0.5)*H
692 x=[ [X1,Z0], [X2,Z0] ]
693 eta=None
694
695 model=MT2DModelTMMode(domain, omega, x, Z_XY, eta, mu=mu0, tol=1e-9, sigma0=SIGMA, directSolver=True)
696
697 # this is the base line:
698 xx=domain.getX()[0]
699 RHO0=3.*(xx+0.3)
700 args0=model.getArguments(RHO0)
701 d0=model.getDefect(RHO0, *args0)
702 dg0=model.getGradient(RHO0, *args0)
703 self.assertTrue(isinstance(dg0, Data))
704 self.assertTrue(dg0.getShape()==())
705
706 X=Function(domain).getX()
707
708 # test 1
709 p=INC
710 RHO1=RHO0+p
711 args1=model.getArguments(RHO1)
712 d1=model.getDefect(RHO1, *args1)
713 self.assertLess( abs( d1-d0-integrate(dg0*p) ), 1e-2*abs(d1-d0) )
714
715 # test 2
716 p=exp(-length(X-(0.2,0.2))**2/10)*INC
717 RHO1=RHO0+p
718 args1=model.getArguments(RHO1)
719 d1=model.getDefect(RHO1, *args1)
720 self.assertLess( abs( d1-d0-integrate(dg0*p) ), 1e-2*abs(d1-d0) )
721
722 # test 3
723 p=sin(length(X)*3*3.14)*INC
724 RHO1=RHO0+p
725 args1=model.getArguments(RHO1)
726 d1=model.getDefect(RHO1, *args1)
727 self.assertLess( abs( d1-d0-integrate(dg0*p) ), 1e-2*abs(d1-d0) )
728
729 if __name__ == '__main__':
730 run_tests(__name__, exit_on_failure=True)
731

  ViewVC Help
Powered by ViewVC 1.1.26