From cf39a15ef701a3748bfcd5516da02f80b2a644da Mon Sep 17 00:00:00 2001 From: rockfool Date: Thu, 30 Jan 2020 17:24:37 -0500 Subject: [PATCH 01/10] fix a bug in ZernikeRadial --- openmc/polynomial.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/openmc/polynomial.py b/openmc/polynomial.py index abcfd8057e..6257526a5c 100644 --- a/openmc/polynomial.py +++ b/openmc/polynomial.py @@ -75,7 +75,9 @@ class ZernikeRadial(Polynomial): def __call__(self, r): import openmc.lib as lib if isinstance(r, Iterable): - return [np.sum(self._norm_coef * lib.calc_zn_rad(self.order, r_i)) + return [np.sum(self._norm_coef * lib.calc_zn_rad(self.order, r_i / self.radius)) for r_i in r] else: - return np.sum(self._norm_coef * lib.calc_zn_rad(self.order, r)) + return np.sum(self._norm_coef * lib.calc_zn_rad(self.order, r / self.radius)) + + From 49f715fee5e1b57d8d79e7d6da4ca6b2829c87d7 Mon Sep 17 00:00:00 2001 From: rockfool Date: Thu, 30 Jan 2020 17:25:21 -0500 Subject: [PATCH 02/10] add reconstruction for Azimuthal Zernike --- openmc/polynomial.py | 60 ++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 60 insertions(+) diff --git a/openmc/polynomial.py b/openmc/polynomial.py index 6257526a5c..79b914b256 100644 --- a/openmc/polynomial.py +++ b/openmc/polynomial.py @@ -81,3 +81,63 @@ class ZernikeRadial(Polynomial): return np.sum(self._norm_coef * lib.calc_zn_rad(self.order, r / self.radius)) +class Zernike(Polynomial): + """Create Zernike polynomials given coefficients and domain. + The Zernike polynomials are defined as in + :class:`ZernikeFilter`. + Parameters + ---------- + coef : Iterable of float + A list of coefficients of each term in radial only Zernike polynomials + radius : float + Domain of Zernike polynomials to be applied on. Default is 1. + r : float + Position to be evaluated, normalized on radius [0,1] + theta : float + Polar angular position to be evaluated, normalized on :math:`[0, 2*\pi]` + + Attributes + ---------- + order : int + The maximum (even) order of Zernike polynomials. + radius : float + Domain of Zernike polynomials to be applied on. Default is 1. + theta : float + Polar angle of Zernike polynomial to be applied on. Default is 0. + norm_coef : iterable of float + The list of coefficients of each term in the polynomials after + normailization. + """ + def __init__(self, coef, radius=1): + super().__init__(coef) + self._order = int(((np.sqrt(8 * len(self.coef) + 1)) -3) / 2) + self.radius = radius + norm_vec = np.ones((len(self.coef))) + for n in range(self._order + 1): + for m in range(-n,(n+1),2): + j = int((n * (n + 2) + m) / 2) + if m==0: + norm_vec[j] = n + 1.0 + else: + norm_vec[j] = 2 * n + 2 + norm_vec /= (np.pi * radius**2) + self._norm_coef = norm_vec * self.coef + + @property + def order(self): + return self._order + + def __call__(self, r, theta=0.0): + import openmc.lib as lib + if isinstance(r, Iterable) and isinstance(theta, Iterable): + return [[np.sum(self._norm_coef * lib.calc_zn(self.order, r_i/self.radius, theta_i)) + for r_i in r] for theta_i in theta] + elif isinstance(r, Iterable) and not isinstance(theta, Iterable): + return [np.sum(self._norm_coef * lib.calc_zn(self.order, r_i/self.radius, theta)) + for r_i in r] + elif not isinstance(r, Iterable) and isinstance(theta, Iterable): + return [np.sum(self._norm_coef * lib.calc_zn(self.order, r/self.radius, theta_i)) + for theta_i in theta] + else: + return np.sum(self._norm_coef * lib.calc_zn(self.order, r/self.radius, theta)) + \ No newline at end of file From 4a44e6a21d7375df0fd58a58e86d494ae2ce5daa Mon Sep 17 00:00:00 2001 From: rockfool Date: Thu, 30 Jan 2020 17:32:10 -0500 Subject: [PATCH 03/10] minor fix --- openmc/polynomial.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/openmc/polynomial.py b/openmc/polynomial.py index 79b914b256..b6acf58457 100644 --- a/openmc/polynomial.py +++ b/openmc/polynomial.py @@ -117,7 +117,7 @@ class Zernike(Polynomial): for m in range(-n,(n+1),2): j = int((n * (n + 2) + m) / 2) if m==0: - norm_vec[j] = n + 1.0 + norm_vec[j] = n + 1 else: norm_vec[j] = 2 * n + 2 norm_vec /= (np.pi * radius**2) @@ -140,4 +140,3 @@ class Zernike(Polynomial): for theta_i in theta] else: return np.sum(self._norm_coef * lib.calc_zn(self.order, r/self.radius, theta)) - \ No newline at end of file From df809d619170852cc9bc9a1eefc52900b9e10f91 Mon Sep 17 00:00:00 2001 From: rockfool Date: Thu, 30 Jan 2020 21:04:49 -0500 Subject: [PATCH 04/10] fix a typo --- openmc/polynomial.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/openmc/polynomial.py b/openmc/polynomial.py index b6acf58457..235800b408 100644 --- a/openmc/polynomial.py +++ b/openmc/polynomial.py @@ -94,7 +94,7 @@ class Zernike(Polynomial): r : float Position to be evaluated, normalized on radius [0,1] theta : float - Polar angular position to be evaluated, normalized on :math:`[0, 2*\pi]` + Azimuthal to be evaluated, normalized on :math:`[0, 2*\pi]` Attributes ---------- @@ -103,7 +103,7 @@ class Zernike(Polynomial): radius : float Domain of Zernike polynomials to be applied on. Default is 1. theta : float - Polar angle of Zernike polynomial to be applied on. Default is 0. + Azimuthal of Zernike polynomial to be applied on. Default is 0. norm_coef : iterable of float The list of coefficients of each term in the polynomials after normailization. From a740ae6d30b17fc6bdc1cf49e05a97201153e351 Mon Sep 17 00:00:00 2001 From: rockfool Date: Thu, 30 Jan 2020 21:08:07 -0500 Subject: [PATCH 05/10] modify radius in unittest for zernikeradial --- tests/unit_tests/test_polynomials.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/unit_tests/test_polynomials.py b/tests/unit_tests/test_polynomials.py index f6ac2711fd..dfcdb1f58a 100644 --- a/tests/unit_tests/test_polynomials.py +++ b/tests/unit_tests/test_polynomials.py @@ -24,7 +24,7 @@ def test_zernike_radial(): ref_vals = np.sum(norm_coeff * raw_zn) - test_vals = zn_rad(rho) + test_vals = zn_rad(rho * zn_rad.radius) assert ref_vals == test_vals @@ -39,6 +39,6 @@ def test_zernike_radial(): ref_vals = [np.sum(norm_coeff * raw_zn1), np.sum(norm_coeff * raw_zn2)] - test_vals = zn_rad(rho) + test_vals = zn_rad([i * zn_rad.radius for i in rho]) assert np.allclose(ref_vals, test_vals) From 93cb6120d16a53e36e8dcdf15ac7cebea4c41619 Mon Sep 17 00:00:00 2001 From: rockfool Date: Wed, 5 Feb 2020 23:25:50 -0500 Subject: [PATCH 06/10] add space next to math sign --- openmc/polynomial.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/openmc/polynomial.py b/openmc/polynomial.py index 235800b408..7f9f0e7e55 100644 --- a/openmc/polynomial.py +++ b/openmc/polynomial.py @@ -110,11 +110,11 @@ class Zernike(Polynomial): """ def __init__(self, coef, radius=1): super().__init__(coef) - self._order = int(((np.sqrt(8 * len(self.coef) + 1)) -3) / 2) + self._order = int(((np.sqrt(8 * len(self.coef) + 1)) - 3) / 2) self.radius = radius norm_vec = np.ones((len(self.coef))) for n in range(self._order + 1): - for m in range(-n,(n+1),2): + for m in range(-n, n + 1, 2): j = int((n * (n + 2) + m) / 2) if m==0: norm_vec[j] = n + 1 @@ -130,13 +130,13 @@ class Zernike(Polynomial): def __call__(self, r, theta=0.0): import openmc.lib as lib if isinstance(r, Iterable) and isinstance(theta, Iterable): - return [[np.sum(self._norm_coef * lib.calc_zn(self.order, r_i/self.radius, theta_i)) + return [[np.sum(self._norm_coef * lib.calc_zn(self.order, r_i / self.radius, theta_i)) for r_i in r] for theta_i in theta] elif isinstance(r, Iterable) and not isinstance(theta, Iterable): - return [np.sum(self._norm_coef * lib.calc_zn(self.order, r_i/self.radius, theta)) + return [np.sum(self._norm_coef * lib.calc_zn(self.order, r_i / self.radius, theta)) for r_i in r] elif not isinstance(r, Iterable) and isinstance(theta, Iterable): - return [np.sum(self._norm_coef * lib.calc_zn(self.order, r/self.radius, theta_i)) + return [np.sum(self._norm_coef * lib.calc_zn(self.order, r / self.radius, theta_i)) for theta_i in theta] else: - return np.sum(self._norm_coef * lib.calc_zn(self.order, r/self.radius, theta)) + return np.sum(self._norm_coef * lib.calc_zn(self.order, r / self.radius, theta)) From d2102397528fc3789d12bfabbf27d5c1385df88c Mon Sep 17 00:00:00 2001 From: rockfool Date: Fri, 21 Feb 2020 09:34:51 -0500 Subject: [PATCH 07/10] address comments from reviewer --- openmc/polynomial.py | 31 +++++++++++++++---------------- 1 file changed, 15 insertions(+), 16 deletions(-) diff --git a/openmc/polynomial.py b/openmc/polynomial.py index 7f9f0e7e55..24742b8e35 100644 --- a/openmc/polynomial.py +++ b/openmc/polynomial.py @@ -1,3 +1,4 @@ +import math import numpy as np import openmc from collections.abc import Iterable @@ -47,8 +48,6 @@ class ZernikeRadial(Polynomial): A list of coefficients of each term in radial only Zernike polynomials radius : float Domain of Zernike polynomials to be applied on. Default is 1. - r : float - Position to be evaluated, normalized on radius [0,1] Attributes ---------- @@ -65,7 +64,7 @@ class ZernikeRadial(Polynomial): super().__init__(coef) self._order = 2 * (len(self.coef) - 1) self.radius = radius - norm_vec = (2 * np.arange(len(self.coef)) + 1) / (np.pi * radius**2) + norm_vec = (2 * np.arange(len(self.coef)) + 1) / (math.pi * radius**2) self._norm_coef = norm_vec * self.coef @property @@ -82,19 +81,16 @@ class ZernikeRadial(Polynomial): class Zernike(Polynomial): - """Create Zernike polynomials given coefficients and domain. - The Zernike polynomials are defined as in - :class:`ZernikeFilter`. + r"""Create Zernike polynomials given coefficients and domain. + + The azimuthal Zernike polynomials are defined as in :class:`ZernikeFilter`. + Parameters ---------- coef : Iterable of float A list of coefficients of each term in radial only Zernike polynomials radius : float Domain of Zernike polynomials to be applied on. Default is 1. - r : float - Position to be evaluated, normalized on radius [0,1] - theta : float - Azimuthal to be evaluated, normalized on :math:`[0, 2*\pi]` Attributes ---------- @@ -110,17 +106,20 @@ class Zernike(Polynomial): """ def __init__(self, coef, radius=1): super().__init__(coef) - self._order = int(((np.sqrt(8 * len(self.coef) + 1)) - 3) / 2) + r"""Solve order from number of coefficients + N = (order + 1)(order + 2) / 2 + """ + self._order = int((math.sqrt(8 * len(self.coef) + 1) - 3) / 2) self.radius = radius - norm_vec = np.ones((len(self.coef))) + norm_vec = np.ones(len(self.coef)) for n in range(self._order + 1): for m in range(-n, n + 1, 2): - j = int((n * (n + 2) + m) / 2) - if m==0: + j = int((n*(n + 2) + m)/2) + if m == 0: norm_vec[j] = n + 1 else: - norm_vec[j] = 2 * n + 2 - norm_vec /= (np.pi * radius**2) + norm_vec[j] = 2*n + 2 + norm_vec /= (math.pi * radius**2) self._norm_coef = norm_vec * self.coef @property From 8d83b391b6481f0eadf482c4ddf60ee7bad42e03 Mon Sep 17 00:00:00 2001 From: rockfool Date: Mon, 24 Feb 2020 15:10:41 -0500 Subject: [PATCH 08/10] add unit test for azimuthal zernike --- tests/unit_tests/test_polynomials.py | 52 ++++++++++++++++++++++++++++ 1 file changed, 52 insertions(+) diff --git a/tests/unit_tests/test_polynomials.py b/tests/unit_tests/test_polynomials.py index dfcdb1f58a..bad8fafd98 100644 --- a/tests/unit_tests/test_polynomials.py +++ b/tests/unit_tests/test_polynomials.py @@ -42,3 +42,55 @@ def test_zernike_radial(): test_vals = zn_rad([i * zn_rad.radius for i in rho]) assert np.allclose(ref_vals, test_vals) + + +def test_zernike(): + import openmc.lib as lib + + coeff = np.asarray([1.1e-1, -3.2e2, 5.3, 7.4, -9.5, 0.005]) + zn_azimuthal = openmc.Zernike(coeff) + assert zn_azimuthal.order == 2 + assert zn_azimuthal.radius == 1 + + coeff = np.asarray([1.5, -3.6, 9.7e-1, -6.8e-1, 0.11, 0.33e2, 0.002, 13.75, \ + 3.1, -7.3, 7.8e-1, -1.1e-1, 2.56, 5.25e3, 0.123]) + zn_azimuthal = openmc.Zernike(coeff, 0.392) + assert zn_azimuthal.order == 4 + assert zn_azimuthal.radius == 0.392 + norm_vec = np.array([1, 4, 4, 6, 3, 6, 8, 8, 8, 8, + 10, 10, 5, 10, 10]) / (np.pi * 0.392 ** 2) + norm_coeff = norm_vec * coeff + + rho = 0.5 + + theta = np.radians(45) + # Reference solution from running the C API for calc_zn + raw_zn = lib.calc_zn(zn_azimuthal.order, rho, theta) + + ref_vals = np.sum(norm_coeff * raw_zn) + + test_vals = zn_azimuthal(rho * zn_azimuthal.radius, theta) + + assert ref_vals == test_vals + + rho = [0.2, 0.5] + + theta = np.radians(30) + #Reference solution from running the C API for calc_zn + raw_zn1 = lib.calc_zn(zn_azimuthal.order, rho[0], theta) + + raw_zn2 = lib.calc_zn(zn_azimuthal.order, rho[1], theta) + + ref_vals = [np.sum(norm_coeff * raw_zn1), np.sum(norm_coeff * raw_zn2)] + + test_vals = zn_azimuthal([i * zn_azimuthal.radius for i in rho], theta) + + assert np.allclose(ref_vals, test_vals) + + rho = 0.2 + + theta = np.radians([30, 60]) + #Reference solution from running the C API for calc_zn + raw_zn1 = lib.calc_zn(zn_azimuthal.order, rho, theta[0]) + + raw_zn2 = lib.calc_zn(zn_azimuthal.order, rho, theta[1]) From 5566f75feebda172c9c11ece28b86a2b5c254465 Mon Sep 17 00:00:00 2001 From: rockfool Date: Mon, 24 Feb 2020 15:21:00 -0500 Subject: [PATCH 09/10] bug fix in unit test --- tests/unit_tests/test_polynomials.py | 29 ++++++++++++++++++++++++++++ 1 file changed, 29 insertions(+) diff --git a/tests/unit_tests/test_polynomials.py b/tests/unit_tests/test_polynomials.py index bad8fafd98..414670a289 100644 --- a/tests/unit_tests/test_polynomials.py +++ b/tests/unit_tests/test_polynomials.py @@ -94,3 +94,32 @@ def test_zernike(): raw_zn1 = lib.calc_zn(zn_azimuthal.order, rho, theta[0]) raw_zn2 = lib.calc_zn(zn_azimuthal.order, rho, theta[1]) + + ref_vals = [np.sum(norm_coeff * raw_zn1), np.sum(norm_coeff * raw_zn2)] + + test_vals = zn_azimuthal(rho * zn_azimuthal.radius, [j for j in theta]) + + assert np.allclose(ref_vals, test_vals) + + rho = [0.2, 0.5] + + theta = np.radians([30, 60]) + #Reference solution from running the C API for calc_zn + raw_zn1 = lib.calc_zn(zn_azimuthal.order, rho[0], theta[0]) + + raw_zn2 = lib.calc_zn(zn_azimuthal.order, rho[1], theta[0]) + + raw_zn3 = lib.calc_zn(zn_azimuthal.order, rho[0], theta[1]) + + raw_zn4 = lib.calc_zn(zn_azimuthal.order, rho[1], theta[1]) + + ref_vals = [np.sum(norm_coeff * raw_zn1), np.sum(norm_coeff * raw_zn2), + np.sum(norm_coeff * raw_zn3), np.sum(norm_coeff * raw_zn4)] + + test_vals = zn_azimuthal([i * zn_azimuthal.radius for i in rho], [j for j in theta]) + + test_vals = np.ravel(test_vals) + + assert np.allclose(ref_vals, test_vals) + + From d409e99eb054466a8056896c9bb27980a7ebf96e Mon Sep 17 00:00:00 2001 From: rockfool Date: Tue, 17 Mar 2020 15:41:17 -0400 Subject: [PATCH 10/10] address comments from reviewers --- openmc/polynomial.py | 5 ++--- tests/unit_tests/test_polynomials.py | 2 +- 2 files changed, 3 insertions(+), 4 deletions(-) diff --git a/openmc/polynomial.py b/openmc/polynomial.py index 24742b8e35..82a6756c82 100644 --- a/openmc/polynomial.py +++ b/openmc/polynomial.py @@ -106,9 +106,8 @@ class Zernike(Polynomial): """ def __init__(self, coef, radius=1): super().__init__(coef) - r"""Solve order from number of coefficients - N = (order + 1)(order + 2) / 2 - """ + # Solve order from number of coefficients + # N = (order + 1)(order + 2) / 2 self._order = int((math.sqrt(8 * len(self.coef) + 1) - 3) / 2) self.radius = radius norm_vec = np.ones(len(self.coef)) diff --git a/tests/unit_tests/test_polynomials.py b/tests/unit_tests/test_polynomials.py index 414670a289..03f1fd5a46 100644 --- a/tests/unit_tests/test_polynomials.py +++ b/tests/unit_tests/test_polynomials.py @@ -52,7 +52,7 @@ def test_zernike(): assert zn_azimuthal.order == 2 assert zn_azimuthal.radius == 1 - coeff = np.asarray([1.5, -3.6, 9.7e-1, -6.8e-1, 0.11, 0.33e2, 0.002, 13.75, \ + coeff = np.asarray([1.5, -3.6, 9.7e-1, -6.8e-1, 0.11, 0.33e2, 0.002, 13.75, 3.1, -7.3, 7.8e-1, -1.1e-1, 2.56, 5.25e3, 0.123]) zn_azimuthal = openmc.Zernike(coeff, 0.392) assert zn_azimuthal.order == 4