Merge branch 'fix-triangle-division-by-zero'

tests/unit_tests.py

@@ -164,6 +164,47 @@ def test(a, b, c):
 
         for (a,b,c) in rand(3, 3):
             test(a,b,c)
+
+    def test_potential_triangle_close(self):
+        v0 = np.array([0.0, 0.0, 0.0])
+        v1 = np.array([1.0, 0.0, 0.0])
+        v2 = np.array([0.0, 1.0, 0.0])
+
+        for x in range(-20, 21):
+            target = np.array([x, 0.0, 0.5])
+            assert np.isclose(B.potential_triangle(v0, v1, v2, target), potential_exact_integrated(v0, v1, v2, target), atol=0.0, rtol=1e-9)
+
+        for x in range(-20, 21):
+            target = np.array([x, 0.0, -0.5])
+            assert np.isclose(B.potential_triangle(v0, v1, v2, target), potential_exact_integrated(v0, v1, v2, target), atol=0.0, rtol=1e-9)
+
+        for y in range(-20, 21):
+            target = np.array([0.0, y, 0.5])
+            assert np.isclose(B.potential_triangle(v0, v1, v2, target), potential_exact_integrated(v0, v1, v2, target), atol=0.0, rtol=1e-9)
+
+        for y in range(-20, 21):
+            target = np.array([0.0, y, -0.5])
+            assert np.isclose(B.potential_triangle(v0, v1, v2, target), potential_exact_integrated(v0, v1, v2, target), atol=0.0, rtol=1e-9)
+
+        for z in range(-20, 21):
+            target = np.array([1.0, 1.0, z])
+            assert np.isclose(B.potential_triangle(v0, v1, v2, target), potential_exact_integrated(v0, v1, v2, target), atol=0.0, rtol=1e-9)
+
+        for z in range(-20, 21):
+            target = np.array([-1.0, -1.0, z])
+            assert np.isclose(B.potential_triangle(v0, v1, v2, target), potential_exact_integrated(v0, v1, v2, target), atol=0.0, rtol=1e-9)
+
+        for z in range(-20, 21):
+            target = np.array([-0.5, 0.5, z])
+            assert np.isclose(B.potential_triangle(v0, v1, v2, target), potential_exact_integrated(v0, v1, v2, target), atol=0.0, rtol=1e-9)
+
+        for z in range(-20, 21):
+            target = np.array([0.5, -0.5, z])
+            assert np.isclose(B.potential_triangle(v0, v1, v2, target), potential_exact_integrated(v0, v1, v2, target), atol=0.0, rtol=1e-9)
+
+        for k in range(-20, 21):
+            target = np.array([k, k, 0.5])
+            assert np.isclose(B.potential_triangle(v0, v1, v2, target), potential_exact_integrated(v0, v1, v2, target), atol=0.0, rtol=1e-9)
 
     def test_derivative_x_quadrants(self):
         def test(a, b, c, z0):
@@ -313,6 +354,54 @@ def test_potential(self):
             correct = potential_exact_integrated(v0, v1, v2, target) 
             approx = B.potential_triangle(v0, v1, v2, target)
             assert np.allclose(correct, approx)
+
+    def test_flux_triangle_close(self):
+        v0 = np.array([0.0, 0.0, 0.0])
+        v1 = np.array([1.0, 0.0, 0.0])
+        v2 = np.array([0.0, 1.0, 0.0])
+
+        normals = [
+            np.array([1.0, 0.0, 0.0]),
+            np.array([0.0, 1.0, 0.0]),
+            np.array([0.0, 0.0, 1.0])
+        ]
+
+        for normal in normals:
+            for x in range(-20, 21):
+                target = np.array([x, 0.0, 0.5])
+                assert np.isclose(B.flux_triangle(v0, v1, v2, target, normal), flux_exact_integrated(v0, v1, v2, target, normal), atol=0.0, rtol=1e-8)
+
+            for x in range(-20, 21):
+                target = np.array([x, 0.0, -0.5])
+                assert np.isclose(B.flux_triangle(v0, v1, v2, target, normal), flux_exact_integrated(v0, v1, v2, target, normal), atol=0.0, rtol=1e-8)
+
+            for y in range(-20, 21):
+                target = np.array([0.0, y, 0.5])
+                assert np.isclose(B.flux_triangle(v0, v1, v2, target, normal), flux_exact_integrated(v0, v1, v2, target, normal), atol=0.0, rtol=1e-8)
+
+            for y in range(-20, 21):
+                target = np.array([0.0, y, -0.5])
+                assert np.isclose(B.flux_triangle(v0, v1, v2, target, normal), flux_exact_integrated(v0, v1, v2, target, normal), atol=0.0, rtol=1e-8)
+
+            for z in range(-20, 21):
+                target = np.array([1.0, 1.0, z])
+                assert np.isclose(B.flux_triangle(v0, v1, v2, target, normal), flux_exact_integrated(v0, v1, v2, target, normal), atol=0.0, rtol=1e-8)
+
+            for z in range(-20, 21):
+                target = np.array([-1.0, -1.0, z])
+                assert np.isclose(B.flux_triangle(v0, v1, v2, target, normal), flux_exact_integrated(v0, v1, v2, target, normal), atol=0.0, rtol=1e-8)
+
+            for z in range(-20, 21):
+                target = np.array([-0.5, 0.5, z])
+                assert np.isclose(B.flux_triangle(v0, v1, v2, target, normal), flux_exact_integrated(v0, v1, v2, target, normal), atol=0.0, rtol=1e-8)
+
+            for z in range(-20, 21):
+                target = np.array([0.5, -0.5, z])
+                assert np.isclose(B.flux_triangle(v0, v1, v2, target, normal), flux_exact_integrated(v0, v1, v2, target, normal), atol=0.0, rtol=1e-8)
+
+            for k in range(-20, 21):
+                target = np.array([k, k, 0.5])
+                assert np.isclose(B.flux_triangle(v0, v1, v2, target, normal), flux_exact_integrated(v0, v1, v2, target, normal), atol=0.0, rtol=1e-8)
 
     def test_flux(self):
         for _ in range(10):

traceon/backend/triangle_contribution.c

@@ -40,6 +40,13 @@ double _potential_integrand(double y, void *args_p) {
 	double xmax = args.x0 + args.a + y/args.c*(args.b-args.a);
 
 	double denom = sqrt((y+args.y0)*(y+args.y0) + args.z*args.z);
+
+	if(denom < 1e-12) {
+		// The asinh(xmax/denom) - asinh(xmin/denom) is numerical 
+		// unstable when denom is small. Taking the taylor expansion
+		// of denom -> 0 we find
+		return log(xmax) - log(xmin);
+	}
     return asinh(xmax/denom) - asinh(xmin/denom);
 }
 
@@ -139,11 +146,23 @@ double _flux_integrand(double y, void *args_p) {
 	double xmin2 = xmin*xmin;
 	double yy02 = (y+y0)*(y+y0);
 	double xmax2 = xmax*xmax;
+	double r2 = z2 + yy02;
 
 	double flux[3];
-    flux[0] = 1/sqrt(z2+yy02+xmax2) - 1/sqrt(z2+yy02+xmin2);
-    flux[1] = -(xmax*(y+y0))/((z2+yy02)*sqrt(z2+yy02+xmax2)) + (xmin*(y+y0))/((z2+yy02)*sqrt(z2+yy02+xmin2));
-    flux[2] = (xmax*z)/((z2+yy02)*sqrt(z2+yy02+xmax2)) - (xmin*z)/((z2+yy02)*sqrt(z2+yy02+xmin2));
+
+	flux[0] = 1 / sqrt(r2 + xmax2) - 1 / sqrt(r2 + xmin2);
+
+	// Singularity when r2 is small...
+	if (fabs(r2) < 1e-9) {
+		flux[1] = ((xmin2 - xmax2) * y0 + (xmin2 - xmax2) * y) / (2.0 * xmax2 * xmin2);
+		flux[2] = -((xmin2 - xmax2) * z) / (2.0 * xmax2 * xmin2);
+	} else {
+		double denom_max = r2 * sqrt(r2 + xmax2);
+		double denom_min = r2 * sqrt(r2 + xmin2);
+
+		flux[1] = -((xmax * (y + y0)) / denom_max) + (xmin * (y + y0)) / denom_min;
+		flux[2] = (xmax * z) / denom_max - (xmin * z) / denom_min;
+	}
 
 	return dot_3d(args.normal, flux);
 }