/[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 5534 - (hide annotations)
Thu Mar 12 06:01:50 2015 UTC (4 years, 1 month ago) by caltinay
File MIME type: text/x-python
File size: 20369 byte(s)
A bit of no-op cleanup in forwardmodels before adding things.

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 sshaw 5288 @unittest.skipIf(not HAVE_RIPLEY, "Ripley module not available")
236 gross 4979 class TestMT2DModelTEMode(unittest.TestCase):
237     def test_API(self):
238 sshaw 5288 domain=ripRectangle(19, 19, d1=mpisize)
239 gross 4979 omega=2.
240     x=[ [0.2,0.5], [0.3,0.5] ]
241     Z_XY=[ complex(1.2,1.5), complex(1.3,2.5) ]
242     eta=1.
243     w0=1.
244     E_x0=1.
245     # now we do a real one
246     acw=MT2DModelTEMode(domain, omega, x, Z_XY, eta, w0=w0, E_x0=E_x0)
247     self.assertEqual(acw.getDomain(), domain)
248     pde=acw.setUpPDE()
249     self.assertIsInstance(pde, LinearPDE)
250     self.assertEqual(pde.getNumEquations(), 2)
251     self.assertEqual(pde.getNumSolutions(), 2)
252     self.assertEqual(pde.getDomain(), domain)
253    
254 caltinay 5534 # other things that should work
255 gross 4979 acw=MT2DModelTEMode(domain, omega, x, Z_XY, eta=[1.,1.], w0=[2.,3.], E_x0=complex(4.5,6) )
256    
257     # these shouldn't work
258     self.assertRaises(ValueError, MT2DModelTEMode, domain, omega, x, [3.], eta=[1.,1.], w0=[2.,3.], E_x0=complex(4.5,6) )
259     self.assertRaises(ValueError, MT2DModelTEMode, domain, omega, x, Z_XY, eta=[1.], w0=[2.,3.], E_x0=complex(4.5,6) )
260     self.assertRaises(ValueError, MT2DModelTEMode, domain, omega, [(6.7,5)], Z_XY, eta=[1.,1.], w0=[2.,3.], E_x0=complex(4.5,6) )
261    
262     def test_PDE(self):
263     omega=2.
264     mu0=0.123
265     SIGMA=15.
266     k=cmath.sqrt(1j*omega*mu0*SIGMA) # Ex=exp(k*z)
267    
268 sshaw 5288 domain=ripRectangle(199,199, d1=mpisize)
269 gross 4979
270 gross 5057 IMP=cmath.sqrt(1j*omega*mu0/SIGMA)
271 gross 4980 Z_XY=[ IMP, IMP ]
272 gross 4979 x=[ [0.3,0.5], [0.6,0.5] ]
273 gross 4980 eta=0.005
274 gross 4979 z=domain.getX()[1]
275     Ex0_ex=exp(-k.real*z)*cos(-k.imag*z)
276     Ex0_ex_z=(-k.real*cos(-k.imag*z)+k.imag*sin(-k.imag*z)) * exp(-k.real*z)
277     Ex1_ex=exp(-k.real*z)*sin(-k.imag*z)
278     Ex1_ex_z=(-k.real*sin(-k.imag*z)-k.imag*cos(-k.imag*z)) * exp(-k.real*z)
279    
280 caltinay 5235 acw=MT2DModelTEMode(domain, omega, x, Z_XY, eta, mu=mu0, fixAtBottom=True, E_x0=Ex0_ex*[1.,0]+ Ex1_ex*[0,1.], tol=1e-9)
281 gross 4979
282     args=acw.getArguments(SIGMA)
283     Ex=args[0]
284 caltinay 5534 Exz=args[1]
285 gross 4979 self.assertTrue(Lsup(Ex[0]-Ex0_ex) <= 1e-4 * Lsup(Ex0_ex))
286     self.assertTrue(Lsup(Ex[1]-Ex1_ex) <= 1e-4 * Lsup(Ex1_ex))
287     self.assertTrue(Lsup(Exz[0]-Ex0_ex_z) <= 1e-2 * Lsup(Ex0_ex_z))
288     self.assertTrue(Lsup(Exz[1]-Ex1_ex_z) <= 1e-2 * Lsup(Ex1_ex_z))
289    
290 gross 4980 argsr=acw.getArguments(0.)
291     ref=acw.getDefect(0., *argsr)
292 caltinay 5534
293 gross 4980 # this should be almost zero:
294 gross 4979 args=acw.getArguments(SIGMA)
295     d=acw.getDefect(SIGMA, *args)
296 caltinay 5534
297 gross 4979 self.assertTrue( d > 0.)
298     self.assertTrue( ref > 0.)
299 caltinay 5534 self.assertTrue( d <= 1e-4 * ref ) # d should be zero (some sort of)
300 gross 4980
301 caltinay 5534 # and this should be zero
302 gross 4980 d0=acw.getDefect(SIGMA, Ex0_ex*[1.,0]+ Ex1_ex*[0,1.], Ex0_ex_z*[1.,0]+ Ex1_ex_z*[0,1.])
303 caltinay 5534 self.assertTrue( d0 <= 1e-8 * ref ) # d should be zero (some sort of)
304 gross 4980
305     # and this too
306     dg=acw.getGradient(SIGMA, Ex0_ex*[1.,0]+ Ex1_ex*[0,1.], Ex0_ex_z*[1.,0]+ Ex1_ex_z*[0,1.])
307     self.assertTrue(isinstance(dg, Data))
308 caltinay 5534 self.assertTrue(dg.getShape()==())
309 gross 4980 self.assertTrue(Lsup(dg) < 1e-10)
310 caltinay 5534
311 gross 4981 def test_Differential(self):
312 caltinay 5534 INC=0.01
313 gross 4980 omega=2.
314     mu0=0.123
315     SIGMA=15.
316     k=cmath.sqrt(1j*omega*mu0*SIGMA) # Ex=exp(k*z)
317    
318 sshaw 5288 domain=ripRectangle(99,99, d1=mpisize)
319 gross 4980
320     IMP=-cmath.sqrt(1j*omega*mu0/SIGMA)
321     Z_XY=[ IMP, IMP ]
322     x=[ [0.3,0.5], [0.6,0.5] ]
323     eta=0.005
324     z=domain.getX()[1]
325     Ex0_ex=exp(-k.real*z)*cos(-k.imag*z)
326     Ex0_ex_z=(-k.real*cos(-k.imag*z)+k.imag*sin(-k.imag*z)) * exp(-k.real*z)
327     Ex1_ex=exp(-k.real*z)*sin(-k.imag*z)
328     Ex1_ex_z=(-k.real*sin(-k.imag*z)-k.imag*cos(-k.imag*z)) * exp(-k.real*z)
329    
330 caltinay 5230 acw=MT2DModelTEMode(domain, omega, x, Z_XY, eta, mu=mu0, fixAtBottom=True, E_x0=Ex0_ex*[1.,0]+ Ex1_ex*[0,1.], tol=1e-9 )
331 caltinay 5534
332 gross 4981 # this is the base line:
333 caltinay 5534 SIGMA0=10.
334 gross 4980 args0=acw.getArguments(SIGMA0)
335     d0=acw.getDefect(SIGMA0, *args0)
336 gross 4981 dg0=acw.getGradient(SIGMA0, *args0)
337     self.assertTrue(isinstance(dg0, Data))
338     self.assertTrue(dg0.getShape()==())
339 caltinay 5534
340 gross 4980 X=Function(domain).getX()
341 gross 4981
342     # test 1
343     p=INC
344     SIGMA1=SIGMA0+p
345 gross 4980 args1=acw.getArguments(SIGMA1)
346     d1=acw.getDefect(SIGMA1, *args1)
347 gross 4981 self.assertTrue( abs( d1-d0-integrate(dg0*p) ) < 1e-2 * abs(d1-d0) )
348 caltinay 5534
349 gross 4981 # test 2
350     p=exp(-length(X-(0.2,0.2))**2/10)*INC
351     SIGMA1=SIGMA0+p
352     args1=acw.getArguments(SIGMA1)
353     d1=acw.getDefect(SIGMA1, *args1)
354     self.assertTrue( abs( d1-d0-integrate(dg0*p) ) < 1e-2 * abs(d1-d0) )
355 gross 4980
356 gross 4981 # test 3
357     p=sin(length(X)*3*3.14)*INC
358     SIGMA1=SIGMA0+p
359     args1=acw.getArguments(SIGMA1)
360     d1=acw.getDefect(SIGMA1, *args1)
361     self.assertTrue( abs( d1-d0-integrate(dg0*p) ) < 1e-2 * abs(d1-d0) )
362 gross 5163
363 sshaw 5288 @unittest.skipIf(not HAVE_RIPLEY, "Ripley module not available")
364 gross 5163 class TestSubsidence(unittest.TestCase):
365     def test_PDE(self):
366     lam=2.
367     mu=1.
368    
369 sshaw 5288 domain=ripBrick(20,20,19, d2=mpisize)
370 caltinay 5534
371 gross 5163 xb=FunctionOnBoundary(domain).getX()
372     m=whereZero(xb[2]-1)
373     w=m*[0,0,1]
374     d=m*2.5
375     acw=Subsidence(domain, w,d, lam, mu )
376 caltinay 5534
377     P0=10.
378 gross 5163 args0=acw.getArguments(P0)
379     u=args0[0]
380     self.assertTrue(Lsup(u[0]) < 1.e-8)
381     self.assertTrue(Lsup(u[1]) < 1.e-8)
382     self.assertTrue(Lsup(u[2]-2.5*domain.getX()[2]) < 1.e-8)
383 caltinay 5534
384 gross 5163 dd=acw.getDefect(P0, *args0)
385 caltinay 5534
386 gross 5163 self.assertTrue( dd >= 0.)
387     self.assertTrue( dd <= 1e-7 * 2.5 )
388 caltinay 5534
389 gross 5163 def test_Differential(self):
390     lam=2.
391     mu=1.
392 caltinay 5534 INC=0.01
393 sshaw 5288 domain=ripBrick(20,20,20*mpisize-1 , d2=mpisize)
394 caltinay 5534
395 gross 5163 xb=FunctionOnBoundary(domain).getX()
396     m=whereZero(xb[2]-1)
397     w=m*[0,0,1]
398     d=m*2.5
399     acw=Subsidence(domain, w,d, lam, mu )
400 caltinay 5534
401 gross 5163 x=Function(domain).getX()
402     P0=x[0]*x[1]
403     args0=acw.getArguments(P0)
404     d0=acw.getDefect(P0, *args0)
405     grad_d=acw.getGradient(P0, *args0)
406    
407     dP=exp(-(length(x-[0.5,0.5,0.5])/0.06)**2)
408     P1=P0+INC*dP
409     args1=acw.getArguments(P1)
410     d1=acw.getDefect(P1, *args1)
411     ref=abs((d1-d0)/INC)
412 caltinay 5534 self.assertTrue(abs((d1-d0)/INC-integrate(grad_d* dP)) < ref * 1.e-5)
413 gross 5163
414     dP=exp(-(length(x-[0.3,0.3,0.5])/0.06)**2)
415     P2=P0-INC*dP
416     args2=acw.getArguments(P2)
417     d2=acw.getDefect(P2, *args2)
418     ref=abs((d2-d0)/INC)
419 caltinay 5534 self.assertTrue(abs((d2-d0)/INC+integrate(grad_d* dP)) < ref * 1.e-5)
420 gross 5163
421 sshaw 5288 @unittest.skipIf(not HAVE_FINLEY, "Finley module not available")
422 gross 5262 class TestDCResistivity(unittest.TestCase):
423 gross 5163
424 gross 5262 def test_PDE2D(self):
425 caltinay 5534 dx_tests=0.1
426 gross 5262 sigma0=1.
427     electrodes=[(0.5-2*dx_tests,1.), (0.5-dx_tests,1.), (0.5+dx_tests,1.), (0.5+2*dx_tests,1.)]
428 sshaw 5288 domain=finRectangle(20,20, d1=mpisize, diracPoints=electrodes, diracTags=["sl0", "sl1", "sr0", "sr1"] )
429 gross 5262 loc=Locator(domain,electrodes[2:])
430    
431     # this creates some reference Data:
432     x=domain.getX()
433     q=whereZero(x[0]-inf(x[0]))+whereZero(x[0]-sup(x[0]))+whereZero(x[1]-inf(x[1]))
434     ppde=LinearPDE(domain, numEquations=1)
435     s=Scalar(0.,DiracDeltaFunctions(domain))
436     s.setTaggedValue("sl0" ,1.)
437     s.setTaggedValue("sl1",-1.)
438     ppde.setValue(A=kronecker(2)*sigma0, q=q, y_dirac=s)
439     pp=ppde.getSolution()
440     uu=loc(pp)
441 caltinay 5534
442 gross 5262 # arguments for DcRes
443     current = 10.
444     sourceInfo = [ "sl0", "sl1" ]
445     sampleTags = [ ("sr0", "sr1") ]
446    
447     sigmaPrimary=7.
448     phiPrimary=pp*current*sigma0/sigmaPrimary
449    
450     uuscale=1-current*sigma0/sigmaPrimary
451     delphi_in = [ (uu[1]-uu[0]) * uuscale]
452 caltinay 5534
453 gross 5262 acw=DcRes(domain, loc, delphi_in, sampleTags, phiPrimary, sigmaPrimary)
454    
455     self.assertTrue(Lsup(phiPrimary-acw.getPrimaryPotential()) < 1.e-10 * Lsup(acw.getPrimaryPotential()))
456    
457 caltinay 5534 SIGMA=10. # matches current
458 gross 5262 args0=acw.getArguments(SIGMA)
459     p=args0[0]
460     u=args0[1]
461    
462     # true secondary potential
463     pps=pp-phiPrimary
464     self.assertTrue(Lsup(p-pps) < 1.e-6 * Lsup(pps))
465    
466    
467     # test return values at electrodes:
468     self.assertTrue(abs(u[0]-uu[0]*uuscale) < 1.e-6 * abs(uu[0]*uuscale))
469     self.assertTrue(abs(u[1]-uu[1]*uuscale) < 1.e-6 * abs(uu[1]*uuscale))
470    
471     # this sould be zero
472     dd=acw.getDefect(SIGMA, *args0)
473     self.assertTrue( dd >= 0.)
474     self.assertTrue( dd <= 1e-7 )
475    
476     def test_Differential2D(self):
477     INC=0.001
478     sigma0=1.
479 caltinay 5534 dx_tests=0.1
480 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.)]
481 sshaw 5288 domain=finRectangle(20,20, d1=mpisize, diracPoints=electrodes, diracTags=["sl0", "sl1", "sr0", "sr1"] )
482 gross 5262 loc=Locator(domain,electrodes[2:])
483    
484     # arguments for DcRes
485     #current = 10.
486     sampleTags = [ ("sr0", "sr1") ]
487    
488     delphi_in = [ 0.05 ]
489    
490     sigmaPrimary=1
491     x=domain.getX()
492     phiPrimary=(x[0]-inf(x[0]))*(x[1]-inf(x[1]))*(x[0]-sup(x[0]))
493    
494     acw=DcRes(domain, loc, delphi_in, sampleTags, phiPrimary, sigmaPrimary)
495    
496 caltinay 5534 #=====================================================================
497 gross 5262 x=Function(domain).getX()
498     SIGMA0=x[0]*x[1]+1
499     args0=acw.getArguments(SIGMA0)
500     d0=acw.getDefect(SIGMA0, *args0)
501     grad_d=acw.getGradient(SIGMA0, *args0)
502 caltinay 5534
503 gross 5262 dS=exp(-(length(x-[0.5,0.5])/0.2)**2)
504     SIGMA1=SIGMA0+INC*dS
505     args1=acw.getArguments(SIGMA1)
506     d1=acw.getDefect(SIGMA1, *args1)
507     ref=abs((d1-d0)/INC)
508 caltinay 5534 self.assertTrue(abs((d1-d0)/INC-integrate(grad_d* dS)) < ref * 1.e-3)
509 gross 5262
510     dS=-exp(-(length(x-[0.5,0.5])/0.2)**2)
511     SIGMA2=SIGMA0+INC*dS
512     args2=acw.getArguments(SIGMA2)
513     d2=acw.getDefect(SIGMA2, *args2)
514     ref=abs((d2-d0)/INC)
515 caltinay 5534 self.assertTrue(abs((d2-d0)/INC-integrate(grad_d* dS)) < ref * 1.e-3)
516 gross 5262
517     dS=-1
518     SIGMA3=SIGMA0+INC*dS
519     args3=acw.getArguments(SIGMA3)
520     d3=acw.getDefect(SIGMA3, *args3)
521     ref=abs((d3-d0)/INC)
522 caltinay 5534 self.assertTrue(abs((d3-d0)/INC-integrate(grad_d* dS)) < ref * 1.e-3)
523 gross 5262
524     dS=1
525     SIGMA4=SIGMA0+INC*dS
526     args4=acw.getArguments(SIGMA4)
527     d4=acw.getDefect(SIGMA4, *args4)
528     ref=abs((d4-d0)/INC)
529 caltinay 5534 self.assertTrue(abs((d4-d0)/INC-integrate(grad_d* dS)) < ref * 1.e-3)
530 gross 5262
531 gross 5163 class TestIsostaticPressure(unittest.TestCase):
532 sshaw 5288 @unittest.skipIf(not HAVE_RIPLEY, "Ripley module not available")
533 gross 5169 def test_all(self):
534 sshaw 5288 domain=ripBrick(50,50,20*mpisize-1, d2=mpisize)
535 gross 5169
536     ps=IsostaticPressure(domain, level0=1., coordinates=None)
537 gross 5163 g=Vector(0., Function(domain))
538 gross 5169 rho=Scalar(100, Function(domain))
539 gross 5163 p0=ps.getPressure(g, rho)
540 gross 5169 p_ref=-(1.-domain.getX()[2])*981.
541     self.assertTrue(Lsup(p0-p_ref) < 1e-6 * Lsup(p_ref))
542    
543     g=Vector([0,0,-10], Function(domain))
544     rho=Scalar(0, Function(domain))
545     p0=ps.getPressure(g, rho)
546     p_ref=-(1.-domain.getX()[2])*26700
547     self.assertTrue(Lsup(p0-p_ref) < 1e-6 * Lsup(p_ref))
548    
549     g=Vector([0,0,-10], Function(domain))
550     rho=Scalar(100, Function(domain))
551     p0=ps.getPressure(g, rho)
552     p_ref=-(1.-domain.getX()[2])*(981.+26700+1000)
553     self.assertTrue(Lsup(p0-p_ref) < 1e-6 * Lsup(p_ref))
554 caltinay 5534
555 sshaw 4984 if __name__ == '__main__':
556     run_tests(__name__, exit_on_failure=True)
557 caltinay 5534

  ViewVC Help
Powered by ViewVC 1.1.26