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

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 caltinay 5558 self.assertLess(Lsup(dg), 1e-10)
132 gross 4688
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 caltinay 5558 self.assertLess(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 caltinay 5558 self.assertLess(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 caltinay 5558 self.assertLess( abs( d1-d0-integrate(dg0[0]*p) ), 1e-2*abs(d1-d0) )
170 gross 4688
171     sigma2=sigma0+p*[0,1]
172     args=acw.getArguments(sigma2)
173     d2=acw.getDefect(sigma2, *args)
174 caltinay 5558 self.assertLess( 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 caltinay 5558 self.assertLess(Lsup(d), 1e-10)
197 jfenwick 4719 #self.assertTrue(d >= 0)
198 caltinay 5534
199 gross 4688 dg=acw.getGradient(sigma_comps, *args)
200    
201     self.assertTrue(isinstance(dg, Data))
202     self.assertTrue(dg.getShape()==(2,))
203 caltinay 5534 self.assertTrue(dg.getFunctionSpace()==Solution(domain))
204 caltinay 5558 self.assertLess(Lsup(dg), 5e-10)
205 gross 4688 # 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 caltinay 5534
213 gross 4688 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 caltinay 5534 # test the gradient numerically:
219 gross 4688 h=0.001
220     X=Function(domain).getX()
221     p=h*sin(length(X)*np.pi)*Lsup(length(sigma0))
222 caltinay 5534
223 gross 4688 sigma1=sigma0+p*[1,0]
224     args=acw.getArguments(sigma1)
225     d1=acw.getDefect(sigma1, *args)
226    
227 caltinay 5558 self.assertLess( abs( d1-d0-integrate(dg0[0]*p) ), 1e-2*abs(d1-d0) )
228 gross 4688
229     sigma2=sigma0+p*[0,1]
230     args=acw.getArguments(sigma2)
231     d2=acw.getDefect(sigma2, *args)
232 caltinay 5558 self.assertLess( abs(d2-d0-integrate(dg0[1]*p)), 1e-2*abs(d2-d0) )
233 gross 4688
234 gross 4979
235 sshaw 5288 @unittest.skipIf(not HAVE_RIPLEY, "Ripley module not available")
236 gross 5163 class TestSubsidence(unittest.TestCase):
237     def test_PDE(self):
238     lam=2.
239     mu=1.
240    
241 sshaw 5288 domain=ripBrick(20,20,19, d2=mpisize)
242 caltinay 5534
243 gross 5163 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 caltinay 5534
249     P0=10.
250 gross 5163 args0=acw.getArguments(P0)
251     u=args0[0]
252 caltinay 5558 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 caltinay 5534
256 gross 5163 dd=acw.getDefect(P0, *args0)
257 caltinay 5534
258 gross 5163 self.assertTrue( dd >= 0.)
259     self.assertTrue( dd <= 1e-7 * 2.5 )
260 caltinay 5534
261 gross 5163 def test_Differential(self):
262     lam=2.
263     mu=1.
264 caltinay 5534 INC=0.01
265 sshaw 5288 domain=ripBrick(20,20,20*mpisize-1 , d2=mpisize)
266 caltinay 5534
267 gross 5163 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 caltinay 5534
273 gross 5163 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 caltinay 5558 self.assertLess(abs((d1-d0)/INC-integrate(grad_d* dP)), ref * 1.e-5)
285 gross 5163
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 caltinay 5558 self.assertLess(abs((d2-d0)/INC+integrate(grad_d* dP)), ref * 1.e-5)
292 gross 5163
293 sshaw 5288 @unittest.skipIf(not HAVE_FINLEY, "Finley module not available")
294 gross 5262 class TestDCResistivity(unittest.TestCase):
295 gross 5163
296 gross 5262 def test_PDE2D(self):
297 caltinay 5534 dx_tests=0.1
298 gross 5262 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 sshaw 5288 domain=finRectangle(20,20, d1=mpisize, diracPoints=electrodes, diracTags=["sl0", "sl1", "sr0", "sr1"] )
301 gross 5262 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 caltinay 5534
314 gross 5262 # 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 caltinay 5534
325 gross 5262 acw=DcRes(domain, loc, delphi_in, sampleTags, phiPrimary, sigmaPrimary)
326    
327 caltinay 5558 self.assertLess(Lsup(phiPrimary-acw.getPrimaryPotential()), 1.e-10 * Lsup(acw.getPrimaryPotential()))
328 gross 5262
329 caltinay 5534 SIGMA=10. # matches current
330 gross 5262 args0=acw.getArguments(SIGMA)
331     p=args0[0]
332     u=args0[1]
333    
334     # true secondary potential
335     pps=pp-phiPrimary
336 caltinay 5558 self.assertLess(Lsup(p-pps), 1.e-6*Lsup(pps))
337 gross 5262
338     # test return values at electrodes:
339 caltinay 5558 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 gross 5262
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 caltinay 5534 dx_tests=0.1
351 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.)]
352 sshaw 5288 domain=finRectangle(20,20, d1=mpisize, diracPoints=electrodes, diracTags=["sl0", "sl1", "sr0", "sr1"] )
353 gross 5262 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 caltinay 5534 #=====================================================================
368 gross 5262 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 caltinay 5534
374 gross 5262 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 caltinay 5558 self.assertLess(abs((d1-d0)/INC-integrate(grad_d* dS)), ref * 1.e-3)
380 gross 5262
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 caltinay 5558 self.assertLess(abs((d2-d0)/INC-integrate(grad_d* dS)), ref * 1.e-3)
387 gross 5262
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 caltinay 5558 self.assertLess(abs((d3-d0)/INC-integrate(grad_d* dS)), ref * 1.e-3)
394 gross 5262
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 caltinay 5558 self.assertLess(abs((d4-d0)/INC-integrate(grad_d* dS)), ref * 1.e-3)
401 gross 5262
402 gross 5163 class TestIsostaticPressure(unittest.TestCase):
403 sshaw 5288 @unittest.skipIf(not HAVE_RIPLEY, "Ripley module not available")
404 gross 5169 def test_all(self):
405 sshaw 5288 domain=ripBrick(50,50,20*mpisize-1, d2=mpisize)
406 gross 5169
407     ps=IsostaticPressure(domain, level0=1., coordinates=None)
408 gross 5163 g=Vector(0., Function(domain))
409 gross 5169 rho=Scalar(100, Function(domain))
410 gross 5163 p0=ps.getPressure(g, rho)
411 gross 5169 p_ref=-(1.-domain.getX()[2])*981.
412 caltinay 5558 self.assertLess(Lsup(p0-p_ref), 1e-6 * Lsup(p_ref))
413 gross 5169
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 caltinay 5558 self.assertLess(Lsup(p0-p_ref), 1e-6 * Lsup(p_ref))
419 gross 5169
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 caltinay 5558 self.assertLess(Lsup(p0-p_ref), 1e-6 * Lsup(p_ref))
425 caltinay 5534
426 gross 5546 @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 caltinay 5558 model=MT2DModelTEMode(domain, omega, x, Z_XY, eta, w0=w0, Ex_top=Ex0)
438     self.assertEqual(model.getDomain(), domain)
439     pde=model.setUpPDE()
440 gross 5546 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 caltinay 5558 model=MT2DModelTEMode(domain, omega, x, Z_XY, eta=None, w0=[2.,3.], Ex_top=complex(4.5,6) )
447 gross 5546
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 caltinay 5558 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 gross 5546
480 caltinay 5558 args=model.getArguments(SIGMA)
481 gross 5546 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 caltinay 5558 argsr=model.getArguments(0.)
489     ref=model.getDefect(0., *argsr)
490 gross 5546
491     # this should be almost zero:
492 caltinay 5558 args=model.getArguments(SIGMA)
493     d=model.getDefect(SIGMA, *args)
494 gross 5546 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 caltinay 5558 d0=model.getDefect(SIGMA, Ex0_ex*[1.,0]+ Ex1_ex*[0,1.], Ex0_ex_z*[1.,0]+ Ex1_ex_z*[0,1.])
505 gross 5546 self.assertTrue( d0 <= 1e-8 * ref ) # d should be zero (some sort of)
506    
507     # and this too
508 caltinay 5558 dg=model.getGradient(SIGMA, Ex0_ex*[1.,0]+ Ex1_ex*[0,1.], Ex0_ex_z*[1.,0]+ Ex1_ex_z*[0,1.])
509 gross 5546 self.assertTrue(isinstance(dg, Data))
510     self.assertTrue(dg.getShape()==())
511 caltinay 5558 self.assertLess(Lsup(dg), 1e-10)
512 gross 5546
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 caltinay 5558 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 gross 5546
540     # this is the base line:
541 gross 5556 xx=domain.getX()[0]
542     SIGMA0=3.*(xx+0.3)
543 caltinay 5558 args0=model.getArguments(SIGMA0)
544     d0=model.getDefect(SIGMA0, *args0)
545     dg0=model.getGradient(SIGMA0, *args0)
546 gross 5546 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 caltinay 5558 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 gross 5546
558     # test 2
559     p=exp(-length(X-(0.2,0.2))**2/10)*INC
560     SIGMA1=SIGMA0+p
561 caltinay 5558 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 gross 5546
565     # test 3
566     p=sin(length(X)*3*3.14)*INC
567     SIGMA1=SIGMA0+p
568 caltinay 5558 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 gross 5546
572 gross 5554
573     @unittest.skipIf(not HAVE_RIPLEY, "Ripley module not available")
574     class TestMT2DModelTMMode(unittest.TestCase):
575     def test_API(self):
576 caltinay 5559 domain=ripRectangle(25, 25, d0=mpisize)
577 gross 5554 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 caltinay 5558 model=MT2DModelTMMode(domain, omega, x, Z_XY, eta, w0=w0, Hx_bottom=Hx0)
585     self.assertEqual(model.getDomain(), domain)
586     pde=model.setUpPDE()
587 gross 5554 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 caltinay 5558 model=MT2DModelTMMode(domain, omega, x, Z_XY, eta=None, w0=[2.,3.], Hx_bottom=complex(4.5,6) )
594 gross 5554
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 caltinay 5555 omega=10.
602 gross 5554 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 caltinay 5559 domain=ripRectangle(NE,NE, d0=mpisize)
608 gross 5554
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 caltinay 5558 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 gross 5554
628 caltinay 5558 args=model.getArguments(RHO)
629 gross 5554 Hx=args[0]
630 caltinay 5558 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 gross 5554
637 caltinay 5558 argsr=model.getArguments(1.)
638     ref=model.getDefect(1., *argsr)
639 gross 5554
640     # this should be almost zero:
641 caltinay 5558 args=model.getArguments(RHO)
642     d=model.getDefect(RHO, *args)
643 gross 5554 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 caltinay 5558 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 gross 5554 # and this should be zero
656 caltinay 5558 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 gross 5554
659     # and this too
660 caltinay 5558 dg=model.getGradient(RHO, Hx0_ex*[1.,0]+Hx1_ex*[0,1.], g_Hx)
661 gross 5554 self.assertTrue(isinstance(dg, Data))
662     self.assertTrue(dg.getShape()==())
663 caltinay 5558 self.assertLess(Lsup(dg), 1e-10)
664 gross 5554
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 caltinay 5559 domain=ripRectangle(NE,NE, d0=mpisize)
675 gross 5554
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 caltinay 5558 model=MT2DModelTMMode(domain, omega, x, Z_XY, eta, mu=mu0, tol=1e-9, directSolver=True)
686 gross 5554
687     # this is the base line:
688 gross 5556 xx=domain.getX()[0]
689     RHO0=3.*(xx+0.3)
690 caltinay 5558 args0=model.getArguments(RHO0)
691     d0=model.getDefect(RHO0, *args0)
692     dg0=model.getGradient(RHO0, *args0)
693 gross 5554 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 caltinay 5558 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 gross 5554
705     # test 2
706     p=exp(-length(X-(0.2,0.2))**2/10)*INC
707     RHO1=RHO0+p
708 caltinay 5558 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 gross 5554
712     # test 3
713     p=sin(length(X)*3*3.14)*INC
714     RHO1=RHO0+p
715 caltinay 5558 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 sshaw 4984 if __name__ == '__main__':
721     run_tests(__name__, exit_on_failure=True)
722 caltinay 5534

  ViewVC Help
Powered by ViewVC 1.1.26