add a neutrino cooling unit test

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,168 @@
+#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_screen.H>
+#include <AMReX_buildInfo.H>
+
+#include <network.H>
+#include <eos.H>
+#include <screen.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();
+
+    screening_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);
+
+      screen_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_screening." + screen_name;
+
+    // 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;
+
+}