/[escript]/release/5.2/speckley/test/python/run_customAssemblersOnSpeckley.py
ViewVC logotype

Contents of /release/5.2/speckley/test/python/run_customAssemblersOnSpeckley.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6692 - (show annotations)
Mon Jun 25 02:31:06 2018 UTC (3 years, 2 months ago) by jfenwick
File MIME type: text/x-python
File size: 10448 byte(s)
Fix

1
2 ##############################################################################
3 #
4 # Copyright (c) 2003-2018 by The University of Queensland
5 # http://www.uq.edu.au
6 #
7 # Primary Business: Queensland, Australia
8 # Licensed under the Apache License, version 2.0
9 # http://www.apache.org/licenses/LICENSE-2.0
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 from __future__ import print_function, division
18
19 __copyright__="""Copyright (c) 2003-2018 by The University of Queensland
20 http://www.uq.edu.au
21 Primary Business: Queensland, Australia"""
22 __license__="""Licensed under the Apache License, version 2.0
23 http://www.apache.org/licenses/LICENSE-2.0"""
24 __url__="https://launchpad.net/escript-finley"
25
26 import esys.escriptcore.utestselect as unittest
27 from esys.escriptcore.testing import *
28 from esys.escript import *
29 from esys.speckley import Rectangle, Brick
30 from esys.escript.linearPDEs import LameEquation, LinearPDESystem, WavePDE
31 from esys.downunder import HTIWave, VTIWave, Ricker
32
33 EXPANDED, SCALAR, CONSTANT = range(3)
34
35 class SpeckleyWaveAssemblerTestBase(unittest.TestCase):
36 def generate_fast_HTI_PDE_solution(self, domain):
37 pde = WavePDE(domain, [("c11", self.c11),
38 ("c23", self.c23), ("c13", self.c13), ("c33", self.c33),
39 ("c44", self.c44), ("c66", self.c66)])
40 pde.getSolverOptions().setSolverMethod(SolverOptions.HRZ_LUMPING)
41 pde.setSymmetryOn()
42 dim = pde.getDim()
43 X = domain.getX()
44 y = Vector([2.,3.,4.][:dim], DiracDeltaFunctions(domain))
45 du = grad(X*X)
46 D = 2500.*kronecker(dim)
47
48 pde.setValue(D=D, y_dirac=y, du=du)
49 return pde.getSolution()
50
51 def generate_slow_HTI_PDE_solution(self, domain):
52 pde = LinearPDESystem(domain)
53 pde.getSolverOptions().setSolverMethod(SolverOptions.HRZ_LUMPING)
54 pde.setSymmetryOn()
55
56 dim = pde.getDim()
57 X = domain.getX()
58 y = Vector([2.,3.,4.][:dim], DiracDeltaFunctions(domain))
59 du = grad(X*X)
60 D = 2500.*kronecker(dim)
61 pde.setValue(X=pde.createCoefficient('X'))
62 sigma = pde.getCoefficient('X')
63 if dim == 3:
64 e11=du[0,0]
65 e22=du[1,1]
66 e33=du[2,2]
67
68 sigma[0,0]=self.c11*e11+self.c13*(e22+e33)
69 sigma[1,1]=self.c13*e11+self.c33*e22+self.c23*e33
70 sigma[2,2]=self.c13*e11+self.c23*e22+self.c33*e33
71
72 s=self.c44*(du[2,1]+du[1,2])
73 sigma[1,2]=s
74 sigma[2,1]=s
75
76 s=self.c66*(du[2,0]+du[0,2])
77 sigma[0,2]=s
78 sigma[2,0]=s
79
80 s=self.c66*(du[0,1]+du[1,0])
81 sigma[0,1]=s
82 sigma[1,0]=s
83
84 else:
85 e11=du[0,0]
86 e22=du[1,1]
87 sigma[0,0]=self.c11*e11+self.c13*e22
88 sigma[1,1]=self.c13*e11+self.c33*e22
89
90 s=self.c66*(du[1,0]+du[0,1])
91 sigma[0,1]=s
92 sigma[1,0]=s
93
94 pde.setValue(D=D, X=-sigma, y_dirac=y)
95 return pde.getSolution()
96
97 def generate_fast_VTI_PDE_solution(self, domain):
98 pde = WavePDE(domain, [("c11", self.c11),
99 ("c12", self.c12), ("c13", self.c13), ("c33", self.c33),
100 ("c44", self.c44), ("c66", self.c66)])
101 pde.getSolverOptions().setSolverMethod(SolverOptions.HRZ_LUMPING)
102 pde.setSymmetryOn()
103 dim = pde.getDim()
104
105 X = domain.getX()
106 y = Vector([2.,3.,4.][:dim], DiracDeltaFunctions(domain))
107 du = grad(X*X)
108 D = 2500.*kronecker(dim)
109
110 pde.setValue(D=D, y_dirac=y, du=du)
111 return pde.getSolution()
112
113 def generate_slow_VTI_PDE_solution(self, domain):
114 pde = LinearPDESystem(domain)
115 pde.getSolverOptions().setSolverMethod(SolverOptions.HRZ_LUMPING)
116 pde.setSymmetryOn()
117 dim = pde.getDim()
118
119 dim = pde.getDim()
120 X = domain.getX()
121 y = Vector([2.,3.,4.][:dim], DiracDeltaFunctions(domain))
122 du = grad(X*X)
123 D = 2500.*kronecker(dim)
124 pde.setValue(X=pde.createCoefficient('X'))
125 sigma = pde.getCoefficient('X')
126
127 if dim == 3:
128 e11=du[0,0]
129 e22=du[1,1]
130 e33=du[2,2]
131 sigma[0,0]=self.c11*e11+self.c12*e22+self.c13*e33
132 sigma[1,1]=self.c12*e11+self.c11*e22+self.c13*e33
133 sigma[2,2]=self.c13*(e11+e22)+self.c33*e33
134
135 s=self.c44*(du[2,1]+du[1,2])
136 sigma[1,2]=s
137 sigma[2,1]=s
138
139 s=self.c44*(du[2,0]+du[0,2])
140 sigma[0,2]=s
141 sigma[2,0]=s
142
143 s=self.c66*(du[0,1]+du[1,0])
144 sigma[0,1]=s
145 sigma[1,0]=s
146 else:
147 e11=du[0,0]
148 e22=du[1,1]
149 sigma[0,0]=self.c11*e11+self.c13*e22
150 sigma[1,1]=self.c13*e11+self.c33*e22
151 s=self.c44*(du[1,0]+du[0,1])
152 sigma[0,1]=s
153 sigma[1,0]=s
154
155 pde.setValue(D=D, X=-sigma, y_dirac=y)
156 return pde.getSolution()
157
158 def run_HTI_assembly(self, domain):
159 model = HTIWave(domain, self.V_p, self.V_s, self.wavelet, "source",
160 source_vector=[0,0,1], eps=0., gamma=0., delta=0.,
161 rho=2000., absorption_zone=None, lumping=True,
162 disable_fast_assemblers=True)
163 self.assertFalse(model.fastAssembler) #ensure the arg is obeyed
164
165 model = HTIWave(domain, self.V_p, self.V_s, self.wavelet, "source",
166 source_vector=[0,0,1], eps=0., gamma=0., delta=0.,
167 rho=2000., absorption_zone=None, lumping=True)
168 self.assertTrue(model.fastAssembler) #ensure fast is actually used
169
170 slow = self.generate_slow_HTI_PDE_solution(domain)
171 fast = self.generate_fast_HTI_PDE_solution(domain)
172
173 self.assertLess(Lsup(fast - slow), 1e-12*Lsup(slow)) #comparison between them
174
175 def run_VTI_assembly(self, domain):
176 model = VTIWave(domain, self.V_p, self.V_s, self.wavelet, "source",
177 source_vector=[0,0,1], eps=0., gamma=0., delta=0.,
178 rho=2000., absorption_zone=None, lumping=True,
179 disable_fast_assemblers=True)
180 self.assertFalse(model.fastAssembler) #ensure the arg is obeyed
181
182 model = VTIWave(domain, self.V_p, self.V_s, self.wavelet, "source",
183 source_vector=[0,0,1], eps=0., gamma=0., delta=0.,
184 rho=2000., absorption_zone=None, lumping=True)
185 self.assertTrue(model.fastAssembler) #ensure fast is actually used
186
187 slow = self.generate_slow_VTI_PDE_solution(domain)
188 fast = self.generate_fast_VTI_PDE_solution(domain)
189
190 self.assertLess(Lsup(fast - slow), 1e-12*Lsup(slow)) #comparison between them
191
192 def test_Function_params(self):
193 for domain in self.domains:
194 self.V_p = Data(2500., (), Function(domain))
195 self.V_s = Data(1250., (), Function(domain))
196 self.c11 = Data(11., (), Function(domain))
197 self.c12 = Data(12., (), Function(domain))
198 self.c13 = Data(13., (), Function(domain))
199 self.c23 = Data(23., (), Function(domain))
200 self.c33 = Data(33., (), Function(domain))
201 self.c44 = Data(44., (), Function(domain))
202 self.c66 = Data(66., (), Function(domain))
203 for i in [self.V_p, self.V_s, self.c11, self.c12, self.c13, self.c23,
204 self.c33, self.c44, self.c66]:
205 i.expand()
206 with self.assertRaises(RuntimeError) as e:
207 self.run_HTI_assembly(domain)
208 self.assertTrue("C tensor elements must be reduced" in str(e.exception))
209 with self.assertRaises(RuntimeError) as e:
210 self.run_VTI_assembly(domain)
211 self.assertTrue("C tensor elements must be reduced" in str(e.exception))
212
213 def test_ReducedFunction_params(self):
214 for domain in self.domains:
215 self.V_p = Data(2500., (), ReducedFunction(domain))
216 self.V_s = Data(1250., (), ReducedFunction(domain))
217 self.c11 = Data(11., (), ReducedFunction(domain))
218 self.c12 = Data(12., (), ReducedFunction(domain))
219 self.c13 = Data(13., (), ReducedFunction(domain))
220 self.c23 = Data(23., (), ReducedFunction(domain))
221 self.c33 = Data(33., (), ReducedFunction(domain))
222 self.c44 = Data(44., (), ReducedFunction(domain))
223 self.c66 = Data(66., (), ReducedFunction(domain))
224 for i in [self.V_p, self.V_s, self.c11, self.c12, self.c13, self.c23,
225 self.c33, self.c44, self.c66]:
226 i.expand()
227 self.run_HTI_assembly(domain)
228 self.run_VTI_assembly(domain)
229
230 def test_Constant_params(self):
231 for domain in self.domains:
232 self.V_p = Scalar(2500., ReducedFunction(domain))
233 self.V_s = Scalar(1250., ReducedFunction(domain))
234 self.c11 = Scalar(11., ReducedFunction(domain))
235 self.c12 = Scalar(12., ReducedFunction(domain))
236 self.c13 = Scalar(13., ReducedFunction(domain))
237 self.c23 = Scalar(23., ReducedFunction(domain))
238 self.c33 = Scalar(33., ReducedFunction(domain))
239 self.c44 = Scalar(44., ReducedFunction(domain))
240 self.c66 = Scalar(66., ReducedFunction(domain))
241 self.run_HTI_assembly(domain)
242 self.run_VTI_assembly(domain)
243
244 class Test_SpeckleyWaveAssembler2D(SpeckleyWaveAssemblerTestBase):
245 def setUp(self):
246 self.domains = []
247 for order in range(2,11):
248 self.domains.append(Rectangle(order,10,10,l0=100,l1=100,diracTags=["source"],
249 diracPoints=[(0,0)]))
250 self.wavelet = Ricker(100.)
251
252 def tearDown(self):
253 del self.domains
254
255 class Test_SpeckleyWaveAssembler3D(SpeckleyWaveAssemblerTestBase):
256 def setUp(self):
257 self.domains = []
258 for order in range(2,11):
259 self.domains.append(Brick(order, 4,4,4,
260 diracTags=["source"], diracPoints=[(0,0,0)]))
261 self.wavelet = Ricker(100.)
262
263 def tearDown(self):
264 del self.domains
265
266
267 if __name__ == '__main__':
268 run_tests(__name__, exit_on_failure=True)
269

  ViewVC Help
Powered by ViewVC 1.1.26