From 19d976115a48f641b896939740284c2ad29dce0c Mon Sep 17 00:00:00 2001 From: Hao Gao Date: Sat, 20 Oct 2018 20:41:48 -0500 Subject: [PATCH 1/9] Raise exception on fp error for urchin --- meshmode/mesh/generation.py | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/meshmode/mesh/generation.py b/meshmode/mesh/generation.py index 5c91895f..f2195c79 100644 --- a/meshmode/mesh/generation.py +++ b/meshmode/mesh/generation.py @@ -574,11 +574,18 @@ def refine_mesh_and_get_urchin_warper(order, m, n, est_rel_interp_tolerance, for i in range(uniform_refinement_rounds): refiner.refine_uniformly() + # Calculating spherical harmonics can cause floating-point errors like overflow + # or NaN. Raise FloatingPointError when this occurs. + old_settings = np.geterr() + np.seterr(all='raise') + nodes_sph = sph_harm(m, n, unwarped_mesh.groups[0].nodes).real lo = np.min(nodes_sph) hi = np.max(nodes_sph) del nodes_sph + np.seterr(**old_settings) + from functools import partial unwarped_mesh = warp_and_refine_until_resolved( refiner, -- GitLab From dba0e63c2978be31b88cc522c037dbe00cc62ba8 Mon Sep 17 00:00:00 2001 From: Hao Gao Date: Sun, 21 Oct 2018 22:36:50 -0500 Subject: [PATCH 2/9] Revert "Raise exception on fp error for urchin" This reverts commit 19d976115a48f641b896939740284c2ad29dce0c. --- meshmode/mesh/generation.py | 7 ------- 1 file changed, 7 deletions(-) diff --git a/meshmode/mesh/generation.py b/meshmode/mesh/generation.py index f2195c79..5c91895f 100644 --- a/meshmode/mesh/generation.py +++ b/meshmode/mesh/generation.py @@ -574,18 +574,11 @@ def refine_mesh_and_get_urchin_warper(order, m, n, est_rel_interp_tolerance, for i in range(uniform_refinement_rounds): refiner.refine_uniformly() - # Calculating spherical harmonics can cause floating-point errors like overflow - # or NaN. Raise FloatingPointError when this occurs. - old_settings = np.geterr() - np.seterr(all='raise') - nodes_sph = sph_harm(m, n, unwarped_mesh.groups[0].nodes).real lo = np.min(nodes_sph) hi = np.max(nodes_sph) del nodes_sph - np.seterr(**old_settings) - from functools import partial unwarped_mesh = warp_and_refine_until_resolved( refiner, -- GitLab From 0c017330cc628b65fce84dba92fade62f81b1bc4 Mon Sep 17 00:00:00 2001 From: Hao Gao Date: Sun, 21 Oct 2018 23:16:11 -0500 Subject: [PATCH 3/9] Raise FloatingPointError when vertices are invalid in warped mesh --- meshmode/mesh/generation.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/meshmode/mesh/generation.py b/meshmode/mesh/generation.py index 5c91895f..6e7fb655 100644 --- a/meshmode/mesh/generation.py +++ b/meshmode/mesh/generation.py @@ -819,6 +819,10 @@ def warp_and_refine_until_resolved( warped_mesh = warp_callable(unwarped_mesh) + # test whether there are invalid values in warped mesh + if np.isnan(warped_mesh.vertices).any(): + raise FloatingPointError("Invalid vertices in warped mesh") + for egrp in warped_mesh.groups: dim, nunit_nodes = egrp.unit_nodes.shape -- GitLab From a086804ba1437d270a84d8dc7f2a2832d476c706 Mon Sep 17 00:00:00 2001 From: Hao Gao Date: Sun, 21 Oct 2018 23:40:40 -0500 Subject: [PATCH 4/9] Attempting to use pyfmmlib for sperical harmonics results in mem error [ci skip] --- meshmode/mesh/generation.py | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/meshmode/mesh/generation.py b/meshmode/mesh/generation.py index 6e7fb655..45aebdbf 100644 --- a/meshmode/mesh/generation.py +++ b/meshmode/mesh/generation.py @@ -542,8 +542,13 @@ def refine_mesh_and_get_urchin_warper(order, m, n, est_rel_interp_tolerance, theta = np.arccos(z/r) phi = np.arctan2(y, x) - import scipy.special as sps - return sps.sph_harm(m, n, phi, theta) + try: + from pyfmmlib import ylgndr_vec + legendre = ylgndr_vec(n, np.cos(phi))[n][m].reshape(phi.shape) + return legendre * np.sqrt(1/(4 * np.pi)) * np.exp(m * theta * 1j) + except ImportError: + import scipy.special as sps + return sps.sph_harm(m, n, phi, theta) def map_coords(pts): r = np.sqrt(np.sum(pts**2, axis=0)) -- GitLab From 0242ace1cdd3e18cf7f56048b544542aca19ee74 Mon Sep 17 00:00:00 2001 From: Hao Gao Date: Mon, 22 Oct 2018 09:03:54 -0500 Subject: [PATCH 5/9] Make error message more concrete [ci skip] --- meshmode/mesh/generation.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/meshmode/mesh/generation.py b/meshmode/mesh/generation.py index 45aebdbf..51229568 100644 --- a/meshmode/mesh/generation.py +++ b/meshmode/mesh/generation.py @@ -825,8 +825,9 @@ def warp_and_refine_until_resolved( warped_mesh = warp_callable(unwarped_mesh) # test whether there are invalid values in warped mesh - if np.isnan(warped_mesh.vertices).any(): - raise FloatingPointError("Invalid vertices in warped mesh") + if not np.isfinite(warped_mesh.vertices).all(): + raise FloatingPointError("Warped mesh contains non-finite vertices " + "(NaN or Inf)") for egrp in warped_mesh.groups: dim, nunit_nodes = egrp.unit_nodes.shape -- GitLab From 12730d086be48256e636518330d4948a7bb0988f Mon Sep 17 00:00:00 2001 From: Hao Gao Date: Mon, 22 Oct 2018 09:12:57 -0500 Subject: [PATCH 6/9] Add finiteness check for nodes in each group of the warped mesh [ci skip] --- meshmode/mesh/generation.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/meshmode/mesh/generation.py b/meshmode/mesh/generation.py index 51229568..4017cd99 100644 --- a/meshmode/mesh/generation.py +++ b/meshmode/mesh/generation.py @@ -829,6 +829,11 @@ def warp_and_refine_until_resolved( raise FloatingPointError("Warped mesh contains non-finite vertices " "(NaN or Inf)") + for group in warped_mesh.groups: + if not np.isfinite(group.nodes).all(): + raise FloatingPointError("Warped mesh contains non-finite nodes " + "(NaN or Inf)") + for egrp in warped_mesh.groups: dim, nunit_nodes = egrp.unit_nodes.shape -- GitLab From 749c3640df860743870d0fecbfaa9056d49a0f1b Mon Sep 17 00:00:00 2001 From: Hao Gao Date: Sun, 28 Oct 2018 10:46:22 -0500 Subject: [PATCH 7/9] Revert back to pure scipy implementation --- meshmode/mesh/generation.py | 10 +++------- 1 file changed, 3 insertions(+), 7 deletions(-) diff --git a/meshmode/mesh/generation.py b/meshmode/mesh/generation.py index 4017cd99..ed806a0f 100644 --- a/meshmode/mesh/generation.py +++ b/meshmode/mesh/generation.py @@ -542,13 +542,9 @@ def refine_mesh_and_get_urchin_warper(order, m, n, est_rel_interp_tolerance, theta = np.arccos(z/r) phi = np.arctan2(y, x) - try: - from pyfmmlib import ylgndr_vec - legendre = ylgndr_vec(n, np.cos(phi))[n][m].reshape(phi.shape) - return legendre * np.sqrt(1/(4 * np.pi)) * np.exp(m * theta * 1j) - except ImportError: - import scipy.special as sps - return sps.sph_harm(m, n, phi, theta) + import scipy.special as sps + # Note phi and theta are interchanged + return sps.sph_harm(m, n, phi, theta) def map_coords(pts): r = np.sqrt(np.sum(pts**2, axis=0)) -- GitLab From 7c86e956d52ecb52760eeebe921f55c2d75d02c1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andreas=20Kl=C3=B6ckner?= Date: Mon, 29 Oct 2018 11:42:02 -0400 Subject: [PATCH 8/9] Update meshmode/mesh/generation.py --- meshmode/mesh/generation.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/meshmode/mesh/generation.py b/meshmode/mesh/generation.py index ed806a0f..5d09c651 100644 --- a/meshmode/mesh/generation.py +++ b/meshmode/mesh/generation.py @@ -543,7 +543,9 @@ def refine_mesh_and_get_urchin_warper(order, m, n, est_rel_interp_tolerance, phi = np.arctan2(y, x) import scipy.special as sps - # Note phi and theta are interchanged + # Note: numpy arguments are called (theta, phi). + # This matches the notation in the QBX3D paper? + ### FIXME return sps.sph_harm(m, n, phi, theta) def map_coords(pts): -- GitLab From 2527164a1e544e6ed5a2b714f559a67dd314f6e9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andreas=20Kl=C3=B6ckner?= Date: Tue, 30 Oct 2018 13:23:13 -0400 Subject: [PATCH 9/9] Update comment on spherical harmonic convention in urchin generation --- meshmode/mesh/generation.py | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/meshmode/mesh/generation.py b/meshmode/mesh/generation.py index 5d09c651..ddecc69f 100644 --- a/meshmode/mesh/generation.py +++ b/meshmode/mesh/generation.py @@ -543,9 +543,13 @@ def refine_mesh_and_get_urchin_warper(order, m, n, est_rel_interp_tolerance, phi = np.arctan2(y, x) import scipy.special as sps - # Note: numpy arguments are called (theta, phi). - # This matches the notation in the QBX3D paper? - ### FIXME + # Note: This matches the spherical harmonic + # convention in the QBX3D paper: + # https://arxiv.org/abs/1805.06106 + # + # Numpy takes arguments in the order (theta, phi) + # *and* swaps their meanings, so passing the + # arguments swapped maintains the intended meaning. return sps.sph_harm(m, n, phi, theta) def map_coords(pts): -- GitLab