final preps for a PR

This commit is contained in:
Adam G Nelson 2018-05-13 19:21:36 -04:00
parent a6e75c3509
commit d7b20fb51e
4 changed files with 23 additions and 87 deletions

View file

@ -1,39 +1,38 @@
from ctypes import (c_int, c_double, POINTER, CFUNCTYPE)
from ctypes import (c_int, c_double, POINTER)
import numpy as np
from numpy.ctypeslib import ndpointer
from . import _dll
# _dll.normal_percentile.restype = c_double
# _dll.normal_percentile.argtypes = [POINTER(c_double)]
_dll.t_percentile.restype = c_double
_dll.t_percentile.argtypes = [POINTER(c_double), POINTER(c_int)]
_dll.calc_pn.restype = None
_dll.calc_pn.argtypes = [POINTER(c_int), 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.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.rotate_angle.restype = None
_dll.rotate_angle.argtypes = [ndpointer(c_double), POINTER(c_double),
ndpointer(c_double), POINTER(c_double)]
_dll.maxwell_spectrum.restype = c_double
_dll.maxwell_spectrum.argtypes = [POINTER(c_double)]
_dll.watt_spectrum.restype = c_double
_dll.watt_spectrum.argtypes = [POINTER(c_double), POINTER(c_double)]
# COMPLEX DATA TYPES?
# _dll.faddeeva.restype = c_double
# _dll.faddeeva.argtypes = [POINTER(c_double)]
# _dll.w_derivative.restype = c_double
# _dll.w_derivative.argtypes = [POINTER(c_double)]
_dll.broaden_wmp_polynomials.restype = None
_dll.broaden_wmp_polynomials.argtypes = [POINTER(c_double), POINTER(c_double),
POINTER(c_int), ndpointer(c_double)]

View file

@ -2,6 +2,16 @@
namespace openmc {
//==============================================================================
// Module constants.
//==============================================================================
// TODO: cmath::M_PI has 3 more digits precision than the Fortran constant we
// use so for now we will reuse the Fortran constant until we are OK with
// modifying test results
extern "C" const double PI {3.1415926535898};
extern "C" const double SQRT_PI {std::sqrt(PI)};
//==============================================================================
// NORMAL_PERCENTILE calculates the percentile of the standard normal
@ -748,59 +758,6 @@ double watt_spectrum_c(const double a, const double b) {
return E_out;
}
//==============================================================================
// FADDEEVA the Faddeeva function, using Stephen Johnson's implementation
//==============================================================================
// 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:
// // w(z) = I/Pi * Integrate[Exp[-t^2]/(z-t), {t, -Infinity, Infinity}]
// // as shown in Equation 63 from Hwang, R. N. "A rigorous pole
// // representation of multilevel cross sections and its practical
// // applications." Nuclear Science and Engineering 96.3 (1987): 192-209.
// //
// // The MIT Faddeeva function evaluates w(z) = exp(-z^2)erfc(-iz). These
// // two forms of the Faddeeva function are related by a transformation.
// //
// // If we call the integral form w_int, and the function form w_fun:
// // For imag(z) > 0, w_int(z) = w_fun(z)
// // For imag(z) < 0, w_int(z) = -conjg(w_fun(conjg(z)))
// // Note that faddeeva_w will interpret zero as machine epsilon
// relerr = 0.;
// if (cimag(z) > 0.) {
// wv = Faddeeva_w(z, relerr);
// } else {
// wv = -conj(Faddeeva_w(conj(z), relerr));
// }
// return wv;
// }
// double complex w_derivative_c(double complex z, int order){
// double complex wv; // The resultant w(z) value
// const double complex twoi_sqrtpi = 2.0 / SQRT_PI * I;
// switch(order) {
// case 0:
// wv = faddeeva_c(z);
// break;
// case 1:
// wv = -2. * z * faddeeva_c(z) + twoi_sqrtpi;
// break;
// default:
// wv = -2. * z * w_derivative_c(z, order - 1) - 2. * (order - 1) *
// w_derivative_c(z, order - 2);
// }
// return wv;
// }
//==============================================================================
// 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) ...

View file

@ -2,21 +2,22 @@
#define MATH_FUNCTIONS_H
#include <cmath>
#include <complex.h>
#include <iostream>
#include "random_lcg.h"
// #include "faddeeva/Faddeeva.h"
namespace openmc {
//==============================================================================
// Module constants.
//==============================================================================
// TODO: cmath::M_PI has 3 more digits precision than the Fortran constant we
// use so for now we will reuse the Fortran constant until we are OK with
// modifying test results
const double PI = 3.1415926535898;
extern "C" const double PI;
const double SQRT_PI = std::sqrt(PI);
extern "C" const double SQRT_PI;
//==============================================================================
// NORMAL_PERCENTILE calculates the percentile of the standard normal
@ -94,14 +95,6 @@ extern "C" double maxwell_spectrum_c(const double T);
extern "C" double watt_spectrum_c(const double a, const double b);
//==============================================================================
// FADDEEVA the Faddeeva function, using Stephen Johnson's implementation
//==============================================================================
// extern "C" double complex faddeeva_c(double complex z) __attribute__ ((const));
// 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) ...

View file

@ -5,19 +5,6 @@ import openmc
import openmc.capi
# def test_normal_percentile():
# # normal_percentile has three branches to consider:
# # p < 0.02425; 0.02425 <= p <= 0.97575; and p > 0.97575
# test_ps = [0.02, 0.4, 0.5, 0.6, 0.98]
# # The reference solutions come from Scipy
# ref_zs = [sp.stats.norm.ppf(p) for p in test_ps]
# test_zs = [openmc.capi.math.normal_percentile(p) for p in test_ps]
# assert np.allclose(ref_zs, test_zs)
def test_t_percentile():
# Permutations include 1 DoF, 2 DoF, and > 2 DoF
# We will test 5 p-values at 3-DoF values