/[escript]/trunk/downunder/test/python/run_forward.py
ViewVC logotype

Contents of /trunk/downunder/test/python/run_forward.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 5448 - (show annotations)
Fri Feb 6 05:31:37 2015 UTC (4 years, 2 months ago) by jfenwick
File MIME type: text/x-python
File size: 20980 byte(s)
Updating all the dates
1 from __future__ import print_function
2 ##############################################################################
3 #
4 # Copyright (c) 2003-2015 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-2015 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 import esys.escriptcore.utestselect as unittest
26 from esys.escriptcore.testing import *
27 import numpy as np
28 import os
29 import sys
30 import cmath
31 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 from esys.escript import getEscriptParamInt
37 from esys.escript.pdetools import Locator
38
39 try:
40 from esys.ripley import Rectangle as ripRectangle, Brick as ripBrick
41 HAVE_RIPLEY = True
42 except ImportError as e:
43 HAVE_RIPLEY = False
44
45 try:
46 from esys.finley import Rectangle as finRectangle
47 HAVE_FINLEY = True
48 except ImportError as e:
49 HAVE_FINLEY = False
50
51 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
65
66 have_direct=getEscriptParamInt("PASO_DIRECT")
67
68 @unittest.skipIf(not HAVE_RIPLEY, "Ripley module not available")
69 @unittest.skipIf(mpisize>1 or have_direct!=1, "more than 1 MPI rank or missing direct solver")
70 class TestAcousticInversion(unittest.TestCase):
71 def test_API(self):
72
73 domain=ripRectangle(20,20, diracPoints=[(0.5,1.)], diracTags=['sss'])
74 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 domain=ripRectangle(100,100, diracPoints=[(0.5,1.)], diracTags=['sss'])
100 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 self.assertLess(abs(d), 1e-10)
130
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 domain=ripRectangle(10,20, diracPoints=[(0.5,1.)], diracTags=['sss'])
184 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 self.assertTrue(Lsup(d) < 1e-10)
204 #self.assertTrue(d >= 0)
205 #self.assertTrue(d < 1e-10)
206
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 @unittest.skipIf(not HAVE_RIPLEY, "Ripley module not available")
243 class TestMT2DModelTEMode(unittest.TestCase):
244 def test_API(self):
245 domain=ripRectangle(19, 19, d1=mpisize)
246 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 domain=ripRectangle(199,199, d1=mpisize)
276
277
278 IMP=cmath.sqrt(1j*omega*mu0/SIGMA)
279 Z_XY=[ IMP, IMP ]
280 x=[ [0.3,0.5], [0.6,0.5] ]
281 eta=0.005
282 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 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
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 argsr=acw.getArguments(0.)
299 ref=acw.getDefect(0., *argsr)
300
301 # this should be almost zero:
302 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
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
313
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
320 def test_Differential(self):
321
322 INC=0.01
323
324 omega=2.
325 mu0=0.123
326 SIGMA=15.
327 k=cmath.sqrt(1j*omega*mu0*SIGMA) # Ex=exp(k*z)
328
329 domain=ripRectangle(99,99, d1=mpisize)
330
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 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
344 # this is the base line:
345 SIGMA0=10.
346 args0=acw.getArguments(SIGMA0)
347 d0=acw.getDefect(SIGMA0, *args0)
348
349 dg0=acw.getGradient(SIGMA0, *args0)
350 self.assertTrue(isinstance(dg0, Data))
351 self.assertTrue(dg0.getShape()==())
352
353
354 X=Function(domain).getX()
355
356 # test 1
357 p=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
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
370 # 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
377 @unittest.skipIf(not HAVE_RIPLEY, "Ripley module not available")
378 class TestSubsidence(unittest.TestCase):
379 def test_PDE(self):
380
381 lam=2.
382 mu=1.
383
384 domain=ripBrick(20,20,19, d2=mpisize)
385
386 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 domain=ripBrick(20,20,20*mpisize-1 , d2=mpisize)
410
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 @unittest.skipIf(not HAVE_FINLEY, "Finley module not available")
440 class TestDCResistivity(unittest.TestCase):
441
442 def test_PDE2D(self):
443
444 dx_tests=0.1
445
446 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 domain=finRectangle(20,20, d1=mpisize, diracPoints=electrodes, diracTags=["sl0", "sl1", "sr0", "sr1"] )
449 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 domain=finRectangle(20,20, d1=mpisize, diracPoints=electrodes, diracTags=["sl0", "sl1", "sr0", "sr1"] )
504 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 class TestIsostaticPressure(unittest.TestCase):
554 @unittest.skipIf(not HAVE_RIPLEY, "Ripley module not available")
555 def test_all(self):
556 domain=ripBrick(50,50,20*mpisize-1, d2=mpisize)
557
558 ps=IsostaticPressure(domain, level0=1., coordinates=None)
559
560 g=Vector(0., Function(domain))
561 rho=Scalar(100, Function(domain))
562 p0=ps.getPressure(g, rho)
563 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
578 if __name__ == '__main__':
579 run_tests(__name__, exit_on_failure=True)
580

  ViewVC Help
Powered by ViewVC 1.1.26