/[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 5559 - (show annotations)
Wed Mar 25 05:32:01 2015 UTC (3 years, 9 months ago) by caltinay
File MIME type: text/x-python
File size: 27086 byte(s)
changed axis of subdivision to avoid element stretching effects.

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.assertLess(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.assertLess(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.assertLess(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.assertLess( 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.assertLess( 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.assertLess(Lsup(d), 1e-10)
197 #self.assertTrue(d >= 0)
198
199 dg=acw.getGradient(sigma_comps, *args)
200
201 self.assertTrue(isinstance(dg, Data))
202 self.assertTrue(dg.getShape()==(2,))
203 self.assertTrue(dg.getFunctionSpace()==Solution(domain))
204 self.assertLess(Lsup(dg), 5e-10)
205 # this shouldn't be zero:
206 sigma0=Data([2*sigma.real, sigma.imag/2], Function(domain) )
207 args=acw.getArguments(sigma0)
208 d0=acw.getDefect(sigma0, *args)
209 self.assertTrue(isinstance(d0, float))
210 self.assertTrue(d0 >= 0)
211 self.assertTrue(d0 > 1e-10)
212
213 dg0=acw.getGradient(sigma0, *args)
214 self.assertTrue(isinstance(dg0, Data))
215 self.assertTrue(dg0.getShape()==(2,))
216 self.assertTrue(dg0.getFunctionSpace()==Solution(domain))
217 self.assertTrue(Lsup(dg0) > 1e-10)
218 # test the gradient numerically:
219 h=0.001
220 X=Function(domain).getX()
221 p=h*sin(length(X)*np.pi)*Lsup(length(sigma0))
222
223 sigma1=sigma0+p*[1,0]
224 args=acw.getArguments(sigma1)
225 d1=acw.getDefect(sigma1, *args)
226
227 self.assertLess( abs( d1-d0-integrate(dg0[0]*p) ), 1e-2*abs(d1-d0) )
228
229 sigma2=sigma0+p*[0,1]
230 args=acw.getArguments(sigma2)
231 d2=acw.getDefect(sigma2, *args)
232 self.assertLess( abs(d2-d0-integrate(dg0[1]*p)), 1e-2*abs(d2-d0) )
233
234
235 @unittest.skipIf(not HAVE_RIPLEY, "Ripley module not available")
236 class TestSubsidence(unittest.TestCase):
237 def test_PDE(self):
238 lam=2.
239 mu=1.
240
241 domain=ripBrick(20,20,19, d2=mpisize)
242
243 xb=FunctionOnBoundary(domain).getX()
244 m=whereZero(xb[2]-1)
245 w=m*[0,0,1]
246 d=m*2.5
247 acw=Subsidence(domain, w,d, lam, mu )
248
249 P0=10.
250 args0=acw.getArguments(P0)
251 u=args0[0]
252 self.assertLess(Lsup(u[0]), 1.e-8)
253 self.assertLess(Lsup(u[1]), 1.e-8)
254 self.assertLess(Lsup(u[2]-2.5*domain.getX()[2]), 1.e-8)
255
256 dd=acw.getDefect(P0, *args0)
257
258 self.assertTrue( dd >= 0.)
259 self.assertTrue( dd <= 1e-7 * 2.5 )
260
261 def test_Differential(self):
262 lam=2.
263 mu=1.
264 INC=0.01
265 domain=ripBrick(20,20,20*mpisize-1 , d2=mpisize)
266
267 xb=FunctionOnBoundary(domain).getX()
268 m=whereZero(xb[2]-1)
269 w=m*[0,0,1]
270 d=m*2.5
271 acw=Subsidence(domain, w,d, lam, mu )
272
273 x=Function(domain).getX()
274 P0=x[0]*x[1]
275 args0=acw.getArguments(P0)
276 d0=acw.getDefect(P0, *args0)
277 grad_d=acw.getGradient(P0, *args0)
278
279 dP=exp(-(length(x-[0.5,0.5,0.5])/0.06)**2)
280 P1=P0+INC*dP
281 args1=acw.getArguments(P1)
282 d1=acw.getDefect(P1, *args1)
283 ref=abs((d1-d0)/INC)
284 self.assertLess(abs((d1-d0)/INC-integrate(grad_d* dP)), ref * 1.e-5)
285
286 dP=exp(-(length(x-[0.3,0.3,0.5])/0.06)**2)
287 P2=P0-INC*dP
288 args2=acw.getArguments(P2)
289 d2=acw.getDefect(P2, *args2)
290 ref=abs((d2-d0)/INC)
291 self.assertLess(abs((d2-d0)/INC+integrate(grad_d* dP)), ref * 1.e-5)
292
293 @unittest.skipIf(not HAVE_FINLEY, "Finley module not available")
294 class TestDCResistivity(unittest.TestCase):
295
296 def test_PDE2D(self):
297 dx_tests=0.1
298 sigma0=1.
299 electrodes=[(0.5-2*dx_tests,1.), (0.5-dx_tests,1.), (0.5+dx_tests,1.), (0.5+2*dx_tests,1.)]
300 domain=finRectangle(20,20, d1=mpisize, diracPoints=electrodes, diracTags=["sl0", "sl1", "sr0", "sr1"] )
301 loc=Locator(domain,electrodes[2:])
302
303 # this creates some reference Data:
304 x=domain.getX()
305 q=whereZero(x[0]-inf(x[0]))+whereZero(x[0]-sup(x[0]))+whereZero(x[1]-inf(x[1]))
306 ppde=LinearPDE(domain, numEquations=1)
307 s=Scalar(0.,DiracDeltaFunctions(domain))
308 s.setTaggedValue("sl0" ,1.)
309 s.setTaggedValue("sl1",-1.)
310 ppde.setValue(A=kronecker(2)*sigma0, q=q, y_dirac=s)
311 pp=ppde.getSolution()
312 uu=loc(pp)
313
314 # arguments for DcRes
315 current = 10.
316 sourceInfo = [ "sl0", "sl1" ]
317 sampleTags = [ ("sr0", "sr1") ]
318
319 sigmaPrimary=7.
320 phiPrimary=pp*current*sigma0/sigmaPrimary
321
322 uuscale=1-current*sigma0/sigmaPrimary
323 delphi_in = [ (uu[1]-uu[0]) * uuscale]
324
325 acw=DcRes(domain, loc, delphi_in, sampleTags, phiPrimary, sigmaPrimary)
326
327 self.assertLess(Lsup(phiPrimary-acw.getPrimaryPotential()), 1.e-10 * Lsup(acw.getPrimaryPotential()))
328
329 SIGMA=10. # matches current
330 args0=acw.getArguments(SIGMA)
331 p=args0[0]
332 u=args0[1]
333
334 # true secondary potential
335 pps=pp-phiPrimary
336 self.assertLess(Lsup(p-pps), 1.e-6*Lsup(pps))
337
338 # test return values at electrodes:
339 self.assertLess(abs(u[0]-uu[0]*uuscale), 1.e-6 * abs(uu[0]*uuscale))
340 self.assertLess(abs(u[1]-uu[1]*uuscale), 1.e-6 * abs(uu[1]*uuscale))
341
342 # this sould be zero
343 dd=acw.getDefect(SIGMA, *args0)
344 self.assertTrue( dd >= 0.)
345 self.assertTrue( dd <= 1e-7 )
346
347 def test_Differential2D(self):
348 INC=0.001
349 sigma0=1.
350 dx_tests=0.1
351 electrodes=[(0.5-2*dx_tests,1.), (0.5-dx_tests,1.), (0.5+dx_tests,1.), (0.5+2*dx_tests,1.)]
352 domain=finRectangle(20,20, d1=mpisize, diracPoints=electrodes, diracTags=["sl0", "sl1", "sr0", "sr1"] )
353 loc=Locator(domain,electrodes[2:])
354
355 # arguments for DcRes
356 #current = 10.
357 sampleTags = [ ("sr0", "sr1") ]
358
359 delphi_in = [ 0.05 ]
360
361 sigmaPrimary=1
362 x=domain.getX()
363 phiPrimary=(x[0]-inf(x[0]))*(x[1]-inf(x[1]))*(x[0]-sup(x[0]))
364
365 acw=DcRes(domain, loc, delphi_in, sampleTags, phiPrimary, sigmaPrimary)
366
367 #=====================================================================
368 x=Function(domain).getX()
369 SIGMA0=x[0]*x[1]+1
370 args0=acw.getArguments(SIGMA0)
371 d0=acw.getDefect(SIGMA0, *args0)
372 grad_d=acw.getGradient(SIGMA0, *args0)
373
374 dS=exp(-(length(x-[0.5,0.5])/0.2)**2)
375 SIGMA1=SIGMA0+INC*dS
376 args1=acw.getArguments(SIGMA1)
377 d1=acw.getDefect(SIGMA1, *args1)
378 ref=abs((d1-d0)/INC)
379 self.assertLess(abs((d1-d0)/INC-integrate(grad_d* dS)), ref * 1.e-3)
380
381 dS=-exp(-(length(x-[0.5,0.5])/0.2)**2)
382 SIGMA2=SIGMA0+INC*dS
383 args2=acw.getArguments(SIGMA2)
384 d2=acw.getDefect(SIGMA2, *args2)
385 ref=abs((d2-d0)/INC)
386 self.assertLess(abs((d2-d0)/INC-integrate(grad_d* dS)), ref * 1.e-3)
387
388 dS=-1
389 SIGMA3=SIGMA0+INC*dS
390 args3=acw.getArguments(SIGMA3)
391 d3=acw.getDefect(SIGMA3, *args3)
392 ref=abs((d3-d0)/INC)
393 self.assertLess(abs((d3-d0)/INC-integrate(grad_d* dS)), ref * 1.e-3)
394
395 dS=1
396 SIGMA4=SIGMA0+INC*dS
397 args4=acw.getArguments(SIGMA4)
398 d4=acw.getDefect(SIGMA4, *args4)
399 ref=abs((d4-d0)/INC)
400 self.assertLess(abs((d4-d0)/INC-integrate(grad_d* dS)), ref * 1.e-3)
401
402 class TestIsostaticPressure(unittest.TestCase):
403 @unittest.skipIf(not HAVE_RIPLEY, "Ripley module not available")
404 def test_all(self):
405 domain=ripBrick(50,50,20*mpisize-1, d2=mpisize)
406
407 ps=IsostaticPressure(domain, level0=1., coordinates=None)
408 g=Vector(0., Function(domain))
409 rho=Scalar(100, Function(domain))
410 p0=ps.getPressure(g, rho)
411 p_ref=-(1.-domain.getX()[2])*981.
412 self.assertLess(Lsup(p0-p_ref), 1e-6 * Lsup(p_ref))
413
414 g=Vector([0,0,-10], Function(domain))
415 rho=Scalar(0, Function(domain))
416 p0=ps.getPressure(g, rho)
417 p_ref=-(1.-domain.getX()[2])*26700
418 self.assertLess(Lsup(p0-p_ref), 1e-6 * Lsup(p_ref))
419
420 g=Vector([0,0,-10], Function(domain))
421 rho=Scalar(100, Function(domain))
422 p0=ps.getPressure(g, rho)
423 p_ref=-(1.-domain.getX()[2])*(981.+26700+1000)
424 self.assertLess(Lsup(p0-p_ref), 1e-6 * Lsup(p_ref))
425
426 @unittest.skipIf(not HAVE_RIPLEY, "Ripley module not available")
427 class TestMT2DModelTEMode(unittest.TestCase):
428 def test_API(self):
429 domain=ripRectangle(25, 25, d1=mpisize)
430 omega=2.
431 x=[ [0.2,0.5], [0.3,0.5] ]
432 Z_XY=[ complex(1.2,1.5), complex(1.3,2.5) ]
433 eta=1.
434 w0=1.
435 Ex0=1.
436 # now we do a real one
437 model=MT2DModelTEMode(domain, omega, x, Z_XY, eta, w0=w0, Ex_top=Ex0)
438 self.assertEqual(model.getDomain(), domain)
439 pde=model.setUpPDE()
440 self.assertIsInstance(pde, LinearPDE)
441 self.assertEqual(pde.getNumEquations(), 2)
442 self.assertEqual(pde.getNumSolutions(), 2)
443 self.assertEqual(pde.getDomain(), domain)
444
445 # other things that should work
446 model=MT2DModelTEMode(domain, omega, x, Z_XY, eta=None, w0=[2.,3.], Ex_top=complex(4.5,6) )
447
448 # these shouldn't work
449 self.assertRaises(ValueError, MT2DModelTEMode, domain, omega, x, [3.], eta=[1.,1.], w0=[2.,3.], Ex_top=complex(4.5,6) )
450 self.assertRaises(ValueError, MT2DModelTEMode, domain, omega, x, Z_XY, eta=[1.], w0=[2.,3.], Ex_top=complex(4.5,6) )
451 self.assertRaises(ValueError, MT2DModelTEMode, domain, omega, [(6.7,5)], Z_XY, eta=[1.,1.], w0=[2.,3.], Ex_top=complex(4.5,6) )
452
453 def test_PDE(self):
454 omega=1.
455 mu0=0.123
456 SIGMA=15.
457 k=cmath.sqrt(1j*omega*mu0*SIGMA) # Ex=exp(k*z)
458 NE=101
459 domain=ripRectangle(NE,NE, d1=mpisize)
460
461 Z0=0.5
462 H=1./NE
463 X1=(int(0.3/H)+0.5)*H
464 X2=(int(0.6/H)+0.5)*H
465
466 IMP=-(1j*omega*mu0)/k*(cmath.exp(k*Z0)-cmath.exp(-k*Z0))/(cmath.exp(k*Z0)+cmath.exp(-k*Z0))
467 Z_XY=[ IMP, IMP ]
468
469 x=[ [X1,Z0], [X2,Z0] ]
470 eta=None
471
472 z=domain.getX()[1]
473 Ex0_ex=cos(k.imag*z)*(exp(k.real*z)-exp(-k.real*z))
474 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
475 Ex1_ex=sin(k.imag*z)*(exp(k.real*z)+exp(-k.real*z))
476 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
477
478 model=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)
479
480 args=model.getArguments(SIGMA)
481 Ex=args[0]
482 Exz=args[1]
483 self.assertTrue(Lsup(Ex[0]-Ex0_ex) <= 1e-4 * Lsup(Ex0_ex))
484 self.assertTrue(Lsup(Ex[1]-Ex1_ex) <= 1e-4 * Lsup(Ex1_ex))
485 self.assertTrue(Lsup(Exz[0]-Ex0_ex_z) <= 1e-2 * Lsup(Ex0_ex_z))
486 self.assertTrue(Lsup(Exz[1]-Ex1_ex_z) <= 1e-2 * Lsup(Ex1_ex_z))
487
488 argsr=model.getArguments(0.)
489 ref=model.getDefect(0., *argsr)
490
491 # this should be almost zero:
492 args=model.getArguments(SIGMA)
493 d=model.getDefect(SIGMA, *args)
494 self.assertTrue( d > 0.)
495 self.assertTrue( ref > 0.)
496 self.assertTrue( d <= 3e-3 * ref ) # d should be zero (some sort of)
497
498 z=ReducedFunction(domain).getX()[1]
499 Ex0_ex=cos(k.imag*z)*(exp(k.real*z)-exp(-k.real*z))
500 Ex1_ex=sin(k.imag*z)*(exp(k.real*z)+exp(-k.real*z))
501 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
502 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
503 # and this should be zero
504 d0=model.getDefect(SIGMA, Ex0_ex*[1.,0]+ Ex1_ex*[0,1.], Ex0_ex_z*[1.,0]+ Ex1_ex_z*[0,1.])
505 self.assertTrue( d0 <= 1e-8 * ref ) # d should be zero (some sort of)
506
507 # and this too
508 dg=model.getGradient(SIGMA, Ex0_ex*[1.,0]+ Ex1_ex*[0,1.], Ex0_ex_z*[1.,0]+ Ex1_ex_z*[0,1.])
509 self.assertTrue(isinstance(dg, Data))
510 self.assertTrue(dg.getShape()==())
511 self.assertLess(Lsup(dg), 1e-10)
512
513 def test_Differential(self):
514 INC=0.001
515 omega=5.
516 mu0=0.123
517 SIGMA=15.
518 k=cmath.sqrt(1j*omega*mu0*SIGMA) # Ex=exp(k*z)
519
520 NE=101
521 domain=ripRectangle(NE,NE, d1=mpisize)
522
523 Z0=0.5
524 IMP=-(1j*omega*mu0)/k*(cmath.exp(k*Z0)-cmath.exp(-k*Z0))/(cmath.exp(k*Z0)+cmath.exp(-k*Z0))
525 Z_XY=[ IMP, IMP ]
526 H=1./NE
527 X1=(int(0.3/H)+0.5)*H
528 X2=(int(0.6/H)+0.5)*H
529 x=[ [X1,Z0], [X2,Z0] ]
530 eta=None
531
532 z=domain.getX()[1]
533 Ex0_ex=cos(k.imag*z)*(exp(k.real*z)-exp(-k.real*z))
534 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
535 Ex1_ex=sin(k.imag*z)*(exp(k.real*z)+exp(-k.real*z))
536 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
537
538 model=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)
539
540 # this is the base line:
541 xx=domain.getX()[0]
542 SIGMA0=3.*(xx+0.3)
543 args0=model.getArguments(SIGMA0)
544 d0=model.getDefect(SIGMA0, *args0)
545 dg0=model.getGradient(SIGMA0, *args0)
546 self.assertTrue(isinstance(dg0, Data))
547 self.assertTrue(dg0.getShape()==())
548
549 X=Function(domain).getX()
550
551 # test 1
552 p=INC
553 SIGMA1=SIGMA0+p
554 args1=model.getArguments(SIGMA1)
555 d1=model.getDefect(SIGMA1, *args1)
556 self.assertLess( abs( d1-d0-integrate(dg0*p) ), 1e-2*abs(d1-d0) )
557
558 # test 2
559 p=exp(-length(X-(0.2,0.2))**2/10)*INC
560 SIGMA1=SIGMA0+p
561 args1=model.getArguments(SIGMA1)
562 d1=model.getDefect(SIGMA1, *args1)
563 self.assertLess( abs( d1-d0-integrate(dg0*p) ), 1e-2*abs(d1-d0) )
564
565 # test 3
566 p=sin(length(X)*3*3.14)*INC
567 SIGMA1=SIGMA0+p
568 args1=model.getArguments(SIGMA1)
569 d1=model.getDefect(SIGMA1, *args1)
570 self.assertLess( abs( d1-d0-integrate(dg0*p) ), 1e-2*abs(d1-d0) )
571
572
573 @unittest.skipIf(not HAVE_RIPLEY, "Ripley module not available")
574 class TestMT2DModelTMMode(unittest.TestCase):
575 def test_API(self):
576 domain=ripRectangle(25, 25, d0=mpisize)
577 omega=2.
578 x=[ [0.2,0.5], [0.3,0.5] ]
579 Z_XY=[ complex(1.2,1.5), complex(1.3,2.5) ]
580 eta=1.
581 w0=1.
582 Hx0=1.
583 # now we do a real one
584 model=MT2DModelTMMode(domain, omega, x, Z_XY, eta, w0=w0, Hx_bottom=Hx0)
585 self.assertEqual(model.getDomain(), domain)
586 pde=model.setUpPDE()
587 self.assertIsInstance(pde, LinearPDE)
588 self.assertEqual(pde.getNumEquations(), 2)
589 self.assertEqual(pde.getNumSolutions(), 2)
590 self.assertEqual(pde.getDomain(), domain)
591
592 # other things that should work
593 model=MT2DModelTMMode(domain, omega, x, Z_XY, eta=None, w0=[2.,3.], Hx_bottom=complex(4.5,6) )
594
595 # these shouldn't work
596 self.assertRaises(ValueError, MT2DModelTMMode, domain, omega, x, [3.], eta=[1.,1.], w0=[2.,3.], Hx_bottom=complex(4.5,6) )
597 self.assertRaises(ValueError, MT2DModelTMMode, domain, omega, x, Z_XY, eta=[1.], w0=[2.,3.], Hx_bottom=complex(4.5,6) )
598 self.assertRaises(ValueError, MT2DModelTMMode, domain, omega, [(6.7,5)], Z_XY, eta=[1.,1.], w0=[2.,3.], Hx_bottom=complex(4.5,6) )
599
600 def test_PDE(self):
601 omega=10.
602 mu0=0.123
603 RHO=0.15
604 k=cmath.sqrt(1j*omega*mu0/RHO) # Hx=exp(k*z)
605 NE=101
606 L=1
607 domain=ripRectangle(NE,NE, d0=mpisize)
608
609 Z0=0.5
610 H=1./NE
611 X1=(int(0.3/H)+0.5)*H
612 X2=(int(0.6/H)+0.5)*H
613
614 IMP=RHO*k*(cmath.exp(k*(Z0-L))-cmath.exp(-k*(Z0-L)))/(cmath.exp(k*(Z0-L))+cmath.exp(-k*(Z0-L)))
615 Z_XY=[ IMP, IMP ]
616
617 x=[ [X1,Z0], [X2,Z0] ]
618 eta=None
619
620 z=domain.getX()[1]
621 Hx0_ex=cos(k.imag*(z-L))*(exp(k.real*(z-L))+exp(-k.real*(z-L)))/2
622 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
623 Hx1_ex=sin(k.imag*(z-L))*(exp(k.real*(z-L))-exp(-k.real*(z-L)))/2
624 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
625
626 model=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)
627
628 args=model.getArguments(RHO)
629 Hx=args[0]
630 g_Hx=args[1]
631 Hxz=g_Hx[:,1]
632 self.assertLess(Lsup(Hx[0]-Hx0_ex), 1e-4 * Lsup(Hx0_ex))
633 self.assertLess(Lsup(Hx[1]-Hx1_ex), 1e-4 * Lsup(Hx1_ex))
634 self.assertLess(Lsup(Hxz[0]-Hx0_ex_z), 1e-2 * Lsup(Hx0_ex_z))
635 self.assertLess(Lsup(Hxz[1]-Hx1_ex_z), 1e-2 * Lsup(Hx1_ex_z))
636
637 argsr=model.getArguments(1.)
638 ref=model.getDefect(1., *argsr)
639
640 # this should be almost zero:
641 args=model.getArguments(RHO)
642 d=model.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 g_Hx = Data(0, (2,2), Hx0_ex_z.getFunctionSpace())
653 g_Hx[0,1] = Hx0_ex_z
654 g_Hx[1,1] = Hx1_ex_z
655 # and this should be zero
656 d0=model.getDefect(RHO, Hx0_ex*[1.,0]+ Hx1_ex*[0,1.], g_Hx)
657 self.assertLess( d0, 1e-8 * ref ) # d should be zero (some sort of)
658
659 # and this too
660 dg=model.getGradient(RHO, Hx0_ex*[1.,0]+Hx1_ex*[0,1.], g_Hx)
661 self.assertTrue(isinstance(dg, Data))
662 self.assertTrue(dg.getShape()==())
663 self.assertLess(Lsup(dg), 1e-10)
664
665 def test_Differential(self):
666 INC=0.001
667 omega=5.
668 mu0=0.123
669 RHO=0.15
670 k=cmath.sqrt(1j*omega*mu0*1/RHO) # Hx=exp(k*z)
671
672 L=1
673 NE=101
674 domain=ripRectangle(NE,NE, d0=mpisize)
675
676 Z0=0.5
677 IMP=RHO*k*(cmath.exp(k*(Z0-L))-cmath.exp(-k*(Z0-L)))/(cmath.exp(k*(Z0-L))+cmath.exp(-k*(Z0-L)))
678 Z_XY=[ IMP, IMP ]
679 H=1./NE
680 X1=(int(0.3/H)+0.5)*H
681 X2=(int(0.6/H)+0.5)*H
682 x=[ [X1,Z0], [X2,Z0] ]
683 eta=None
684
685 model=MT2DModelTMMode(domain, omega, x, Z_XY, eta, mu=mu0, tol=1e-9, directSolver=True)
686
687 # this is the base line:
688 xx=domain.getX()[0]
689 RHO0=3.*(xx+0.3)
690 args0=model.getArguments(RHO0)
691 d0=model.getDefect(RHO0, *args0)
692 dg0=model.getGradient(RHO0, *args0)
693 self.assertTrue(isinstance(dg0, Data))
694 self.assertTrue(dg0.getShape()==())
695
696 X=Function(domain).getX()
697
698 # test 1
699 p=INC
700 RHO1=RHO0+p
701 args1=model.getArguments(RHO1)
702 d1=model.getDefect(RHO1, *args1)
703 self.assertLess( abs( d1-d0-integrate(dg0*p) ), 1e-2*abs(d1-d0) )
704
705 # test 2
706 p=exp(-length(X-(0.2,0.2))**2/10)*INC
707 RHO1=RHO0+p
708 args1=model.getArguments(RHO1)
709 d1=model.getDefect(RHO1, *args1)
710 self.assertLess( abs( d1-d0-integrate(dg0*p) ), 1e-2*abs(d1-d0) )
711
712 # test 3
713 p=sin(length(X)*3*3.14)*INC
714 RHO1=RHO0+p
715 args1=model.getArguments(RHO1)
716 d1=model.getDefect(RHO1, *args1)
717 self.assertLess( abs( d1-d0-integrate(dg0*p) ), 1e-2*abs(d1-d0) )
718
719
720 if __name__ == '__main__':
721 run_tests(__name__, exit_on_failure=True)
722

  ViewVC Help
Powered by ViewVC 1.1.26