From 35a571258d05a3142893a89064624603ad0e6b7e Mon Sep 17 00:00:00 2001 From: TaigoFr Date: Sat, 19 Jun 2021 00:51:37 +0100 Subject: [PATCH] Add Cubic Spline related classes Cubic Spline is an interpolator that is quick and useful for 1D data. The Lagrange class could be used for interpolating 1D data, but Cubic Spline allows easily for non-uniform spacing 1D grid. This then can be used to interpolate data that exist in params.txt via the FileInterpolator1D. For instance, I have used this to solve 1D constraints in spherical symmetry elsewhere and then importing the data to the 3D grid. It requires a matrix inversion, and that is done by a LU solver with the TriDiagonalMatrix class. Two tests added. CubicSplineInterpolatorTest, to test the cubic spline interpolation, and FileInterpolator1DTest, to test the importing of file data to interpolate. The later can be also useful as an example of how to use these classes. If "N", "x" and "y" are set in a separate file 'fileWithData.txt', one easy way to import them to the 'params.txt' without changing the file is running with the command: 'Main.ex params.txt FILE=fileWithData.txt' --- .../CubicSplineInterpolator.cpp | 181 ++++++++++++++++++ .../CubicSplineInterpolator.hpp | 60 ++++++ .../FileInterpolator1D.hpp | 85 ++++++++ .../TriDiagonalMatrix.cpp | 81 ++++++++ .../TriDiagonalMatrix.hpp | 43 +++++ .../CubicSplineInterpolatorTest.cpp | 75 ++++++++ .../CubicSplineInterpolatorTest.inputs | 8 + Tests/CubicSplineInterpolatorTest/GNUmakefile | 21 ++ .../SimulationParameters.hpp | 19 ++ .../UserVariables.hpp | 26 +++ .../FileInterpolator1DTest.cpp | 79 ++++++++ .../FileInterpolator1DTest.inputs | 4 + Tests/FileInterpolator1DTest/GNUmakefile | 21 ++ .../SimulationParameters.hpp | 19 ++ .../FileInterpolator1DTest/UserVariables.hpp | 26 +++ 15 files changed, 748 insertions(+) create mode 100644 Source/utils/CubicSpline1DInterpolation/CubicSplineInterpolator.cpp create mode 100644 Source/utils/CubicSpline1DInterpolation/CubicSplineInterpolator.hpp create mode 100644 Source/utils/CubicSpline1DInterpolation/FileInterpolator1D.hpp create mode 100644 Source/utils/CubicSpline1DInterpolation/TriDiagonalMatrix.cpp create mode 100644 Source/utils/CubicSpline1DInterpolation/TriDiagonalMatrix.hpp create mode 100644 Tests/CubicSplineInterpolatorTest/CubicSplineInterpolatorTest.cpp create mode 100644 Tests/CubicSplineInterpolatorTest/CubicSplineInterpolatorTest.inputs create mode 100644 Tests/CubicSplineInterpolatorTest/GNUmakefile create mode 100644 Tests/CubicSplineInterpolatorTest/SimulationParameters.hpp create mode 100644 Tests/CubicSplineInterpolatorTest/UserVariables.hpp create mode 100644 Tests/FileInterpolator1DTest/FileInterpolator1DTest.cpp create mode 100644 Tests/FileInterpolator1DTest/FileInterpolator1DTest.inputs create mode 100644 Tests/FileInterpolator1DTest/GNUmakefile create mode 100644 Tests/FileInterpolator1DTest/SimulationParameters.hpp create mode 100644 Tests/FileInterpolator1DTest/UserVariables.hpp diff --git a/Source/utils/CubicSpline1DInterpolation/CubicSplineInterpolator.cpp b/Source/utils/CubicSpline1DInterpolation/CubicSplineInterpolator.cpp new file mode 100644 index 000000000..c0fe74b74 --- /dev/null +++ b/Source/utils/CubicSpline1DInterpolation/CubicSplineInterpolator.cpp @@ -0,0 +1,181 @@ +/* GRChombo + * Copyright 2012 The GRChombo collaboration. + * Please refer to LICENSE in GRChombo's root directory. + */ + +// Chombo includes +#include "parstream.H" //Gives us pout() + +#include "CubicSplineInterpolator.hpp" +#include "TriDiagonalMatrix.hpp" + +#include + +// Chombo namespace +#include "UsingNamespace.H" + +CubicSplineInterpolator::CubicSplineInterpolator(const std::vector &a_x, + const std::vector &a_y) + : solved(false), points_set(false), m_bc_left(second_deriv), + m_bc_right(second_deriv), m_bc_left_value(0.), m_bc_right_value(0.) +{ + allow_extrapolation(false); + set_points(a_x, a_y); +} + +void CubicSplineInterpolator::allow_extrapolation(bool allow) +{ + m_allow_extrapolation = allow; +} + +void CubicSplineInterpolator::set_points(const std::vector &a_x, + const std::vector &a_y, + bool do_solve) +{ + if (a_x.size() > 0) + { + assert(a_x.size() == a_y.size()); + m_x = a_x; + m_y = a_y; + points_set = true; + if (do_solve) + solve(); + } +} + +void CubicSplineInterpolator::set_boundary_conditions(BoundaryType a_left, + double a_left_value, + BoundaryType a_right, + double a_right_value, + bool solve_if_points_set) +{ + m_bc_left = a_left; + m_bc_right = a_right; + m_bc_left_value = a_left_value; + m_bc_right_value = a_right_value; + + if (points_set && solve_if_points_set) + solve(); +} + +void CubicSplineInterpolator::solve() +{ + unsigned n = size(); + + TriDiagonalMatrix mat(n); + std::vector rhs(n); + + for (unsigned i = 1; i < n - 1; i++) + { + mat.upper(i) = m_x[i] - m_x[i + 1]; + mat.diag(i) = 2. * (m_x[i - 1] - m_x[i + 1]); + mat.lower(i - 1) = m_x[i - 1] - m_x[i]; + rhs[i] = 6. * ((m_y[i - 1] - m_y[i]) / (m_x[i - 1] - m_x[i]) - + (m_y[i] - m_y[i + 1]) / (m_x[i] - m_x[i + 1])); + } + + // left boundary conditions + if (m_bc_left == second_deriv) + { + // K[0] = f'' + mat.diag(0) = 1.; + mat.upper(0) = 0.; + rhs[0] = m_bc_left_value; + } + else if (m_bc_left == first_deriv) + { + // f', needs to be re-expressed in terms of K: + // (2K[0]+K[1])(x[0]-x[1]) = 6 ( f' - (y[0]-y[1])/(x[0]-x[1]) ) + mat.diag(0) = 2. * (m_x[0] - m_x[1]); + mat.upper(0) = (m_x[0] - m_x[1]); + rhs[0] = 6. * (m_bc_left_value - (m_y[0] - m_y[1]) / (m_x[0] - m_x[1])); + } + else + { + assert(false); + } + + // right boundary conditions + if (m_bc_right == second_deriv) + { + // K[n-1] = f'' + mat.diag(n - 1) = 1.; + mat.lower(n - 2) = 0.; + rhs[n - 1] = m_bc_right_value; + } + else if (m_bc_right == first_deriv) + { + // f', needs to be re-expressed in terms of K: + // (K[n-2]+2K[n-1])(x[n-2]-x[n-1]) = + // = 6 ( (y[n-2]-y[n-1])/(x[n-2]-x[n-1]) - f' ) + mat.diag(n - 1) = 2. * (m_x[n - 2] - m_x[n - 1]); + mat.lower(n - 2) = (m_x[n - 2] - m_x[n - 1]); + rhs[0] = 6. * ((m_y[n - 2] - m_y[n - 1]) / (m_x[n - 2] - m_x[n - 1]) - + m_bc_right_value); + } + else + { + assert(false); + } + + m_K = mat.lu_solve(rhs); + solved = true; +} + +double CubicSplineInterpolator::interpolate(double x, unsigned derivative) const +{ + assert(solved); + + unsigned n = size(); + + unsigned i; + for (i = 0; i < n; ++i) + if (x < m_x[i]) + break; + + if (i == 0 || i == n) + { + std::string msg = + "CubicSplineInterpolator::interpolate - extrapolating at " + + std::to_string(x) + " outside of valid range [" + + std::to_string(m_x[0]) + "," + std::to_string(m_x[n - 1]) + "]"; + pout() << msg << std::endl; + assert(m_allow_extrapolation); + + if (i == 0) // use segment 0 + ++i; + else // use segment n-1 + --i; + } + + double d_left = (x - m_x[i - 1]); + double d_right = (x - m_x[i]); + double delta = (m_x[i - 1] - m_x[i]); + + double answer; + switch (derivative) + { + case 0: + answer = + m_K[i - 1] / 6. * + (d_right * d_right * d_right / delta - d_right * delta) - + m_K[i] / 6. * (d_left * d_left * d_left / delta - d_left * delta) + + (m_y[i - 1] * d_right - m_y[i] * d_left) / delta; + break; + case 1: + answer = m_K[i - 1] / 6. * (3. * d_right * d_right / delta - delta) - + m_K[i] / 6. * (3 * d_left * d_left / delta - delta) + + (m_y[i - 1] - m_y[i]) / delta; + break; + case 2: + answer = (m_K[i - 1] * d_right - m_K[i] * d_left) / delta; + break; + case 3: + answer = (m_K[i - 1] - m_K[i]) / delta; + break; + default: // all others + answer = 0.; + } + + return answer; +} diff --git a/Source/utils/CubicSpline1DInterpolation/CubicSplineInterpolator.hpp b/Source/utils/CubicSpline1DInterpolation/CubicSplineInterpolator.hpp new file mode 100644 index 000000000..c84be8849 --- /dev/null +++ b/Source/utils/CubicSpline1DInterpolation/CubicSplineInterpolator.hpp @@ -0,0 +1,60 @@ +/* GRChombo + * Copyright 2012 The GRChombo collaboration. + * Please refer to LICENSE in GRChombo's root directory. + */ + +#ifndef CUBICSPLINEINTERPOLATOR_HPP_ +#define CUBICSPLINEINTERPOLATOR_HPP_ + +#include +#include + +//! Calculate 'x's, 'y's and 'K's in Mathematica, import them to here +//! and this class interpolates with cubic spline +class CubicSplineInterpolator +{ + public: + enum BoundaryType + { + first_deriv = 1, + second_deriv + }; + + CubicSplineInterpolator( + const std::vector &a_x = std::vector(), + const std::vector &a_y = std::vector()); + + void allow_extrapolation(bool allow); + void set_points(const std::vector &a_x, + const std::vector &a_y, bool do_solve = true); + // optional, default is 2nd derivative = 0 + void set_boundary_conditions(BoundaryType a_left, double a_left_value, + BoundaryType a_right, double a_right_value, + bool solve_if_points_set = true); + + double interpolate(double x, unsigned derivative = 0) const; + + void solve(); + + inline unsigned size() const { return m_x.size(); } + inline std::vector &x() { return m_x; } + inline std::vector x() const { return m_x; } + inline std::vector &y() { return m_y; } + inline std::vector y() const { return m_y; } + + // return 2nd derivatives + inline std::vector &K() { return m_K; } + inline std::vector K() const { return m_K; } + + private: + std::vector m_x, m_y; + bool m_allow_extrapolation; + + bool solved, points_set; + BoundaryType m_bc_left, m_bc_right; + double m_bc_left_value, m_bc_right_value; + + std::vector m_K; // storing 2nd derivatives on each point +}; + +#endif /* CUBICSPLINEINTERPOLATOR_HPP_ */ diff --git a/Source/utils/CubicSpline1DInterpolation/FileInterpolator1D.hpp b/Source/utils/CubicSpline1DInterpolation/FileInterpolator1D.hpp new file mode 100644 index 000000000..85ce1d16a --- /dev/null +++ b/Source/utils/CubicSpline1DInterpolation/FileInterpolator1D.hpp @@ -0,0 +1,85 @@ +/* GRChombo + * Copyright 2012 The GRChombo collaboration. + * Please refer to LICENSE in GRChombo's root directory. + */ + +#ifndef FILEINTERPOLATOR1D_HPP_ +#define FILEINTERPOLATOR1D_HPP_ + +#include "CubicSplineInterpolator.hpp" +#include "GRParmParse.hpp" + +#include +#include + +//! Calculate 'x's, 'y's elsewhere, import them to here +//! and this class interpolates with 4th order Lagrange +//! (left as a separate class to keep 'CubicSplineInterpolator' free from any +//! Chombo/GRChombo code) +class FileInterpolator1D +{ + public: + FileInterpolator1D(const std::string &N_label, const std::string &x_label, + const std::string &y_label, + double y_default_if_missing_file = 0.) + { + GRParmParse pp; + + int n = 0; + std::vector xs, ys; + if (pp.contains(N_label.c_str())) + { + pp.load(N_label.c_str(), n, 0); + pp.load(x_label.c_str(), xs, n, {0.}); + pp.load(y_label.c_str(), ys, n, {0.}); + } + + assert(n >= 0); + if (n >= 2) + spline.set_points(xs, ys); + else if (n == 1) + { + spline.allow_extrapolation(true); + // *2 such that it avoid numerical erros if 'x' is very big (+1 + // would do little to x=1e20) +1 to still work for x=0 + spline.set_points({xs[0], xs[0] * 2. + 1.}, {ys[0], ys[0]}); + } + else + { + pout() + << "Parameter: " << N_label + << " not found in parameter file. FileInterpolator1D will set " + << y_label + << " to its default value = " << y_default_if_missing_file + << " for all " << x_label << std::endl; + + spline.allow_extrapolation(true); + spline.set_points({-1.e5, 1.e5}, {y_default_if_missing_file, + y_default_if_missing_file}); + } + } + + inline void allow_extrapolation(bool allow) + { + if (spline.size() >= 2) + spline.allow_extrapolation(allow); + } + + // optional, default is 2nd derivative = 0 + inline void set_boundary_conditions( + CubicSplineInterpolator::BoundaryType a_left, double a_left_value, + CubicSplineInterpolator::BoundaryType a_right, double a_right_value) + { + spline.set_boundary_conditions(a_left, a_left_value, a_right, + a_right_value); + } + inline double interpolate(double x, int derivative = 0) + { + return spline.interpolate(x, derivative); + } + + private: + CubicSplineInterpolator spline; +}; + +#endif /* FILEINTERPOLATOR1D_HPP_ */ diff --git a/Source/utils/CubicSpline1DInterpolation/TriDiagonalMatrix.cpp b/Source/utils/CubicSpline1DInterpolation/TriDiagonalMatrix.cpp new file mode 100644 index 000000000..3233840f7 --- /dev/null +++ b/Source/utils/CubicSpline1DInterpolation/TriDiagonalMatrix.cpp @@ -0,0 +1,81 @@ +/* GRChombo + * Copyright 2012 The GRChombo collaboration. + * Please refer to LICENSE in GRChombo's root directory. + */ + +#include "TriDiagonalMatrix.hpp" +#include + +void TriDiagonalMatrix::resize(unsigned a_size) +{ + assert(a_size > 0); + m_upper.resize(a_size - 1); + m_diag.resize(a_size); + m_lower.resize(a_size - 1); +} + +double &TriDiagonalMatrix::upper(unsigned i) +{ + assert(i < size() - 1); + return m_upper[i]; +} +double TriDiagonalMatrix::upper(unsigned i) const +{ + assert(i < size() - 1); + return m_upper[i]; +} + +double &TriDiagonalMatrix::diag(unsigned i) +{ + assert(i < size()); + return m_diag[i]; +} +double TriDiagonalMatrix::diag(unsigned i) const +{ + assert(i < size()); + return m_diag[i]; +} + +double &TriDiagonalMatrix::lower(unsigned i) +{ + assert(i < size() - 1); + return m_lower[i]; +} +double TriDiagonalMatrix::lower(unsigned i) const +{ + assert(i < size() - 1); + return m_lower[i]; +} + +std::vector TriDiagonalMatrix::lu_solve(const std::vector &b) +{ + lu_decomposition(); + return back_sub(forward_sub(b)); +} + +void TriDiagonalMatrix::lu_decomposition() +{ + for (unsigned i = 1; i < size(); i++) + { + lower(i - 1) /= diag(i - 1); + diag(i) -= lower(i - 1) * upper(i - 1); + } +} + +std::vector TriDiagonalMatrix::back_sub(const std::vector &b) +{ + assert(size() == b.size()); + std::vector out(b); + out[size() - 1] /= diag(size() - 1); + for (int i = size() - 2; i >= 0; i--) + out[i] = (out[i] - upper(i) * out[i + 1]) / diag(i); + return out; +} +std::vector TriDiagonalMatrix::forward_sub(const std::vector &d) +{ + assert(size() == d.size()); + std::vector out(d); + for (unsigned i = 1; i < size(); i++) + out[i] -= lower(i - 1) * out[i - 1]; + return out; +} diff --git a/Source/utils/CubicSpline1DInterpolation/TriDiagonalMatrix.hpp b/Source/utils/CubicSpline1DInterpolation/TriDiagonalMatrix.hpp new file mode 100644 index 000000000..42d1180e6 --- /dev/null +++ b/Source/utils/CubicSpline1DInterpolation/TriDiagonalMatrix.hpp @@ -0,0 +1,43 @@ +/* GRChombo + * Copyright 2012 The GRChombo collaboration. + * Please refer to LICENSE in GRChombo's root directory. + */ + +#ifndef TRIDIAGONALMATRIX_HPP_ +#define TRIDIAGONALMATRIX_HPP_ + +#include + +class TriDiagonalMatrix +{ + public: + TriDiagonalMatrix() {} // default constructor + TriDiagonalMatrix(unsigned a_size) { resize(a_size); } + + inline unsigned size() const { return m_diag.size(); } + + void resize(unsigned a_size); + + double &upper(unsigned i); + double upper(unsigned i) const; + + double &diag(unsigned i); + double diag(unsigned i) const; + + double &lower(unsigned i); + double lower(unsigned i) const; + + std::vector lu_solve(const std::vector &b); + + private: + void lu_decomposition(); + std::vector back_sub(const std::vector &); + std::vector forward_sub(const std::vector &); + + private: + std::vector m_upper; // upper band + std::vector m_diag; // diagonal + std::vector m_lower; // lower band +}; + +#endif /* TRIDIAGONALMATRIX_HPP_ */ diff --git a/Tests/CubicSplineInterpolatorTest/CubicSplineInterpolatorTest.cpp b/Tests/CubicSplineInterpolatorTest/CubicSplineInterpolatorTest.cpp new file mode 100644 index 000000000..e2136bbc6 --- /dev/null +++ b/Tests/CubicSplineInterpolatorTest/CubicSplineInterpolatorTest.cpp @@ -0,0 +1,75 @@ +/* GRChombo + * Copyright 2012 The GRChombo collaboration. + * Please refer to LICENSE in GRChombo's root directory. + */ + +#ifdef CH_LANG_CC +/* + * _______ __ + * / ___/ / ___ __ _ / / ___ + * / /__/ _ \/ _ \/ V \/ _ \/ _ \ + * \___/_//_/\___/_/_/_/_.__/\___/ + * Please refer to LICENSE, in Chombo's root directory. + */ +#endif + +// Chombo includes +#include "parstream.H" //Gives us pout() + +#include "FileInterpolator1D.hpp" +#include "SetupFunctions.hpp" +#include + +// Chombo namespace +#include "UsingNamespace.H" + +int runCubicSplineInterpolatorTest(int argc, char *argv[]) +{ + // Load the parameter file and construct the SimulationParameter class + std::string in_string = argv[argc - 1]; + pout() << in_string << std::endl; + char const *in_file = argv[argc - 1]; + GRParmParse pp(0, argv + argc, NULL, in_file); + + FileInterpolator1D file("N", "x", "y"); + file.allow_extrapolation(false); + + // using sin(2 * pi * x) as a test function + // notice it already satisfies the default of: f''(0) = 0 = f''(1) + // file.set_boundary_conditions(CubicSplineInterpolator::second_deriv, 0., + // CubicSplineInterpolator::second_deriv, 0.); + + double point = 1. / M_PI; // some irrational number + double exact = sin(2.); + double result = file.interpolate(point); + + pout() << result << std::endl; + pout() << exact << std::endl; + + double err = abs(result / exact - 1.) * 100.; + + pout() << std::setprecision(9); + pout() << "Exact value is: " << exact << std::endl; + pout() << "Cubic Spline Interpolator gave: " << result << " (" << err + << "% error)" << std::endl; + + bool wrong = (err > 1e-3); + + return wrong; +} + +int main(int argc, char *argv[]) +{ + mainSetup(argc, argv); + + int status = runCubicSplineInterpolatorTest(argc, argv); + + if (status == 0) + std::cout << "CubicSplineInterpolator test passed." << endl; + else + std::cout << "CubicSplineInterpolator test failed with return code " + << status << endl; + + mainFinalize(); + return status; +} diff --git a/Tests/CubicSplineInterpolatorTest/CubicSplineInterpolatorTest.inputs b/Tests/CubicSplineInterpolatorTest/CubicSplineInterpolatorTest.inputs new file mode 100644 index 000000000..c36cb9255 --- /dev/null +++ b/Tests/CubicSplineInterpolatorTest/CubicSplineInterpolatorTest.inputs @@ -0,0 +1,8 @@ + +# points from function sin(2*pi*x) in the interval [0, 1] + +N = 100 + +x = 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 0.11 0.12 0.13 0.14 0.15 0.16 0.17 0.18 0.19 0.2 0.21 0.22 0.23 0.24 0.25 0.26 0.27 0.28 0.29 0.3 0.31 0.32 0.33 0.34 0.35 0.36 0.37 0.38 0.39 0.4 0.41 0.42 0.43 0.44 0.45 0.46 0.47 0.48 0.49 0.5 0.51 0.52 0.53 0.54 0.55 0.56 0.57 0.58 0.59 0.6 0.61 0.62 0.63 0.64 0.65 0.66 0.67 0.68 0.69 0.7 0.71 0.72 0.73 0.74 0.75 0.76 0.77 0.78 0.79 0.8 0.81 0.82 0.83 0.84 0.85 0.86 0.87 0.88 0.89 0.9 0.91 0.92 0.93 0.94 0.95 0.96 0.97 0.98 0.99 1. + +y = 0 0.0627905 0.125333 0.187381 0.24869 0.309017 0.368125 0.425779 0.481754 0.535827 0.587785 0.637424 0.684547 0.728969 0.770513 0.809017 0.844328 0.876307 0.904827 0.929776 0.951057 0.968583 0.982287 0.992115 0.998027 1. 0.998027 0.992115 0.982287 0.968583 0.951057 0.929776 0.904827 0.876307 0.844328 0.809017 0.770513 0.728969 0.684547 0.637424 0.587785 0.535827 0.481754 0.425779 0.368125 0.309017 0.24869 0.187381 0.125333 0.0627905 0 -0.0627905 -0.125333 -0.187381 -0.24869 -0.309017 -0.368125 -0.425779 -0.481754 -0.535827 -0.587785 -0.637424 -0.684547 -0.728969 -0.770513 -0.809017 -0.844328 -0.876307 -0.904827 -0.929776 -0.951057 -0.968583 -0.982287 -0.992115 -0.998027 -1. -0.998027 -0.992115 -0.982287 -0.968583 -0.951057 -0.929776 -0.904827 -0.876307 -0.844328 -0.809017 -0.770513 -0.728969 -0.684547 -0.637424 -0.587785 -0.535827 -0.481754 -0.425779 -0.368125 -0.309017 -0.24869 -0.187381 -0.125333 -0.0627905 0 \ No newline at end of file diff --git a/Tests/CubicSplineInterpolatorTest/GNUmakefile b/Tests/CubicSplineInterpolatorTest/GNUmakefile new file mode 100644 index 000000000..9f27d3e60 --- /dev/null +++ b/Tests/CubicSplineInterpolatorTest/GNUmakefile @@ -0,0 +1,21 @@ +# -*- Mode: Makefile -*- + +### This makefile produces an executable for each name in the `ebase' +### variable using the libraries named in the `LibNames' variable. + +# Included makefiles need an absolute path to the Chombo installation +# CHOMBO_HOME := Please set the CHOMBO_HOME locally (e.g. export CHOMBO_HOME=... in bash) + +GRCHOMBO_SOURCE = $(shell pwd)/../../Source + +ebase := CubicSplineInterpolatorTest + +LibNames := AMRTimeDependent AMRTools BoxTools + +src_dirs := $(GRCHOMBO_SOURCE)/GRChomboCore \ + $(GRCHOMBO_SOURCE)/utils \ + $(GRCHOMBO_SOURCE)/simd \ + $(GRCHOMBO_SOURCE)/AMRInterpolator \ + $(GRCHOMBO_SOURCE)/utils/CubicSpline1DInterpolation + +include $(CHOMBO_HOME)/mk/Make.test diff --git a/Tests/CubicSplineInterpolatorTest/SimulationParameters.hpp b/Tests/CubicSplineInterpolatorTest/SimulationParameters.hpp new file mode 100644 index 000000000..5e00a7de6 --- /dev/null +++ b/Tests/CubicSplineInterpolatorTest/SimulationParameters.hpp @@ -0,0 +1,19 @@ +/* GRChombo + * Copyright 2012 The GRChombo collaboration. + * Please refer to LICENSE in GRChombo's root directory. + */ + +#ifndef SIMULATIONPARAMETERS_HPP_ +#define SIMULATIONPARAMETERS_HPP_ + +// General includes +#include "ChomboParameters.hpp" +#include "GRParmParse.hpp" + +class SimulationParameters : public ChomboParameters +{ + public: + SimulationParameters(GRParmParse &pp) : ChomboParameters(pp) {} +}; + +#endif /* SIMULATIONPARAMETERS_HPP_ */ diff --git a/Tests/CubicSplineInterpolatorTest/UserVariables.hpp b/Tests/CubicSplineInterpolatorTest/UserVariables.hpp new file mode 100644 index 000000000..f7642667d --- /dev/null +++ b/Tests/CubicSplineInterpolatorTest/UserVariables.hpp @@ -0,0 +1,26 @@ +/* GRChombo + * Copyright 2012 The GRChombo collaboration. + * Please refer to LICENSE in GRChombo's root directory. + */ + +#ifndef USERVARIABLES_HPP +#define USERVARIABLES_HPP + +#include "EmptyDiagnosticVariables.hpp" +#include +#include + +// assign enum to each variable +enum +{ + NUM_VARS +}; + +namespace UserVariables +{ +static const std::array variable_names = {}; +} + +#include "UserVariables.inc.hpp" + +#endif /* USERVARIABLES_HPP */ diff --git a/Tests/FileInterpolator1DTest/FileInterpolator1DTest.cpp b/Tests/FileInterpolator1DTest/FileInterpolator1DTest.cpp new file mode 100644 index 000000000..35dcee0e0 --- /dev/null +++ b/Tests/FileInterpolator1DTest/FileInterpolator1DTest.cpp @@ -0,0 +1,79 @@ +/* GRChombo + * Copyright 2012 The GRChombo collaboration. + * Please refer to LICENSE in GRChombo's root directory. + */ + +#ifdef CH_LANG_CC +/* + * _______ __ + * / ___/ / ___ __ _ / / ___ + * / /__/ _ \/ _ \/ V \/ _ \/ _ \ + * \___/_//_/\___/_/_/_/_.__/\___/ + * Please refer to LICENSE, in Chombo's root directory. + */ +#endif + +// Chombo includes: +#include "parstream.H" //Gives us pout() + +// Other includes +#include "FileInterpolator1D.hpp" +#include "GRAMR.hpp" +#include "GRParmParse.hpp" +#include "SetupFunctions.hpp" + +#include + +// Chombo namespace +#include "UsingNamespace.H" + +double test_function(double x) { return sin(x * M_PI); } +double test_function_1st_der(double x) { return M_PI * cos(x * M_PI); } + +int runFileInterpolator1D(int argc, char *argv[]) +{ + std::string in_string = argv[argc - 1]; + pout() << in_string << std::endl; + char const *in_file = argv[argc - 1]; + GRParmParse pp(0, argv + argc, NULL, in_file); + + FileInterpolator1D file("N", "x", "y"); + + double test_point = 1. / M_PI; // value surely not in the points in the file + double exact = test_function(test_point); + double val = file.interpolate(test_point); + + double exact_der = test_function_1st_der(test_point); + double val_der = file.interpolate(test_point, 1); + + double err = abs(val / exact - 1.) * 100.; + double err_der = abs(val_der / exact_der - 1.) * 100.; + + pout() << std::setprecision(9); + pout() << "Exact value is: " << exact << std::endl; + pout() << "FileInterpolator1D gave: " << val << " (" << err << "% error)" + << std::endl; + pout() << "Exact 1st derivative value is: " << exact_der << std::endl; + pout() << "FileInterpolator1D gave: " << val_der << " (" << err_der + << "% error)" << std::endl; + + bool wrong = (err + err_der > 1e-2); + + return wrong; +} + +int main(int argc, char *argv[]) +{ + mainSetup(argc, argv); + + int status = runFileInterpolator1D(argc, argv); + + if (status == 0) + std::cout << "FileInterpolator1D test passed." << endl; + else + std::cout << "FileInterpolator1D test failed with return code " + << status << endl; + + mainFinalize(); + return status; +} diff --git a/Tests/FileInterpolator1DTest/FileInterpolator1DTest.inputs b/Tests/FileInterpolator1DTest/FileInterpolator1DTest.inputs new file mode 100644 index 000000000..4052da909 --- /dev/null +++ b/Tests/FileInterpolator1DTest/FileInterpolator1DTest.inputs @@ -0,0 +1,4 @@ +N = 21 +x = 0. 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1. +# sin(x * M_PI) +y = 0. 0.156434 0.309017 0.45399 0.587785 0.707107 0.809017 0.891007 0.951057 0.987688 1. 0.987688 0.951057 0.891007 0.809017 0.707107 0.587785 0.45399 0.309017 0.156434 0. diff --git a/Tests/FileInterpolator1DTest/GNUmakefile b/Tests/FileInterpolator1DTest/GNUmakefile new file mode 100644 index 000000000..f8bb30a52 --- /dev/null +++ b/Tests/FileInterpolator1DTest/GNUmakefile @@ -0,0 +1,21 @@ +# -*- Mode: Makefile -*- + +### This makefile produces an executable for each name in the `ebase' +### variable using the libraries named in the `LibNames' variable. + +# Included makefiles need an absolute path to the Chombo installation +# CHOMBO_HOME := Please set the CHOMBO_HOME locally (e.g. export CHOMBO_HOME=... in bash) + +GRCHOMBO_SOURCE = $(shell pwd)/../../Source + +ebase := FileInterpolator1DTest + +LibNames := AMRTimeDependent AMRTools BoxTools + +src_dirs := $(GRCHOMBO_SOURCE)/GRChomboCore \ + $(GRCHOMBO_SOURCE)/utils \ + $(GRCHOMBO_SOURCE)/simd \ + $(GRCHOMBO_SOURCE)/AMRInterpolator \ + $(GRCHOMBO_SOURCE)/utils/CubicSpline1DInterpolation + +include $(CHOMBO_HOME)/mk/Make.test diff --git a/Tests/FileInterpolator1DTest/SimulationParameters.hpp b/Tests/FileInterpolator1DTest/SimulationParameters.hpp new file mode 100644 index 000000000..5e00a7de6 --- /dev/null +++ b/Tests/FileInterpolator1DTest/SimulationParameters.hpp @@ -0,0 +1,19 @@ +/* GRChombo + * Copyright 2012 The GRChombo collaboration. + * Please refer to LICENSE in GRChombo's root directory. + */ + +#ifndef SIMULATIONPARAMETERS_HPP_ +#define SIMULATIONPARAMETERS_HPP_ + +// General includes +#include "ChomboParameters.hpp" +#include "GRParmParse.hpp" + +class SimulationParameters : public ChomboParameters +{ + public: + SimulationParameters(GRParmParse &pp) : ChomboParameters(pp) {} +}; + +#endif /* SIMULATIONPARAMETERS_HPP_ */ diff --git a/Tests/FileInterpolator1DTest/UserVariables.hpp b/Tests/FileInterpolator1DTest/UserVariables.hpp new file mode 100644 index 000000000..f7642667d --- /dev/null +++ b/Tests/FileInterpolator1DTest/UserVariables.hpp @@ -0,0 +1,26 @@ +/* GRChombo + * Copyright 2012 The GRChombo collaboration. + * Please refer to LICENSE in GRChombo's root directory. + */ + +#ifndef USERVARIABLES_HPP +#define USERVARIABLES_HPP + +#include "EmptyDiagnosticVariables.hpp" +#include +#include + +// assign enum to each variable +enum +{ + NUM_VARS +}; + +namespace UserVariables +{ +static const std::array variable_names = {}; +} + +#include "UserVariables.inc.hpp" + +#endif /* USERVARIABLES_HPP */