/[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 6651 - (hide annotations)
Wed Feb 7 02:12:08 2018 UTC (14 months, 1 week ago) by jfenwick
File MIME type: text/x-python
File size: 27149 byte(s)
Make everyone sad by touching all the files

Copyright dates update

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

  ViewVC Help
Powered by ViewVC 1.1.26