add a neutrino cooling unit test (#1331)

unit_test/test_neutrino_cooling/GNUmakefile

@@ -0,0 +1,43 @@
+PRECISION  = DOUBLE
+PROFILE    = FALSE
+
+DEBUG      = FALSE
+
+DIM        = 3
+
+COMP	   = gnu
+
+USE_MPI    = FALSE
+USE_OMP    = FALSE
+
+USE_REACT = TRUE
+USE_CONDUCTIVITY = FALSE
+
+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 := aprox21
+
+# This isn't actually used but we need VODE to compile with CUDA
+INTEGRATOR_DIR := VODE
+
+CONDUCTIVITY_DIR := stellar
+
+SCREEN_METHOD := screen5
+
+EXTERN_SEARCH += .
+
+Bpack   := ./Make.package
+Blocs   := .
+
+include $(MICROPHYSICS_HOME)/unit_test/Make.unit_test
+
+

unit_test/test_neutrino_cooling/Make.package

@@ -0,0 +1,7 @@
+CEXE_sources += main.cpp
+
+CEXE_headers += test_sneut.H
+
+CEXE_sources += variables.cpp
+CEXE_sources += neutrino_util.cpp
+CEXE_headers += variables.H

unit_test/test_neutrino_cooling/README.md

@@ -0,0 +1,4 @@
+# test_neutrino_cooling
+
+Test the neutrino cooling routine, sneut5
+

unit_test/test_neutrino_cooling/_parameters

@@ -0,0 +1,11 @@
+@namespace: unit_test
+
+dens_min      real       1.d6
+dens_max      real       1.d9
+temp_min      real       1.d6
+temp_max      real       1.d12
+
+metalicity_max  real     0.1d0
+
+small_temp    real        1.d4
+small_dens    real        1.d-4

unit_test/test_neutrino_cooling/inputs

@@ -0,0 +1,9 @@
+n_cell = 16
+max_grid_size = 32
+
+unit_test.dens_min = 10.e0
+unit_test.dens_max = 5.e9
+unit_test.temp_min = 1.e6
+unit_test.temp_max = 1.e10
+
+unit_test.metalicity_max = 0.5e0

unit_test/test_neutrino_cooling/main.cpp

@@ -0,0 +1,166 @@
+#include <AMReX_PlotFileUtil.H>
+#include <AMReX_ParmParse.H>
+#include <AMReX_Print.H>
+
+#include <AMReX_Geometry.H>
+#include <AMReX_MultiFab.H>
+#include <AMReX_BCRec.H>
+
+
+using namespace amrex;
+
+#include <test_sneut.H>
+#include <AMReX_buildInfo.H>
+
+#include <network.H>
+#include <eos.H>
+#include <sneut5.H>
+
+#include <variables.H>
+
+#include <cmath>
+#include <unit_test.H>
+
+using namespace unit_test_rp;
+
+int main (int argc, char* argv[])
+{
+    amrex::Initialize(argc, argv);
+
+    main_main();
+
+    amrex::Finalize();
+    return 0;
+}
+
+void main_main ()
+{
+
+    // AMREX_SPACEDIM: number of dimensions
+    int n_cell, max_grid_size;
+
+    // inputs parameters
+    {
+        // ParmParse is way of reading inputs from the inputs file
+        ParmParse pp;
+
+        // We need to get n_cell from the inputs file - this is the
+        // number of cells on each side of a square (or cubic) domain.
+        pp.get("n_cell", n_cell);
+
+        // The domain is broken into boxes of size max_grid_size
+        max_grid_size = 32;
+        pp.query("max_grid_size", max_grid_size);
+
+    }
+
+
+    Vector<int> is_periodic(AMREX_SPACEDIM,0);
+    for (int idim=0; idim < AMREX_SPACEDIM; ++idim) {
+      is_periodic[idim] = 1;
+    }
+
+    // make BoxArray and Geometry
+    BoxArray ba;
+    Geometry geom;
+    {
+        IntVect dom_lo(AMREX_D_DECL(       0,        0,        0));
+        IntVect dom_hi(AMREX_D_DECL(n_cell-1, n_cell-1, n_cell-1));
+        Box domain(dom_lo, dom_hi);
+
+        // Initialize the boxarray "ba" from the single box "bx"
+        ba.define(domain);
+
+        // Break up boxarray "ba" into chunks no larger than
+        // "max_grid_size" along a direction
+        ba.maxSize(max_grid_size);
+
+        // This defines the physical box, [0, 1] in each direction.
+        RealBox real_box({AMREX_D_DECL(0.0, 0.0, 0.0)},
+                         {AMREX_D_DECL(1.0, 1.0, 1.0)});
+
+        // This defines a Geometry object
+        geom.define(domain, &real_box,
+                    CoordSys::cartesian, is_periodic.data());
+    }
+
+    // Nghost = number of ghost cells for each array
+    int Nghost = 0;
+
+    // do the runtime parameter initializations and microphysics inits
+    if (ParallelDescriptor::IOProcessor()) {
+      std::cout << "reading extern runtime parameters ..." << std::endl;
+    }
+
+    ParmParse ppa("amr");
+
+    init_unit_test();
+
+    eos_init(small_temp, small_dens);
+
+    network_init();
+
+    plot_t vars;
+
+    vars = init_variables();
+
+    amrex::Vector<std::string> names;
+    get_varnames(vars, names);
+
+    // time = starting time in the simulation
+    Real time = 0.0_rt;
+
+    // How Boxes are distributed among MPI processes
+    DistributionMapping dm(ba);
+
+    // we allocate our main multifabs
+    MultiFab state(ba, dm, vars.n_plot_comps, Nghost);
+
+    // Initialize the state to zero; we will fill
+    // it in below in do_eos.
+    state.setVal(0.0_rt);
+
+    // What time is it now?  We'll use this to compute total run time.
+    Real strt_time = ParallelDescriptor::second();
+
+
+    Real dlogrho = 0.0e0_rt;
+    Real dlogT   = 0.0e0_rt;
+    Real dmetal  = 0.0e0_rt;
+
+    if (n_cell > 1) {
+        dlogrho = (std::log10(dens_max) - std::log10(dens_min)) / static_cast<Real>(n_cell - 1);
+        dlogT   = (std::log10(temp_max) - std::log10(temp_min))/ static_cast<Real>(n_cell - 1);
+        dmetal  = (metalicity_max  - 0.0_rt)/ static_cast<Real>(n_cell - 1);
+    }
+
+    // Initialize the state and compute the different thermodynamics
+    // by inverting the EOS
+    for ( MFIter mfi(state); mfi.isValid(); ++mfi )
+    {
+      const Box& bx = mfi.validbox();
+
+      Array4<Real> const sp = state.array(mfi);
+
+      neut_test_C(bx, dlogrho, dlogT, dmetal, vars, sp);
+
+    }
+
+    // Call the timer again and compute the maximum difference between
+    // the start time and stop time over all processors
+    Real stop_time = ParallelDescriptor::second() - strt_time;
+    const int IOProc = ParallelDescriptor::IOProcessorNumber();
+    ParallelDescriptor::ReduceRealMax(stop_time, IOProc);
+
+
+    std::string name = "test_sneut5";
+
+    // Write a plotfile
+    WriteSingleLevelPlotfile(name, state, names, geom, time, 0);
+
+    write_job_info(name);
+
+    // Tell the I/O Processor to write out the "run time"
+    amrex::Print() << "Run time = " << stop_time << std::endl;
+
+}

unit_test/test_neutrino_cooling/neutrino_util.cpp

@@ -0,0 +1,93 @@
+#include <AMReX_PlotFileUtil.H>
+#include <AMReX_ParmParse.H>
+#include <AMReX_Print.H>
+
+#include <AMReX_Geometry.H>
+#include <AMReX_MultiFab.H>
+#include <AMReX_BCRec.H>
+
+#include <variables.H>
+#include <network.H>
+#include <eos.H>
+
+#include <sneut5.H>
+
+#include <cmath>
+
+using namespace amrex;
+using namespace unit_test_rp;
+
+void neut_test_C(const Box& bx,
+                 const Real dlogrho, const Real dlogT, const Real dmetal,
+                 const plot_t& vars,
+                 Array4<Real> const sp) {
+
+  using namespace Species;
+
+  const int ih1 = network_spec_index("hydrogen-1");
+  if (ih1 < 0) amrex::Error("Error: ih1 not found");
+
+  const int ihe4 = network_spec_index("helium-4");
+  if (ihe4 < 0) amrex::Error("Error: ihe4 not found");
+
+
+  amrex::ParallelFor(bx,
+  [=] AMREX_GPU_DEVICE (int i, int j, int k)
+  {
+
+    // set the composition -- approximately solar
+    Real metalicity = 0.0 + static_cast<Real> (k) * dmetal;
+
+    Real xn[NumSpec];
+
+    // for now... the screening using 1-based indexing
+    Array1D<Real, 1, NumSpec> ymass;
+
+    for (int n = 0; n < NumSpec; n++) {
+      xn[n] = metalicity / static_cast<Real>(NumSpec - 2);
+    }
+    xn[ih1] = 0.75_rt - 0.5_rt * metalicity;
+    xn[ihe4] = 0.25_rt - 0.5_rt * metalicity;
+
+    Real temp_zone = std::pow(10.0, std::log10(temp_min) + static_cast<Real>(j)*dlogT);
+
+    Real dens_zone = std::pow(10.0, std::log10(dens_min) + static_cast<Real>(i)*dlogrho);
+
+    // store default state
+    sp(i, j, k, vars.irho) = dens_zone;
+    sp(i, j, k, vars.itemp) = temp_zone;
+    for (int n = 0; n < NumSpec; n++) {
+      sp(i, j, k, vars.ispec+n) = xn[n];
+    }
+
+    // compute abar and zbar
+
+    Real abar = 0.0;
+    for (int n = 0; n < NumSpec; n++) {
+      abar += xn[n] / aion[n];
+    }
+    abar = 1.0_rt / abar;
+
+    Real zbar = 0.0;
+    for (int n = 0; n < NumSpec; n++) {
+      zbar += zion[n] * xn[n] / aion[n];
+    }
+    zbar *= abar;
+
+    Real snu;
+    Real dsnudt;
+    Real dsnudd;
+    Real dsnuda;
+    Real dsnudz;
+
+    sneut5(temp_zone, dens_zone, abar, zbar,
+           snu, dsnudt, dsnudd, dsnuda, dsnudz);
+
+    sp(i, j, k, vars.isneut) = snu;
+    sp(i, j, k, vars.isneutdt) = dsnudt;
+    sp(i, j, k, vars.isneutda) = dsnuda;
+    sp(i, j, k, vars.isneutdz) = dsnudz;
+
+  });
+
+}

unit_test/test_neutrino_cooling/test_sneut.H

@@ -0,0 +1,14 @@
+#ifndef TEST_SCREEN_H
+#define TEST_SCREEN_H
+
+#include <extern_parameters.H>
+#include <variables.H>
+
+void main_main();
+
+void neut_test_C(const Box& bx,
+                 const Real dlogrho, const Real dlogT, const Real dmetal,
+                 const plot_t& vars,
+                 Array4<Real> const sp);
+
+#endif

unit_test/test_neutrino_cooling/variables.H

@@ -0,0 +1,37 @@
+#include <vector>
+#include <string>
+
+#ifndef VARIABLES_H
+#define VARIABLES_H
+
+#include <AMReX_Vector.H>
+
+class plot_t {
+
+public:
+
+
+  int irho = -1;
+  int itemp = -1;
+  int ispec = -1;
+  int isneut = -1;
+  int isneutdt = -1;
+  int isneutda = -1;
+  int isneutdz = -1;
+
+  int n_plot_comps = 0;
+
+  int next_index(const int num) {
+    int next = n_plot_comps;
+    n_plot_comps += num;
+    return next;
+  }
+
+};
+
+plot_t init_variables();
+
+void get_varnames(const plot_t& p, amrex::Vector<std::string>& names);
+
+
+#endif