/[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 5290 - (hide annotations)
Wed Dec 3 01:15:33 2014 UTC (4 years, 4 months ago) by sshaw
File MIME type: text/x-python
File size: 20978 byte(s)
fixing some broken imports from earlier commit
1 sshaw 5288 from __future__ import print_function
2 gross 4688 ##############################################################################
3     #
4     # Copyright (c) 2003-2014 by University of Queensland
5     # 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     __copyright__="""Copyright (c) 2003-2014 by University of Queensland
18     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     from esys.ripley import Rectangle, Brick as ripRectangle, ripBrick
41     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 jfenwick 4713 have_direct=getEscriptParamInt("PASO_DIRECT")
67 gross 4688
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    
76    
77     data=Data([1,2], FunctionOnBoundary(domain))
78     F=Data([2,3], Function(domain))
79     w=1.
80     self.assertRaises(ValueError, AcousticWaveForm, domain, omega, w, data, 1.) # F is a scalar
81     self.assertRaises(ValueError, AcousticWaveForm, domain, omega, w, [1,2], F) # data is not Data
82     self.assertRaises(ValueError, AcousticWaveForm, domain, omega, w, Data([1,2], Function(domain)), F) # data is not on boundary
83     self.assertRaises(ValueError, AcousticWaveForm, domain, omega, w, Scalar(1, Function(domain)), F) # data is not of shape (2,)
84     self.assertRaises(ValueError, AcousticWaveForm, domain, omega, [1,2], data, F) # w is not a scalar
85     self.assertRaises(ValueError, AcousticWaveForm, domain, omega, Scalar(1, Function(domain)), data, F) # w is not a scalar
86    
87     # now we do a real one
88     acw=AcousticWaveForm(domain, omega, w, data, F)
89     self.assertEqual(acw.getDomain(), domain)
90     pde=acw.setUpPDE()
91     self.assertIsInstance(pde, LinearPDE)
92     self.assertEqual(pde.getNumEquations(), 2)
93     self.assertEqual(pde.getNumSolutions(), 2)
94     self.assertEqual(pde.getDomain(), domain)
95    
96    
97     def test_numeric2DscaleF(self):
98    
99 sshaw 5288 domain=ripRectangle(100,100, diracPoints=[(0.5,1.)], diracTags=['sss'])
100 gross 4688 omega=2.
101    
102     # test solution is u = a * z where a is complex
103     a=complex(3.45, 0.56)
104     sigma=complex(1e-3, 0.056)
105    
106    
107     data=Data([a.real, a.imag], FunctionOnBoundary(domain))
108     mydata=data.copy()
109    
110     z=FunctionOnBoundary(domain).getX()[1]
111     w=whereZero(z-1.)
112     # source:
113     F=Data( [1,0],Function(domain))
114     #
115     acw=AcousticWaveForm(domain, omega, w, data, F, coordinates=None, fixAtBottom=False, tol=1e-8, saveMemory=True, scaleF=True)
116     # check rescaled data
117     surv=acw.getSurvey()
118     self.assertAlmostEqual( integrate(length(surv[0])**2 * surv[1]), 1.)
119    
120     mydata_scale=sqrt( integrate(w*length(mydata)**2) )
121     self.assertAlmostEqual( acw.getSourceScaling(z*[1, 0.]) , a/mydata_scale )
122     self.assertAlmostEqual( acw.getSourceScaling(mydata) , 1./mydata_scale )
123    
124     # this should be zero:
125     sigma_comps=[sigma.real, sigma.imag]
126     args=acw.getArguments(sigma_comps)
127     d=acw.getDefect(sigma_comps, *args)
128     self.assertTrue(isinstance(d, float))
129 caltinay 5148 self.assertLess(abs(d), 1e-10)
130 gross 4688
131     dg=acw.getGradient(sigma_comps, *args)
132     self.assertTrue(isinstance(dg, Data))
133     self.assertTrue(dg.getShape()==(2,))
134     self.assertTrue(dg.getFunctionSpace()==Solution(domain))
135     self.assertTrue(Lsup(dg) < 1e-10)
136    
137     # this shuld be zero' too
138     sigma_comps=[2*sigma.real, sigma.imag/2.]
139     args=acw.getArguments(sigma_comps)
140     d=acw.getDefect(sigma_comps, *args)
141     self.assertTrue(isinstance(d, float))
142     self.assertTrue(abs(d)< 1e-10)
143    
144     dg=acw.getGradient(sigma_comps, *args)
145     self.assertTrue(isinstance(dg, Data))
146     self.assertTrue(dg.getShape()==(2,))
147     self.assertTrue(dg.getFunctionSpace()==Solution(domain))
148     self.assertTrue(Lsup(dg) < 1e-10)
149    
150     # this shouldn't be zero:
151     sigma0=[2*sigma.real, 10*a.imag]*(27*Function(domain).getX()[0]-Function(domain).getX()[1])
152     args=acw.getArguments(sigma0)
153     d0=acw.getDefect(sigma0, *args)
154     self.assertTrue(isinstance(d0, float))
155     self.assertTrue(d0 >= 0)
156     self.assertTrue(d0 > 1e-10)
157    
158     dg0=acw.getGradient(sigma0, *args)
159     self.assertTrue(isinstance(dg0, Data))
160     self.assertTrue(dg0.getShape()==(2,))
161     self.assertTrue(dg0.getFunctionSpace()==Solution(domain))
162     self.assertTrue(Lsup(dg0) > 1e-10)
163    
164     # test the gradient numerrically:
165     h=0.002
166     X=Function(domain).getX()
167     # .. increment:
168     p=h*exp(-(length(X-[0.6,0.6])/10)**2)*Lsup(length(sigma0))
169    
170    
171     sigma1=sigma0+p*[1,0]
172     args=acw.getArguments(sigma1)
173     d1=acw.getDefect(sigma1, *args)
174     self.assertTrue( abs( d1-d0-integrate(dg0[0]*p) ) < 1e-2 * abs(d1-d0) )
175    
176     sigma2=sigma0+p*[0,1]
177     args=acw.getArguments(sigma2)
178     d2=acw.getDefect(sigma2, *args)
179     self.assertTrue( abs(d2-d0-integrate(dg0[1]*p)) < 1e-2 * abs(d2-d0) )
180    
181     def test_numeric2DnoscaleF(self):
182    
183 sshaw 5288 domain=ripRectangle(10,20, diracPoints=[(0.5,1.)], diracTags=['sss'])
184 gross 4688 omega=1.5
185    
186     # test solution is u = a * z where a is complex
187     a=complex(3.45, 0.56)
188     sigma=complex(1e-3, 0.056)
189    
190    
191     data=Data([a.real, a.imag], FunctionOnBoundary(domain))
192     z=FunctionOnBoundary(domain).getX()[1]
193     w=whereZero(z-1.)
194     # F = - a*omega* sigma
195     F=Data( [-(a*omega**2*sigma).real, -(a*omega**2*sigma).imag ],Function(domain))
196    
197     acw=AcousticWaveForm(domain, omega, w, data, F, coordinates=None, fixAtBottom=False, tol=1e-8, saveMemory=True, scaleF=False)
198     # this should be zero:
199     sigma_comps=[sigma.real, sigma.imag]
200     args=acw.getArguments(sigma_comps)
201     d=acw.getDefect(sigma_comps, *args)
202     self.assertTrue(isinstance(d, float))
203 jfenwick 4719 self.assertTrue(Lsup(d) < 1e-10)
204     #self.assertTrue(d >= 0)
205     #self.assertTrue(d < 1e-10)
206 gross 4688
207     dg=acw.getGradient(sigma_comps, *args)
208    
209     self.assertTrue(isinstance(dg, Data))
210     self.assertTrue(dg.getShape()==(2,))
211     self.assertTrue(dg.getFunctionSpace()==Solution(domain))
212     self.assertTrue(Lsup(dg) < 5e-10)
213     # this shouldn't be zero:
214     sigma0=Data([2*sigma.real, sigma.imag/2], Function(domain) )
215     args=acw.getArguments(sigma0)
216     d0=acw.getDefect(sigma0, *args)
217     self.assertTrue(isinstance(d0, float))
218     self.assertTrue(d0 >= 0)
219     self.assertTrue(d0 > 1e-10)
220    
221     dg0=acw.getGradient(sigma0, *args)
222     self.assertTrue(isinstance(dg0, Data))
223     self.assertTrue(dg0.getShape()==(2,))
224     self.assertTrue(dg0.getFunctionSpace()==Solution(domain))
225     self.assertTrue(Lsup(dg0) > 1e-10)
226     # test the gradient numerrically:
227     h=0.001
228     X=Function(domain).getX()
229     p=h*sin(length(X)*np.pi)*Lsup(length(sigma0))
230    
231     sigma1=sigma0+p*[1,0]
232     args=acw.getArguments(sigma1)
233     d1=acw.getDefect(sigma1, *args)
234    
235     self.assertTrue( abs( d1-d0-integrate(dg0[0]*p) ) < 1e-2 * abs(d1-d0) )
236    
237     sigma2=sigma0+p*[0,1]
238     args=acw.getArguments(sigma2)
239     d2=acw.getDefect(sigma2, *args)
240     self.assertTrue( abs(d2-d0-integrate(dg0[1]*p)) < 1e-2 * abs(d2-d0) )
241    
242 sshaw 5288 @unittest.skipIf(not HAVE_RIPLEY, "Ripley module not available")
243 gross 4979 class TestMT2DModelTEMode(unittest.TestCase):
244     def test_API(self):
245 sshaw 5288 domain=ripRectangle(19, 19, d1=mpisize)
246 gross 4979 omega=2.
247     x=[ [0.2,0.5], [0.3,0.5] ]
248     Z_XY=[ complex(1.2,1.5), complex(1.3,2.5) ]
249     eta=1.
250     w0=1.
251     E_x0=1.
252     # now we do a real one
253     acw=MT2DModelTEMode(domain, omega, x, Z_XY, eta, w0=w0, E_x0=E_x0)
254     self.assertEqual(acw.getDomain(), domain)
255     pde=acw.setUpPDE()
256     self.assertIsInstance(pde, LinearPDE)
257     self.assertEqual(pde.getNumEquations(), 2)
258     self.assertEqual(pde.getNumSolutions(), 2)
259     self.assertEqual(pde.getDomain(), domain)
260    
261     # other things that should work
262     acw=MT2DModelTEMode(domain, omega, x, Z_XY, eta=[1.,1.], w0=[2.,3.], E_x0=complex(4.5,6) )
263    
264     # these shouldn't work
265     self.assertRaises(ValueError, MT2DModelTEMode, domain, omega, x, [3.], eta=[1.,1.], w0=[2.,3.], E_x0=complex(4.5,6) )
266     self.assertRaises(ValueError, MT2DModelTEMode, domain, omega, x, Z_XY, eta=[1.], w0=[2.,3.], E_x0=complex(4.5,6) )
267     self.assertRaises(ValueError, MT2DModelTEMode, domain, omega, [(6.7,5)], Z_XY, eta=[1.,1.], w0=[2.,3.], E_x0=complex(4.5,6) )
268    
269     def test_PDE(self):
270     omega=2.
271     mu0=0.123
272     SIGMA=15.
273     k=cmath.sqrt(1j*omega*mu0*SIGMA) # Ex=exp(k*z)
274    
275 sshaw 5288 domain=ripRectangle(199,199, d1=mpisize)
276 gross 4979
277    
278 gross 5057 IMP=cmath.sqrt(1j*omega*mu0/SIGMA)
279 gross 4980 Z_XY=[ IMP, IMP ]
280 gross 4979 x=[ [0.3,0.5], [0.6,0.5] ]
281 gross 4980 eta=0.005
282 gross 4979 z=domain.getX()[1]
283     Ex0_ex=exp(-k.real*z)*cos(-k.imag*z)
284     Ex0_ex_z=(-k.real*cos(-k.imag*z)+k.imag*sin(-k.imag*z)) * exp(-k.real*z)
285     Ex1_ex=exp(-k.real*z)*sin(-k.imag*z)
286     Ex1_ex_z=(-k.real*sin(-k.imag*z)-k.imag*cos(-k.imag*z)) * exp(-k.real*z)
287    
288 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)
289 gross 4979
290     args=acw.getArguments(SIGMA)
291     Ex=args[0]
292     Exz=args[1]
293     self.assertTrue(Lsup(Ex[0]-Ex0_ex) <= 1e-4 * Lsup(Ex0_ex))
294     self.assertTrue(Lsup(Ex[1]-Ex1_ex) <= 1e-4 * Lsup(Ex1_ex))
295     self.assertTrue(Lsup(Exz[0]-Ex0_ex_z) <= 1e-2 * Lsup(Ex0_ex_z))
296     self.assertTrue(Lsup(Exz[1]-Ex1_ex_z) <= 1e-2 * Lsup(Ex1_ex_z))
297    
298 gross 4980 argsr=acw.getArguments(0.)
299     ref=acw.getDefect(0., *argsr)
300    
301     # this should be almost zero:
302 gross 4979 args=acw.getArguments(SIGMA)
303     d=acw.getDefect(SIGMA, *args)
304    
305     self.assertTrue( d > 0.)
306     self.assertTrue( ref > 0.)
307     self.assertTrue( d <= 1e-4 * ref ) # d should be zero (some sort of)
308 gross 4980
309     # and this should be zero
310     d0=acw.getDefect(SIGMA, Ex0_ex*[1.,0]+ Ex1_ex*[0,1.], Ex0_ex_z*[1.,0]+ Ex1_ex_z*[0,1.])
311     self.assertTrue( d0 <= 1e-8 * ref ) # d should be zero (some sort of)
312 gross 4979
313 gross 4980
314     # and this too
315     dg=acw.getGradient(SIGMA, Ex0_ex*[1.,0]+ Ex1_ex*[0,1.], Ex0_ex_z*[1.,0]+ Ex1_ex_z*[0,1.])
316     self.assertTrue(isinstance(dg, Data))
317     self.assertTrue(dg.getShape()==())
318     self.assertTrue(Lsup(dg) < 1e-10)
319 gross 4979
320 gross 4981 def test_Differential(self):
321 gross 4980
322 gross 4981 INC=0.01
323 gross 4980
324     omega=2.
325     mu0=0.123
326     SIGMA=15.
327     k=cmath.sqrt(1j*omega*mu0*SIGMA) # Ex=exp(k*z)
328    
329 sshaw 5288 domain=ripRectangle(99,99, d1=mpisize)
330 gross 4980
331    
332     IMP=-cmath.sqrt(1j*omega*mu0/SIGMA)
333     Z_XY=[ IMP, IMP ]
334     x=[ [0.3,0.5], [0.6,0.5] ]
335     eta=0.005
336     z=domain.getX()[1]
337     Ex0_ex=exp(-k.real*z)*cos(-k.imag*z)
338     Ex0_ex_z=(-k.real*cos(-k.imag*z)+k.imag*sin(-k.imag*z)) * exp(-k.real*z)
339     Ex1_ex=exp(-k.real*z)*sin(-k.imag*z)
340     Ex1_ex_z=(-k.real*sin(-k.imag*z)-k.imag*cos(-k.imag*z)) * exp(-k.real*z)
341    
342 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 )
343 gross 4981
344     # this is the base line:
345     SIGMA0=10.
346 gross 4980 args0=acw.getArguments(SIGMA0)
347     d0=acw.getDefect(SIGMA0, *args0)
348    
349 gross 4981 dg0=acw.getGradient(SIGMA0, *args0)
350     self.assertTrue(isinstance(dg0, Data))
351     self.assertTrue(dg0.getShape()==())
352 gross 4980
353 gross 4981
354 gross 4980 X=Function(domain).getX()
355 gross 4981
356     # test 1
357     p=INC
358     SIGMA1=SIGMA0+p
359 gross 4980 args1=acw.getArguments(SIGMA1)
360     d1=acw.getDefect(SIGMA1, *args1)
361 gross 4981 self.assertTrue( abs( d1-d0-integrate(dg0*p) ) < 1e-2 * abs(d1-d0) )
362    
363     # test 2
364     p=exp(-length(X-(0.2,0.2))**2/10)*INC
365     SIGMA1=SIGMA0+p
366     args1=acw.getArguments(SIGMA1)
367     d1=acw.getDefect(SIGMA1, *args1)
368     self.assertTrue( abs( d1-d0-integrate(dg0*p) ) < 1e-2 * abs(d1-d0) )
369 gross 4980
370 gross 4981 # test 3
371     p=sin(length(X)*3*3.14)*INC
372     SIGMA1=SIGMA0+p
373     args1=acw.getArguments(SIGMA1)
374     d1=acw.getDefect(SIGMA1, *args1)
375     self.assertTrue( abs( d1-d0-integrate(dg0*p) ) < 1e-2 * abs(d1-d0) )
376 gross 5163
377 sshaw 5288 @unittest.skipIf(not HAVE_RIPLEY, "Ripley module not available")
378 gross 5163 class TestSubsidence(unittest.TestCase):
379     def test_PDE(self):
380    
381     lam=2.
382     mu=1.
383    
384 sshaw 5288 domain=ripBrick(20,20,19, d2=mpisize)
385 gross 4979
386 gross 5163 xb=FunctionOnBoundary(domain).getX()
387     m=whereZero(xb[2]-1)
388     w=m*[0,0,1]
389     d=m*2.5
390     acw=Subsidence(domain, w,d, lam, mu )
391    
392     P0=10.
393     args0=acw.getArguments(P0)
394     u=args0[0]
395     self.assertTrue(Lsup(u[0]) < 1.e-8)
396     self.assertTrue(Lsup(u[1]) < 1.e-8)
397     self.assertTrue(Lsup(u[2]-2.5*domain.getX()[2]) < 1.e-8)
398    
399     dd=acw.getDefect(P0, *args0)
400    
401     self.assertTrue( dd >= 0.)
402     self.assertTrue( dd <= 1e-7 * 2.5 )
403     def test_Differential(self):
404    
405     lam=2.
406     mu=1.
407    
408     INC=0.01
409 sshaw 5288 domain=ripBrick(20,20,20*mpisize-1 , d2=mpisize)
410 gross 5163
411     xb=FunctionOnBoundary(domain).getX()
412     m=whereZero(xb[2]-1)
413     w=m*[0,0,1]
414     d=m*2.5
415     acw=Subsidence(domain, w,d, lam, mu )
416    
417    
418     x=Function(domain).getX()
419     P0=x[0]*x[1]
420     args0=acw.getArguments(P0)
421     d0=acw.getDefect(P0, *args0)
422     grad_d=acw.getGradient(P0, *args0)
423    
424    
425     dP=exp(-(length(x-[0.5,0.5,0.5])/0.06)**2)
426     P1=P0+INC*dP
427     args1=acw.getArguments(P1)
428     d1=acw.getDefect(P1, *args1)
429     ref=abs((d1-d0)/INC)
430     self.assertTrue(abs((d1-d0)/INC-integrate(grad_d* dP)) < ref * 1.e-5)
431    
432     dP=exp(-(length(x-[0.3,0.3,0.5])/0.06)**2)
433     P2=P0-INC*dP
434     args2=acw.getArguments(P2)
435     d2=acw.getDefect(P2, *args2)
436     ref=abs((d2-d0)/INC)
437     self.assertTrue(abs((d2-d0)/INC+integrate(grad_d* dP)) < ref * 1.e-5)
438    
439 sshaw 5288 @unittest.skipIf(not HAVE_FINLEY, "Finley module not available")
440 gross 5262 class TestDCResistivity(unittest.TestCase):
441 gross 5163
442 gross 5262 def test_PDE2D(self):
443    
444     dx_tests=0.1
445 gross 5163
446 gross 5262 sigma0=1.
447     electrodes=[(0.5-2*dx_tests,1.), (0.5-dx_tests,1.), (0.5+dx_tests,1.), (0.5+2*dx_tests,1.)]
448 sshaw 5288 domain=finRectangle(20,20, d1=mpisize, diracPoints=electrodes, diracTags=["sl0", "sl1", "sr0", "sr1"] )
449 gross 5262 loc=Locator(domain,electrodes[2:])
450    
451     # this creates some reference Data:
452     x=domain.getX()
453     q=whereZero(x[0]-inf(x[0]))+whereZero(x[0]-sup(x[0]))+whereZero(x[1]-inf(x[1]))
454     ppde=LinearPDE(domain, numEquations=1)
455     s=Scalar(0.,DiracDeltaFunctions(domain))
456     s.setTaggedValue("sl0" ,1.)
457     s.setTaggedValue("sl1",-1.)
458     ppde.setValue(A=kronecker(2)*sigma0, q=q, y_dirac=s)
459     pp=ppde.getSolution()
460     uu=loc(pp)
461    
462     # arguments for DcRes
463     current = 10.
464     sourceInfo = [ "sl0", "sl1" ]
465     sampleTags = [ ("sr0", "sr1") ]
466    
467     sigmaPrimary=7.
468     phiPrimary=pp*current*sigma0/sigmaPrimary
469    
470     uuscale=1-current*sigma0/sigmaPrimary
471     delphi_in = [ (uu[1]-uu[0]) * uuscale]
472    
473     acw=DcRes(domain, loc, delphi_in, sampleTags, phiPrimary, sigmaPrimary)
474    
475     self.assertTrue(Lsup(phiPrimary-acw.getPrimaryPotential()) < 1.e-10 * Lsup(acw.getPrimaryPotential()))
476    
477     SIGMA=10. # matches current
478     args0=acw.getArguments(SIGMA)
479     p=args0[0]
480     u=args0[1]
481    
482     # true secondary potential
483     pps=pp-phiPrimary
484     self.assertTrue(Lsup(p-pps) < 1.e-6 * Lsup(pps))
485    
486    
487     # test return values at electrodes:
488     self.assertTrue(abs(u[0]-uu[0]*uuscale) < 1.e-6 * abs(uu[0]*uuscale))
489     self.assertTrue(abs(u[1]-uu[1]*uuscale) < 1.e-6 * abs(uu[1]*uuscale))
490    
491     # this sould be zero
492     dd=acw.getDefect(SIGMA, *args0)
493     self.assertTrue( dd >= 0.)
494     self.assertTrue( dd <= 1e-7 )
495    
496     def test_Differential2D(self):
497    
498     INC=0.001
499    
500     sigma0=1.
501     dx_tests=0.1
502     electrodes=[(0.5-2*dx_tests,1.), (0.5-dx_tests,1.), (0.5+dx_tests,1.), (0.5+2*dx_tests,1.)]
503 sshaw 5288 domain=finRectangle(20,20, d1=mpisize, diracPoints=electrodes, diracTags=["sl0", "sl1", "sr0", "sr1"] )
504 gross 5262 loc=Locator(domain,electrodes[2:])
505    
506     # arguments for DcRes
507     #current = 10.
508     sampleTags = [ ("sr0", "sr1") ]
509    
510     delphi_in = [ 0.05 ]
511    
512     sigmaPrimary=1
513     x=domain.getX()
514     phiPrimary=(x[0]-inf(x[0]))*(x[1]-inf(x[1]))*(x[0]-sup(x[0]))
515    
516     acw=DcRes(domain, loc, delphi_in, sampleTags, phiPrimary, sigmaPrimary)
517    
518     #===========================================================================
519     x=Function(domain).getX()
520     SIGMA0=x[0]*x[1]+1
521     args0=acw.getArguments(SIGMA0)
522     d0=acw.getDefect(SIGMA0, *args0)
523     grad_d=acw.getGradient(SIGMA0, *args0)
524    
525     dS=exp(-(length(x-[0.5,0.5])/0.2)**2)
526     SIGMA1=SIGMA0+INC*dS
527     args1=acw.getArguments(SIGMA1)
528     d1=acw.getDefect(SIGMA1, *args1)
529     ref=abs((d1-d0)/INC)
530     self.assertTrue(abs((d1-d0)/INC-integrate(grad_d* dS)) < ref * 1.e-3)
531    
532     dS=-exp(-(length(x-[0.5,0.5])/0.2)**2)
533     SIGMA2=SIGMA0+INC*dS
534     args2=acw.getArguments(SIGMA2)
535     d2=acw.getDefect(SIGMA2, *args2)
536     ref=abs((d2-d0)/INC)
537     self.assertTrue(abs((d2-d0)/INC-integrate(grad_d* dS)) < ref * 1.e-3)
538    
539     dS=-1
540     SIGMA3=SIGMA0+INC*dS
541     args3=acw.getArguments(SIGMA3)
542     d3=acw.getDefect(SIGMA3, *args3)
543     ref=abs((d3-d0)/INC)
544     self.assertTrue(abs((d3-d0)/INC-integrate(grad_d* dS)) < ref * 1.e-3)
545    
546     dS=1
547     SIGMA4=SIGMA0+INC*dS
548     args4=acw.getArguments(SIGMA4)
549     d4=acw.getDefect(SIGMA4, *args4)
550     ref=abs((d4-d0)/INC)
551     self.assertTrue(abs((d4-d0)/INC-integrate(grad_d* dS)) < ref * 1.e-3)
552    
553 gross 5163 class TestIsostaticPressure(unittest.TestCase):
554 sshaw 5288 @unittest.skipIf(not HAVE_RIPLEY, "Ripley module not available")
555 gross 5169 def test_all(self):
556 sshaw 5288 domain=ripBrick(50,50,20*mpisize-1, d2=mpisize)
557 gross 5169
558     ps=IsostaticPressure(domain, level0=1., coordinates=None)
559 gross 5163
560     g=Vector(0., Function(domain))
561 gross 5169 rho=Scalar(100, Function(domain))
562 gross 5163 p0=ps.getPressure(g, rho)
563 gross 5169 p_ref=-(1.-domain.getX()[2])*981.
564     self.assertTrue(Lsup(p0-p_ref) < 1e-6 * Lsup(p_ref))
565    
566     g=Vector([0,0,-10], Function(domain))
567     rho=Scalar(0, Function(domain))
568     p0=ps.getPressure(g, rho)
569     p_ref=-(1.-domain.getX()[2])*26700
570     self.assertTrue(Lsup(p0-p_ref) < 1e-6 * Lsup(p_ref))
571    
572     g=Vector([0,0,-10], Function(domain))
573     rho=Scalar(100, Function(domain))
574     p0=ps.getPressure(g, rho)
575     p_ref=-(1.-domain.getX()[2])*(981.+26700+1000)
576     self.assertTrue(Lsup(p0-p_ref) < 1e-6 * Lsup(p_ref))
577 gross 5163
578 sshaw 4984 if __name__ == '__main__':
579     run_tests(__name__, exit_on_failure=True)
580 gross 4688

  ViewVC Help
Powered by ViewVC 1.1.26