/[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 5555 - (hide annotations)
Tue Mar 24 06:31:43 2015 UTC (4 years, 1 month ago) by caltinay
File MIME type: text/x-python
File size: 26965 byte(s)
Correction to TM mode gradient

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     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 gross 5554
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 caltinay 5555 omega=10.
603 gross 5554 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 caltinay 5555 RHO0=2. # was 100.15
686 gross 5554 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 caltinay 5555 self.assertLess( abs( d1-d0-integrate(dg0*p) ), 1e-2 * abs(d1-d0) )
700 gross 5554
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 sshaw 4984 if __name__ == '__main__':
715     run_tests(__name__, exit_on_failure=True)
716 caltinay 5534

  ViewVC Help
Powered by ViewVC 1.1.26