[yt-svn] commit/yt: ngoldbaum: Merged in jzuhone/yt (pull request #2466)

[commits-noreply at bitbucket.org]

1 new commit in yt:

https://bitbucket.org/yt_analysis/yt/commits/ad6d4f04b990/
Changeset:   ad6d4f04b990
Branch:      yt
User:        ngoldbaum
Date:        2016-12-15 17:57:34+00:00
Summary:     Merged in jzuhone/yt (pull request #2466)

Do exact cell volume computation for curvilinear coordinates
Affected #:  5 files

diff -r 208722c07b452c9121fae826a53f42b5ce81d802 -r ad6d4f04b990edfeae34479a220c6629df074070 yt/geometry/coordinates/cylindrical_coordinates.py
--- a/yt/geometry/coordinates/cylindrical_coordinates.py
+++ b/yt/geometry/coordinates/cylindrical_coordinates.py
@@ -68,11 +68,13 @@
                            units = "")
 
         def _CylindricalVolume(field, data):
-            return data["index", "dtheta"] \
-                 * data["index", "r"] \
-                 * data["index", "dr"] \
-                 * data["index", "dz"]
-        registry.add_field(("index", "cell_volume"), sampling_type="cell", 
+            r = data["index", "r"]
+            dr = data["index", "dr"]
+            vol = 0.5*((r+0.5*dr)**2-(r-0.5*dr)**2)
+            vol *= data["index", "dtheta"]
+            vol *= data["index", "dz"]
+            return vol
+        registry.add_field(("index", "cell_volume"), sampling_type="cell",
                  function=_CylindricalVolume,
                  units = "code_length**3")
 

diff -r 208722c07b452c9121fae826a53f42b5ce81d802 -r ad6d4f04b990edfeae34479a220c6629df074070 yt/geometry/coordinates/geographic_coordinates.py
--- a/yt/geometry/coordinates/geographic_coordinates.py
+++ b/yt/geometry/coordinates/geographic_coordinates.py
@@ -68,15 +68,17 @@
                            units = "code_length")
 
         def _SphericalVolume(field, data):
-            # r**2 sin theta dr dtheta dphi
             # We can use the transformed coordinates here.
-            vol = data["index", "r"]**2.0
-            vol *= data["index", "dr"]
-            vol *= np.sin(data["index", "theta"])
-            vol *= data["index", "dtheta"]
+            # Here we compute the spherical volume element exactly
+            r = data["index", "r"]
+            dr = data["index", "dr"]
+            theta = data["index", "theta"]
+            dtheta = data["index", "dtheta"]
+            vol = ((r+0.5*dr)**3-(r-0.5*dr)**3) / 3.0
+            vol *= np.cos(theta-0.5*dtheta)-np.cos(theta+0.5*dtheta)
             vol *= data["index", "dphi"]
             return vol
-        registry.add_field(("index", "cell_volume"), sampling_type="cell", 
+        registry.add_field(("index", "cell_volume"), sampling_type="cell",
                  function=_SphericalVolume,
                  units = "code_length**3")
 

diff -r 208722c07b452c9121fae826a53f42b5ce81d802 -r ad6d4f04b990edfeae34479a220c6629df074070 yt/geometry/coordinates/spherical_coordinates.py
--- a/yt/geometry/coordinates/spherical_coordinates.py
+++ b/yt/geometry/coordinates/spherical_coordinates.py
@@ -66,11 +66,13 @@
                            units = "")
 
         def _SphericalVolume(field, data):
-            # r**2 sin theta dr dtheta dphi
-            vol = data["index", "r"]**2.0
-            vol *= data["index", "dr"]
-            vol *= np.sin(data["index", "theta"])
-            vol *= data["index", "dtheta"]
+            # Here we compute the spherical volume element exactly
+            r = data["index", "r"]
+            dr = data["index", "dr"]
+            theta = data["index", "theta"]
+            dtheta = data["index", "dtheta"]
+            vol = ((r+0.5*dr)**3-(r-0.5*dr)**3)/3.0
+            vol *= np.cos(theta-0.5*dtheta)-np.cos(theta+0.5*dtheta)
             vol *= data["index", "dphi"]
             return vol
         registry.add_field(("index", "cell_volume"), sampling_type="cell", 

diff -r 208722c07b452c9121fae826a53f42b5ce81d802 -r ad6d4f04b990edfeae34479a220c6629df074070 yt/geometry/coordinates/tests/test_geographic_coordinates.py
--- a/yt/geometry/coordinates/tests/test_geographic_coordinates.py
+++ b/yt/geometry/coordinates/tests/test_geographic_coordinates.py
@@ -80,10 +80,8 @@
                         dd["index","dlatitude"]*np.pi/180.0
     yield assert_equal, dd["index","dphi"], \
                         dd["index","dlongitude"]*np.pi/180.0
-    # Note our terrible agreement here.
     yield assert_rel_equal, dd["cell_volume"].sum(dtype="float64"), \
-                        (4.0/3.0) * np.pi * (outer_r**3 - inner_r**3), \
-                        3
+                        (4.0/3.0) * np.pi * (outer_r**3 - inner_r**3), 10
     yield assert_equal, dd["index", "path_element_depth"], \
                         dd["index", "ddepth"]
     yield assert_equal, dd["index", "path_element_depth"], \

diff -r 208722c07b452c9121fae826a53f42b5ce81d802 -r ad6d4f04b990edfeae34479a220c6629df074070 yt/geometry/coordinates/tests/test_spherical_coordinates.py
--- a/yt/geometry/coordinates/tests/test_spherical_coordinates.py
+++ b/yt/geometry/coordinates/tests/test_spherical_coordinates.py
@@ -29,8 +29,7 @@
     # don't think it is avoidable as of right now.  Real datasets will almost
     # certainly be correct, if this is correct to 3 decimel places.
     yield assert_almost_equal, dd["cell_volume"].sum(dtype="float64"), \
-                        (4.0/3.0) * np.pi*ds.domain_width[0]**3, \
-                        3 # Allow for GENEROUS error on the volume ...
+                        (4.0/3.0) * np.pi*ds.domain_width[0]**3
     yield assert_equal, dd["index", "path_element_r"], dd["index", "dr"]
     yield assert_equal, dd["index", "path_element_theta"], \
                         dd["index", "r"] * dd["index", "dtheta"]

Repository URL: https://bitbucket.org/yt_analysis/yt/

--

This is a commit notification from bitbucket.org. You are receiving
this because you have the service enabled, addressing the recipient of
this email.