crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-28 06:05:58 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Added in Zernike polynomials to C++ code

Browse source
This commit is contained in:
Adam G Nelson 2018-05-06 20:14:54 -04:00
parent feecf0134e
commit 468126f8c4
6 changed files with 211 additions and 161 deletions
Download patch file Download diff file Expand all files Collapse all files

8
openmc/capi/math.py
View file

@ -11,15 +11,15 @@ _dll.t_percentile.restype = c_double
_dll.t_percentile.argtypes = [POINTER(c_double), POINTER(c_int)]
_dll.calc_pn.restype = c_double
_dll.calc_pn.argtypes = [POINTER(c_int), POINTER(c_double)]
_dll.evaluate_legendre.restype = c_double
_dll.evaluate_legendre.argtypes = [POINTER(c_int), POINTER(c_double),
POINTER(c_double)]
_dll.calc_rn.restype = None
_dll.calc_rn.argtypes = [POINTER(c_int), ndpointer(c_double),
ndpointer(c_double)]
_dll.calc_zn.restype = None
_dll.calc_zn.argtypes = [POINTER(c_int), POINTER(c_double), POINTER(c_double),
ndpointer(c_double)]
_dll.evaluate_legendre.restype = c_double
_dll.evaluate_legendre.argtypes = [POINTER(c_int), POINTER(c_double),
POINTER(c_double)]
_dll.rotate_angle.restype = None
_dll.rotate_angle.argtypes = [ndpointer(c_double), POINTER(c_double),
ndpointer(c_double), POINTER(c_double)]
@ -169,7 +169,7 @@ def calc_zn(n, rho, phi):
"""
num_bins = ((n + 1) * (n + 2)) / 2
num_bins = ((n + 1) * (n + 2)) // 2
zn = np.zeros(num_bins, dtype=np.float64)
_dll.calc_zn(c_int(n), c_double(rho), c_double(phi), zn)
return zn

8
src/faddeeva/Faddeeva.h
View file

@ -1,5 +1,5 @@
/* Copyright (c) 2012 Massachusetts Institute of Technology
*
*
* Permission is hereby granted, free of charge, to any person obtaining
* a copy of this software and associated documentation files (the
* "Software"), to deal in the Software without restriction, including
@ -7,17 +7,17 @@
* distribute, sublicense, and/or sell copies of the Software, and to
* permit persons to whom the Software is furnished to do so, subject to
* the following conditions:
*
*
* The above copyright notice and this permission notice shall be
* included in all copies or substantial portions of the Software.
*
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
* MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
* LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
* OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
* WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
* WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/
/* Available at: http://ab-initio.mit.edu/Faddeeva

162
src/math.F90
View file

@ -43,14 +43,6 @@ module math
real(C_DOUBLE) :: pnx
end function calc_pn_cc
subroutine calc_rn_cc(n, uvw, rn) bind(C, name='calc_rn_c')
use ISO_C_BINDING
implicit none
integer(C_INT), value, intent(in) :: n
real(C_DOUBLE), intent(in) :: uvw(3)
real(C_DOUBLE), intent(in) :: rn(2 * n + 1)
end subroutine calc_rn_cc
pure function evaluate_legendre_cc(n, data, x) &
bind(C, name='evaluate_legendre_c') result(val)
use ISO_C_BINDING
@ -61,6 +53,23 @@ module math
real(C_DOUBLE) :: val
end function evaluate_legendre_cc
subroutine calc_rn_cc(n, uvw, rn) bind(C, name='calc_rn_c')
use ISO_C_BINDING
implicit none
integer(C_INT), value, intent(in) :: n
real(C_DOUBLE), intent(in) :: uvw(3)
real(C_DOUBLE), intent(out) :: rn(2 * n + 1)
end subroutine calc_rn_cc
subroutine calc_zn_cc(n, rho, phi, zn) bind(C, name='calc_zn_c')
use ISO_C_BINDING
implicit none
integer(C_INT), value, intent(in) :: n
real(C_DOUBLE), value, intent(in) :: rho
real(C_DOUBLE), value, intent(in) :: phi
real(C_DOUBLE), intent(out) :: zn(((n + 1) * (n + 2)) / 2)
end subroutine calc_zn_cc
subroutine rotate_angle_cc(uvw, mu, phi) bind(C, name='rotate_angle_c')
use ISO_C_BINDING
implicit none
@ -200,102 +209,13 @@ contains
!===============================================================================
subroutine calc_zn(n, rho, phi, zn) bind(C)
! This procedure uses the modified Kintner's method for calculating Zernike
! polynomials as outlined in Chong, C. W., Raveendran, P., & Mukundan,
! R. (2003). A comparative analysis of algorithms for fast computation of
! Zernike moments. Pattern Recognition, 36(3), 731-742.
integer(C_INT), intent(in) :: n ! Maximum order
real(C_DOUBLE), intent(in) :: rho ! Radial location in the unit disk
real(C_DOUBLE), intent(in) :: phi ! Theta (radians) location in the unit disk
real(C_DOUBLE), intent(out) :: zn(:) ! The resulting list of coefficients
! The resulting list of coefficients
real(C_DOUBLE), intent(out) :: zn(((n + 1) * (n + 2)) / 2)
real(C_DOUBLE) :: sin_phi, cos_phi ! Sine and Cosine of phi
real(C_DOUBLE) :: sin_phi_vec(n+1) ! Contains sin(n*phi)
real(C_DOUBLE) :: cos_phi_vec(n+1) ! Contains cos(n*phi)
real(C_DOUBLE) :: zn_mat(n+1, n+1) ! Matrix form of the coefficients which is
! easier to work with
real(C_DOUBLE) :: k1, k2, k3, k4 ! Variables for R_m_n calculation
integer(C_INT) :: i,p,q ! Loop counters
real(C_DOUBLE), parameter :: SQRT_N_1(0:10) = [&
sqrt(1.0_8), sqrt(2.0_8), sqrt(3.0_8), sqrt(4.0_8), &
sqrt(5.0_8), sqrt(6.0_8), sqrt(7.0_8), sqrt(8.0_8), &
sqrt(9.0_8), sqrt(10.0_8), sqrt(11.0_8)]
real(C_DOUBLE), parameter :: SQRT_2N_2(0:10) = SQRT_N_1*sqrt(2.0_8)
! n == radial degree
! m == azimuthal frequency
! ==========================================================================
! Determine vector of sin(n*phi) and cos(n*phi). This takes advantage of the
! following recurrence relations so that only a single sin/cos have to be
! evaluated (http://mathworld.wolfram.com/Multiple-AngleFormulas.html)
!
! sin(nx) = 2 cos(x) sin((n-1)x) - sin((n-2)x)
! cos(nx) = 2 cos(x) cos((n-1)x) - cos((n-2)x)
sin_phi = sin(phi)
cos_phi = cos(phi)
sin_phi_vec(1) = 1.0_8
cos_phi_vec(1) = 1.0_8
sin_phi_vec(2) = 2.0_8 * cos_phi
cos_phi_vec(2) = cos_phi
do i = 3, n+1
sin_phi_vec(i) = 2.0_8 * cos_phi * sin_phi_vec(i-1) - sin_phi_vec(i-2)
cos_phi_vec(i) = 2.0_8 * cos_phi * cos_phi_vec(i-1) - cos_phi_vec(i-2)
end do
do i = 1, n+1
sin_phi_vec(i) = sin_phi_vec(i) * sin_phi
end do
! ==========================================================================
! Calculate R_pq(rho)
! Fill the main diagonal first (Eq. 3.9 in Chong)
do p = 0, n
zn_mat(p+1, p+1) = rho**p
end do
! Fill in the second diagonal (Eq. 3.10 in Chong)
do q = 0, n-2
zn_mat(q+2+1, q+1) = (q+2) * zn_mat(q+2+1, q+2+1) - (q+1) * zn_mat(q+1, q+1)
end do
! Fill in the rest of the values using the original results (Eq. 3.8 in Chong)
do p = 4, n
k2 = 2 * p * (p - 1) * (p - 2)
do q = p-4, 0, -2
k1 = (p + q) * (p - q) * (p - 2) / 2
k3 = -q**2*(p - 1) - p * (p - 1) * (p - 2)
k4 = -p * (p + q - 2) * (p - q - 2) / 2
zn_mat(p+1, q+1) = ((k2 * rho**2 + k3) * zn_mat(p-2+1, q+1) + k4 * zn_mat(p-4+1, q+1)) / k1
end do
end do
! Roll into a single vector for easier computation later
! The vector is ordered (0,0), (1,-1), (1,1), (2,-2), (2,0),
! (2, 2), .... in (n,m) indices
! Note that the cos and sin vectors are offset by one
! sin_phi_vec = [sin(x), sin(2x), sin(3x) ...]
! cos_phi_vec = [1.0, cos(x), cos(2x)... ]
i = 1
do p = 0, n
do q = -p, p, 2
if (q < 0) then
zn(i) = zn_mat(p+1, abs(q)+1) * sin_phi_vec(abs(q)) * SQRT_2N_2(p)
else if (q == 0) then
zn(i) = zn_mat(p+1, q+1) * SQRT_N_1(p)
else
zn(i) = zn_mat(p+1, q+1) * cos_phi_vec(abs(q)+1) * SQRT_2N_2(p)
end if
i = i + 1
end do
end do
call calc_zn_cc(n, rho, phi, zn)
end subroutine calc_zn
!===============================================================================
@ -310,44 +230,11 @@ contains
real(C_DOUBLE), intent(out) :: uvw(3) ! rotated directional cosine
real(C_DOUBLE), optional :: phi ! azimuthal angle
real(C_DOUBLE) :: phi_ ! azimuthal angle
real(C_DOUBLE) :: sinphi ! sine of azimuthal angle
real(C_DOUBLE) :: cosphi ! cosine of azimuthal angle
real(C_DOUBLE) :: a ! sqrt(1 - mu^2)
real(C_DOUBLE) :: b ! sqrt(1 - w^2)
real(C_DOUBLE) :: u0 ! original cosine in x direction
real(C_DOUBLE) :: v0 ! original cosine in y direction
real(C_DOUBLE) :: w0 ! original cosine in z direction
! Copy original directional cosines
u0 = uvw0(1)
v0 = uvw0(2)
w0 = uvw0(3)
! Sample azimuthal angle in [0,2pi) if none provided
uvw = uvw0
if (present(phi)) then
phi_ = phi
call rotate_angle_cc(uvw, mu, phi)
else
phi_ = TWO * PI * prn()
end if
! Precompute factors to save flops
sinphi = sin(phi_)
cosphi = cos(phi_)
a = sqrt(max(ZERO, ONE - mu*mu))
b = sqrt(max(ZERO, ONE - w0*w0))
! Need to treat special case where sqrt(1 - w**2) is close to zero by
! expanding about the v component rather than the w component
if (b > 1e-10) then
uvw(1) = mu*u0 + a*(u0*w0*cosphi - v0*sinphi)/b
uvw(2) = mu*v0 + a*(v0*w0*cosphi + u0*sinphi)/b
uvw(3) = mu*w0 - a*b*cosphi
else
b = sqrt(ONE - v0*v0)
uvw(1) = mu*u0 + a*(u0*v0*cosphi + w0*sinphi)/b
uvw(2) = mu*v0 - a*b*cosphi
uvw(3) = mu*w0 + a*(v0*w0*cosphi - u0*sinphi)/b
call rotate_angle_cc(uvw, mu, -10._8)
end if
end subroutine rotate_angle
@ -492,6 +379,9 @@ contains
factors(i+3) = factors(i+1)*(E + (ONE + TWO * i) * half_inv_dopp2)
end if
end do
! call broaden_wmp_polynomials_cc(E, dopp, n, factors)
end subroutine broaden_wmp_polynomials
end module math

135
src/math_functions.cpp
View file

@ -169,9 +169,10 @@ double __attribute__ ((const)) calc_pn_c(int n, double x) {
double __attribute__ ((const)) evaluate_legendre_c(int n, double data[],
double x) {
double val;
int l;
val = 0.5 * data[0];
for (int l = 1; l < n; l++) {
for (l = 1; l < n; l++) {
val += (static_cast<double>(l) + 0.5) * data[l] * calc_pn_c(l, x);
}
return val;
@ -577,6 +578,122 @@ void calc_rn_c(int n, double uvw[3], double rn[]){
}
}
//==============================================================================
// CALC_ZN calculates the n-th order modified Zernike polynomial moment for a
// given angle (rho, theta) location in the unit disk. The normalization of the
// polynomials is such tha the integral of Z_pq*Z_pq over the unit disk is
// exactly pi
//==============================================================================
void calc_zn_c(int n, double rho, double phi, double zn[]) {
// This procedure uses the modified Kintner's method for calculating Zernike
// polynomials as outlined in Chong, C. W., Raveendran, P., & Mukundan,
// R. (2003). A comparative analysis of algorithms for fast computation of
// Zernike moments. Pattern Recognition, 36(3), 731-742.
double sin_phi; // Cosine of phi
double cos_phi; // Sine of phi
double sin_phi_vec[n + 1]; // Sin[n * phi]
double cos_phi_vec[n + 1]; // Cos[n * phi]
double zn_mat[n + 1][n + 1]; // Matrix forms of the coefficients which are
// easier to work with
// Variables for R_m_n calculation
double k1;
double k2;
double k3;
double k4;
// Loop counters
int i;
int p;
int q;
const double SQRT_N_1[11] =
{std::sqrt(1.), std::sqrt(2.), std::sqrt(3.), std::sqrt(4.),
std::sqrt(5.), std::sqrt(6.), std::sqrt(7.), std::sqrt(8.),
std::sqrt(9.), std::sqrt(10.), std::sqrt(11.)};
const double SQRT_2N_2[11] =
{std::sqrt(2.) * std::sqrt(1.), std::sqrt(2.) * std::sqrt(2.),
std::sqrt(2.) * std::sqrt(3.), std::sqrt(2.) * std::sqrt(4.),
std::sqrt(2.) * std::sqrt(5.), std::sqrt(2.) * std::sqrt(6.),
std::sqrt(2.) * std::sqrt(7.), std::sqrt(2.) * std::sqrt(8.),
std::sqrt(2.) * std::sqrt(9.), std::sqrt(2.) * std::sqrt(10.),
std::sqrt(2.) * std::sqrt(11.)};
// n == radial degree
// m == azimuthal frequency
// ===========================================================================
// Determine vector of sin(n*phi) and cos(n*phi). This takes advantage of the
// following recurrence relations so that only a single sin/cos have to be
// evaluated (http://mathworld.wolfram.com/Multiple-AngleFormulas.html)
//
// sin(nx) = 2 cos(x) sin((n-1)x) - sin((n-2)x)
// cos(nx) = 2 cos(x) cos((n-1)x) - cos((n-2)x)
sin_phi = std::sin(phi);
cos_phi = std::cos(phi);
sin_phi_vec[0] = 1.0;
cos_phi_vec[0] = 1.0;
sin_phi_vec[1] = 2.0 * cos_phi;
cos_phi_vec[1] = cos_phi;
for (i = 2; i <= n; i++) {
sin_phi_vec[i] = 2. * cos_phi * sin_phi_vec[i - 1] - sin_phi_vec[i - 2];
cos_phi_vec[i] = 2. * cos_phi * cos_phi_vec[i - 1] - cos_phi_vec[i - 2];
}
for (i = 0; i <= n; i++) {
sin_phi_vec[i] *= sin_phi;
}
// ===========================================================================
// Calculate R_pq(rho)
// Fill the main diagonal first (Eq 3.9 in Chong)
for (p = 0; p <= n; p++) {
zn_mat[p][p] = std::pow(rho, p);
}
// Fill the 2nd diagonal (Eq 3.10 in Chong)
for (q = 0; q <= n - 2; q++) {
zn_mat[q][q+2] = (q + 2) * zn_mat[q+2][q+2] - (q + 1) * zn_mat[q][q];
}
// Fill in the rest of the values using the original results (Eq. 3.8 in Chong)
for (p = 4; p <= n; p++) {
k2 = static_cast<double> (2 * p * (p - 1) * (p - 2));
for (q = p - 4; q >= 0; q -= 2) {
k1 = static_cast<double>((p + q) * (p - q) * (p - 2)) / 2.;
k3 = static_cast<double>(-q * q * (p - 1) - p * (p - 1) * (p - 2));
k4 = static_cast<double>(-p * (p + q - 2) * (p - q - 2)) / 2.;
zn_mat[q][p] =
((k2 * rho * rho + k3) * zn_mat[q][p-2] + k4 * zn_mat[q][p-4]) / k1;
}
}
// Roll into a single vector for easier computation later
// The vector is ordered (0,0), (1,-1), (1,1), (2,-2), (2,0),
// (2, 2), .... in (n,m) indices
// Note that the cos and sin vectors are offset by one
// sin_phi_vec = [sin(x), sin(2x), sin(3x) ...]
// cos_phi_vec = [1.0, cos(x), cos(2x)... ]
i = 0;
for (p = 0; p <= n; p++) {
for (q = -p; q <= p; q += 2) {
if (q < 0) {
zn[i] = zn_mat[std::abs(q)][p] * sin_phi_vec[std::abs(q) - 1] * SQRT_2N_2[p];
} else if (q == 0) {
zn[i] = zn_mat[q][p] * SQRT_N_1[p];
} else {
zn[i] = zn_mat[q][p] * cos_phi_vec[q] * SQRT_2N_2[p];
}
i++;
}
}
}
//==============================================================================
// ROTATE_ANGLE rotates direction std::cosines through a polar angle whose
@ -678,8 +795,8 @@ double watt_spectrum_c(double a, double b) {
// FADDEEVA the Faddeeva function, using Stephen Johnson's implementation
//==============================================================================
// std::complex<double> faddeeva_c(std::complex<double> z) {
// std::complex<double> wv; // The resultant w(z) value
// double complex __attribute__ ((const)) faddeeva_c(double complex z) {
// double complex wv; // The resultant w(z) value
// double relerr; // Target relative error in the inner loop of MIT Faddeeva
// // Technically, the value we want is given by the equation:
@ -698,19 +815,19 @@ double watt_spectrum_c(double a, double b) {
// // Note that faddeeva_w will interpret zero as machine epsilon
// relerr = 0.;
// if (z.imag() > 0.) {
// wv = Faddeeva::w(z, relerr);
// if (cimag(z) > 0.) {
// wv = Faddeeva_w(z, relerr);
// } else {
// wv = -std::conj(Faddeeva::w(std::conj(z), relerr));
// wv = -conj(Faddeeva_w(conj(z), relerr));
// }
// return wv;
// }
// std::complex<double> w_derivative_c(std::complex<double> z, int order){
// std::complex<double> wv; // The resultant w(z) value
// double complex w_derivative_c(double complex z, int order){
// double complex wv; // The resultant w(z) value
// const std::complex<double> twoi_sqrtpi(0.0, 2.0 / std::sqrt(PI));
// const double complex twoi_sqrtpi = 2.0 / std::sqrt(PI) * I;
// switch(order) {
// case 0:

20
src/math_functions.h
View file

@ -2,11 +2,11 @@
#define MATH_FUNCTIONS_H
#include <cmath>
#include <complex>
#include <complex.h>
#include <iostream>
#include "random_lcg.h"
// #include "faddeeva/Faddeeva.hh"
// #include "faddeeva/Faddeeva.h"
namespace openmc {
@ -51,6 +51,15 @@ extern "C" double evaluate_legendre_c(int n, double data[], double x)
extern "C" void calc_rn_c(int n, double uvw[3], double rn[]);
//==============================================================================
// CALC_ZN calculates the n-th order modified Zernike polynomial moment for a
// given angle (rho, theta) location in the unit disk. The normalization of the
// polynomials is such tha the integral of Z_pq*Z_pq over the unit disk is
// exactly pi
//==============================================================================
extern "C" void calc_zn_c(int n, double rho, double phi, double zn[]);
//==============================================================================
// ROTATE_ANGLE rotates direction cosines through a polar angle whose cosine is
// mu and through an azimuthal angle sampled uniformly. Note that this is done
@ -83,16 +92,17 @@ extern "C" double watt_spectrum_c(double a, double b);
// FADDEEVA the Faddeeva function, using Stephen Johnson's implementation
//==============================================================================
// extern "C" std::complex<double> faddeeva_c(std::complex<double> z);
// extern "C" double complex faddeeva_c(double complex z) __attribute__ ((const));
// extern "C" std::complex<double> w_derivative_c(std::complex<double> z, int order);
// extern "C" double complex w_derivative_c(double complex z, int order);
//==============================================================================
// BROADEN_WMP_POLYNOMIALS Doppler broadens the windowed multipole curvefit.
// The curvefit is a polynomial of the form a/E + b/sqrt(E) + c + d sqrt(E) ...
//==============================================================================
extern "C" void broaden_wmp_polynomials_c(double E, double dopp, int n, double factors[]);
extern "C" void broaden_wmp_polynomials_c(double E, double dopp, int n,
double factors[]);
} // namespace openmc
#endif // MATH_FUNCTIONS_H

39
tests/unit_tests/test_math.py
View file

@ -115,7 +115,38 @@ def test_calc_rn():
def test_calc_zn():
pass
n = 10
rho = 0.5
phi = 0.5
# Reference solution from running the Fortran implementation
ref_vals = np.array(
[1.00000000e+00, 4.79425539e-01, 8.77582562e-01,
5.15293637e-01, -8.66025404e-01, 3.30866239e-01,
3.52667735e-01, -8.47512624e-01, -1.55136145e+00,
2.50093775e-02, 1.79715684e-01, -1.33048245e+00,
-2.79508497e-01, -8.54292956e-01, -8.22482403e-02,
6.47865100e-02, -1.18780200e+00, 5.18993370e-01,
9.50010991e-01, -8.42327938e-02, -8.67263404e-02,
8.25035501e-03, -7.44248626e-01, 1.52505281e+00,
1.15751620e+00, 9.79225149e-01, 3.40611006e-01,
-5.78783240e-02, -1.09619759e-02, -3.17938327e-01,
1.90147482e+00, 2.84658914e-01, 5.21064646e-01,
1.34842791e-01, 4.25607546e-01, -2.92642715e-02,
-1.25423479e-02, -4.67751162e-02, 1.50696182e+00,
-6.13603897e-01, -8.67187500e-01, -3.93990531e-01,
-6.89672461e-01, 3.28139254e-01, -1.08327149e-02,
-8.53837419e-03, 7.04712042e-02, 7.73660979e-01,
-1.63799891e+00, -8.12396290e-01, -1.48708143e+00,
-1.16158437e-01, -1.03565982e+00, 1.88131088e-01,
-1.84122553e-03, -4.39233743e-03, 9.01295675e-02,
1.32511582e-01, -1.83260987e+00, -7.93994967e-01,
-2.97978009e-01, -5.09818305e-01, 8.38707753e-01,
-9.29602211e-01, 7.78441102e-02, 1.29931014e-03])
test_vals = openmc.capi.math.calc_zn(n, rho, phi)
assert np.allclose(ref_vals, test_vals)
def test_rotate_angle():
@ -128,7 +159,7 @@ def test_rotate_angle():
test_uvw = openmc.capi.math.rotate_angle(uvw0, mu, phi)
assert np.allclose(ref_uvw, test_uvw)
assert np.array_equal(ref_uvw, test_uvw)
# Repeat for mu = 1 (no change)
mu = 1.
@ -136,7 +167,7 @@ def test_rotate_angle():
test_uvw = openmc.capi.math.rotate_angle(uvw0, mu, phi)
assert np.allclose(ref_uvw, test_uvw)
assert np.array_equal(ref_uvw, test_uvw)
# Need to test phi=None somehow...
@ -180,3 +211,5 @@ def test_broaden_wmp_polynomials():
test_val = openmc.capi.math.broaden_wmp_polynomials(test_E, test_dopp, n)
assert np.allclose(ref_val, test_val)
test_calc_zn()