add a unit test for the cubic interpolation in pynucastro nets (#1891)

Author: zingale

Docs/source/unit_tests.rst

@@ -117,7 +117,7 @@ by options in the input file.
 One-zone tests
 ==============
 
-.. index:: burn_cell, burn_cell_primordial_chem, burn_cell_sdc, eos_cell, jac_cell, nse_table_cell, nse_net_cell, part_func_cell
+.. index:: burn_cell, burn_cell_primordial_chem, burn_cell_sdc, eos_cell, jac_cell, nse_table_cell, nse_net_cell, part_func_cell, interp_cell
 
 * ``burn_cell`` :
 
@@ -138,6 +138,11 @@ One-zone tests
   given a $\rho$, $T$, and $X_k$, call the equation of state and print out
   the thermodynamic information.  See :ref:`sec:eos_cell` for more information.
 
+* ``interp_cell`` :
+
+  This tests the cubic interpolant used in pynucastro networks for interpolating
+  a rate that is given as pairs of $(T, N_A \langle \sigma v \rangle)$.
+
 * ``jac_cell`` :
 
   for a single thermodynamic state, compute the analytic Jacobian

unit_test/interp_cell/GNUmakefile

@@ -0,0 +1,39 @@
+PRECISION  = DOUBLE
+PROFILE    = FALSE
+
+DEBUG      = FALSE
+
+DIM        = 3
+
+COMP	   = gnu
+
+USE_MPI    = FALSE
+USE_OMP    = FALSE
+
+USE_REACT = TRUE
+
+EBASE = main
+
+BL_NO_FORT = TRUE
+
+# define the location of the Microphysics top directory
+MICROPHYSICS_HOME  := ../..
+
+# This sets the EOS directory
+EOS_DIR     := helmholtz
+
+# This sets the network directory
+NETWORK_DIR := partition_test
+
+CONDUCTIVITY_DIR := stellar
+
+INTEGRATOR_DIR =  VODE
+
+EXTERN_SEARCH += .
+
+Bpack   := ./Make.package
+Blocs   := .
+
+include $(MICROPHYSICS_HOME)/unit_test/Make.unit_test
+
+

unit_test/interp_cell/Make.package

@@ -0,0 +1,2 @@
+CEXE_sources += main.cpp
+CEXE_headers += interp_cell.H

unit_test/interp_cell/README.md

@@ -0,0 +1,4 @@
+# `interp_cell`
+
+Test the interpolation functions in `interp_tools.H` in the pynucastro
+networks.

unit_test/interp_cell/_parameters

@@ -0,0 +1,3 @@
+@namespace: unit_test
+
+

unit_test/interp_cell/interp_cell.H

@@ -0,0 +1,52 @@
+#ifndef INTERP_CELL_H
+#define INTERP_CELL_H
+
+#include <AMReX_Array.H>
+
+#include <extern_parameters.H>
+
+#include <interp_tools.H>
+
+amrex::Real test_function(amrex::Real x) {
+    amrex::Real x0{-2.5};
+    return (2.0_rt/3.0_rt) * amrex::Math::powi<3>(x - x0) - 1.5_rt * amrex::Math::powi<2>(x - x0) + 7.0_rt * (x - x0) - 10.0_rt;
+}
+
+AMREX_INLINE
+void interp_cell_c()
+{
+
+    // create some fake data
+    // our test function is a cubic so we should fit it
+    // exactly with out interpolation
+    amrex::Array1D<amrex::Real, 1, 10> x_array{-3.2, -2.5, -1.23, -0.6, 0.01, 0.54, 0.98, 1.35, 1.87, 2.5};
+    amrex::Array1D<amrex::Real, 1, 10> f_array{};
+    std::cout << "initial data: " << std::endl;
+    for (int i = x_array.lo(); i <= x_array.hi(); ++i) {
+        f_array(i) = test_function(x_array(i));
+        std::cout << x_array(i) << ", " << f_array(i) << std::endl;
+    }
+    std::cout << std::endl;
+
+    // now we see if we recover things well
+    {
+        amrex::Real x0{-1.1};
+        constexpr int do_derivative = 1;
+        const auto [f, fprime] = interp_net::cubic_interp_uneven<do_derivative>(x0, x_array, f_array);
+
+        std::cout << f << " " << test_function(x0) << std::endl;
+    }
+
+    {
+        amrex::Real x0{1.42};
+        constexpr int do_derivative = 1;
+        const auto [f, fprime] = interp_net::cubic_interp_uneven<do_derivative>(x0, x_array, f_array);
+
+        std::cout << f << " " << test_function(x0) << std::endl;
+    }
+
+
+
+
+}
+#endif

unit_test/interp_cell/main.cpp

@@ -0,0 +1,24 @@
+#include <iostream>
+
+#include <extern_parameters.H>
+#include <eos.H>
+#include <network.H>
+#include <interp_cell.H>
+#include <unit_test.H>
+
+int main(int argc, char *argv[]) {
+
+  amrex::Initialize(argc, argv);
+
+  init_unit_test();
+
+  // C++ EOS initialization (must be done after init_extern_parameters)
+  eos_init(unit_test_rp::small_temp, unit_test_rp::small_dens);
+
+  // C++ Network, RHS, screening, rates initialization
+  network_init();
+
+  interp_cell_c();
+
+  amrex::Finalize();
+}