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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4433 - (show annotations)
Fri May 31 12:09:58 2013 UTC (6 years, 3 months ago) by gross
File MIME type: text/x-python
File size: 12312 byte(s)
some clarifications on geodetic coordinates. 
order of background magnetic flux density component has been corrected: input is now B_east, B_north, B_vertical.


1
2 ##############################################################################
3 #
4 # Copyright (c) 2003-2013 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 since 2012 by School of Earth Sciences
13 #
14 ##############################################################################
15
16 __copyright__="""Copyright (c) 2003-2013 by University of Queensland
17 http://www.uq.edu.au
18 Primary Business: Queensland, Australia"""
19 __license__="""Licensed under the Open Software License version 3.0
20 http://www.opensource.org/licenses/osl-3.0.php"""
21 __url__="https://launchpad.net/escript-finley"
22
23
24
25 import logging
26 import unittest
27 from esys.downunder.coordinates import *
28 from esys.ripley import Brick, Rectangle
29 import esys.escript
30 from esys.escript import Lsup, sin, cos
31 import esys.escript.unitsSI as U
32 import os
33 import sys
34 from math import pi
35 # this is mainly to avoid warning messages
36 logging.basicConfig(format='%(name)s: %(message)s', level=logging.INFO)
37
38 try:
39 WORKDIR=os.environ['DOWNUNDER_WORKDIR']
40 except KeyError:
41 WORKDIR='.'
42
43 NE=20
44 RTOL=1e-8
45
46 class TestCoordinates(unittest.TestCase):
47 def test_CartesianReferenceSystem(self):
48 cs=CartesianReferenceSystem()
49
50 self.assertEqual(cs.getName(), "CARTESIAN", "wrong name")
51 self.assertTrue(str(cs).startswith("CARTESIAN (id"), "wrong str")
52 self.assertTrue(cs == cs, "wrong cs == cs")
53 self.assertFalse(cs != cs, "wrong cs == cs")
54 self.assertFalse(cs == None, "wrong cs == None")
55 self.assertTrue(cs != None, "wrong cs == None")
56
57 cs2=CartesianReferenceSystem()
58 self.assertTrue(cs == cs2, "wrong cs == cs2")
59 self.assertFalse(cs != cs2, "wrong cs == cs2")
60 self.assertTrue(cs2 == cs, "wrong cs2 == cs")
61 self.assertFalse(cs2 != cs, "wrong cs2 == cs")
62
63 def test_GeodeticReferenceSystem(self):
64
65 cs = GeodeticReferenceSystem(a= 4*U.m, f=0.75, angular_unit=0.1, name="A1")
66 self.assertEqual(cs.getName(), "A1", "wrong name")
67 self.assertEqual(cs.getAngularUnit(), 0.1, "wrong angular_unit")
68 self.assertEqual(cs.getFlattening(), 0.75, "wrong flattening")
69 self.assertEqual(cs.getSemiMajorAxis(), 4., "wrong SemiMajorAxis")
70 self.assertEqual(cs.getSemiMinorAxis(), 1., "wrong SemiMinorAxis")
71 self.assertTrue(cs == cs, "wrong cs == cs")
72 self.assertFalse(cs != cs, "wrong cs == cs")
73 self.assertFalse(cs == None, "wrong cs == None")
74 self.assertTrue(cs != None, "wrong cs == None")
75
76 # lets try a few stupid things:
77 self.assertRaises(ValueError, GeodeticReferenceSystem, a= -4)
78 self.assertRaises(ValueError, GeodeticReferenceSystem, f=-0.75)
79 self.assertRaises(ValueError, GeodeticReferenceSystem, f=1.)
80 self.assertRaises(ValueError, GeodeticReferenceSystem, angular_unit=-0.1)
81
82 # cs == CartesianReferenceSystem() ?
83 cs2=CartesianReferenceSystem()
84 self.assertFalse(cs == cs2, "wrong cs == cs2")
85 self.assertTrue(cs != cs2, "wrong cs == cs2")
86 self.assertFalse(cs2 == cs, "wrong cs2 == cs")
87 self.assertTrue(cs2 != cs, "wrong cs2 == cs")
88
89 cs2 = GeodeticReferenceSystem(a= 4*U.m, f=0.75, angular_unit=0.1, name="A2")
90 self.assertTrue(cs == cs2, "wrong cs == cs2")
91 self.assertFalse(cs != cs2, "wrong cs == cs2")
92 self.assertTrue(cs2 == cs, "wrong cs2 == cs")
93 self.assertFalse(cs2 != cs, "wrong cs2 == cs")
94
95 # different semi majr axis?
96 cs2 = GeodeticReferenceSystem(a= 2*U.m, f=0.75, angular_unit=0.1, name="A2")
97 self.assertFalse(cs == cs2, "wrong cs == cs2")
98 self.assertTrue(cs != cs2, "wrong cs == cs2")
99 self.assertFalse(cs2 == cs, "wrong cs2 == cs")
100 self.assertTrue(cs2 != cs, "wrong cs2 == cs")
101
102 # different flattening
103 cs2 = GeodeticReferenceSystem(a= 4*U.m, f=0.7, angular_unit=0.1, name="A2")
104 self.assertFalse(cs == cs2, "wrong cs == cs2")
105 self.assertTrue(cs != cs2, "wrong cs == cs2")
106 self.assertFalse(cs2 == cs, "wrong cs2 == cs")
107 self.assertTrue(cs2 != cs, "wrong cs2 == cs")
108
109 # different angular_unit
110 cs2 = GeodeticReferenceSystem(a= 4*U.m, f=0.75, angular_unit=0.2, name="A2")
111 self.assertFalse(cs == cs2, "wrong cs == cs2")
112 self.assertTrue(cs != cs2, "wrong cs == cs2")
113 self.assertFalse(cs2 == cs, "wrong cs2 == cs")
114 self.assertTrue(cs2 != cs, "wrong cs2 == cs")
115
116 sp=SphericalReferenceSystem(R=1*U.km)
117 self.assertEqual(sp.getName(), "SPHERE", "wrong name")
118 self.assertAlmostEqual(sp.getAngularUnit(), 1*U.DEG, msg="wrong angular_unit")
119 self.assertAlmostEqual(sp.getFlattening(), 0., msg="wrong flattening")
120 self.assertAlmostEqual(sp.getSemiMajorAxis(), 1000.,msg= "wrong SemiMajorAxis")
121
122 sp=SphericalReferenceSystem()
123 self.assertEqual(sp.getName(), "SPHERE", msg= "wrong name")
124 self.assertAlmostEqual(sp.getAngularUnit(), 1*U.DEG, msg="wrong angular_unit")
125 self.assertAlmostEqual(sp.getFlattening(), 0., msg="wrong flattening")
126 self.assertAlmostEqual(sp.getSemiMajorAxis(), 6378137, msg="wrong SemiMajorAxis")
127
128
129 sp=WGS84ReferenceSystem()
130 self.assertEqual(sp.getName(), "WGS84", msg= "wrong name")
131 self.assertAlmostEqual(sp.getAngularUnit(), 1*U.DEG, msg="wrong angular_unit")
132 self.assertAlmostEqual(sp.getFlattening(), 0.0033528106647474805, msg="wrong flattening")
133 self.assertAlmostEqual(sp.getSemiMajorAxis(), 6378137, msg="wrong SemiMajorAxis")
134
135 sp=GRS80ReferenceSystem()
136 self.assertEqual(sp.getName(), "GRS80", msg= "wrong name")
137 self.assertAlmostEqual(sp.getAngularUnit(), 1*U.DEG, msg="wrong angular_unit")
138 self.assertAlmostEqual(sp.getFlattening(), 0.003352810681182319, msg="wrong flattening")
139 self.assertAlmostEqual(sp.getSemiMajorAxis(), 6378137, msg="wrong SemiMajorAxis")
140
141 def test_CartesianTransformation3D(self):
142
143 dom=Brick(NE,NE,NE, l0=10, l1=10, l2=10)
144
145 cs=CartesianReferenceSystem()
146 tf=cs.createTransformation(dom)
147
148 self.assertEqual(tf.getReferenceSystem(),cs , "wrong reference")
149 self.assertEqual(tf.getDomain(),dom , "wrong reference")
150 self.assertTrue(tf.isCartesian(), "wrong isCartesian check")
151
152
153 v=tf.getVolumeFactor()
154 self.assertTrue(isinstance(v, esys.escript.Data), "wrong volume factor type")
155 self.assertEqual(v.getFunctionSpace(), esys.escript.Function(dom), "wrong volume factor type")
156 error=Lsup(v-1.)
157 self.assertTrue(error<=RTOL, "volume factor")
158
159 s=tf.getScalingFactors()
160 self.assertTrue(isinstance(s, esys.escript.Data), "scaling factor type")
161 self.assertEqual(s.getShape(), (dom.getDim(),), "scaling factor length")
162 self.assertEqual(s.getFunctionSpace(), esys.escript.Function(dom), "wrong 0-th scaling factor type")
163
164 error=Lsup(s[0]-1.)
165 self.assertTrue(error<=RTOL, "0-th scaling factor")
166 error=Lsup(s[1]-1.)
167 self.assertTrue(error<=RTOL, "1-th scaling factor")
168 error=Lsup(s[2]-1.)
169 self.assertTrue(error<=RTOL, "2-th scaling factor")
170
171 def test_CartesianTransformation2D(self):
172
173 dom=Rectangle(NE,NE, l0=10, l1=10)
174
175 cs=CartesianReferenceSystem()
176 tf=cs.createTransformation(dom)
177
178 self.assertEqual(tf.getReferenceSystem(),cs , "wrong reference")
179 self.assertEqual(tf.getDomain(),dom , "wrong reference")
180 self.assertTrue(tf.isCartesian(), "wrong isCartesian check")
181
182
183 v=tf.getVolumeFactor()
184 self.assertTrue(isinstance(v, esys.escript.Data), "wrong volume factor type")
185 self.assertEqual(v.getFunctionSpace(), esys.escript.Function(dom), "wrong volume factor type")
186 error=Lsup(v-1.)
187 self.assertTrue(error<=RTOL, "volume factor")
188
189 s=tf.getScalingFactors()
190 self.assertTrue(isinstance(s, esys.escript.Data), "scaling factor type")
191 self.assertEqual(s.getShape(), (dom.getDim(),), "scaling factor length")
192 self.assertEqual(s.getFunctionSpace(), esys.escript.Function(dom), "wrong 0-th scaling factor type")
193
194 error=Lsup(s[0]-1.)
195 self.assertTrue(error<=RTOL, "0-th scaling factor")
196
197 error=Lsup(s[1]-1.)
198 self.assertTrue(error<=RTOL, "1-th scaling factor")
199
200
201
202 def test_SpericalTransformation3D(self):
203
204 dom=Brick(NE,NE,NE, l0=90, l1=45, l2=10.)
205
206 cs=SphericalReferenceSystem()
207 tf=cs.createTransformation(dom)
208
209 self.assertEqual(tf.getReferenceSystem(),cs , "wrong reference")
210 self.assertEqual(tf.getDomain(),dom , "wrong reference")
211 self.assertFalse(tf.isCartesian(), "wrong isCartesian check")
212
213 R=6378137.0
214 x=esys.escript.Function(dom).getX()
215 phi=(90.-x[1])/180.*pi
216 lam=x[0]/180.*pi
217 h=x[2]*1000.
218 r=h+R
219
220 v=tf.getVolumeFactor()
221 self.assertTrue(isinstance(v, esys.escript.Data), "wrong volume factor type")
222 self.assertEqual(v.getFunctionSpace(), esys.escript.Function(dom), "wrong volume factor type")
223 error=Lsup(v- r**2*sin(phi)*(pi/180.)**2*1000. )
224 self.assertTrue(error<=RTOL * R*R*(pi/180.)**2*1000., "volume factor")
225
226 s=tf.getScalingFactors()
227 self.assertTrue(isinstance(s, esys.escript.Data), "scaling factor type")
228 self.assertEqual(s.getShape(), (dom.getDim(),), "scaling factor length")
229 self.assertEqual(s.getFunctionSpace(), esys.escript.Function(dom), "wrong 0-th scaling factor type")
230
231 error=Lsup(s[1]-1/r/pi*180.)
232 self.assertTrue(error<=RTOL/R/pi*180., "0-th scaling factor")
233
234 error=Lsup(s[0]-1/(r*sin(phi))/pi*180.)
235 self.assertTrue(error<=RTOL/R/pi*180., "1-th scaling factor")
236
237 error=Lsup(s[2]-1./1000.)
238 self.assertTrue(error<=RTOL/1000., "2-th scaling factor")
239
240 def test_SpericalTransformation2D(self):
241
242 dom=Rectangle(NE,NE, l0=45., l1=10.)
243
244 cs=SphericalReferenceSystem()
245 tf=cs.createTransformation(dom)
246
247 self.assertEqual(tf.getReferenceSystem(),cs , "wrong reference")
248 self.assertEqual(tf.getDomain(),dom , "wrong reference")
249 self.assertFalse(tf.isCartesian(), "wrong isCartesian check")
250
251 R=6378137.0
252 x=esys.escript.Function(dom).getX()
253 phi=(90.-x[0])/180.*pi
254 h=x[1]*1000.
255 r=h+R
256
257 v=tf.getVolumeFactor()
258 self.assertTrue(isinstance(v, esys.escript.Data), "wrong volume factor type")
259 self.assertEqual(v.getFunctionSpace(), esys.escript.Function(dom), "wrong volume factor type")
260 error=Lsup(v-r*pi/180.*1000.)
261 self.assertTrue(error<=RTOL*R, "volume factor")
262
263 s=tf.getScalingFactors()
264 self.assertTrue(isinstance(s, esys.escript.Data), "scaling factor type")
265 self.assertEqual(s.getShape(), (dom.getDim(),), "scaling factor length")
266 self.assertEqual(s.getFunctionSpace(), esys.escript.Function(dom), "wrong 0-th scaling factor type")
267
268 error=Lsup(s[0]-1./r/pi*180.)
269 self.assertTrue(error<=RTOL/R/pi*180., "0-th scaling factor")
270 error=Lsup(s[1]-1./1000)
271 self.assertTrue(error<=RTOL/1000., "1-th scaling factor")
272
273
274 if __name__ == "__main__":
275 suite = unittest.TestSuite()
276 suite.addTest(unittest.makeSuite(TestCoordinates))
277 s=unittest.TextTestRunner(verbosity=2).run(suite)
278 if not s.wasSuccessful(): sys.exit(1)
279

  ViewVC Help
Powered by ViewVC 1.1.26