/[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 5554 - (show annotations)
Mon Mar 23 06:54:53 2015 UTC (4 years, 1 month ago) by gross
File MIME type: text/x-python
File size: 26965 byte(s)
fixes and test for TM mode 
1 from __future__ import print_function
2 ##############################################################################
3 #
4 # Copyright (c) 2003-2015 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-2015 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 import logging
25 import esys.escriptcore.utestselect as unittest
26 from esys.escriptcore.testing import *
27 import numpy as np
28 import os
29 import sys
30 import cmath
31 from esys.downunder import *
32 from esys.escript import unitsSI as U
33 from esys.escript import *
34 from esys.weipa import saveSilo
35 from esys.escript.linearPDEs import LinearSinglePDE, LinearPDE
36 from esys.escript import getEscriptParamInt
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 try:
56 TEST_DATA_ROOT=os.environ['DOWNUNDER_TEST_DATA_ROOT']
57 except KeyError:
58 TEST_DATA_ROOT='ref_data'
59
60 try:
61 WORKDIR=os.environ['DOWNUNDER_WORKDIR']
62 except KeyError:
63 WORKDIR='.'
64
65
66 have_direct=getEscriptParamInt("PASO_DIRECT")
67
68 @unittest.skipIf(not HAVE_RIPLEY, "Ripley module not available")
69 @unittest.skipIf(mpisize>1 or have_direct!=1, "more than 1 MPI rank or missing direct solver")
70 class TestAcousticInversion(unittest.TestCase):
71 def test_API(self):
72
73 domain=ripRectangle(20,20, diracPoints=[(0.5,1.)], diracTags=['sss'])
74 omega=2.
75
76 data=Data([1,2], FunctionOnBoundary(domain))
77 F=Data([2,3], Function(domain))
78 w=1.
79 self.assertRaises(ValueError, AcousticWaveForm, domain, omega, w, data, 1.) # F is a scalar
80 self.assertRaises(ValueError, AcousticWaveForm, domain, omega, w, [1,2], F) # data is not Data
81 self.assertRaises(ValueError, AcousticWaveForm, domain, omega, w, Data([1,2], Function(domain)), F) # data is not on boundary
82 self.assertRaises(ValueError, AcousticWaveForm, domain, omega, w, Scalar(1, Function(domain)), F) # data is not of shape (2,)
83 self.assertRaises(ValueError, AcousticWaveForm, domain, omega, [1,2], data, F) # w is not a scalar
84 self.assertRaises(ValueError, AcousticWaveForm, domain, omega, Scalar(1, Function(domain)), data, F) # w is not a scalar
85
86 # now we do a real one
87 acw=AcousticWaveForm(domain, omega, w, data, F)
88 self.assertEqual(acw.getDomain(), domain)
89 pde=acw.setUpPDE()
90 self.assertIsInstance(pde, LinearPDE)
91 self.assertEqual(pde.getNumEquations(), 2)
92 self.assertEqual(pde.getNumSolutions(), 2)
93 self.assertEqual(pde.getDomain(), domain)
94
95 def test_numeric2DscaleF(self):
96 domain=ripRectangle(100,100, diracPoints=[(0.5,1.)], diracTags=['sss'])
97 omega=2.
98
99 # test solution is u = a * z where a is complex
100 a=complex(3.45, 0.56)
101 sigma=complex(1e-3, 0.056)
102
103 data=Data([a.real, a.imag], FunctionOnBoundary(domain))
104 mydata=data.copy()
105
106 z=FunctionOnBoundary(domain).getX()[1]
107 w=whereZero(z-1.)
108 # source:
109 F=Data( [1,0],Function(domain))
110
111 acw=AcousticWaveForm(domain, omega, w, data, F, coordinates=None, fixAtBottom=False, tol=1e-8, saveMemory=True, scaleF=True)
112 # check rescaled data
113 surv=acw.getSurvey()
114 self.assertAlmostEqual( integrate(length(surv[0])**2 * surv[1]), 1.)
115
116 mydata_scale=sqrt( integrate(w*length(mydata)**2) )
117 self.assertAlmostEqual( acw.getSourceScaling(z*[1, 0.]) , a/mydata_scale )
118 self.assertAlmostEqual( acw.getSourceScaling(mydata) , 1./mydata_scale )
119
120 # this should be zero:
121 sigma_comps=[sigma.real, sigma.imag]
122 args=acw.getArguments(sigma_comps)
123 d=acw.getDefect(sigma_comps, *args)
124 self.assertTrue(isinstance(d, float))
125 self.assertLess(abs(d), 1e-10)
126
127 dg=acw.getGradient(sigma_comps, *args)
128 self.assertTrue(isinstance(dg, Data))
129 self.assertTrue(dg.getShape()==(2,))
130 self.assertTrue(dg.getFunctionSpace()==Solution(domain))
131 self.assertTrue(Lsup(dg) < 1e-10)
132
133 # this shuld be zero' too
134 sigma_comps=[2*sigma.real, sigma.imag/2.]
135 args=acw.getArguments(sigma_comps)
136 d=acw.getDefect(sigma_comps, *args)
137 self.assertTrue(isinstance(d, float))
138 self.assertTrue(abs(d)< 1e-10)
139
140 dg=acw.getGradient(sigma_comps, *args)
141 self.assertTrue(isinstance(dg, Data))
142 self.assertTrue(dg.getShape()==(2,))
143 self.assertTrue(dg.getFunctionSpace()==Solution(domain))
144 self.assertTrue(Lsup(dg) < 1e-10)
145
146 # this shouldn't be zero:
147 sigma0=[2*sigma.real, 10*a.imag]*(27*Function(domain).getX()[0]-Function(domain).getX()[1])
148 args=acw.getArguments(sigma0)
149 d0=acw.getDefect(sigma0, *args)
150 self.assertTrue(isinstance(d0, float))
151 self.assertTrue(d0 >= 0)
152 self.assertTrue(d0 > 1e-10)
153
154 dg0=acw.getGradient(sigma0, *args)
155 self.assertTrue(isinstance(dg0, Data))
156 self.assertTrue(dg0.getShape()==(2,))
157 self.assertTrue(dg0.getFunctionSpace()==Solution(domain))
158 self.assertTrue(Lsup(dg0) > 1e-10)
159
160 # test the gradient numerrically:
161 h=0.002
162 X=Function(domain).getX()
163 # .. increment:
164 p=h*exp(-(length(X-[0.6,0.6])/10)**2)*Lsup(length(sigma0))
165
166 sigma1=sigma0+p*[1,0]
167 args=acw.getArguments(sigma1)
168 d1=acw.getDefect(sigma1, *args)
169 self.assertTrue( abs( d1-d0-integrate(dg0[0]*p) ) < 1e-2 * abs(d1-d0) )
170
171 sigma2=sigma0+p*[0,1]
172 args=acw.getArguments(sigma2)
173 d2=acw.getDefect(sigma2, *args)
174 self.assertTrue( abs(d2-d0-integrate(dg0[1]*p)) < 1e-2 * abs(d2-d0) )
175
176 def test_numeric2DnoscaleF(self):
177 domain=ripRectangle(10,20, diracPoints=[(0.5,1.)], diracTags=['sss'])
178 omega=1.5
179
180 # test solution is u = a * z where a is complex
181 a=complex(3.45, 0.56)
182 sigma=complex(1e-3, 0.056)
183
184 data=Data([a.real, a.imag], FunctionOnBoundary(domain))
185 z=FunctionOnBoundary(domain).getX()[1]
186 w=whereZero(z-1.)
187 # F = - a*omega* sigma
188 F=Data( [-(a*omega**2*sigma).real, -(a*omega**2*sigma).imag ],Function(domain))
189
190 acw=AcousticWaveForm(domain, omega, w, data, F, coordinates=None, fixAtBottom=False, tol=1e-8, saveMemory=True, scaleF=False)
191 # this should be zero:
192 sigma_comps=[sigma.real, sigma.imag]
193 args=acw.getArguments(sigma_comps)
194 d=acw.getDefect(sigma_comps, *args)
195 self.assertTrue(isinstance(d, float))
196 self.assertTrue(Lsup(d) < 1e-10)
197 #self.assertTrue(d >= 0)
198 #self.assertTrue(d < 1e-10)
199
200 dg=acw.getGradient(sigma_comps, *args)
201
202 self.assertTrue(isinstance(dg, Data))
203 self.assertTrue(dg.getShape()==(2,))
204 self.assertTrue(dg.getFunctionSpace()==Solution(domain))
205 self.assertTrue(Lsup(dg) < 5e-10)
206 # this shouldn't be zero:
207 sigma0=Data([2*sigma.real, sigma.imag/2], Function(domain) )
208 args=acw.getArguments(sigma0)
209 d0=acw.getDefect(sigma0, *args)
210 self.assertTrue(isinstance(d0, float))
211 self.assertTrue(d0 >= 0)
212 self.assertTrue(d0 > 1e-10)
213
214 dg0=acw.getGradient(sigma0, *args)
215 self.assertTrue(isinstance(dg0, Data))
216 self.assertTrue(dg0.getShape()==(2,))
217 self.assertTrue(dg0.getFunctionSpace()==Solution(domain))
218 self.assertTrue(Lsup(dg0) > 1e-10)
219 # test the gradient numerically:
220 h=0.001
221 X=Function(domain).getX()
222 p=h*sin(length(X)*np.pi)*Lsup(length(sigma0))
223
224 sigma1=sigma0+p*[1,0]
225 args=acw.getArguments(sigma1)
226 d1=acw.getDefect(sigma1, *args)
227
228 self.assertTrue( abs( d1-d0-integrate(dg0[0]*p) ) < 1e-2 * abs(d1-d0) )
229
230 sigma2=sigma0+p*[0,1]
231 args=acw.getArguments(sigma2)
232 d2=acw.getDefect(sigma2, *args)
233 self.assertTrue( abs(d2-d0-integrate(dg0[1]*p)) < 1e-2 * abs(d2-d0) )
234
235
236 @unittest.skipIf(not HAVE_RIPLEY, "Ripley module not available")
237 class TestSubsidence(unittest.TestCase):
238 def test_PDE(self):
239 lam=2.
240 mu=1.
241
242 domain=ripBrick(20,20,19, d2=mpisize)
243
244 xb=FunctionOnBoundary(domain).getX()
245 m=whereZero(xb[2]-1)
246 w=m*[0,0,1]
247 d=m*2.5
248 acw=Subsidence(domain, w,d, lam, mu )
249
250 P0=10.
251 args0=acw.getArguments(P0)
252 u=args0[0]
253 self.assertTrue(Lsup(u[0]) < 1.e-8)
254 self.assertTrue(Lsup(u[1]) < 1.e-8)
255 self.assertTrue(Lsup(u[2]-2.5*domain.getX()[2]) < 1.e-8)
256
257 dd=acw.getDefect(P0, *args0)
258
259 self.assertTrue( dd >= 0.)
260 self.assertTrue( dd <= 1e-7 * 2.5 )
261
262 def test_Differential(self):
263 lam=2.
264 mu=1.
265 INC=0.01
266 domain=ripBrick(20,20,20*mpisize-1 , d2=mpisize)
267
268 xb=FunctionOnBoundary(domain).getX()
269 m=whereZero(xb[2]-1)
270 w=m*[0,0,1]
271 d=m*2.5
272 acw=Subsidence(domain, w,d, lam, mu )
273
274 x=Function(domain).getX()
275 P0=x[0]*x[1]
276 args0=acw.getArguments(P0)
277 d0=acw.getDefect(P0, *args0)
278 grad_d=acw.getGradient(P0, *args0)
279
280 dP=exp(-(length(x-[0.5,0.5,0.5])/0.06)**2)
281 P1=P0+INC*dP
282 args1=acw.getArguments(P1)
283 d1=acw.getDefect(P1, *args1)
284 ref=abs((d1-d0)/INC)
285 self.assertTrue(abs((d1-d0)/INC-integrate(grad_d* dP)) < ref * 1.e-5)
286
287 dP=exp(-(length(x-[0.3,0.3,0.5])/0.06)**2)
288 P2=P0-INC*dP
289 args2=acw.getArguments(P2)
290 d2=acw.getDefect(P2, *args2)
291 ref=abs((d2-d0)/INC)
292 self.assertTrue(abs((d2-d0)/INC+integrate(grad_d* dP)) < ref * 1.e-5)
293
294 @unittest.skipIf(not HAVE_FINLEY, "Finley module not available")
295 class TestDCResistivity(unittest.TestCase):
296
297 def test_PDE2D(self):
298 dx_tests=0.1
299 sigma0=1.
300 electrodes=[(0.5-2*dx_tests,1.), (0.5-dx_tests,1.), (0.5+dx_tests,1.), (0.5+2*dx_tests,1.)]
301 domain=finRectangle(20,20, d1=mpisize, diracPoints=electrodes, diracTags=["sl0", "sl1", "sr0", "sr1"] )
302 loc=Locator(domain,electrodes[2:])
303
304 # this creates some reference Data:
305 x=domain.getX()
306 q=whereZero(x[0]-inf(x[0]))+whereZero(x[0]-sup(x[0]))+whereZero(x[1]-inf(x[1]))
307 ppde=LinearPDE(domain, numEquations=1)
308 s=Scalar(0.,DiracDeltaFunctions(domain))
309 s.setTaggedValue("sl0" ,1.)
310 s.setTaggedValue("sl1",-1.)
311 ppde.setValue(A=kronecker(2)*sigma0, q=q, y_dirac=s)
312 pp=ppde.getSolution()
313 uu=loc(pp)
314
315 # arguments for DcRes
316 current = 10.
317 sourceInfo = [ "sl0", "sl1" ]
318 sampleTags = [ ("sr0", "sr1") ]
319
320 sigmaPrimary=7.
321 phiPrimary=pp*current*sigma0/sigmaPrimary
322
323 uuscale=1-current*sigma0/sigmaPrimary
324 delphi_in = [ (uu[1]-uu[0]) * uuscale]
325
326 acw=DcRes(domain, loc, delphi_in, sampleTags, phiPrimary, sigmaPrimary)
327
328 self.assertTrue(Lsup(phiPrimary-acw.getPrimaryPotential()) < 1.e-10 * Lsup(acw.getPrimaryPotential()))
329
330 SIGMA=10. # matches current
331 args0=acw.getArguments(SIGMA)
332 p=args0[0]
333 u=args0[1]
334
335 # true secondary potential
336 pps=pp-phiPrimary
337 self.assertTrue(Lsup(p-pps) < 1.e-6 * Lsup(pps))
338
339
340 # test return values at electrodes:
341 self.assertTrue(abs(u[0]-uu[0]*uuscale) < 1.e-6 * abs(uu[0]*uuscale))
342 self.assertTrue(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.assertTrue(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.assertTrue(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.assertTrue(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.assertTrue(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.assertTrue(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.assertTrue(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.assertTrue(Lsup(p0-p_ref) < 1e-6 * Lsup(p_ref))
427
428 @unittest.skipIf(not HAVE_RIPLEY, "Ripley module not available")
429 class TestMT2DModelTEMode(unittest.TestCase):
430 def test_API(self):
431 domain=ripRectangle(25, 25, d1=mpisize)
432 omega=2.
433 x=[ [0.2,0.5], [0.3,0.5] ]
434 Z_XY=[ complex(1.2,1.5), complex(1.3,2.5) ]
435 eta=1.
436 w0=1.
437 Ex0=1.
438 # now we do a real one
439 acw=MT2DModelTEMode(domain, omega, x, Z_XY, eta, w0=w0, Ex_top=Ex0)
440 self.assertEqual(acw.getDomain(), domain)
441 pde=acw.setUpPDE()
442 self.assertIsInstance(pde, LinearPDE)
443 self.assertEqual(pde.getNumEquations(), 2)
444 self.assertEqual(pde.getNumSolutions(), 2)
445 self.assertEqual(pde.getDomain(), domain)
446
447 # other things that should work
448 acw=MT2DModelTEMode(domain, omega, x, Z_XY, eta=None, w0=[2.,3.], Ex_top=complex(4.5,6) )
449
450 # these shouldn't work
451 self.assertRaises(ValueError, MT2DModelTEMode, domain, omega, x, [3.], eta=[1.,1.], w0=[2.,3.], Ex_top=complex(4.5,6) )
452 self.assertRaises(ValueError, MT2DModelTEMode, domain, omega, x, Z_XY, eta=[1.], w0=[2.,3.], Ex_top=complex(4.5,6) )
453 self.assertRaises(ValueError, MT2DModelTEMode, domain, omega, [(6.7,5)], Z_XY, eta=[1.,1.], w0=[2.,3.], Ex_top=complex(4.5,6) )
454
455 def test_PDE(self):
456 omega=1.
457 mu0=0.123
458 SIGMA=15.
459 k=cmath.sqrt(1j*omega*mu0*SIGMA) # Ex=exp(k*z)
460 NE=101
461 domain=ripRectangle(NE,NE, d1=mpisize)
462
463 Z0=0.5
464 H=1./NE
465 X1=(int(0.3/H)+0.5)*H
466 X2=(int(0.6/H)+0.5)*H
467
468 IMP=-(1j*omega*mu0)/k*(cmath.exp(k*Z0)-cmath.exp(-k*Z0))/(cmath.exp(k*Z0)+cmath.exp(-k*Z0))
469 Z_XY=[ IMP, IMP ]
470
471 x=[ [X1,Z0], [X2,Z0] ]
472 eta=None
473
474 z=domain.getX()[1]
475 Ex0_ex=cos(k.imag*z)*(exp(k.real*z)-exp(-k.real*z))
476 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
477 Ex1_ex=sin(k.imag*z)*(exp(k.real*z)+exp(-k.real*z))
478 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
479
480 acw=MT2DModelTEMode(domain, omega, x, Z_XY, eta, mu=mu0, fixAtTop=True, Ex_top=Ex0_ex*[1.,0]+ Ex1_ex*[0,1.], tol=1e-9, directSolver=True)
481
482 args=acw.getArguments(SIGMA)
483 Ex=args[0]
484 Exz=args[1]
485 self.assertTrue(Lsup(Ex[0]-Ex0_ex) <= 1e-4 * Lsup(Ex0_ex))
486 self.assertTrue(Lsup(Ex[1]-Ex1_ex) <= 1e-4 * Lsup(Ex1_ex))
487 self.assertTrue(Lsup(Exz[0]-Ex0_ex_z) <= 1e-2 * Lsup(Ex0_ex_z))
488 self.assertTrue(Lsup(Exz[1]-Ex1_ex_z) <= 1e-2 * Lsup(Ex1_ex_z))
489
490 argsr=acw.getArguments(0.)
491 ref=acw.getDefect(0., *argsr)
492
493 # this should be almost zero:
494 args=acw.getArguments(SIGMA)
495 d=acw.getDefect(SIGMA, *args)
496 self.assertTrue( d > 0.)
497 self.assertTrue( ref > 0.)
498 self.assertTrue( d <= 3e-3 * ref ) # d should be zero (some sort of)
499
500 z=ReducedFunction(domain).getX()[1]
501 Ex0_ex=cos(k.imag*z)*(exp(k.real*z)-exp(-k.real*z))
502 Ex1_ex=sin(k.imag*z)*(exp(k.real*z)+exp(-k.real*z))
503 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
504 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
505 # and this should be zero
506 d0=acw.getDefect(SIGMA, Ex0_ex*[1.,0]+ Ex1_ex*[0,1.], Ex0_ex_z*[1.,0]+ Ex1_ex_z*[0,1.])
507 self.assertTrue( d0 <= 1e-8 * ref ) # d should be zero (some sort of)
508
509 # and this too
510 dg=acw.getGradient(SIGMA, Ex0_ex*[1.,0]+ Ex1_ex*[0,1.], Ex0_ex_z*[1.,0]+ Ex1_ex_z*[0,1.])
511 self.assertTrue(isinstance(dg, Data))
512 self.assertTrue(dg.getShape()==())
513 self.assertTrue(Lsup(dg) < 1e-10)
514
515 def test_Differential(self):
516 INC=0.001
517 omega=5.
518 mu0=0.123
519 SIGMA=15.
520 k=cmath.sqrt(1j*omega*mu0*SIGMA) # Ex=exp(k*z)
521
522 NE=101
523 domain=ripRectangle(NE,NE, d1=mpisize)
524
525 Z0=0.5
526 IMP=-(1j*omega*mu0)/k*(cmath.exp(k*Z0)-cmath.exp(-k*Z0))/(cmath.exp(k*Z0)+cmath.exp(-k*Z0))
527 Z_XY=[ IMP, IMP ]
528 H=1./NE
529 X1=(int(0.3/H)+0.5)*H
530 X2=(int(0.6/H)+0.5)*H
531 x=[ [X1,Z0], [X2,Z0] ]
532 eta=None
533
534 z=domain.getX()[1]
535 Ex0_ex=cos(k.imag*z)*(exp(k.real*z)-exp(-k.real*z))
536 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
537 Ex1_ex=sin(k.imag*z)*(exp(k.real*z)+exp(-k.real*z))
538 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
539
540 acw=MT2DModelTEMode(domain, omega, x, Z_XY, eta, mu=mu0, fixAtTop=True, Ex_top=Ex0_ex*[1.,0]+ Ex1_ex*[0,1.], tol=1e-9, directSolver=True)
541
542 # this is the base line:
543 SIGMA0=1.
544 args0=acw.getArguments(SIGMA0)
545 d0=acw.getDefect(SIGMA0, *args0)
546 dg0=acw.getGradient(SIGMA0, *args0)
547 self.assertTrue(isinstance(dg0, Data))
548 self.assertTrue(dg0.getShape()==())
549
550 X=Function(domain).getX()
551
552 # test 1
553 p=INC
554 SIGMA1=SIGMA0+p
555 args1=acw.getArguments(SIGMA1)
556 d1=acw.getDefect(SIGMA1, *args1)
557 self.assertTrue( abs( d1-d0-integrate(dg0*p) ) < 1e-2 * abs(d1-d0) )
558
559 # test 2
560 p=exp(-length(X-(0.2,0.2))**2/10)*INC
561 SIGMA1=SIGMA0+p
562 args1=acw.getArguments(SIGMA1)
563 d1=acw.getDefect(SIGMA1, *args1)
564 self.assertTrue( abs( d1-d0-integrate(dg0*p) ) < 1e-2 * abs(d1-d0) )
565
566 # test 3
567 p=sin(length(X)*3*3.14)*INC
568 SIGMA1=SIGMA0+p
569 args1=acw.getArguments(SIGMA1)
570 d1=acw.getDefect(SIGMA1, *args1)
571 self.assertTrue( abs( d1-d0-integrate(dg0*p) ) < 1e-2 * abs(d1-d0) )
572
573
574 @unittest.skipIf(not HAVE_RIPLEY, "Ripley module not available")
575 class TestMT2DModelTMMode(unittest.TestCase):
576 def test_API(self):
577 domain=ripRectangle(25, 25, d1=mpisize)
578 omega=2.
579 x=[ [0.2,0.5], [0.3,0.5] ]
580 Z_XY=[ complex(1.2,1.5), complex(1.3,2.5) ]
581 eta=1.
582 w0=1.
583 Hx0=1.
584 # now we do a real one
585 acw=MT2DModelTMMode(domain, omega, x, Z_XY, eta, w0=w0, Hx_bottom=Hx0)
586 self.assertEqual(acw.getDomain(), domain)
587 pde=acw.setUpPDE()
588 self.assertIsInstance(pde, LinearPDE)
589 self.assertEqual(pde.getNumEquations(), 2)
590 self.assertEqual(pde.getNumSolutions(), 2)
591 self.assertEqual(pde.getDomain(), domain)
592
593 # other things that should work
594 acw=MT2DModelTMMode(domain, omega, x, Z_XY, eta=None, w0=[2.,3.], Hx_bottom=complex(4.5,6) )
595
596 # these shouldn't work
597 self.assertRaises(ValueError, MT2DModelTMMode, domain, omega, x, [3.], eta=[1.,1.], w0=[2.,3.], Hx_bottom=complex(4.5,6) )
598 self.assertRaises(ValueError, MT2DModelTMMode, domain, omega, x, Z_XY, eta=[1.], w0=[2.,3.], Hx_bottom=complex(4.5,6) )
599 self.assertRaises(ValueError, MT2DModelTMMode, domain, omega, [(6.7,5)], Z_XY, eta=[1.,1.], w0=[2.,3.], Hx_bottom=complex(4.5,6) )
600
601 def test_PDE(self):
602 omega=1.
603 mu0=0.123
604 RHO=0.15
605 k=cmath.sqrt(1j*omega*mu0/RHO) # Hx=exp(k*z)
606 NE=101
607 L=1
608 domain=ripRectangle(NE,NE, d1=mpisize)
609
610 Z0=0.5
611 H=1./NE
612 X1=(int(0.3/H)+0.5)*H
613 X2=(int(0.6/H)+0.5)*H
614
615 IMP=RHO*k*(cmath.exp(k*(Z0-L))-cmath.exp(-k*(Z0-L)))/(cmath.exp(k*(Z0-L))+cmath.exp(-k*(Z0-L)))
616 Z_XY=[ IMP, IMP ]
617
618 x=[ [X1,Z0], [X2,Z0] ]
619 eta=None
620
621 z=domain.getX()[1]
622 Hx0_ex=cos(k.imag*(z-L))*(exp(k.real*(z-L))+exp(-k.real*(z-L)))/2
623 Hx0_ex_z=(-sin(k.imag*(z-L))*k.imag*(exp(k.real*(z-L))+exp(-k.real*(z-L)))+cos(k.imag*(z-L))*(exp(k.real*(z-L))-exp(-k.real*(z-L)))*k.real)/2
624 Hx1_ex=sin(k.imag*(z-L))*(exp(k.real*(z-L))-exp(-k.real*(z-L)))/2
625 Hx1_ex_z=(cos(k.imag*(z-L))*k.imag*(exp(k.real*(z-L))-exp(-k.real*(z-L)))+sin(k.imag*(z-L))*(exp(k.real*(z-L))+exp(-k.real*(z-L)))*k.real)/2
626
627 acw=MT2DModelTMMode(domain, omega, x, Z_XY, eta, mu=mu0, fixAtBottom=True, Hx_bottom=Hx0_ex*[1.,0]+ Hx1_ex*[0,1.], tol=1e-9, directSolver=True)
628
629 args=acw.getArguments(RHO)
630 Hx=args[0]
631 Hxz=args[1]
632 self.assertTrue(Lsup(Hx[0]-Hx0_ex) <= 1e-4 * Lsup(Hx0_ex))
633 self.assertTrue(Lsup(Hx[1]-Hx1_ex) <= 1e-4 * Lsup(Hx1_ex))
634 self.assertTrue(Lsup(Hxz[0]-Hx0_ex_z) <= 1e-2 * Lsup(Hx0_ex_z))
635 self.assertTrue(Lsup(Hxz[1]-Hx1_ex_z) <= 1e-2 * Lsup(Hx1_ex_z))
636
637 argsr=acw.getArguments(1.)
638 ref=acw.getDefect(1., *argsr)
639
640 # this should be almost zero:
641 args=acw.getArguments(RHO)
642 d=acw.getDefect(RHO, *args)
643 self.assertTrue( d > 0.)
644 self.assertTrue( ref > 0.)
645 self.assertTrue( d <= 3e-3 * ref ) # d should be zero (some sort of)
646
647 z=ReducedFunction(domain).getX()[1]
648 Hx0_ex=cos(k.imag*(z-L))*(exp(k.real*(z-L))+exp(-k.real*(z-L)))/2
649 Hx0_ex_z=(-sin(k.imag*(z-L))*k.imag*(exp(k.real*(z-L))+exp(-k.real*(z-L)))+cos(k.imag*(z-L))*(exp(k.real*(z-L))-exp(-k.real*(z-L)))*k.real)/2
650 Hx1_ex=sin(k.imag*(z-L))*(exp(k.real*(z-L))-exp(-k.real*(z-L)))/2
651 Hx1_ex_z=(cos(k.imag*(z-L))*k.imag*(exp(k.real*(z-L))-exp(-k.real*(z-L)))+sin(k.imag*(z-L))*(exp(k.real*(z-L))+exp(-k.real*(z-L)))*k.real)/2
652 # and this should be zero
653 d0=acw.getDefect(RHO, Hx0_ex*[1.,0]+ Hx1_ex*[0,1.], Hx0_ex_z*[1.,0]+ Hx1_ex_z*[0,1.])
654 self.assertTrue( d0 <= 1e-8 * ref ) # d should be zero (some sort of)
655
656 # and this too
657 dg=acw.getGradient(RHO, Hx0_ex*[1.,0]+ Hx1_ex*[0,1.], Hx0_ex_z*[1.,0]+ Hx1_ex_z*[0,1.])
658 self.assertTrue(isinstance(dg, Data))
659 self.assertTrue(dg.getShape()==())
660 self.assertTrue(Lsup(dg) < 1e-10)
661
662 def test_Differential(self):
663 INC=0.001
664 omega=5.
665 mu0=0.123
666 RHO=0.15
667 k=cmath.sqrt(1j*omega*mu0*1/RHO) # Hx=exp(k*z)
668
669 L=1
670 NE=101
671 domain=ripRectangle(NE,NE, d1=mpisize)
672
673 Z0=0.5
674 IMP=RHO*k*(cmath.exp(k*(Z0-L))-cmath.exp(-k*(Z0-L)))/(cmath.exp(k*(Z0-L))+cmath.exp(-k*(Z0-L)))
675 Z_XY=[ IMP, IMP ]
676 H=1./NE
677 X1=(int(0.3/H)+0.5)*H
678 X2=(int(0.6/H)+0.5)*H
679 x=[ [X1,Z0], [X2,Z0] ]
680 eta=None
681
682 acw=MT2DModelTMMode(domain, omega, x, Z_XY, eta, mu=mu0, tol=1e-9, directSolver=True)
683
684 # this is the base line:
685 RHO0=1. # was 100.15
686 args0=acw.getArguments(RHO0)
687 d0=acw.getDefect(RHO0, *args0)
688 dg0=acw.getGradient(RHO0, *args0)
689 self.assertTrue(isinstance(dg0, Data))
690 self.assertTrue(dg0.getShape()==())
691
692 X=Function(domain).getX()
693
694 # test 1
695 p=INC
696 RHO1=RHO0+p
697 args1=acw.getArguments(RHO1)
698 d1=acw.getDefect(RHO1, *args1)
699 self.assertTrue( abs( d1-d0-integrate(dg0*p) ) < 1e-2 * abs(d1-d0) )
700
701 # test 2
702 p=exp(-length(X-(0.2,0.2))**2/10)*INC
703 RHO1=RHO0+p
704 args1=acw.getArguments(RHO1)
705 d1=acw.getDefect(RHO1, *args1)
706 self.assertTrue( abs( d1-d0-integrate(dg0*p) ) < 1e-2 * abs(d1-d0) )
707
708 # test 3
709 p=sin(length(X)*3*3.14)*INC
710 RHO1=RHO0+p
711 args1=acw.getArguments(RHO1)
712 d1=acw.getDefect(RHO1, *args1)
713 self.assertTrue( abs( d1-d0-integrate(dg0*p) ) < 1e-2 * abs(d1-d0) )
714 if __name__ == '__main__':
715 run_tests(__name__, exit_on_failure=True)
716

  ViewVC Help
Powered by ViewVC 1.1.26