diff --git a/Jenkinsfile b/Jenkinsfile index be2611aa..04fafc5b 100644 --- a/Jenkinsfile +++ b/Jenkinsfile @@ -1,6 +1,6 @@ pipeline { triggers { pollSCM('') } // Run tests whenever a new commit is detected. - agent { dockerfile {args '--gpus all'}} // Use the Dockerfile defined in the root Flash-X directory + agent { dockerfile {args '--gpus all -v /mnt/scratch/tables:/tables:ro'}} // Use the Dockerfile defined in the root Flash-X directory environment { // Get rid of Read -1, expected , errno =1 error // See https://github.com/open-mpi/ompi/issues/4948 @@ -16,7 +16,7 @@ pipeline { sh 'nvidia-smi' sh 'nvcc -V' sh 'git submodule update --init' - sh 'cp makefiles/GNUmakefile_jenkins Exec/GNUmakefile' + sh 'cp makefiles/GNUmakefile_jenkins_HDF5_CUDA Exec/GNUmakefile' dir('Exec'){ sh 'make generate; make -j' } @@ -66,8 +66,8 @@ pipeline { dir('Exec'){ sh 'python ../Scripts/initial_conditions/st4_linear_moment_ffi.py' sh 'mpirun -np 4 ./main3d.gnu.TPROF.MPI.CUDA.ex ../sample_inputs/inputs_1d_fiducial' - sh 'python ../Scripts/data_reduction/reduce_data.py' sh 'python ../Scripts/data_reduction/reduce_data_fft.py' + sh 'python ../Scripts/data_reduction/reduce_data.py' sh 'python ../Scripts/data_reduction/combine_files.py plt _reduced_data.h5' sh 'python ../Scripts/data_reduction/combine_files.py plt _reduced_data_fft_power.h5' sh 'python ../Scripts/babysitting/avgfee.py' @@ -85,14 +85,14 @@ pipeline { sh 'make realclean; make generate; make -j' sh 'python ../Scripts/initial_conditions/st4_linear_moment_ffi_3F.py' sh 'mpirun -np 4 ./main3d.gnu.TPROF.MPI.ex ../sample_inputs/inputs_1d_fiducial' - sh 'python3 ../Scripts/babysitting/avgfee_HDF5.py' + /*sh 'python3 ../Scripts/babysitting/avgfee_HDF5.py'*/ sh 'rm -rf plt*' } }} stage('Collisions flavor instability'){ steps{ dir('Exec'){ - sh 'cp ../makefiles/GNUmakefile_jenkins GNUmakefile' + sh 'cp ../makefiles/GNUmakefile_jenkins_HDF5_CUDA GNUmakefile' sh 'make realclean; make generate NUM_FLAVORS=2; make -j NUM_FLAVORS=2' sh 'python ../Scripts/initial_conditions/st8_coll_inst_test.py' sh 'mpirun -np 4 ./main3d.gnu.TPROF.MPI.CUDA.ex ../sample_inputs/inputs_collisional_instability_test' @@ -104,7 +104,7 @@ pipeline { stage('Collisions to equilibrium'){ steps{ dir('Exec'){ - sh 'cp ../makefiles/GNUmakefile_jenkins GNUmakefile' + sh 'cp ../makefiles/GNUmakefile_jenkins_HDF5_CUDA GNUmakefile' sh 'make realclean; make generate NUM_FLAVORS=3; make -j NUM_FLAVORS=3' sh 'python ../Scripts/initial_conditions/st7_empty_particles.py' sh 'mpirun -np 4 ./main3d.gnu.TPROF.MPI.CUDA.ex ../sample_inputs/inputs_coll_equi_test' diff --git a/Scripts/tests/coll_inst_test.py b/Scripts/tests/coll_inst_test.py index cc3ba35a..4d22f404 100644 --- a/Scripts/tests/coll_inst_test.py +++ b/Scripts/tests/coll_inst_test.py @@ -273,4 +273,4 @@ def QKE(t,y): plt.yscale('log') plt.legend() plt.savefig('EMU_Julien_LucasLSA_Neu.pdf') - plt.close() \ No newline at end of file + plt.close() diff --git a/Source/DataReducer.H b/Source/DataReducer.H index 008ad833..c3f2209e 100644 --- a/Source/DataReducer.H +++ b/Source/DataReducer.H @@ -7,6 +7,11 @@ #include #include #include + +//We use the AMReX binary format to write data for now +//That's because the HDF5 write format gives errors when running with CUDA. +#undef AMREX_USE_HDF5 + #ifdef AMREX_USE_HDF5 #include <../submodules/HighFive/include/highfive/H5File.hpp> #include <../submodules/HighFive/include/highfive/H5DataSpace.hpp> diff --git a/Source/DataReducer.cpp b/Source/DataReducer.cpp index 4c3c15ab..b90799e1 100644 --- a/Source/DataReducer.cpp +++ b/Source/DataReducer.cpp @@ -4,6 +4,11 @@ #include "ArithmeticArray.H" #include #include + +//We use the AMReX binary format to write data for now +//That's because the HDF5 write format gives errors when running with CUDA. +#undef AMREX_USE_HDF5 + #ifdef AMREX_USE_HDF5 #include <../submodules/HighFive/include/highfive/H5File.hpp> #include <../submodules/HighFive/include/highfive/H5DataSpace.hpp> @@ -107,6 +112,8 @@ void DataReducer::InitializeFiles() outfile << j << ":sumTrN\t"; j++; outfile << j << ":sumTrHN\t"; + j++; + outfile << j << ":Vphase\t"; outfile << std::endl; outfile.close(); @@ -127,18 +134,23 @@ DataReducer::WriteReducedData0D(const amrex::Geometry& geom, // Do reductions over the particles // //==================================// using PType = typename FlavoredNeutrinoContainer::ParticleType; - amrex::ReduceOps reduce_ops; - auto particleResult = amrex::ParticleReduce< ReduceData >(neutrinos, - [=] AMREX_GPU_DEVICE(const PType& p) noexcept -> amrex::GpuTuple { + amrex::ReduceOps reduce_ops; + auto particleResult = amrex::ParticleReduce< ReduceData >(neutrinos, + [=] AMREX_GPU_DEVICE(const PType& p) noexcept -> amrex::GpuTuple { Real TrHN = p.rdata(PIdx::TrHN); - Real TrN = 0; + Real TrN = 0; + Real Vphase = p.rdata(PIdx::Vphase); #include "generated_files/DataReducer.cpp_fill_particles" - return GpuTuple{TrN,TrHN}; + return GpuTuple{TrN,TrHN, Vphase}; }, reduce_ops); Real TrN = amrex::get<0>(particleResult); Real TrHN = amrex::get<1>(particleResult); + Real Vphase = amrex::get<2>(particleResult); ParallelDescriptor::ReduceRealSum(TrN); ParallelDescriptor::ReduceRealSum(TrHN); + ParallelDescriptor::ReduceRealSum(Vphase); + + printf("TrN=%g, TrHN=%g, Vphase=%g\n", TrN, TrHN, Vphase); //=============================// // Do reductions over the grid // @@ -326,6 +338,7 @@ DataReducer::WriteReducedData0D(const amrex::Geometry& geom, outfile << N_offdiag_mag << "\t"; outfile << TrN << "\t"; outfile << TrHN << "\t"; + outfile << Vphase << "\t"; outfile << std::endl; outfile.close(); #endif diff --git a/Source/EosTable.H b/Source/EosTable.H new file mode 100644 index 00000000..7f32acfe --- /dev/null +++ b/Source/EosTable.H @@ -0,0 +1,71 @@ +#ifndef EOS_TABLE_H +#define EOS_TABLE_H + +void ReadEosTable(const std::string nuceos_table_name); + +// OLD TODO: remove hard coded constants +// OLD TODO: introduce defines for table index of variables +#define HAVEGR 1 +#define MAX2(x, y) (((x) > (y)) ? (x) : (y)) +#define MIN2(x, y) (((x) < (y)) ? (x) : (y)) +#define NTABLES 19 +#define LENGTHGF 6.77269222552442e-06 +#define TIMEGF 2.03040204956746e05 +#define RHOGF 1.61887093132742e-18 +#define PRESSGF 1.80123683248503e-39 +#define EPSGF 1.11265005605362e-21 +#define INVRHOGF 6.17714470405638e17 +#define INVEPSGF 8.98755178736818e20 +#define INVPRESSGF 5.55174079257738e38 +#define CLIGHT 1 + +// table key +// 0 logpress +// 1 logenergy +// 2 entropy +// 3 munu +// 4 cs2 +// 5 dedt +// 6 dpdrhoe +// 7 dpderho +// 8 muhat +// 9 mu_e +// 10 mu_p +// 11 mu_n +// 12 Xa +// 13 Xh +// 14 Xn +// 15 Xp +// 16 Abar +// 17 Zbar +// 18 Gamma +// enum eos_var {i_logpress=0, i_logenergy, i_entropy, i_munu, i_cs2, i_dedt, + //i_dpdrhoe, i_dpderho, i_muhat, i_mu_e, i_mu_p, i_mu_n, i_Xa, + //i_Xh, i_Xn, i_Xp, i_Abar, i_Zbar, i_Gamma}; +//} + +namespace nuc_eos_private { + +// table data + extern double *alltables; + extern double *epstable; + extern double *logrho; + extern double *logtemp; + extern double *yes; + +#define EOSVAR(VAR) helperVarsReal[helperVarsEnumReal::VAR] +//WARNING: If the number of variables is increased here, then memory allocation for helperVarsReal/helperVarsInt pointer in GRHydroX_ReadEOSTable.cxx will also need to be increased. + enum helperVarsEnumReal { eos_rhomin, eos_rhomax, eos_tempmin, eos_tempmax, eos_yemin, eos_yemax, + energy_shift, dlintemp, dlintempi, drholintempi, dlintempyei, drholintempyei, + dtemp, dtempi, drho, drhoi, dye, dyei, drhotempi, drhoyei, dtempyei, + drhotempyei, eos_epsmin, eos_epsmax }; + extern double *helperVarsReal; + +#define EOSVAR_INT(VAR) helperVarsInt[helperVarsEnumInt::VAR] + enum helperVarsEnumInt { nrho, ntemp, nye } ; + extern int *helperVarsInt; + + +} + +#endif // EOS_TABLE_H diff --git a/Source/EosTableFunctions.H b/Source/EosTableFunctions.H new file mode 100644 index 00000000..fdd6aba9 --- /dev/null +++ b/Source/EosTableFunctions.H @@ -0,0 +1,42 @@ +#ifndef EOSTABLEFUCNTIONS_HXX +#define EOSTABLEFUCNTIONS_HXX + + +#include "EosTableHelpers.H" + +struct EOS_tabulated { + double *alltables, *epstable, *logrho, *logtemp, *yes, *helperVarsReal; + int *helperVarsInt; + + //constructor for Tabulated EOS + AMREX_GPU_DEVICE AMREX_GPU_HOST EOS_tabulated() = default;//Need to keep it + AMREX_GPU_DEVICE AMREX_GPU_HOST EOS_tabulated(double *alltables, double *epstable, double *logrho, + double *logtemp, double *yes, double *helperVarsReal, int *helperVarsInt): + alltables(alltables), epstable(epstable), logrho(logrho), logtemp(logtemp), + yes(yes), helperVarsReal(helperVarsReal), helperVarsInt(helperVarsInt) {} + + //--------------- get helperVarsReal pointer --------------------------------------------- + AMREX_GPU_DEVICE AMREX_GPU_HOST double *get_helperVarsReal() const { + return helperVarsReal; + };//get_helperVarsReal + + + //--------------- get mu_e and muhat for tabulated EOS --------------------------------------------- + AMREX_GPU_DEVICE AMREX_GPU_HOST void get_mue_muhat(double rho, double temperature, double Ye, + double &mue, double &muhat, + int &keyerr, int &anyerr) const { + + nuc_eos_mue_muhat(rho, temperature, Ye, mue, muhat, keyerr, anyerr, + alltables, logrho, logtemp, yes, helperVarsReal, helperVarsInt); + return; + + /*Actual steps: + (1) check bounds + (2) get interp spots + (3) do interpolation to get mu_e and muhat + */ + };//get_mue_muhat + +}; //struct EOS_tabulated + +#endif // EOSTABLEFUCNTIONS_HXX \ No newline at end of file diff --git a/Source/EosTableHelpers.H b/Source/EosTableHelpers.H new file mode 100644 index 00000000..70ac0637 --- /dev/null +++ b/Source/EosTableHelpers.H @@ -0,0 +1,176 @@ +#ifndef EOSTABLEHELPERS_HXX +#define EOSTABLEHELPERS_HXX + +#include "EosTable.H" + +//Check whether input (rho, temperature, Ye) are within the table bounds. +static inline AMREX_GPU_DEVICE AMREX_GPU_HOST int +checkbounds(const double xrho, const double xtemp, const double xye, double *helperVarsReal) { + + using namespace nuc_eos_private; + // keyerr codes: + // 101 -- Y_e too high + // 102 -- Y_e too low + // 103 -- temp too high (if keytemp = 1) + // 104 -- temp too low (if keytemp = 1) + // 105 -- rho too high + // 106 -- rho too low + + if(xrho > EOSVAR(eos_rhomax)) { + return 105; + } + if(xrho < EOSVAR(eos_rhomin)) { + return 106; + } + if(xye > EOSVAR(eos_yemax)) { + return 101; + } + if(xye < EOSVAR(eos_yemin)) { + return 102; + } + if(xtemp > EOSVAR(eos_tempmax)) { + return 103; + } + if(xtemp < EOSVAR(eos_tempmin)) { + return 104; + } + return 0; + +}//function checkbounds + +//Get the indices to prepare for interpolation. +static inline AMREX_GPU_DEVICE AMREX_GPU_HOST void +get_interp_spots(const double x, const double y, const double z, double &delx, + double &dely, double &delz, int *idx, double *logrho, + double *logtemp, double *yes, double *helperVarsReal, int *helperVarsInt) +{ + using namespace nuc_eos_private; + + int ix = 1 + (int)( (x - logrho[0] - 1.0e-10) * EOSVAR(drhoi) ); + int iy = 1 + (int)( (y - logtemp[0] - 1.0e-10) * EOSVAR(dtempi) ); + int iz = 1 + (int)( (z - yes[0] - 1.0e-10) * EOSVAR(dyei) ); + + const int nrho = EOSVAR_INT(nrho); + const int ntemp = EOSVAR_INT(ntemp); + const int nye = EOSVAR_INT(nye); + + ix = MAX2( 1, MIN2( ix, nrho-1 ) ); + iy = MAX2( 1, MIN2( iy, ntemp-1 ) ); + iz = MAX2( 1, MIN2( iz, nye-1 ) ); + + idx[0] = NTABLES*(ix + nrho*(iy + ntemp*iz)); + idx[1] = NTABLES*((ix-1) + nrho*(iy + ntemp*iz)); + idx[2] = NTABLES*(ix + nrho*((iy-1) + ntemp*iz)); + idx[3] = NTABLES*(ix + nrho*(iy + ntemp*(iz-1))); + idx[4] = NTABLES*((ix-1) + nrho*((iy-1) + ntemp*iz)); + idx[5] = NTABLES*((ix-1) + nrho*(iy + ntemp*(iz-1))); + idx[6] = NTABLES*(ix + nrho*((iy-1) + ntemp*(iz-1))); + idx[7] = NTABLES*((ix-1) + nrho*((iy-1) + ntemp*(iz-1))); + + // set up aux vars for interpolation + delx = logrho[ix] - x; + dely = logtemp[iy] - y; + delz = yes[iz] - z; + + return; +} // function get_interp_spots + +//Perform the interpolation +static inline AMREX_GPU_DEVICE AMREX_GPU_HOST void +nuc_eos_C_linterp_one(const int *idx, const double delx, const double dely, + const double delz, double &f, const int iv, + double *alltables, double *helperVarsReal) +{ + using namespace nuc_eos_private; + + // helper variables + double fh[8], a[8]; + + fh[0] = alltables[iv+idx[0]]; + fh[1] = alltables[iv+idx[1]]; + fh[2] = alltables[iv+idx[2]]; + fh[3] = alltables[iv+idx[3]]; + fh[4] = alltables[iv+idx[4]]; + fh[5] = alltables[iv+idx[5]]; + fh[6] = alltables[iv+idx[6]]; + fh[7] = alltables[iv+idx[7]]; + + // set up coeffs of interpolation polynomical and + // evaluate function values + a[0] = fh[0]; + a[1] = EOSVAR(drhoi) * ( fh[1] - fh[0] ); + a[2] = EOSVAR(dtempi) * ( fh[2] - fh[0] ); + a[3] = EOSVAR(dyei) * ( fh[3] - fh[0] ); + a[4] = EOSVAR(drhotempi) * ( fh[4] - fh[1] - fh[2] + fh[0] ); + a[5] = EOSVAR(drhoyei) * ( fh[5] - fh[1] - fh[3] + fh[0] ); + a[6] = EOSVAR(dtempyei) * ( fh[6] - fh[2] - fh[3] + fh[0] ); + a[7] = EOSVAR(drhotempyei) * ( fh[7] - fh[0] + fh[1] + fh[2] + + fh[3] - fh[4] - fh[5] - fh[6] ); + + f = a[0] + a[1] * delx + + a[2] * dely + + a[3] * delz + + a[4] * delx * dely + + a[5] * delx * delz + + a[6] * dely * delz + + a[7] * delx * dely * delz; + + return; +} // function nuc_eos_C_linterp_one + + + +//Main function for mu_e and muhat given (rho, temperature, Ye) +static inline AMREX_GPU_DEVICE AMREX_GPU_HOST void +nuc_eos_mue_muhat(const double rho, const double temperature, + const double Ye, double &mue, double &muhat, + int &keyerr, int &anyerr, double *alltables, + double *logrho, double *logtemp, double *yes, + double *helperVarsReal, int *helperVarsInt) +{ + using namespace nuc_eos_private; + + int anyerr_ = 0; + // check if we are fine + int keyerr_ = checkbounds(rho, temperature, Ye, helperVarsReal); + if(keyerr_ != 0) anyerr_ = 1; + + keyerr = keyerr_ ; anyerr = anyerr_; + // Abort if there is any error in checkbounds. + // This should never happen and the program should abort with + // a fatal error anyway. No point in doing any further EOS calculations. + if(anyerr_){ + printf("Error in checkbounds::%s::%d::%s, keyerr = %d\n", __FILE__, __LINE__, __func__, keyerr); + printf(" rho = %.5e, temperature = %.5e, Ye = %.6f \n", rho, temperature, Ye); + return; + } + + int idx[8]; + double delx,dely,delz; + const double lrho = log(rho); + const double ltemp = log(temperature); + + get_interp_spots(lrho, ltemp, Ye, delx, dely, delz, idx, + logrho, logtemp, yes, helperVarsReal, helperVarsInt); + + double mue_ , muhat_; + { + const int iv = 9; + nuc_eos_C_linterp_one(idx, delx, dely, delz, mue_, iv, alltables, helperVarsReal); + } + + { + const int iv = 8; + nuc_eos_C_linterp_one(idx, delx, dely, delz, muhat_, iv, alltables, helperVarsReal); + } + + // Assign values to reference variables: + mue = mue_; + muhat = muhat_; + + return; +}//nuc_eos_mue_muhat + + + +#endif //EOSTABLEHELPERS_HXX \ No newline at end of file diff --git a/Source/Evolve.cpp b/Source/Evolve.cpp index 62011ae7..31d23957 100644 --- a/Source/Evolve.cpp +++ b/Source/Evolve.cpp @@ -1,8 +1,15 @@ +#include "FlavoredNeutrinoContainer.H" #include "Evolve.H" #include "Constants.H" #include "ParticleInterpolator.H" #include +#include "EosTableFunctions.H" +#include "EosTable.H" + +#include "NuLibTableFunctions.H" +#include "NuLibTable.H" + using namespace amrex; namespace GIdx @@ -149,7 +156,48 @@ void interpolate_rhs_from_mesh(FlavoredNeutrinoContainer& neutrinos_rhs, const M const int shape_factor_order_x = geom.Domain().length(0) > 1 ? SHAPE_FACTOR_ORDER : 0; const int shape_factor_order_y = geom.Domain().length(1) > 1 ? SHAPE_FACTOR_ORDER : 0; const int shape_factor_order_z = geom.Domain().length(2) > 1 ? SHAPE_FACTOR_ORDER : 0; - + + //Create EoS table object + using namespace nuc_eos_private; + EOS_tabulated EOS_tabulated_obj(alltables, epstable, logrho, logtemp, + yes, helperVarsReal, helperVarsInt); + + //Create NuLib table object + using namespace nulib_private; + NuLib_tabulated NuLib_tabulated_obj(alltables_nulib, logrho_nulib, logtemp_nulib, + yes_nulib, helperVarsReal_nulib, helperVarsInt_nulib); + + + + //--------------------------- + const int lev = 0; + +#ifdef _OPENMP +#pragma omp parallel +#endif + for (FNParIter pti(neutrinos_rhs, lev); pti.isValid(); ++pti) + { + const int np = pti.numParticles(); + //printf("(Evolve.cpp) number of particles: np = %d \n", np); + FlavoredNeutrinoContainer::ParticleType* pstruct = &(pti.GetArrayOfStructs()[0]); + + amrex::ParallelFor (np, [=] AMREX_GPU_DEVICE (int i) { + FlavoredNeutrinoContainer::ParticleType& p = pstruct[i]; + + //printf("(Inside Evolve.cpp)i =%d, Vphase = %g \n", i, p.rdata(PIdx::Vphase)); + /*p.pos(0) = p.rdata(PIdx::x); + p.pos(1) = p.rdata(PIdx::y); + p.pos(2) = p.rdata(PIdx::z); + + p.rdata(PIdx::x) = p.pos(0); + p.rdata(PIdx::y) = p.pos(1); + p.rdata(PIdx::z) = p.pos(2);*/ + + }); + } + //-------------------------- + + amrex::MeshToParticle(neutrinos_rhs, state, 0, [=] AMREX_GPU_DEVICE (FlavoredNeutrinoContainer::ParticleType& p, amrex::Array4 const& sarr) @@ -178,18 +226,24 @@ void interpolate_rhs_from_mesh(FlavoredNeutrinoContainer& neutrinos_rhs, const M } // Declare matrices to be used in quantum kinetic equation calculation - Real IMFP_abs[NUM_FLAVORS][NUM_FLAVORS]; // Neutrino inverse mean free path matrix: diag( k_e , k_u , k_t ) - Real IMFP_absbar[NUM_FLAVORS][NUM_FLAVORS]; // Antineutrino inverse mean free path matrix: diag( kbar_e , kbar_u , kbar_t ) + Real IMFP_abs[NUM_FLAVORS][NUM_FLAVORS]; // Neutrino inverse mean free path matrix for nucleon absortion: diag( k_e , k_u , k_t ) + Real IMFP_absbar[NUM_FLAVORS][NUM_FLAVORS]; // Antineutrino inverse mean free path matrix for nucleon absortion: diag( kbar_e , kbar_u , kbar_t ) + Real IMFP_scat[NUM_FLAVORS][NUM_FLAVORS]; // Neutrino inverse mean free path matrix for scatteting: diag( k_e , k_u , k_t ) + Real IMFP_scatbar[NUM_FLAVORS][NUM_FLAVORS]; // Antineutrino inverse mean free path matrix for scatteting: diag( kbar_e , kbar_u , kbar_t ) Real f_eq[NUM_FLAVORS][NUM_FLAVORS]; // Neutrino equilibrium Fermi-dirac distribution matrix: f_eq = diag( f_e , f_u , f_t ) Real f_eqbar[NUM_FLAVORS][NUM_FLAVORS]; // Antineutrino equilibrium Fermi-dirac distribution matrix: f_eq = diag( fbar_e , fbar_u , fbar_t ) - + Real munu[NUM_FLAVORS][NUM_FLAVORS]; // Neutrino chemical potential matrix: munu = diag ( munu_e , munu_x) + Real munubar[NUM_FLAVORS][NUM_FLAVORS]; // Antineutrino chemical potential matrix: munu = diag ( munubar_e , munubar_x) + // Initialize matrices with zeros for (int i=0; iIMFP_abs[0][i]; // Read absorption inverse mean free path from input parameters file. IMFP_absbar[i][i] = parms->IMFP_abs[1][i]; // Read absorption inverse mean free path from input parameters file. + munu[i][i] = parms->munu[0][i]; // Read neutrino chemical potential from input parameters file. + munubar[i][i] = parms->munu[1][i]; // Read antineutrino chemical potential from input parameters file. + + } + } + // If opacity_method is 2, the code interpolate inverse mean free paths from NuLib table and electron neutrino chemical potential from EoS table to compute the collision term. + else if(parms->IMFP_method==2){ + + // Assign temperature, electron fraction, and density at the particle's position to new variables for interpolation of chemical potentials and inverse mean free paths. + Real rho = rho_pp; // Density of background matter at this particle's position g/cm^3 + Real temperature = T_pp; // Temperature of background matter at this particle's position 0.05 //MeV + Real Ye = Ye_pp; // Electron fraction of background matter at this particle's position + + //-------------------- Values from EoS table ------------------------------ + double mue_out, muhat_out; // mue_out : Electron chemical potential. muhat_out : neutron minus proton chemical potential + int keyerr, anyerr; + EOS_tabulated_obj.get_mue_muhat(rho, temperature, Ye, mue_out, muhat_out, keyerr, anyerr); + if (anyerr) assert(0); //If there is an error in interpolation call, stop execution. + +//#define DEBUG_INTERPOLATION_TABLES +#ifdef DEBUG_INTERPOLATION_TABLES + printf("(Evolve.cpp) mu_e interpolated = %f\n", mue_out); + printf("(Evolve.cpp) muhat interpolated = %f\n", muhat_out); +#endif + // munu_val : electron neutrino chemical potential + const double munu_val = mue_out - muhat_out; //munu -> "mu_e" - "muhat" + + munu[0][0] = munu_val; // Save neutrino chemical potential from EOS table in chemical potential matrix + munubar[0][0] = -1.0 * munu_val; // Save antineutrino chemical potential from EOS table in chemical potential matrix + + //--------------------- Values from NuLib table --------------------------- + double *helperVarsReal_nulib = NuLib_tabulated_obj.get_helperVarsReal_nulib(); + int idx_group = NULIBVAR(idx_group); + //FIXME: specify neutrino energy using the following: + // double neutrino_energy = p.rdata(PIdx::pupt); locate energy bin using this. + + //idx_species = {0 for electron neutrino, 1 for electron antineutrino and 2 for all other heavier ones} + //electron neutrino: [0, 0] + int idx_species = 0; + double absorption_opacity, scattering_opacity; + NuLib_tabulated_obj.get_opacities(rho, temperature, Ye, absorption_opacity, scattering_opacity, + keyerr, anyerr, idx_species, idx_group); + if (anyerr) assert(0); + +#ifdef DEBUG_INTERPOLATION_TABLES + printf("(Evolve.cpp) absorption_opacity[e] interpolated = %17.6g\n", absorption_opacity); + printf("(Evolve.cpp) scattering_opacity[e] interpolated = %17.6g\n", scattering_opacity); +#endif + + IMFP_abs[0][0] = absorption_opacity; + IMFP_scat[0][0] = scattering_opacity; + + //electron antineutrino: [1, 0] + idx_species = 1; + NuLib_tabulated_obj.get_opacities(rho, temperature, Ye, absorption_opacity, scattering_opacity, + keyerr, anyerr, idx_species, idx_group); + if (anyerr) assert(0); + +#ifdef DEBUG_INTERPOLATION_TABLES + printf("(Evolve.cpp) absorption_opacity[a] interpolated = %17.6g\n", absorption_opacity); + printf("(Evolve.cpp) scattering_opacity[a] interpolated = %17.6g\n", scattering_opacity); +#endif + + IMFP_absbar[0][0] = absorption_opacity; + IMFP_scatbar[0][0] = scattering_opacity; + + //heavier ones: muon neutrino[0,1], muon antineutruino[1,1], tau neutrino[0,2], tau antineutrino[1,2] + idx_species = 2; + NuLib_tabulated_obj.get_opacities(rho, temperature, Ye, absorption_opacity, scattering_opacity, + keyerr, anyerr, idx_species, idx_group); + if (anyerr) assert(0); + +#ifdef DEBUG_INTERPOLATION_TABLES + printf("(Evolve.cpp) absorption_opacity[x] interpolated = %17.6g\n", absorption_opacity); + printf("(Evolve.cpp) scattering_opacity[x] interpolated = %17.6g\n", scattering_opacity); +#endif + + for (int i=1; ineutrino or 1->antineutrino + // for(int j=1; jelectron, 1->heavy(muon), 2->heavy(tau); all heavy same for current table + IMFP_abs[i][i] = absorption_opacity ; // ... fix it ... + IMFP_absbar[i][i] = absorption_opacity ; // ... fix it ... + IMFP_scat[i][i] = scattering_opacity ; // ... fix it ... + IMFP_scatbar[i][i] = scattering_opacity ; // ... fix it ... + // } + } + //----------------------------------------------------------------------- - // Calculate the Fermi-Dirac distribution for neutrinos and antineutrinos. - f_eq[i][i] = 1. / ( 1. + exp( ( p.rdata( PIdx::pupt ) - parms->munu[0][i] ) / T_pp ) ); - f_eqbar[i][i] = 1. / ( 1. + exp( ( p.rdata( PIdx::pupt ) - parms->munu[1][i] ) / T_pp ) ); + } + else AMREX_ASSERT_WITH_MESSAGE(false, "only available opacity_method is 0, 1 or 2"); - // Include the Pauli blocking term - if (parms->Do_Pauli_blocking == 1){ - IMFP_abs[i][i] = IMFP_abs[i][i] / ( 1 - f_eq[i][i] ) ; // Multiply the absortion inverse mean free path by the Pauli blocking term 1 / (1 - f_eq). - IMFP_absbar[i][i] = IMFP_absbar[i][i] / ( 1 - f_eqbar[i][i] ) ; // Multiply the absortion inverse mean free path by the Pauli blocking term 1 / (1 - f_eq). - } + for (int i=0; iDo_Pauli_blocking == 1){ + IMFP_abs[i][i] = IMFP_abs[i][i] / ( 1 - f_eq[i][i] ) ; // Multiply the absortion inverse mean free path by the Pauli blocking term 1 / (1 - f_eq). + IMFP_absbar[i][i] = IMFP_absbar[i][i] / ( 1 - f_eqbar[i][i] ) ; // Multiply the absortion inverse mean free path by the Pauli blocking term 1 / (1 - f_eq). } + } - else AMREX_ASSERT_WITH_MESSAGE(false, "only available opacity_method is 0 or 1"); + // Compute the time derivative of \( N_{ab} \) using the Quantum Kinetic Equations (QKE). #include "generated_files/Evolve.cpp_dfdt_fill" // set the dfdt values into p.rdata diff --git a/Source/FlavoredNeutrinoContainer.H b/Source/FlavoredNeutrinoContainer.H index 42e13c2a..0dcb4d98 100644 --- a/Source/FlavoredNeutrinoContainer.H +++ b/Source/FlavoredNeutrinoContainer.H @@ -61,6 +61,9 @@ public: } }; +//Create an enum to pass as a template argument to the particle creation function. +enum BoundaryParticleCreationDirection {I_PLUS, I_MINUS, J_PLUS, J_MINUS, K_PLUS, K_MINUS}; + class FlavoredNeutrinoContainer : public amrex::ParticleContainer { @@ -76,6 +79,9 @@ public: void InitParticles(const TestParams* parms); + template + void CreateParticlesAtBoundary(const TestParams* parms, const Real current_dt); + void SyncLocation(int type); void UpdateLocationFrom(FlavoredNeutrinoContainer& Other); @@ -98,4 +104,5 @@ public: static ApplyFlavoredNeutrinoRHS particle_apply_rhs; }; + #endif diff --git a/Source/FlavoredNeutrinoContainerInit.cpp b/Source/FlavoredNeutrinoContainerInit.cpp index 65d68b8f..adbcd81c 100644 --- a/Source/FlavoredNeutrinoContainerInit.cpp +++ b/Source/FlavoredNeutrinoContainerInit.cpp @@ -290,6 +290,7 @@ InitParticles(const TestParams* parms) if(parms->IMFP_method == 1){ p.rdata(PIdx::Vphase) = dx[0]*dx[1]*dx[2]*4*MathConst::pi*(pow(p.rdata(PIdx::pupt)+parms->delta_E/2,3)-pow(p.rdata(PIdx::pupt)-parms->delta_E/2,3))/(3*ndirs_per_loc*parms->nppc[0]*parms->nppc[1]*parms->nppc[2]); + //printf("(Inside FlavoredNeutrinoContainerInit.cpp) Vphase = %g \n", p.rdata(PIdx::Vphase)); } //=====================// @@ -360,3 +361,455 @@ InitParticles(const TestParams* parms) ParallelDescriptor::ReduceRealMin(pupt_min); #include "generated_files/FlavoredNeutrinoContainerInit.cpp_Vvac_fill" } // InitParticles() + + + +//==================================================================================================================// +//========================================= CreateParticlesAtBoundary ==============================================// +//==================================================================================================================// +template +void FlavoredNeutrinoContainer:: +CreateParticlesAtBoundary(const TestParams* parms, const Real current_dt) +{ + BL_PROFILE("FlavoredNeutrinoContainer::CreateParticlesAtBoundary"); + + const int lev = 0; + const auto dx = Geom(lev).CellSizeArray(); + const auto plo = Geom(lev).ProbLoArray(); + const auto& a_bounds = Geom(lev).ProbDomain(); + + //printf("plo = [%f, %f, %f] \n", plo[0], plo[1], plo[2]); + + const int nlocs_per_cell = AMREX_D_TERM( parms->nppc[0], + *parms->nppc[1], + *parms->nppc[2]); + + // array of direction vectors + //TODO: We can use a different custom file to set particle data at boundary points. + Gpu::ManagedVector > particle_data = read_particle_data(parms->particle_data_filename);; + auto* particle_data_p = particle_data.dataPtr(); + + // determine the number of directions per location + int ndirs_per_loc = particle_data.size(); + //amrex::Print() << "Using " << ndirs_per_loc << " directions." << std::endl; + const Real scale_fac = dx[0]*dx[1]*dx[2]/nlocs_per_cell; + + //FIXME: We need to use outflow boundary condition, not periodic boundary condition. Put an assert(!periodic) here. + + + // Loop over multifabs // +#ifdef _OPENMP +#pragma omp parallel +#endif + for (MFIter mfi = MakeMFIter(lev); mfi.isValid(); ++mfi) + { + + //Box tilebox () const noexcept: Return the tile Box at the current index. + const Box& tile_box = mfi.tilebox(); + + const int ncellx = parms->ncell[0]; + const int ncelly = parms->ncell[1]; + const int ncellz = parms->ncell[2]; + + //These actually represent the global indices of the tilebox. + const auto lo = amrex::lbound(tile_box); + const auto hi = amrex::ubound(tile_box); + + //printf("tile_box = [%d, %d, %d] x [%d, %d, %d] \n", lo.x, lo.y, lo.z, hi.x, hi.y, hi.z); + + Gpu::ManagedVector counts(tile_box.numPts(), 0); //PODVector > counts(n, 0) + + unsigned int* pcount = counts.dataPtr(); + + Gpu::ManagedVector offsets(tile_box.numPts()); + unsigned int* poffset = offsets.dataPtr(); + + const int buffer = 0; //TODO: TODO: Set this appropriately. + + // Determine how many particles to add to the particle tile per cell + //This loop runs over all the particles in a given box. + //For each particle, it calculates a unique "cellid". + //It then adds the pcount for that cell by adding ndirs_per_loc value to it (which is number of particles per location emitted). + //From amrex documentation: Tiling is turned off if GPU is enabled so that more parallelism is exposed to GPU kernels. + //Also note that when tiling is off, tilebox returns validbox. + amrex::ParallelFor(tile_box, + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + for (int i_part=0; i_partnppc, i_loc); + + //Real x = plo[0] + (i + r[0])*dx[0]; + //Real x = plo[0] for i_plus direction i.e VC, (y,z) remain same CC. + //Real x = plo[0] + (i+1)*dx[0] for i_minus direction i.e. VC, (y,z) remain same CC. + //Real y = plo[1] + (j + r[1])*dx[1]; + //Real z = plo[2] + (k + r[2])*dx[2]; + Real x, y, z; + + bool create_particle_this_cell = false; + + //Create particles at outer boundary and set face centered coordinates. + //VC=vertex-centered; CC=cell-centered; + switch (DIRECTION) + { + //Create particles in +ve x direction at lower x boundary. + case BoundaryParticleCreationDirection::I_PLUS: + if (i==0+buffer) create_particle_this_cell = true; + x = plo[0]+dx[0]*buffer; //VC, lower x boundary + y = plo[1] + (j + r[1])*dx[1]; //CC + z = plo[2] + (k + r[2])*dx[2]; //CC + break; + + //Create particles in -ve x direction at upper x boundary. + case BoundaryParticleCreationDirection::I_MINUS: + if (i==ncellx-1-buffer) create_particle_this_cell = true; + x = plo[0] + (ncellx-buffer)*dx[0]; //VC, upper x boundary + y = plo[1] + (j + r[1])*dx[1]; //CC + z = plo[2] + (k + r[2])*dx[2]; //CC + break; + + //Create particles in +ve y direction at lower y boundary. + case BoundaryParticleCreationDirection::J_PLUS: + if (j==0+buffer) create_particle_this_cell = true; + y = plo[1]+dx[1]*buffer; //VC, lower y boundary + x = plo[0] + (i + r[0])*dx[0]; //CC + z = plo[2] + (k + r[2])*dx[2]; //CC + break; + + //Create particles in -ve y direction at upper y boundary. + case BoundaryParticleCreationDirection::J_MINUS: + if (j==ncelly-1-buffer) create_particle_this_cell = true; + y = plo[1] + (ncelly-buffer)*dx[1]; //VC, upper y boundary + x = plo[0] + (i + r[0])*dx[0]; //CC + z = plo[2] + (k + r[2])*dx[2]; //CC + break; + + //Create particles in +ve z direction at lower z boundary. + case BoundaryParticleCreationDirection::K_PLUS: + if (k==0+buffer) create_particle_this_cell = true; + z = plo[2]+dx[2]*buffer; //VC, lower z boundary + x = plo[0] + (i + r[0])*dx[0]; //CC + y = plo[1] + (j + r[1])*dx[1]; //CC + break; + + //Create particles in -ve z direction at upper z boundary. + case BoundaryParticleCreationDirection::K_MINUS: + if (k==ncellz-1-buffer) create_particle_this_cell = true; + z = plo[2] + (ncellz-buffer)*dx[2]; //VC, upper z boundary + x = plo[0] + (i + r[0])*dx[0]; //CC + y = plo[1] + (j + r[1])*dx[1]; //CC + break; + + default: + printf("Invalid direction specified. \n"); + assert(0); + break; + } + + //printf("x = %f, y = %f, z = %f \n", x, y, z); + + if (!create_particle_this_cell) continue; + //printf("CREATE PARTRICLE AT: i = %d, j = %d, k = %d \n", i, j, k); + + for(int i_direction=0; i_direction= 0); + AMREX_ASSERT(p.rdata(PIdx::N11_Re ) >= 0); + AMREX_ASSERT(p.rdata(PIdx::N00_Rebar) >= 0); + AMREX_ASSERT(p.rdata(PIdx::N11_Rebar) >= 0); +#if NUM_FLAVORS==3 + AMREX_ASSERT(p.rdata(PIdx::N22_Re ) >= 0); + AMREX_ASSERT(p.rdata(PIdx::N22_Rebar) >= 0); +#endif + + // Set particle position + p.pos(0) = x; + p.pos(1) = y; + p.pos(2) = z; + + // Set particle integrated position + p.rdata(PIdx::x) = x; + p.rdata(PIdx::y) = y; + p.rdata(PIdx::z) = z; + p.rdata(PIdx::time) = 0; + + // scale particle numbers based on number of points per cell and the cell volume + p.rdata(PIdx::N00_Re ) *= scale_fac; + p.rdata(PIdx::N11_Re ) *= scale_fac; + p.rdata(PIdx::N00_Rebar) *= scale_fac; + p.rdata(PIdx::N11_Rebar) *= scale_fac; +#if NUM_FLAVORS==3 + p.rdata(PIdx::N22_Re ) *= scale_fac; + p.rdata(PIdx::N22_Rebar) *= scale_fac; +#endif + + if(parms->IMFP_method == 1){ + const Real V_momentum = 4*MathConst::pi*(pow(p.rdata(PIdx::pupt)+parms->delta_E/2,3)-pow(p.rdata(PIdx::pupt)-parms->delta_E/2,3))/(3*ndirs_per_loc*parms->nppc[0]*parms->nppc[1]*parms->nppc[2]); + + //p.rdata(PIdx::Vphase) = dx[0]*dx[1]*dx[2]*V_momentum; + const Real dt = current_dt; + const Real clight = PhysConst::c; + const Real pupx_ = p.rdata(PIdx::pupx); //TODO: Review this. + const Real pupy_ = p.rdata(PIdx::pupy); //TODO: Review this. + const Real pupz_ = p.rdata(PIdx::pupz); //TODO: Review this. + const Real pupt_ = p.rdata(PIdx::pupt); //TODO: Review this. + + switch (DIRECTION) + { + //Create particles in +ve x direction at lower x boundary. + case BoundaryParticleCreationDirection::I_PLUS: + p.rdata(PIdx::Vphase) = dx[1]*dx[2]*clight*dt*V_momentum*std::abs(pupx_/pupt_); + //printf("(Inside FlavoredNeutrinoContainerInit.cpp:I_PLUS) Vphase = %g \n", p.rdata(PIdx::Vphase)); + break; + + //Create particles in -ve x direction at upper x boundary. + case BoundaryParticleCreationDirection::I_MINUS: + p.rdata(PIdx::Vphase) = dx[1]*dx[2]*clight*dt*V_momentum*std::abs(pupx_/pupt_); + //printf("(Inside FlavoredNeutrinoContainerInit.cpp:I_MINUS) Vphase = %g \n", p.rdata(PIdx::Vphase)); + break; + + //Create particles in +ve y direction at lower y boundary. + case BoundaryParticleCreationDirection::J_PLUS: + p.rdata(PIdx::Vphase) = dx[0]*dx[2]*clight*dt*V_momentum*std::abs(pupy_/pupt_); + //printf("(Inside FlavoredNeutrinoContainerInit.cpp:J_PLUS) Vphase = %g \n", p.rdata(PIdx::Vphase)); + break; + + //Create particles in -ve y direction at upper y boundary. + case BoundaryParticleCreationDirection::J_MINUS: + p.rdata(PIdx::Vphase) = dx[0]*dx[2]*clight*dt*V_momentum*std::abs(pupy_/pupt_); + //printf("(Inside FlavoredNeutrinoContainerInit.cpp:J_MINUS) Vphase = %g \n", p.rdata(PIdx::Vphase)); + break; + + //Create particles in +ve z direction at lower z boundary. + case BoundaryParticleCreationDirection::K_PLUS: + p.rdata(PIdx::Vphase) = dx[0]*dx[1]*clight*dt*V_momentum*std::abs(pupz_/pupt_); + //printf("(Inside FlavoredNeutrinoContainerInit.cpp:K_PLUS) Vphase = %g \n", p.rdata(PIdx::Vphase)); + break; + + //Create particles in -ve z direction at upper z boundary. + case BoundaryParticleCreationDirection::K_MINUS: + p.rdata(PIdx::Vphase) = dx[0]*dx[1]*clight*dt*V_momentum*std::abs(pupz_/pupt_); + //printf("(Inside FlavoredNeutrinoContainerInit.cpp:K_MINUS) Vphase = %g \n", p.rdata(PIdx::Vphase)); + break; + + default: + printf("Invalid direction specified. \n"); + assert(0); + break; + } + } + + //TODO: Do not apply perturbation in this case. + //=====================// + // Apply Perturbations // + //=====================// + if(parms->perturbation_type == 0){ + // random perturbations to the off-diagonals + /*Real rand; + symmetric_uniform(&rand, engine); + p.rdata(PIdx::N01_Re) = parms->perturbation_amplitude*rand * (p.rdata(PIdx::N00_Re ) - p.rdata(PIdx::N11_Re )); + symmetric_uniform(&rand, engine); + p.rdata(PIdx::N01_Im) = parms->perturbation_amplitude*rand * (p.rdata(PIdx::N00_Re ) - p.rdata(PIdx::N11_Re )); + symmetric_uniform(&rand, engine); + p.rdata(PIdx::N01_Rebar) = parms->perturbation_amplitude*rand * (p.rdata(PIdx::N00_Rebar) - p.rdata(PIdx::N11_Rebar)); + symmetric_uniform(&rand, engine); + p.rdata(PIdx::N01_Imbar) = parms->perturbation_amplitude*rand * (p.rdata(PIdx::N00_Rebar) - p.rdata(PIdx::N11_Rebar));*/ + p.rdata(PIdx::N01_Re) = - p.rdata(PIdx::N11_Re); + p.rdata(PIdx::N01_Im) = - p.rdata(PIdx::N11_Re); + p.rdata(PIdx::N01_Rebar) = - p.rdata(PIdx::N11_Rebar); + p.rdata(PIdx::N01_Imbar) = - p.rdata(PIdx::N11_Rebar); +#if NUM_FLAVORS==3 + /*symmetric_uniform(&rand, engine); + p.rdata(PIdx::N02_Re) = parms->perturbation_amplitude*rand * (p.rdata(PIdx::N00_Re ) - p.rdata(PIdx::N22_Re )); + symmetric_uniform(&rand, engine); + p.rdata(PIdx::N02_Im) = parms->perturbation_amplitude*rand * (p.rdata(PIdx::N00_Re ) - p.rdata(PIdx::N22_Re )); + symmetric_uniform(&rand, engine); + p.rdata(PIdx::N12_Re) = parms->perturbation_amplitude*rand * (p.rdata(PIdx::N11_Re ) - p.rdata(PIdx::N22_Re )); + symmetric_uniform(&rand, engine); + p.rdata(PIdx::N12_Im) = parms->perturbation_amplitude*rand * (p.rdata(PIdx::N11_Re ) - p.rdata(PIdx::N22_Re )); + symmetric_uniform(&rand, engine); + p.rdata(PIdx::N02_Rebar) = parms->perturbation_amplitude*rand * (p.rdata(PIdx::N00_Rebar) - p.rdata(PIdx::N22_Rebar)); + symmetric_uniform(&rand, engine); + p.rdata(PIdx::N02_Imbar) = parms->perturbation_amplitude*rand * (p.rdata(PIdx::N00_Rebar) - p.rdata(PIdx::N22_Rebar)); + symmetric_uniform(&rand, engine); + p.rdata(PIdx::N12_Rebar) = parms->perturbation_amplitude*rand * (p.rdata(PIdx::N11_Rebar) - p.rdata(PIdx::N22_Rebar)); + symmetric_uniform(&rand, engine); + p.rdata(PIdx::N12_Imbar) = parms->perturbation_amplitude*rand * (p.rdata(PIdx::N11_Rebar) - p.rdata(PIdx::N22_Rebar));*/ + p.rdata(PIdx::N02_Re) = - p.rdata(PIdx::N22_Re); + p.rdata(PIdx::N02_Im) = - p.rdata(PIdx::N22_Re); + p.rdata(PIdx::N12_Re) = - p.rdata(PIdx::N22_Re); + p.rdata(PIdx::N12_Im) = - p.rdata(PIdx::N22_Re); + p.rdata(PIdx::N02_Rebar) = - p.rdata(PIdx::N22_Rebar); + p.rdata(PIdx::N02_Imbar) = - p.rdata(PIdx::N22_Rebar); + p.rdata(PIdx::N12_Rebar) = - p.rdata(PIdx::N22_Rebar); + p.rdata(PIdx::N12_Imbar) = - p.rdata(PIdx::N22_Rebar); +#endif + } + if(parms->perturbation_type == 1){ + // Perturb real part of e-mu component only sinusoidally in z + /*Real nu_k = (2.*M_PI) / parms->perturbation_wavelength_cm; + p.rdata(PIdx::N01_Re) = parms->perturbation_amplitude*sin(nu_k*p.pos(2)) * (p.rdata(PIdx::N00_Re ) - p.rdata(PIdx::N11_Re )); + p.rdata(PIdx::N01_Rebar) = parms->perturbation_amplitude*sin(nu_k*p.pos(2)) * (p.rdata(PIdx::N00_Rebar) - p.rdata(PIdx::N11_Rebar));*/ + p.rdata(PIdx::N01_Re) = - p.rdata(PIdx::N11_Re); + p.rdata(PIdx::N01_Rebar) = - p.rdata(PIdx::N11_Rebar); + } + + } // loop over direction + } // loop over location + }); // loop over grid cells + + //printf("Finished creating particles at boundary. \n"); + } // loop over multifabs + +} // CreateParticlesAtBoundary() + +//We need to explicitly instantiate the template function for different use cases. +//DIRECTION == BoundaryParticleCreationDirection::I_PLUS (+ve x direction at lower x boundary.) +template void FlavoredNeutrinoContainer::CreateParticlesAtBoundary(const TestParams* parms, const Real current_dt); +//DIRECTION == BoundaryParticleCreationDirection::I_MINUS (-ve x direction at upper x boundary.) +template void FlavoredNeutrinoContainer::CreateParticlesAtBoundary(const TestParams* parms, const Real current_dt); +//DIRECTION == BoundaryParticleCreationDirection::J_PLUS (+ve y direction at lower y boundary.) +template void FlavoredNeutrinoContainer::CreateParticlesAtBoundary(const TestParams* parms, const Real current_dt); +//DIRECTION == BoundaryParticleCreationDirection::J_MINUS (-ve y direction at upper y boundary.) +template void FlavoredNeutrinoContainer::CreateParticlesAtBoundary(const TestParams* parms, const Real current_dt); +//DIRECTION == BoundaryParticleCreationDirection::K_PLUS (+ve z direction at lower z boundary.) +template void FlavoredNeutrinoContainer::CreateParticlesAtBoundary(const TestParams* parms, const Real current_dt); +//DIRECTION == BoundaryParticleCreationDirection::K_MINUS (-ve z direction at upper z boundary.) +template void FlavoredNeutrinoContainer::CreateParticlesAtBoundary(const TestParams* parms, const Real current_dt); + + diff --git a/Source/IO.cpp b/Source/IO.cpp index 8d5bd473..2a839ef8 100644 --- a/Source/IO.cpp +++ b/Source/IO.cpp @@ -2,6 +2,10 @@ #include #include +//We use the AMReX binary format to write data for now +//That's because the HDF5 write format gives errors when running with CUDA. +#undef AMREX_USE_HDF5 + #ifdef AMREX_USE_HDF5 #include "hdf5.h" #endif diff --git a/Source/Make.package b/Source/Make.package index b038507d..4941a84f 100644 --- a/Source/Make.package +++ b/Source/Make.package @@ -4,6 +4,9 @@ CEXE_sources += IO.cpp CEXE_sources += FlavoredNeutrinoContainerInit.cpp CEXE_sources += Evolve.cpp CEXE_sources += FlavoredNeutrinoContainer.cpp +CEXE_sources += ReadEosTable.cpp +CEXE_sources += ReadNuLibTable.cpp +CEXE_sources += ReadInput_RhoTempYe.cpp CEXE_headers += Evolve.H CEXE_headers += FlavoredNeutrinoContainer.H @@ -13,4 +16,11 @@ CEXE_headers += IO.H CEXE_headers += Parameters.H CEXE_headers += ParticleInterpolator.H CEXE_headers += DataReducer.H -CEXE_headers += ArithmeticArray.H \ No newline at end of file +CEXE_headers += ArithmeticArray.H +CEXE_headers += EosTable.H +CEXE_headers += EosTableFunctions.H +CEXE_headers += EosTableHelpers.H +CEXE_headers += NuLibTable.H +CEXE_headers += NuLibTableFunctions.H +CEXE_headers += NuLibTableHelpers.H +CEXE_headers += ReadInput_RhoTempYe.H diff --git a/Source/NuLibTable.H b/Source/NuLibTable.H new file mode 100644 index 00000000..262460ed --- /dev/null +++ b/Source/NuLibTable.H @@ -0,0 +1,37 @@ +#ifndef NULIB_TABLE_H +#define NULIB_TABLE_H + +void ReadNuLibTable(const std::string nulib_table_name); + +#define NTABLES_NULIB 2 + +// table key +// 0 absoprtion_opacity +// 1 scattering opacity + + +namespace nulib_private { + +// table data + extern double *alltables_nulib; + //extern double *epstable; + extern double *logrho_nulib; + extern double *logtemp_nulib; + extern double *yes_nulib; + extern double *species_nulib; + extern double *group_nulib; + +#define NULIBVAR(VAR) helperVarsReal_nulib[helperVarsEnumReal_nulib::VAR] +//WARNING: If the number of variables is increased here, then memory allocation for helperVarsReal/helperVarsInt pointer in GRHydroX_ReadEOSTable.cxx will also need to be increased. + enum helperVarsEnumReal_nulib { eos_rhomin, eos_rhomax, eos_tempmin, eos_tempmax, eos_yemin, eos_yemax, + dlintemp, dlintempi, drholintempi, dlintempyei, drholintempyei, + dtemp, dtempi, drho, drhoi, dye, dyei, drhotempi, drhoyei, dtempyei, + drhotempyei, eos_epsmin, eos_epsmax, idx_group}; + extern double *helperVarsReal_nulib; + +#define NULIBVAR_INT(VAR) helperVarsInt_nulib[helperVarsEnumInt_nulib::VAR] + enum helperVarsEnumInt_nulib { nrho, ntemp, nye, nspecies, ngroup } ; + extern int *helperVarsInt_nulib; +} + +#endif // NULIB_TABLE_H diff --git a/Source/NuLibTableFunctions.H b/Source/NuLibTableFunctions.H new file mode 100644 index 00000000..95f2c198 --- /dev/null +++ b/Source/NuLibTableFunctions.H @@ -0,0 +1,44 @@ +#ifndef NULIBTABLEFUCNTIONS_HXX +#define NULIBTABLEFUCNTIONS_HXX + + +#include "NuLibTableHelpers.H" + +struct NuLib_tabulated { + double *alltables_nulib, *logrho_nulib, *logtemp_nulib, *yes_nulib, *helperVarsReal_nulib; + int *helperVarsInt_nulib; + + //constructor for NuLib_tabulated + AMREX_GPU_DEVICE AMREX_GPU_HOST NuLib_tabulated() = default;//Need to keep it + AMREX_GPU_DEVICE AMREX_GPU_HOST NuLib_tabulated(double *alltables_nulib, double *logrho_nulib, + double *logtemp_nulib, double *yes_nulib, double *helperVarsReal_nulib, int *helperVarsInt_nulib): + alltables_nulib(alltables_nulib), logrho_nulib(logrho_nulib), logtemp_nulib(logtemp_nulib), + yes_nulib(yes_nulib), helperVarsReal_nulib(helperVarsReal_nulib), helperVarsInt_nulib(helperVarsInt_nulib) {} + + //--------------- get helperVarsReal pointer --------------------------------------------- + AMREX_GPU_DEVICE AMREX_GPU_HOST double *get_helperVarsReal_nulib() const { + return helperVarsReal_nulib; + };//get_helperVarsReal_nulib + + + //--------------- get opacaties --------------------------------------------- + AMREX_GPU_DEVICE AMREX_GPU_HOST void get_opacities(double rho, double temperature, double Ye, + double &absorption_opacity, double &scattering_opacity, + int &keyerr, int &anyerr, const int idx_species, const int idx_group) const { + + nulib_opacities(rho, temperature, Ye, absorption_opacity, scattering_opacity, keyerr, anyerr, + alltables_nulib, logrho_nulib, logtemp_nulib, yes_nulib, helperVarsReal_nulib, helperVarsInt_nulib, + idx_species, idx_group); + return; + + /*Actual steps: + (1) check bounds + (2) get interp spots + (3) do interpolation to get absorption and scattering opacity + */ + };//get_opacities + + +}; //struct NuLib_tabulated + +#endif // NULIBTABLEFUCNTIONS_HXX \ No newline at end of file diff --git a/Source/NuLibTableHelpers.H b/Source/NuLibTableHelpers.H new file mode 100644 index 00000000..cc41bc1c --- /dev/null +++ b/Source/NuLibTableHelpers.H @@ -0,0 +1,191 @@ +#ifndef NULIBTABLEHELPERS_HXX +#define NULIBTABLEHELPERS_HXX + +#include "NuLibTable.H" + +//Check whether input (rho, temperature, Ye) are within the table bounds. +static inline AMREX_GPU_DEVICE AMREX_GPU_HOST int +checkbounds_nulib(const double xrho, const double xtemp, const double xye, double *helperVarsReal_nulib) { + + using namespace nulib_private; + // keyerr codes: + // 101 -- Y_e too high + // 102 -- Y_e too low + // 103 -- temp too high (if keytemp = 1) + // 104 -- temp too low (if keytemp = 1) + // 105 -- rho too high + // 106 -- rho too low + + if(xrho > NULIBVAR(eos_rhomax)) { + return 105; + } + if(xrho < NULIBVAR(eos_rhomin)) { + return 106; + } + if(xye > NULIBVAR(eos_yemax)) { + return 101; + } + if(xye < NULIBVAR(eos_yemin)) { + return 102; + } + if(xtemp > NULIBVAR(eos_tempmax)) { + return 103; + } + if(xtemp < NULIBVAR(eos_tempmin)) { + return 104; + } + return 0; + +}//function checkbounds_nulib + +//Get the indices to prepare for interpolation. +static inline AMREX_GPU_DEVICE AMREX_GPU_HOST void +get_interp_spots_nulib(const double x, const double y, const double z, double &delx, + double &dely, double &delz, int *idx, double *logrho_nulib, + double *logtemp_nulib, double *yes_nulib, double *helperVarsReal_nulib, int *helperVarsInt_nulib, + const int idx_species, const int idx_group) +{ + using namespace nulib_private; + + int ix = 1 + (int)( (x - logrho_nulib[0] - 1.0e-10) * NULIBVAR(drhoi) ); + int iy = 1 + (int)( (y - logtemp_nulib[0] - 1.0e-10) * NULIBVAR(dtempi) ); + int iz = 1 + (int)( (z - yes_nulib[0] - 1.0e-10) * NULIBVAR(dyei) ); + + const int nrho = NULIBVAR_INT(nrho); + const int ntemp = NULIBVAR_INT(ntemp); + const int nye = NULIBVAR_INT(nye); + const int nspecies = NULIBVAR_INT(nspecies); + + ix = MAX2( 1, MIN2( ix, nrho-1 ) ); + iy = MAX2( 1, MIN2( iy, ntemp-1 ) ); + iz = MAX2( 1, MIN2( iz, nye-1 ) ); + + //int indnew = iv + NTABLES_NULIB*(i + nrho_*(j + ntemp_*k)); //indexing in EoS table + //int indnew = iv + NTABLES_NULIB*(i + nrho_*(j + ntemp_*(k + nye_*(l + nspecies_*m)))); //indexing in NuLib table + + idx[0] = NTABLES_NULIB*(ix + nrho*(iy + ntemp*(iz + nye*(idx_species + nspecies*idx_group)))); + idx[1] = NTABLES_NULIB*((ix-1) + nrho*(iy + ntemp*(iz + nye*(idx_species + nspecies*idx_group)))); + idx[2] = NTABLES_NULIB*(ix + nrho*((iy-1) + ntemp*(iz + nye*(idx_species + nspecies*idx_group)))); + idx[3] = NTABLES_NULIB*(ix + nrho*(iy + ntemp*((iz-1) + nye*(idx_species + nspecies*idx_group)))); + idx[4] = NTABLES_NULIB*((ix-1) + nrho*((iy-1) + ntemp*(iz + nye*(idx_species + nspecies*idx_group)))); + idx[5] = NTABLES_NULIB*((ix-1) + nrho*(iy + ntemp*((iz-1) + nye*(idx_species + nspecies*idx_group)))); + idx[6] = NTABLES_NULIB*(ix + nrho*((iy-1) + ntemp*((iz-1) + nye*(idx_species + nspecies*idx_group)))); + idx[7] = NTABLES_NULIB*((ix-1) + nrho*((iy-1) + ntemp*((iz-1) + nye*(idx_species + nspecies*idx_group)))); + + //idx[0] = NTABLES_NULIB*(ix + nrho*(iy + ntemp*iz)); + //idx[1] = NTABLES_NULIB*((ix-1) + nrho*(iy + ntemp*iz)); + //idx[2] = NTABLES_NULIB*(ix + nrho*((iy-1) + ntemp*iz)); + //idx[3] = NTABLES_NULIB*(ix + nrho*(iy + ntemp*(iz-1))); + //idx[4] = NTABLES_NULIB*((ix-1) + nrho*((iy-1) + ntemp*iz)); + //idx[5] = NTABLES_NULIB*((ix-1) + nrho*(iy + ntemp*(iz-1))); + //idx[6] = NTABLES_NULIB*(ix + nrho*((iy-1) + ntemp*(iz-1))); + //idx[7] = NTABLES_NULIB*((ix-1) + nrho*((iy-1) + ntemp*(iz-1))); + + // set up aux vars for interpolation + delx = logrho_nulib[ix] - x; + dely = logtemp_nulib[iy] - y; + delz = yes_nulib[iz] - z; + + return; +} // function get_interp_spots_nulib + +//Perform the interpolation +static inline AMREX_GPU_DEVICE AMREX_GPU_HOST void +nuc_eos_C_linterp_one_nulib(const int *idx, const double delx, const double dely, + const double delz, double &f, const int iv, + double *alltables_nulib, double *helperVarsReal_nulib) +{ + using namespace nulib_private; + + // helper variables + double fh[8], a[8]; + + fh[0] = alltables_nulib[iv+idx[0]]; + fh[1] = alltables_nulib[iv+idx[1]]; + fh[2] = alltables_nulib[iv+idx[2]]; + fh[3] = alltables_nulib[iv+idx[3]]; + fh[4] = alltables_nulib[iv+idx[4]]; + fh[5] = alltables_nulib[iv+idx[5]]; + fh[6] = alltables_nulib[iv+idx[6]]; + fh[7] = alltables_nulib[iv+idx[7]]; + + // set up coeffs of interpolation polynomical and + // evaluate function values + a[0] = fh[0]; + a[1] = NULIBVAR(drhoi) * ( fh[1] - fh[0] ); + a[2] = NULIBVAR(dtempi) * ( fh[2] - fh[0] ); + a[3] = NULIBVAR(dyei) * ( fh[3] - fh[0] ); + a[4] = NULIBVAR(drhotempi) * ( fh[4] - fh[1] - fh[2] + fh[0] ); + a[5] = NULIBVAR(drhoyei) * ( fh[5] - fh[1] - fh[3] + fh[0] ); + a[6] = NULIBVAR(dtempyei) * ( fh[6] - fh[2] - fh[3] + fh[0] ); + a[7] = NULIBVAR(drhotempyei) * ( fh[7] - fh[0] + fh[1] + fh[2] + + fh[3] - fh[4] - fh[5] - fh[6] ); + + f = a[0] + a[1] * delx + + a[2] * dely + + a[3] * delz + + a[4] * delx * dely + + a[5] * delx * delz + + a[6] * dely * delz + + a[7] * delx * dely * delz; + + return; +} // function nuc_eos_C_linterp_one_nulib + + + +//Main function for absorption and scattering opacities (rho, temperature, Ye) +static inline AMREX_GPU_DEVICE AMREX_GPU_HOST void +nulib_opacities(const double rho, const double temperature, + const double Ye, double &absorption_opacity, double &scattering_opacity, + int &keyerr, int &anyerr, double *alltables_nulib, + double *logrho_nulib, double *logtemp_nulib, double *yes_nulib, + double *helperVarsReal_nulib, int *helperVarsInt_nulib, + const int idx_species, const int idx_group) +{ + using namespace nulib_private; + + int anyerr_ = 0; + // check if we are fine + int keyerr_ = checkbounds_nulib(rho, temperature, Ye, helperVarsReal_nulib); + if(keyerr_ != 0) anyerr_ = 1; + + keyerr = keyerr_ ; anyerr = anyerr_; + // Abort if there is any error in checkbounds. + // This should never happen and the program should abort with + // a fatal error anyway. No point in doing any further EOS calculations. + if(anyerr_){ + printf("Error in checkbounds::%s::%d::%s, keyerr = %d\n", __FILE__, __LINE__, __func__, keyerr); + printf(" rho = %.5e, temperature = %.5e, Ye = %.6f \n", rho, temperature, Ye); + return; + } + + int idx[8]; + double delx,dely,delz; + const double lrho = log(rho); + const double ltemp = log(temperature); + + get_interp_spots_nulib(lrho, ltemp, Ye, delx, dely, delz, idx, + logrho_nulib, logtemp_nulib, yes_nulib, helperVarsReal_nulib, helperVarsInt_nulib, + idx_species, idx_group); + + double absorption_opacity_ , scattering_opacity_; + { + const int iv = 0; + nuc_eos_C_linterp_one_nulib(idx, delx, dely, delz, absorption_opacity_, iv, alltables_nulib, helperVarsReal_nulib); + } + + { + const int iv = 1; + nuc_eos_C_linterp_one_nulib(idx, delx, dely, delz, scattering_opacity_, iv, alltables_nulib, helperVarsReal_nulib); + } + + // Assign values to reference variables: + absorption_opacity = absorption_opacity_; + scattering_opacity = scattering_opacity_; + + return; +}//nulib_opacities + + +#endif //NULIBTABLEHELPERS_HXX \ No newline at end of file diff --git a/Source/Parameters.H b/Source/Parameters.H index 201e0cd3..3e909954 100644 --- a/Source/Parameters.H +++ b/Source/Parameters.H @@ -54,6 +54,10 @@ struct TestParams : public amrex::Gpu::Managed // attenuation to hamiltonians Real attenuation_hamiltonians; + //HDF5 table names (with full path) for EoS and NuLib tables + std::string nuceos_table_name; + std::string nulib_table_name; + void Initialize(){ ParmParse pp; pp.get("ncell", ncell); @@ -137,12 +141,15 @@ struct TestParams : public amrex::Gpu::Managed munu[0][i] *= 1e6*CGSUnitsConst::eV; munu[1][i] *= 1e6*CGSUnitsConst::eV; } - pp.get("Do_Pauli_blocking", Do_Pauli_blocking); // If 1, it will multiply the opacities by 1 / (1 - f_eq); if 0, do nothing. - pp.get("delta_E", delta_E); + pp.get("Do_Pauli_blocking", Do_Pauli_blocking); // If 1, it will multiply the opacities by 1 / (1 - f_eq); if 0, do nothing. + pp.get("delta_E", delta_E); delta_E*= 1e6*CGSUnitsConst::eV; // erg - } - else AMREX_ASSERT_WITH_MESSAGE(false, "only available opacity_method is 0(do not collisions) and 1(do collisions)"); + } else if (IMFP_method==2){ + //HDF5 table names (with full path) for EoS and NuLib tables + pp.get("nuceos_table_name", nuceos_table_name); + pp.get("nulib_table_name", nulib_table_name); + } else AMREX_ASSERT_WITH_MESSAGE(false, "only available opacity_method is 0(do not collisions) and 1(do collisions)"); } }; diff --git a/Source/ReadEosTable.cpp b/Source/ReadEosTable.cpp new file mode 100644 index 00000000..eba72b95 --- /dev/null +++ b/Source/ReadEosTable.cpp @@ -0,0 +1,325 @@ +#include + +#include + +#define H5_USE_16_API 1 +#include "hdf5.h" + +#include "EosTable.H" + +// mini NoMPI +#define HAVE_CAPABILITY_MPI //FIXME: This should be defined only when USE_MPI = TRUE +#ifdef HAVE_CAPABILITY_MPI +#include +#define BCAST(buffer, size) MPI_Bcast(buffer, size, MPI_BYTE, my_reader_process, MPI_COMM_WORLD) +#else +#define BCAST(buffer, size) do { /* do nothing */ } while(0) +#endif + +// Catch HDF5 errors +#define HDF5_ERROR(fn_call) \ + if(doIO) { \ + int _error_code = fn_call; \ + if (_error_code < 0) { \ + AMREX_ASSERT_WITH_MESSAGE(false, "HDF5 call failed"); \ + } \ + } + +using namespace amrex; + +static int file_is_readable(std::string filename); +static int file_is_readable(std::string filename) +{ + FILE* fp = NULL; + fp = fopen(filename.c_str(), "r"); + if(fp != NULL) + { + fclose(fp); + return 1; + } + return 0; +} + +namespace nuc_eos_private { + double *alltables; + double *epstable; + double *logrho; + double *logtemp; + double *yes; + double *helperVarsReal; + int *helperVarsInt; +} + +//TODO: Pass the /path/to/table here in the function argument +void ReadEosTable(const std::string nuceos_table_name) { + using namespace nuc_eos_private; + + //std::string nuceos_table_name = "/home/sshanka/000_UTK_projects/Emu/Exec/SFHo.h5"; + amrex::Print() << "(ReadEosTable.cpp) Using table: " << nuceos_table_name << std::endl; + + //TODO: + int my_reader_process = 0; //reader_process; + /*if (my_reader_process < 0 || my_reader_process >= CCTK_nProcs(cctkGH)) + { + CCTK_VWarn(CCTK_WARN_COMPLAIN, __LINE__, __FILE__, CCTK_THORNSTRING, + "Requested IO process %d out of range. Reverting to process 0.", my_reader_process); + my_reader_process = 0; + }*/ + + const int read_table_on_single_process = 1; + //const int doIO = !read_table_on_single_process || CCTK_MyProc(cctkGH) == my_reader_process; //TODO: + const int doIO = 1; + + hid_t file; + if (doIO && !file_is_readable(nuceos_table_name)) { + AMREX_ASSERT_WITH_MESSAGE(false, "Could not read nuceos_table_name"); + } + + HDF5_ERROR(file = H5Fopen(nuceos_table_name.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT)); + +// Use these two defines to easily read in a lot of variables in the same way +// The first reads in one variable of a given type completely +#define READ_BCAST_EOS_HDF5(NAME,VAR,TYPE,MEM,NELEMS) \ + do { \ + hid_t dataset; \ + HDF5_ERROR(dataset = H5Dopen(file, NAME)); \ + HDF5_ERROR(H5Dread(dataset, TYPE, MEM, H5S_ALL, H5P_DEFAULT, VAR)); \ + if (read_table_on_single_process) \ + BCAST (VAR, sizeof(*(VAR))*(NELEMS)); \ + HDF5_ERROR(H5Dclose(dataset)); \ + } while (0) +// The second reads a given variable into a hyperslab of the alltables_temp array +#define READ_BCAST_EOSTABLE_HDF5(NAME,OFF,DIMS) \ + do { \ + READ_BCAST_EOS_HDF5(NAME,&alltables_temp[(OFF)*(DIMS)[1]],H5T_NATIVE_DOUBLE,H5S_ALL,(DIMS)[1]); \ + } while (0) + + int nrho_; + int ntemp_; + int nye_; + // Read size of tables + READ_BCAST_EOS_HDF5("pointsrho", &nrho_, H5T_NATIVE_INT, H5S_ALL, 1); + READ_BCAST_EOS_HDF5("pointstemp", &ntemp_, H5T_NATIVE_INT, H5S_ALL, 1); + READ_BCAST_EOS_HDF5("pointsye", &nye_, H5T_NATIVE_INT, H5S_ALL, 1); + + printf("(ReadEosTable.cpp) nrho = %d, ntemp = %d, nye = %d\n", nrho_, ntemp_, nye_); + + //Allocate managed memory arena on unified memory + ManagedArenaAllocator myManagedArena; + ManagedArenaAllocator myManagedArena_Int; + + // Allocate memory for tables + double *alltables_temp; + if (!(alltables_temp = myManagedArena.allocate(nrho_ * ntemp_ * nye_ * NTABLES) )) { + printf("(ReadEosTable.cpp) Cannot allocate memory for EOS table"); + assert(0); + } + if (!(logrho = myManagedArena.allocate(nrho_) )) { + printf("(ReadEosTable.cpp) Cannot allocate memory for EOS table"); + assert(0); + } + if (!(logtemp = myManagedArena.allocate(ntemp_) )) { + printf("(ReadEosTable.cpp) Cannot allocate memory for EOS table"); + assert(0); + } + if (!(yes = myManagedArena.allocate(nye_) )) { + printf("(ReadEosTable.cpp) Cannot allocate memory for EOS table"); + assert(0); + } + + // Prepare HDF5 to read hyperslabs into alltables_temp + hsize_t table_dims[2] = {NTABLES, (hsize_t)nrho_ * ntemp_ * nye_}; + hsize_t var3[2] = { 1, (hsize_t)nrho_ * ntemp_ * nye_}; + hid_t mem3 = H5Screate_simple(2, table_dims, NULL); + + // Read alltables_temp + READ_BCAST_EOSTABLE_HDF5("logpress", 0, table_dims); + READ_BCAST_EOSTABLE_HDF5("logenergy", 1, table_dims); + READ_BCAST_EOSTABLE_HDF5("entropy", 2, table_dims); + READ_BCAST_EOSTABLE_HDF5("munu", 3, table_dims); + READ_BCAST_EOSTABLE_HDF5("cs2", 4, table_dims); + READ_BCAST_EOSTABLE_HDF5("dedt", 5, table_dims); + READ_BCAST_EOSTABLE_HDF5("dpdrhoe", 6, table_dims); + READ_BCAST_EOSTABLE_HDF5("dpderho", 7, table_dims); + // chemical potentials + READ_BCAST_EOSTABLE_HDF5("muhat", 8, table_dims); + READ_BCAST_EOSTABLE_HDF5("mu_e", 9, table_dims); + READ_BCAST_EOSTABLE_HDF5("mu_p", 10, table_dims); + READ_BCAST_EOSTABLE_HDF5("mu_n", 11, table_dims); + // compositions + READ_BCAST_EOSTABLE_HDF5("Xa", 12, table_dims); + READ_BCAST_EOSTABLE_HDF5("Xh", 13, table_dims); + READ_BCAST_EOSTABLE_HDF5("Xn", 14, table_dims); + READ_BCAST_EOSTABLE_HDF5("Xp", 15, table_dims); + // average nucleus + READ_BCAST_EOSTABLE_HDF5("Abar", 16, table_dims); + READ_BCAST_EOSTABLE_HDF5("Zbar", 17, table_dims); + // Gamma + READ_BCAST_EOSTABLE_HDF5("gamma", 18, table_dims); + + double energy_shift_; + // Read additional tables and variables + READ_BCAST_EOS_HDF5("logrho", logrho, H5T_NATIVE_DOUBLE, H5S_ALL, nrho_); + READ_BCAST_EOS_HDF5("logtemp", logtemp, H5T_NATIVE_DOUBLE, H5S_ALL, ntemp_); + READ_BCAST_EOS_HDF5("ye", yes, H5T_NATIVE_DOUBLE, H5S_ALL, nye_); + READ_BCAST_EOS_HDF5("energy_shift", &energy_shift_, H5T_NATIVE_DOUBLE, H5S_ALL, 1); + + HDF5_ERROR(H5Sclose(mem3)); + HDF5_ERROR(H5Fclose(file)); + + + // change ordering of alltables array so that + // the table kind is the fastest changing index + if (!(alltables = myManagedArena.allocate(nrho_ * ntemp_ * nye_ * NTABLES) )) { + printf("(ReadEosTable.cpp) Cannot allocate memory for EOS table"); + assert(0); + } + + for(int iv = 0;iv epsmax) && (epstable[i] < 1.0e150)){ + epsmax = epstable[i]; + } + if (epstable[i] < epsmin){ + epsmin = epstable[i]; + } + } + + //TODO: Is it correct to subtract energy_shift here? + EOSVAR(eos_epsmin) = epsmin - energy_shift_; + EOSVAR(eos_epsmax) = epsmax - energy_shift_; + + printf("(ReadEosTable.cpp) EOS:rhomin = %.5e g/cm^3\n", EOSVAR(eos_rhomin)); + printf("(ReadEosTable.cpp) EOS:rhomax = %.5e g/cm^3\n", EOSVAR(eos_rhomax)); + printf("(ReadEosTable.cpp) EOS:tempmin = %.4f MeV\n", EOSVAR(eos_tempmin)); + printf("(ReadEosTable.cpp) EOS:tempmax = %.4f MeV\n", EOSVAR(eos_tempmax)); + printf("(ReadEosTable.cpp) EOS:yemin = %.4f\n", EOSVAR(eos_yemin)); + printf("(ReadEosTable.cpp) EOS:yemax = %.4f\n", EOSVAR(eos_yemax)); + + printf("(ReadEosTable.cpp) Finished reading EoS table!\n"); + +} // ReadEOSTable + + + diff --git a/Source/ReadInput_RhoTempYe.H b/Source/ReadInput_RhoTempYe.H new file mode 100644 index 00000000..25414e10 --- /dev/null +++ b/Source/ReadInput_RhoTempYe.H @@ -0,0 +1,14 @@ +#ifndef READINPUT_RHOTEMPYE_H +#define READINPUT_RHOTEMPYE_H + +#include +#include +#include +#include +#include +#include + +void set_rho_T_Ye(MultiFab& state, const Geometry& geom); + + +#endif diff --git a/Source/ReadInput_RhoTempYe.cpp b/Source/ReadInput_RhoTempYe.cpp new file mode 100644 index 00000000..413a5aa4 --- /dev/null +++ b/Source/ReadInput_RhoTempYe.cpp @@ -0,0 +1,50 @@ +#include "Evolve.H" +#include "Constants.H" + +#include + + +void set_rho_T_Ye(MultiFab& state, const Geometry& geom) +{ + // Create an alias of the MultiFab so set_rho_T_Ye only sets rho, T and Ye. + int start_comp = GIdx::rho; + int num_comps = 3; //We only want to set GIdx::rho, GIdx::T and GIdx::Ye + MultiFab rho_T_ye_state(state, amrex::make_alias, start_comp, num_comps); + + amrex::GpuArray dx = geom.CellSizeArray(); + //const auto plo = geom.ProbLoArray(); + //const auto dxi = geom.InvCellSizeArray(); + //const Real inv_cell_volume = dxi[0]*dxi[1]*dxi[2]; + + //const int shape_factor_order_x = geom.Domain().length(0) > 1 ? SHAPE_FACTOR_ORDER : 0; + //const int shape_factor_order_y = geom.Domain().length(1) > 1 ? SHAPE_FACTOR_ORDER : 0; + //const int shape_factor_order_z = geom.Domain().length(2) > 1 ? SHAPE_FACTOR_ORDER : 0; + + //always access mf comp index as (GIdx::rho - start_comp) + //Example: Amrex tutorials -> ExampleCodes/MPMD/Case-2/main.cpp. + + for(amrex::MFIter mfi(rho_T_ye_state); mfi.isValid(); ++mfi){ + const amrex::Box& bx = mfi.validbox(); + const amrex::Array4& mf_array = rho_T_ye_state.array(mfi); + + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k){ + + //x, y and z are the coordinates. + //This is not really needed. Cook up some assert statement to make sure we are at the same (x, y, z) that the table is vlaue is referring to. + amrex::Real x = (i+0.5) * dx[0]; + amrex::Real y = (j+0.5) * dx[1]; + amrex::Real z = (k+0.5) * dx[2]; + + //printf("Inside MFIter: x=%f, y=%f, z=%f\n", x, y, z); + + //TODO: Find the global (i, j, k) from the amrex domain. call them (ig, jg, kg). + //TODO: Then get the values from GPU-array for (ig, jg, kg) and set them to corresponding MultiFabs here. + mf_array(i, j, k, GIdx::rho - start_comp) = -404.0; //FIXME: + mf_array(i, j, k, GIdx::T - start_comp) = -404.0; //FIXME: + mf_array(i, j, k, GIdx::Ye - start_comp) = -404.0; //FIXME: + }); + } + +} + + diff --git a/Source/ReadNuLibTable.cpp b/Source/ReadNuLibTable.cpp new file mode 100644 index 00000000..6947ad3f --- /dev/null +++ b/Source/ReadNuLibTable.cpp @@ -0,0 +1,304 @@ +#include + +#include + +#define H5_USE_16_API 1 +#include "hdf5.h" + +#include "NuLibTable.H" + +// mini NoMPI +#define HAVE_CAPABILITY_MPI //FIXME: This should be defined only when USE_MPI = TRUE +#ifdef HAVE_CAPABILITY_MPI +#include +#define BCAST(buffer, size) MPI_Bcast(buffer, size, MPI_BYTE, my_reader_process, MPI_COMM_WORLD) +#else +#define BCAST(buffer, size) do { /* do nothing */ } while(0) +#endif + +// Catch HDF5 errors +#define HDF5_ERROR(fn_call) \ + if(doIO) { \ + int _error_code = fn_call; \ + if (_error_code < 0) { \ + AMREX_ASSERT_WITH_MESSAGE(false, "HDF5 call failed"); \ + } \ + } + +using namespace amrex; + +static int file_is_readable(std::string filename); +static int file_is_readable(std::string filename) +{ + FILE* fp = NULL; + fp = fopen(filename.c_str(), "r"); + if(fp != NULL) + { + fclose(fp); + return 1; + } + return 0; +} + +namespace nulib_private { + double *alltables_nulib; + //double *epstable; + double *logrho_nulib; + double *logtemp_nulib; + double *yes_nulib; + double *species_nulib; //TODO: Get rid of this? + double *group_nulib; + double *helperVarsReal_nulib; + int *helperVarsInt_nulib; +} + +//TODO: Pass the /path/to/table here in the function argument +void ReadNuLibTable(const std::string nulib_table_name) { + using namespace nulib_private; + + //std::string nulib_table_name = "/mnt/scratch/tables/NuLib/NuLib_SFHo.h5"; + amrex::Print() << "(ReadNuLibTable.cpp) Using table: " << nulib_table_name << std::endl; + + //TODO: + int my_reader_process = 0; //reader_process; + /*if (my_reader_process < 0 || my_reader_process >= CCTK_nProcs(cctkGH)) + { + CCTK_VWarn(CCTK_WARN_COMPLAIN, __LINE__, __FILE__, CCTK_THORNSTRING, + "Requested IO process %d out of range. Reverting to process 0.", my_reader_process); + my_reader_process = 0; + }*/ + + const int read_table_on_single_process = 1; + //const int doIO = !read_table_on_single_process || CCTK_MyProc(cctkGH) == my_reader_process; //TODO: + const int doIO = 1; + + hid_t file; + if (doIO && !file_is_readable(nulib_table_name)) { + AMREX_ASSERT_WITH_MESSAGE(false, "Could not read nulib_table_name"); + } + + HDF5_ERROR(file = H5Fopen(nulib_table_name.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT)); + +// Use these two defines to easily read in a lot of variables in the same way +// The first reads in one variable of a given type completely +#define READ_BCAST_EOS_HDF5(NAME,VAR,TYPE,MEM,NELEMS) \ + do { \ + hid_t dataset; \ + HDF5_ERROR(dataset = H5Dopen(file, NAME)); \ + HDF5_ERROR(H5Dread(dataset, TYPE, MEM, H5S_ALL, H5P_DEFAULT, VAR)); \ + if (read_table_on_single_process) \ + BCAST (VAR, sizeof(*(VAR))*(NELEMS)); \ + HDF5_ERROR(H5Dclose(dataset)); \ + } while (0) +// The second reads a given variable into a hyperslab of the alltables_temp array +#define READ_BCAST_EOSTABLE_HDF5(NAME,OFF,DIMS) \ + do { \ + READ_BCAST_EOS_HDF5(NAME,&alltables_temp[(OFF)*(DIMS)[1]],H5T_NATIVE_DOUBLE,H5S_ALL,(DIMS)[1]); \ + } while (0) + + int nrho_; + int ntemp_; + int nye_; + int nspecies_; + int ngroup_; + + // Read size of tables + READ_BCAST_EOS_HDF5("nrho", &nrho_, H5T_NATIVE_INT, H5S_ALL, 1); + READ_BCAST_EOS_HDF5("ntemp", &ntemp_, H5T_NATIVE_INT, H5S_ALL, 1); + READ_BCAST_EOS_HDF5("nye", &nye_, H5T_NATIVE_INT, H5S_ALL, 1); + READ_BCAST_EOS_HDF5("number_species", &nspecies_, H5T_NATIVE_INT, H5S_ALL, 1); + READ_BCAST_EOS_HDF5("number_groups", &ngroup_, H5T_NATIVE_INT, H5S_ALL, 1); + + assert(nspecies_ == 2); //For now, the code only works when NuLib table has (e,a,x) values. + + printf("(ReadNuLibTable.cpp) nrho = %d, ntemp = %d, nye = %d, nspecies=%d, ngroup=%d\n", nrho_, ntemp_, nye_, nspecies_, ngroup_); + + //Allocate managed memory arena on unified memory + ManagedArenaAllocator myManagedArena; + ManagedArenaAllocator myManagedArena_Int; + + // Allocate memory for tables + double *alltables_temp; + if (!(alltables_temp = myManagedArena.allocate(nrho_ * ntemp_ * nye_ * nspecies_ * ngroup_ * NTABLES_NULIB) )) { + printf("(ReadNuLibTable.cpp) Cannot allocate memory for NuLib table"); + assert(0); + } + if (!(logrho_nulib = myManagedArena.allocate(nrho_) )) { + printf("(ReadNuLibTable.cpp) Cannot allocate memory for NuLib table"); + assert(0); + } + if (!(logtemp_nulib = myManagedArena.allocate(ntemp_) )) { + printf("(ReadNuLibTable.cpp) Cannot allocate memory for NuLib table"); + assert(0); + } + if (!(yes_nulib = myManagedArena.allocate(nye_) )) { + printf("(ReadNuLibTable.cpp) Cannot allocate memory for NuLib table"); + assert(0); + } + if (!(species_nulib = myManagedArena.allocate(nspecies_) )) { + printf("(ReadNuLibTable.cpp) Cannot allocate memory for NuLib table"); + assert(0); + } + if (!(group_nulib = myManagedArena.allocate(ngroup_) )) { + printf("(ReadNuLibTable.cpp) Cannot allocate memory for NuLib table"); + assert(0); + } + + //Allocate memory for energy bin determination. + double *energy_bottom; + double *energy_top; + if (!(energy_bottom = myManagedArena.allocate(ngroup_) )) { + printf("(ReadNuLibTable.cpp) Cannot allocate memory for NuLib table"); + assert(0); + } + if (!(energy_top = myManagedArena.allocate(ngroup_) )) { + printf("(ReadNuLibTable.cpp) Cannot allocate memory for NuLib table"); + assert(0); + } + + // Prepare HDF5 to read hyperslabs into alltables_temp + hsize_t table_dims[2] = {NTABLES_NULIB, (hsize_t)nrho_ * ntemp_ * nye_ * nspecies_ * ngroup_}; + //hsize_t var3[2] = { 1, (hsize_t)nrho_ * ntemp_ * nye_ * nspecies_ * ngroup_}; + hid_t mem3 = H5Screate_simple(2, table_dims, NULL); + + // Read alltables_temp + READ_BCAST_EOSTABLE_HDF5("absorption_opacity", 0, table_dims); + READ_BCAST_EOSTABLE_HDF5("scattering_opacity", 1, table_dims); + + // Read additional tables and variables + //This is not log yet. + READ_BCAST_EOS_HDF5("rho_points", logrho_nulib, H5T_NATIVE_DOUBLE, H5S_ALL, nrho_); + READ_BCAST_EOS_HDF5("temp_points", logtemp_nulib, H5T_NATIVE_DOUBLE, H5S_ALL, ntemp_); + READ_BCAST_EOS_HDF5("ye_points", yes_nulib, H5T_NATIVE_DOUBLE, H5S_ALL, nye_); + READ_BCAST_EOS_HDF5("neutrino_energies", group_nulib, H5T_NATIVE_DOUBLE, H5S_ALL, ngroup_); + + READ_BCAST_EOS_HDF5("bin_bottom", energy_bottom, H5T_NATIVE_DOUBLE, H5S_ALL, ngroup_); + READ_BCAST_EOS_HDF5("bin_top", energy_top, H5T_NATIVE_DOUBLE, H5S_ALL, ngroup_); + + HDF5_ERROR(H5Sclose(mem3)); + HDF5_ERROR(H5Fclose(file)); + + // change ordering of alltables_nulib array so that + // the table kind is the fastest changing index + if (!(alltables_nulib = myManagedArena.allocate(nrho_ * ntemp_ * nye_ * nspecies_ * ngroup_ * NTABLES_NULIB) )) { + printf("(ReadNuLibTable.cpp) Cannot allocate memory for NuLib table"); + assert(0); + } + + for(int iv = 0;iv= energy_bottom[i] && given_energy <= energy_top[i]){ + idx_group_ = i; + break; + } + } + + printf("Given neutrino energy = %f, selected bin index = %d\n", given_energy, idx_group); + myManagedArena.deallocate(energy_bottom, ngroup_); + myManagedArena.deallocate(energy_top, ngroup_); + //---------------------------------------------------------------------------------------------- + + // convert any other quantities to log. + /*for(int i=0;i is_periodic(AMREX_SPACEDIM, 1); + + //The BC will be set using parameter file. + //Option 0: use periodic BC + //Option 1: create particles at boundary. + + //FIXME: FIXME: Define this in parameter file. + const int BC_type = 0; //0=periodic, 1=outer. + + int BC_type_val; + enum BC_type_enum {PERIODIC, OUTER}; + + if (BC_type == 0){ + BC_type_val = BC_type_enum::PERIODIC; //use periodic BC + } else if (BC_type == 1){ + BC_type_val = BC_type_enum::OUTER; //use outer BC + } else { + amrex::Abort("BC_type is incorrect."); + } + + int periodic_flag; + if (BC_type_val == BC_type_enum::PERIODIC){ + //1=yes, use periodic + periodic_flag = 1; + } else if (BC_type_val == BC_type_enum::OUTER){ + //2=no, do not use periodic. + periodic_flag = 0; + } else { + amrex::Abort("BC_type is incorrect."); + } + + Vector is_periodic(AMREX_SPACEDIM, periodic_flag); + Vector domain_lo_bc_types(AMREX_SPACEDIM, BCType::int_dir); Vector domain_hi_bc_types(AMREX_SPACEDIM, BCType::int_dir); + //Vector domain_lo_bc_types(AMREX_SPACEDIM, BCType::foextrap); + //Vector domain_hi_bc_types(AMREX_SPACEDIM, BCType::foextrap); + // Define the index space of the domain @@ -72,24 +107,48 @@ void evolve_flavor(const TestParams* parms) const IntVect ngrow(1 + (1+shape_factor_order_vec)/2); for(int i=0; incell[i] >= ngrow[i]); + //printf("ngrow = [%d, %d, %d] \n", ngrow[0], ngrow[1], ngrow[2]); + // We want 1 component (this is one real scalar field on the domain) const int ncomp = GIdx::ncomp; // Create a MultiFab to hold our grid state data and initialize to 0.0 MultiFab state(ba, dm, ncomp, ngrow); + //FIXME: FIXME: Define this in parameter file. + const int read_rho_T_Ye_from_table = 0; + // initialize with NaNs ... state.setVal(0.0); - state.setVal(parms->rho_in,GIdx::rho,1); // g/ccm - state.setVal(parms->Ye_in,GIdx::Ye,1); - state.setVal(parms->kT_in,GIdx::T,1); // erg - state.FillBoundary(geom.periodicity()); + //If reading from table, call function "set_rho_T_Ye". + //Else set rho, T and Ye to constant value throughout the grid using values from parameter file. + if (read_rho_T_Ye_from_table){ + set_rho_T_Ye(state, geom); + } else { + state.setVal(parms->rho_in,GIdx::rho,1); // g/ccm + state.setVal(parms->Ye_in,GIdx::Ye,1); + state.setVal(parms->kT_in,GIdx::T,1); // erg + } + + state.FillBoundary(geom.periodicity()); + // initialize the grid variable names GIdx::Initialize(); + //We only need HDF5 tables if IMFP_method is 2. + if(parms->IMFP_method==2){ + // read the EoS table + amrex::Print() << "Reading EoS table... " << std::endl; + ReadEosTable(parms->nuceos_table_name); + + // read the NuLib table + amrex::Print() << "Reading NuLib table... " << std::endl; + ReadNuLibTable(parms->nulib_table_name); + } + // Initialize particles on the domain - amrex::Print() << "Initializing particles... "; + amrex::Print() << "Initializing particles... " << std::endl; // We store old-time and new-time data FlavoredNeutrinoContainer neutrinos_old(geom, dm, ba); @@ -162,6 +221,22 @@ void evolve_flavor(const TestParams* parms) // Use the latest-time neutrino data auto& neutrinos = neutrinos_new; + const Real current_dt = integrator.get_timestep(); //FIXME: FIXME: Pass this to neutrinos.CreateParticlesAtBoundary. + + //FIXME: Think carefully where to call this function. + //Create particles at outer boundary + if (BC_type_val == BC_type_enum::OUTER){ + neutrinos.CreateParticlesAtBoundary(parms, current_dt); + neutrinos.CreateParticlesAtBoundary(parms, current_dt); + neutrinos.CreateParticlesAtBoundary(parms, current_dt); + neutrinos.CreateParticlesAtBoundary(parms, current_dt); + neutrinos.CreateParticlesAtBoundary(parms, current_dt); + neutrinos.CreateParticlesAtBoundary(parms, current_dt); + } + + //Create particles at inner boundary + //TODO: This needs to be implemented. + // Update the new time particle locations in the domain with their // integrated coordinates. neutrinos.SyncLocation(Sync::CoordinateToPosition); @@ -177,7 +252,9 @@ void evolve_flavor(const TestParams* parms) const int step = integrator.get_step_number(); const Real time = integrator.get_time(); + printf("Writing reduced data to file... \n"); rd.WriteReducedData0D(geom, state, neutrinos, time, step+1); + printf("Done. \n"); run_fom += neutrinos.TotalNumberOfParticles(); @@ -195,8 +272,11 @@ void evolve_flavor(const TestParams* parms) // Note: this won't be the same as the new-time grid data // because the last deposit_to_mesh call was at either the old time (forward Euler) // or the final RK stage, if using Runge-Kutta. + printf("Setting next timestep... \n"); const Real dt = compute_dt(geom, state, neutrinos, parms); integrator.set_timestep(dt); + //printf("current_dt = %g, dt = %g \n", current_dt, dt); + printf("Done. \n"); }; // Attach our RHS and post timestep hooks to the integrator diff --git a/makefiles/GNUmakefile_jenkins_HDF5_CUDA b/makefiles/GNUmakefile_jenkins_HDF5_CUDA new file mode 100644 index 00000000..72f9e71a --- /dev/null +++ b/makefiles/GNUmakefile_jenkins_HDF5_CUDA @@ -0,0 +1,20 @@ +NUM_FLAVORS = 2 +SHAPE_FACTOR_ORDER = 2 +NUM_MOMENTS = 3 + +COMP = gnu + +DEBUG = FALSE + +USE_MPI = TRUE +USE_OMP = FALSE +USE_ACC = FALSE +USE_CUDA = TRUE +AMREX_CUDA_ARCH=60 + +USE_HDF5=TRUE +HDF5_HOME=/usr/lib/x86_64-linux-gnu/hdf5/openmpi + +EMU_HOME = .. +AMREX_HOME = ../submodules/amrex +include ../Make.Emu diff --git a/sample_inputs/inputs_1d_fiducial b/sample_inputs/inputs_1d_fiducial index 0557d522..e15347ca 100644 --- a/sample_inputs/inputs_1d_fiducial +++ b/sample_inputs/inputs_1d_fiducial @@ -85,4 +85,4 @@ deltaCP_degrees = 222 ################# # opacity stuff # ################# -IMFP_method = 0 \ No newline at end of file +IMFP_method = 0 diff --git a/sample_inputs/inputs_bipolar_test b/sample_inputs/inputs_bipolar_test index c36df0fb..9f12b3e1 100644 --- a/sample_inputs/inputs_bipolar_test +++ b/sample_inputs/inputs_bipolar_test @@ -86,4 +86,4 @@ deltaCP_degrees = 222 ################# # opacity stuff # ################# -IMFP_method = 0 \ No newline at end of file +IMFP_method = 0 diff --git a/sample_inputs/inputs_fast_flavor b/sample_inputs/inputs_fast_flavor index c67c0043..1c27ba8f 100644 --- a/sample_inputs/inputs_fast_flavor +++ b/sample_inputs/inputs_fast_flavor @@ -86,4 +86,4 @@ deltaCP_degrees = 222 ################# # opacity stuff # ################# -IMFP_method = 0 \ No newline at end of file +IMFP_method = 0 diff --git a/sample_inputs/inputs_fast_flavor_nonzerok b/sample_inputs/inputs_fast_flavor_nonzerok index 808cd76e..863fe0a4 100644 --- a/sample_inputs/inputs_fast_flavor_nonzerok +++ b/sample_inputs/inputs_fast_flavor_nonzerok @@ -86,4 +86,4 @@ deltaCP_degrees = 222 ################# # opacity stuff # ################# -IMFP_method = 0 \ No newline at end of file +IMFP_method = 0 diff --git a/sample_inputs/inputs_msw_test b/sample_inputs/inputs_msw_test index 163af873..80a3bbea 100644 --- a/sample_inputs/inputs_msw_test +++ b/sample_inputs/inputs_msw_test @@ -86,4 +86,4 @@ deltaCP_degrees = 222 ################# # opacity stuff # ################# -IMFP_method = 0 \ No newline at end of file +IMFP_method = 0