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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 5556 - (hide annotations)
Wed Mar 25 00:11:28 2015 UTC (4 years ago) by gross
File MIME type: text/x-python
File size: 27025 byte(s)
test added to fail TM forward
1 sshaw 5288 from __future__ import print_function
2 gross 4688 ##############################################################################
3     #
4 jfenwick 5448 # Copyright (c) 2003-2015 by University of Queensland
5 gross 4688 # 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 jfenwick 5448 __copyright__="""Copyright (c) 2003-2015 by University of Queensland
18 gross 4688 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 jfenwick 4938 import esys.escriptcore.utestselect as unittest
26 sshaw 4984 from esys.escriptcore.testing import *
27 gross 4688 import numpy as np
28     import os
29     import sys
30 gross 4979 import cmath
31 gross 4688 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 jfenwick 4713 from esys.escript import getEscriptParamInt
37 gross 5262 from esys.escript.pdetools import Locator
38 gross 4688
39 sshaw 5288 try:
40 sshaw 5308 from esys.ripley import Rectangle as ripRectangle, Brick as ripBrick
41 sshaw 5288 HAVE_RIPLEY = True
42     except ImportError as e:
43     HAVE_RIPLEY = False
44    
45     try:
46 sshaw 5290 from esys.finley import Rectangle as finRectangle
47 sshaw 5288 HAVE_FINLEY = True
48     except ImportError as e:
49     HAVE_FINLEY = False
50    
51 gross 4688 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 jfenwick 4713
65 gross 4688
66 caltinay 5534 have_direct=getEscriptParamInt("PASO_DIRECT")
67    
68 sshaw 5288 @unittest.skipIf(not HAVE_RIPLEY, "Ripley module not available")
69 jfenwick 4713 @unittest.skipIf(mpisize>1 or have_direct!=1, "more than 1 MPI rank or missing direct solver")
70 gross 4688 class TestAcousticInversion(unittest.TestCase):
71     def test_API(self):
72 sshaw 5288
73     domain=ripRectangle(20,20, diracPoints=[(0.5,1.)], diracTags=['sss'])
74 gross 4688 omega=2.
75 caltinay 5534
76 gross 4688 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 caltinay 5534
86 gross 4688 # 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 sshaw 5288 domain=ripRectangle(100,100, diracPoints=[(0.5,1.)], diracTags=['sss'])
97 gross 4688 omega=2.
98 caltinay 5534
99 gross 4688 # 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 caltinay 5534
103 gross 4688 data=Data([a.real, a.imag], FunctionOnBoundary(domain))
104     mydata=data.copy()
105 caltinay 5534
106 gross 4688 z=FunctionOnBoundary(domain).getX()[1]
107     w=whereZero(z-1.)
108     # source:
109     F=Data( [1,0],Function(domain))
110 caltinay 5534
111 gross 4688 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 caltinay 5534
120 gross 4688 # 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 caltinay 5148 self.assertLess(abs(d), 1e-10)
126 gross 4688
127     dg=acw.getGradient(sigma_comps, *args)
128     self.assertTrue(isinstance(dg, Data))
129     self.assertTrue(dg.getShape()==(2,))
130 caltinay 5534 self.assertTrue(dg.getFunctionSpace()==Solution(domain))
131 gross 4688 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 caltinay 5534
140 gross 4688 dg=acw.getGradient(sigma_comps, *args)
141     self.assertTrue(isinstance(dg, Data))
142     self.assertTrue(dg.getShape()==(2,))
143 caltinay 5534 self.assertTrue(dg.getFunctionSpace()==Solution(domain))
144 gross 4688 self.assertTrue(Lsup(dg) < 1e-10)
145 caltinay 5534
146 gross 4688 # 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 caltinay 5534
154 gross 4688 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 caltinay 5534
160 gross 4688 # 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 caltinay 5534
176 gross 4688 def test_numeric2DnoscaleF(self):
177 sshaw 5288 domain=ripRectangle(10,20, diracPoints=[(0.5,1.)], diracTags=['sss'])
178 gross 4688 omega=1.5
179 caltinay 5534
180 gross 4688 # 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 caltinay 5534
184 gross 4688 data=Data([a.real, a.imag], FunctionOnBoundary(domain))
185     z=FunctionOnBoundary(domain).getX()[1]
186     w=whereZero(z-1.)
187 caltinay 5534 # F = - a*omega* sigma
188 gross 4688 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 jfenwick 4719 self.assertTrue(Lsup(d) < 1e-10)
197     #self.assertTrue(d >= 0)
198     #self.assertTrue(d < 1e-10)
199 caltinay 5534
200 gross 4688 dg=acw.getGradient(sigma_comps, *args)
201    
202     self.assertTrue(isinstance(dg, Data))
203     self.assertTrue(dg.getShape()==(2,))
204 caltinay 5534 self.assertTrue(dg.getFunctionSpace()==Solution(domain))
205 gross 4688 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 caltinay 5534
214 gross 4688 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 caltinay 5534 # test the gradient numerically:
220 gross 4688 h=0.001
221     X=Function(domain).getX()
222     p=h*sin(length(X)*np.pi)*Lsup(length(sigma0))
223 caltinay 5534
224 gross 4688 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 gross 4979
236 sshaw 5288 @unittest.skipIf(not HAVE_RIPLEY, "Ripley module not available")
237 gross 5163 class TestSubsidence(unittest.TestCase):
238     def test_PDE(self):
239     lam=2.
240     mu=1.
241    
242 sshaw 5288 domain=ripBrick(20,20,19, d2=mpisize)
243 caltinay 5534
244 gross 5163 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 caltinay 5534
250     P0=10.
251 gross 5163 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 caltinay 5534
257 gross 5163 dd=acw.getDefect(P0, *args0)
258 caltinay 5534
259 gross 5163 self.assertTrue( dd >= 0.)
260     self.assertTrue( dd <= 1e-7 * 2.5 )
261 caltinay 5534
262 gross 5163 def test_Differential(self):
263     lam=2.
264     mu=1.
265 caltinay 5534 INC=0.01
266 sshaw 5288 domain=ripBrick(20,20,20*mpisize-1 , d2=mpisize)
267 caltinay 5534
268 gross 5163 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 caltinay 5534
274 gross 5163 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 caltinay 5534 self.assertTrue(abs((d1-d0)/INC-integrate(grad_d* dP)) < ref * 1.e-5)
286 gross 5163
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 caltinay 5534 self.assertTrue(abs((d2-d0)/INC+integrate(grad_d* dP)) < ref * 1.e-5)
293 gross 5163
294 sshaw 5288 @unittest.skipIf(not HAVE_FINLEY, "Finley module not available")
295 gross 5262 class TestDCResistivity(unittest.TestCase):
296 gross 5163
297 gross 5262 def test_PDE2D(self):
298 caltinay 5534 dx_tests=0.1
299 gross 5262 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 sshaw 5288 domain=finRectangle(20,20, d1=mpisize, diracPoints=electrodes, diracTags=["sl0", "sl1", "sr0", "sr1"] )
302 gross 5262 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 caltinay 5534
315 gross 5262 # 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 caltinay 5534
326 gross 5262 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 caltinay 5534 SIGMA=10. # matches current
331 gross 5262 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 caltinay 5534 dx_tests=0.1
353 gross 5262 electrodes=[(0.5-2*dx_tests,1.), (0.5-dx_tests,1.), (0.5+dx_tests,1.), (0.5+2*dx_tests,1.)]
354 sshaw 5288 domain=finRectangle(20,20, d1=mpisize, diracPoints=electrodes, diracTags=["sl0", "sl1", "sr0", "sr1"] )
355 gross 5262 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 caltinay 5534 #=====================================================================
370 gross 5262 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 caltinay 5534
376 gross 5262 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 caltinay 5534 self.assertTrue(abs((d1-d0)/INC-integrate(grad_d* dS)) < ref * 1.e-3)
382 gross 5262
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 caltinay 5534 self.assertTrue(abs((d2-d0)/INC-integrate(grad_d* dS)) < ref * 1.e-3)
389 gross 5262
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 caltinay 5534 self.assertTrue(abs((d3-d0)/INC-integrate(grad_d* dS)) < ref * 1.e-3)
396 gross 5262
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 caltinay 5534 self.assertTrue(abs((d4-d0)/INC-integrate(grad_d* dS)) < ref * 1.e-3)
403 gross 5262
404 gross 5163 class TestIsostaticPressure(unittest.TestCase):
405 sshaw 5288 @unittest.skipIf(not HAVE_RIPLEY, "Ripley module not available")
406 gross 5169 def test_all(self):
407 sshaw 5288 domain=ripBrick(50,50,20*mpisize-1, d2=mpisize)
408 gross 5169
409     ps=IsostaticPressure(domain, level0=1., coordinates=None)
410 gross 5163 g=Vector(0., Function(domain))
411 gross 5169 rho=Scalar(100, Function(domain))
412 gross 5163 p0=ps.getPressure(g, rho)
413 gross 5169 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 caltinay 5534
428 gross 5546 @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 gross 5556 xx=domain.getX()[0]
544     SIGMA0=3.*(xx+0.3)
545 gross 5546 args0=acw.getArguments(SIGMA0)
546     d0=acw.getDefect(SIGMA0, *args0)
547     dg0=acw.getGradient(SIGMA0, *args0)
548     self.assertTrue(isinstance(dg0, Data))
549     self.assertTrue(dg0.getShape()==())
550    
551     X=Function(domain).getX()
552    
553     # test 1
554     p=INC
555     SIGMA1=SIGMA0+p
556     args1=acw.getArguments(SIGMA1)
557     d1=acw.getDefect(SIGMA1, *args1)
558     self.assertTrue( abs( d1-d0-integrate(dg0*p) ) < 1e-2 * abs(d1-d0) )
559    
560     # test 2
561     p=exp(-length(X-(0.2,0.2))**2/10)*INC
562     SIGMA1=SIGMA0+p
563     args1=acw.getArguments(SIGMA1)
564     d1=acw.getDefect(SIGMA1, *args1)
565     self.assertTrue( abs( d1-d0-integrate(dg0*p) ) < 1e-2 * abs(d1-d0) )
566    
567     # test 3
568     p=sin(length(X)*3*3.14)*INC
569     SIGMA1=SIGMA0+p
570     args1=acw.getArguments(SIGMA1)
571     d1=acw.getDefect(SIGMA1, *args1)
572     self.assertTrue( abs( d1-d0-integrate(dg0*p) ) < 1e-2 * abs(d1-d0) )
573    
574 gross 5554
575     @unittest.skipIf(not HAVE_RIPLEY, "Ripley module not available")
576     class TestMT2DModelTMMode(unittest.TestCase):
577     def test_API(self):
578     domain=ripRectangle(25, 25, d1=mpisize)
579     omega=2.
580     x=[ [0.2,0.5], [0.3,0.5] ]
581     Z_XY=[ complex(1.2,1.5), complex(1.3,2.5) ]
582     eta=1.
583     w0=1.
584     Hx0=1.
585     # now we do a real one
586     acw=MT2DModelTMMode(domain, omega, x, Z_XY, eta, w0=w0, Hx_bottom=Hx0)
587     self.assertEqual(acw.getDomain(), domain)
588     pde=acw.setUpPDE()
589     self.assertIsInstance(pde, LinearPDE)
590     self.assertEqual(pde.getNumEquations(), 2)
591     self.assertEqual(pde.getNumSolutions(), 2)
592     self.assertEqual(pde.getDomain(), domain)
593    
594     # other things that should work
595     acw=MT2DModelTMMode(domain, omega, x, Z_XY, eta=None, w0=[2.,3.], Hx_bottom=complex(4.5,6) )
596    
597     # these shouldn't work
598     self.assertRaises(ValueError, MT2DModelTMMode, domain, omega, x, [3.], eta=[1.,1.], w0=[2.,3.], Hx_bottom=complex(4.5,6) )
599     self.assertRaises(ValueError, MT2DModelTMMode, domain, omega, x, Z_XY, eta=[1.], w0=[2.,3.], Hx_bottom=complex(4.5,6) )
600     self.assertRaises(ValueError, MT2DModelTMMode, domain, omega, [(6.7,5)], Z_XY, eta=[1.,1.], w0=[2.,3.], Hx_bottom=complex(4.5,6) )
601    
602     def test_PDE(self):
603 caltinay 5555 omega=10.
604 gross 5554 mu0=0.123
605     RHO=0.15
606     k=cmath.sqrt(1j*omega*mu0/RHO) # Hx=exp(k*z)
607     NE=101
608     L=1
609     domain=ripRectangle(NE,NE, d1=mpisize)
610    
611     Z0=0.5
612     H=1./NE
613     X1=(int(0.3/H)+0.5)*H
614     X2=(int(0.6/H)+0.5)*H
615    
616     IMP=RHO*k*(cmath.exp(k*(Z0-L))-cmath.exp(-k*(Z0-L)))/(cmath.exp(k*(Z0-L))+cmath.exp(-k*(Z0-L)))
617     Z_XY=[ IMP, IMP ]
618    
619     x=[ [X1,Z0], [X2,Z0] ]
620     eta=None
621    
622     z=domain.getX()[1]
623     Hx0_ex=cos(k.imag*(z-L))*(exp(k.real*(z-L))+exp(-k.real*(z-L)))/2
624     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
625     Hx1_ex=sin(k.imag*(z-L))*(exp(k.real*(z-L))-exp(-k.real*(z-L)))/2
626     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
627    
628     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)
629    
630     args=acw.getArguments(RHO)
631     Hx=args[0]
632     Hxz=args[1]
633     self.assertTrue(Lsup(Hx[0]-Hx0_ex) <= 1e-4 * Lsup(Hx0_ex))
634     self.assertTrue(Lsup(Hx[1]-Hx1_ex) <= 1e-4 * Lsup(Hx1_ex))
635     self.assertTrue(Lsup(Hxz[0]-Hx0_ex_z) <= 1e-2 * Lsup(Hx0_ex_z))
636     self.assertTrue(Lsup(Hxz[1]-Hx1_ex_z) <= 1e-2 * Lsup(Hx1_ex_z))
637    
638     argsr=acw.getArguments(1.)
639     ref=acw.getDefect(1., *argsr)
640    
641     # this should be almost zero:
642     args=acw.getArguments(RHO)
643     d=acw.getDefect(RHO, *args)
644     self.assertTrue( d > 0.)
645     self.assertTrue( ref > 0.)
646     self.assertTrue( d <= 3e-3 * ref ) # d should be zero (some sort of)
647    
648     z=ReducedFunction(domain).getX()[1]
649     Hx0_ex=cos(k.imag*(z-L))*(exp(k.real*(z-L))+exp(-k.real*(z-L)))/2
650     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
651     Hx1_ex=sin(k.imag*(z-L))*(exp(k.real*(z-L))-exp(-k.real*(z-L)))/2
652     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
653     # and this should be zero
654     d0=acw.getDefect(RHO, Hx0_ex*[1.,0]+ Hx1_ex*[0,1.], Hx0_ex_z*[1.,0]+ Hx1_ex_z*[0,1.])
655     self.assertTrue( d0 <= 1e-8 * ref ) # d should be zero (some sort of)
656    
657     # and this too
658     dg=acw.getGradient(RHO, Hx0_ex*[1.,0]+ Hx1_ex*[0,1.], Hx0_ex_z*[1.,0]+ Hx1_ex_z*[0,1.])
659     self.assertTrue(isinstance(dg, Data))
660     self.assertTrue(dg.getShape()==())
661     self.assertTrue(Lsup(dg) < 1e-10)
662    
663     def test_Differential(self):
664     INC=0.001
665     omega=5.
666     mu0=0.123
667     RHO=0.15
668     k=cmath.sqrt(1j*omega*mu0*1/RHO) # Hx=exp(k*z)
669    
670     L=1
671     NE=101
672     domain=ripRectangle(NE,NE, d1=mpisize)
673    
674     Z0=0.5
675     IMP=RHO*k*(cmath.exp(k*(Z0-L))-cmath.exp(-k*(Z0-L)))/(cmath.exp(k*(Z0-L))+cmath.exp(-k*(Z0-L)))
676     Z_XY=[ IMP, IMP ]
677     H=1./NE
678     X1=(int(0.3/H)+0.5)*H
679     X2=(int(0.6/H)+0.5)*H
680     x=[ [X1,Z0], [X2,Z0] ]
681     eta=None
682    
683     acw=MT2DModelTMMode(domain, omega, x, Z_XY, eta, mu=mu0, tol=1e-9, directSolver=True)
684    
685     # this is the base line:
686 gross 5556 xx=domain.getX()[0]
687     RHO0=3.*(xx+0.3)
688 gross 5554 args0=acw.getArguments(RHO0)
689     d0=acw.getDefect(RHO0, *args0)
690     dg0=acw.getGradient(RHO0, *args0)
691     self.assertTrue(isinstance(dg0, Data))
692     self.assertTrue(dg0.getShape()==())
693    
694     X=Function(domain).getX()
695    
696     # test 1
697     p=INC
698     RHO1=RHO0+p
699     args1=acw.getArguments(RHO1)
700     d1=acw.getDefect(RHO1, *args1)
701 caltinay 5555 self.assertLess( abs( d1-d0-integrate(dg0*p) ), 1e-2 * abs(d1-d0) )
702 gross 5554
703     # test 2
704     p=exp(-length(X-(0.2,0.2))**2/10)*INC
705     RHO1=RHO0+p
706     args1=acw.getArguments(RHO1)
707     d1=acw.getDefect(RHO1, *args1)
708     self.assertTrue( abs( d1-d0-integrate(dg0*p) ) < 1e-2 * abs(d1-d0) )
709    
710     # test 3
711     p=sin(length(X)*3*3.14)*INC
712     RHO1=RHO0+p
713     args1=acw.getArguments(RHO1)
714     d1=acw.getDefect(RHO1, *args1)
715     self.assertTrue( abs( d1-d0-integrate(dg0*p) ) < 1e-2 * abs(d1-d0) )
716 sshaw 4984 if __name__ == '__main__':
717     run_tests(__name__, exit_on_failure=True)
718 caltinay 5534

  ViewVC Help
Powered by ViewVC 1.1.26