diff --git a/demos/cosmology/CMakeLists.txt b/demos/cosmology/CMakeLists.txt index a6ec824545..e98e6fd13c 100644 --- a/demos/cosmology/CMakeLists.txt +++ b/demos/cosmology/CMakeLists.txt @@ -5,6 +5,13 @@ # ----------------------------------------------------------------------------- message(STATUS "Configuring demos/cosmology/") -add_executable(StructureFormation StructureFormation.cpp) +add_executable(StructureFormation + StructureFormation.cpp + mc-4-Initializer/Cosmology.cpp + mc-4-Initializer/DataBase.cpp + mc-4-Initializer/InputParser.cpp +) + +target_compile_definitions(StructureFormation PRIVATE USENAMESPACE) target_link_libraries(StructureFormation PRIVATE IPPL::ippl) diff --git a/demos/cosmology/GravityLoadBalancer.hpp b/demos/cosmology/GravityLoadBalancer.hpp index a2543315cd..fa51740753 100644 --- a/demos/cosmology/GravityLoadBalancer.hpp +++ b/demos/cosmology/GravityLoadBalancer.hpp @@ -149,26 +149,22 @@ class LoadBalancer { if (ippl::Comm->size() < 2) { return false; } - if (std::strcmp(TestName, "UniformPlasmaTest") == 0) { - return (nstep % loadbalancefreq_m == 0); - } else { - int local = 0; - std::vector res(ippl::Comm->size()); - double equalPart = static_cast(totalP) / ippl::Comm->size(); - double dev = std::abs(static_cast(pc_m->getLocalNum()) - equalPart) / totalP; - if (dev > loadbalancethreshold_m) { - local = 1; - } - MPI_Allgather(&local, 1, MPI_INT, res.data(), 1, MPI_INT, - ippl::Comm->getCommunicator()); - for (unsigned int i = 0; i < res.size(); i++) { - if (res[i] == 1) { - return true; - } + int local = 0; + std::vector res(ippl::Comm->size()); + double equalPart = static_cast(totalP) / ippl::Comm->size(); + double dev = std::abs(static_cast(pc_m->getLocalNum()) - equalPart) / totalP; + if (dev > loadbalancethreshold_m) { + local = 1; + } + MPI_Allgather(&local, 1, MPI_INT, res.data(), 1, MPI_INT, ippl::Comm->getCommunicator()); + + for (unsigned int i = 0; i < res.size(); i++) { + if (res[i] == 1) { + return true; } - return false; } + return false; } private: @@ -182,4 +178,4 @@ class LoadBalancer { ORB orb; ///< ORB for domain partitioning }; -#endif // IPPL_LOAD_BALANCER_H \ No newline at end of file +#endif // IPPL_LOAD_BALANCER_H diff --git a/demos/cosmology/GravityManager.h b/demos/cosmology/GravityManager.h index 4dd9995f86..3b9aa2defa 100644 --- a/demos/cosmology/GravityManager.h +++ b/demos/cosmology/GravityManager.h @@ -12,6 +12,7 @@ #include "Random/InverseTransformSampling.h" #include "Random/NormalDistribution.h" #include "Random/Randn.h" +#include "mc-4-Initializer/InputParser.h" using view_type = typename ippl::detail::ViewType, 1>::view_type; @@ -42,11 +43,13 @@ class GravityManager : public ippl::PicManager * @param lbt_ Load balance threshold. * @param solver_ Solver type. * @param stepMethod_ Time stepping method type. + * @param par_ Inputfile parser */ GravityManager(size_type totalP_, int nt_, Vector_t& nr_, double lbt_, - std::string& solver_, std::string& stepMethod_) + std::string& solver_, std::string& stepMethod_, initializer::InputParser par_) : ippl::PicManager, FieldContainer, LoadBalancer>() + , parser_m(par_) , totalP_m(totalP_) , nt_m(nt_) , nr_m(nr_) @@ -65,6 +68,11 @@ class GravityManager : public ippl::PicManager */ std::string folder; + /** + * @brief Access to the input file with constants and simulation parameters. + */ + initializer::InputParser parser_m; + /** * @brief Get the total number of particles. * @@ -211,8 +219,10 @@ class GravityManager : public ippl::PicManager */ void InitialiseTime() { Inform mes("Inititalise: "); - this->O_m = 0.3; - this->O_L = 0.7; + parser_m.getByName("Omega_m", this->O_m); + parser_m.getByName("Omega_L", this->O_L); + // this->O_m = 0.3; // \todo need to from input file + // this->O_L = 0.7; // \todo need to from input file this->t_L = 2 / (3 * this->Hubble0 * sqrt(this->O_L)); this->a_m = 1 / (1 + this->z_m); this->Dloga = 1. / (this->nt_m) * log((1 + this->z_m) / (1 + this->z_f)); @@ -241,8 +251,8 @@ class GravityManager : public ippl::PicManager * @brief Pre-step method called before each simulation step. */ void pre_step() override { - Inform mes("Pre-step"); - mes << "Done" << endl; + // Inform mes("Pre-step"); + // mes << "Done" << endl; } /** @@ -262,9 +272,9 @@ class GravityManager : public ippl::PicManager // dynamic time step this->dt_m = this->Dloga / this->Hubble_m; - mes << "Finished time step: " << this->it_m << endl; - mes << " time: " << this->time_m << ", timestep: " << this->dt_m << ", z: " << this->z_m - << ", a: " << this->a_m << endl; + mes << "Step: " << this->it_m; + mes << " comological time: " << this->time_m << ", dt: " << this->dt_m + << ", z: " << this->z_m << ", a: " << this->a_m << endl; } // Grid to particle and particle to grid methods @@ -291,10 +301,10 @@ class GravityManager : public ippl::PicManager */ void scatterCIC() { Inform mes("scatter "); - mes << "starting ..." << endl; + this->fcontainer_m->getRho() = 0.0; - ippl::ParticleAttrib* m = &this->pcontainer_m->m; + ippl::ParticleAttrib* m = &this->pcontainer_m->m; typename Base::particle_position_type* R = &this->pcontainer_m->R; Field_t* rho = &this->fcontainer_m->getRho(); Vector_t rmin = rmin_m; @@ -303,7 +313,6 @@ class GravityManager : public ippl::PicManager scatter(*m, *rho, *R); double relError = std::fabs((M_m - (*rho).sum()) / M_m); - mes << "relative error: " << relError << endl; size_type TotalParticles = 0; size_type localParticles = this->pcontainer_m->getLocalNum(); @@ -311,7 +320,8 @@ class GravityManager : public ippl::PicManager ippl::Comm->reduce(localParticles, TotalParticles, 1, std::plus()); if (ippl::Comm->rank() == 0) { - if (TotalParticles != totalP_m || relError > 1e-10) { + if (TotalParticles != totalP_m + || relError > 10. * Kokkos::Experimental::epsilon_v) { mes << "Time step: " << it_m << endl; mes << "Total particles in the sim. " << totalP_m << " " << "after update: " << TotalParticles << endl; @@ -351,8 +361,8 @@ class GravityManager : public ippl::PicManager double Hubble0; ///< Hubble constant today (73.8 km/sec/Mpc). double G; ///< Gravitational constant. [kpc^3/(Msun s^2)] double rho_crit0; ///< Critical density today. [Msun/kpc^3] - double O_m; ///< Matter density parameter. [1] - double O_L; ///< Dark energy density parameter. [1] + float O_m; ///< Matter density parameter. [1] + float O_L; ///< Dark energy density parameter. [1] double t_L; ///< Characteristic time scale. [s] double z_m; ///< Initial redshift. [1] double z_f; ///< Final redshift. [1] @@ -368,4 +378,4 @@ class GravityManager : public ippl::PicManager std::array decomp_m; ///< Decomposition flags for each dimension. [bool] double rhoNorm_m; ///< Normalized density. [1] }; -#endif \ No newline at end of file +#endif diff --git a/demos/cosmology/GravityParticleContainer.hpp b/demos/cosmology/GravityParticleContainer.hpp index e670da5779..778c4b3663 100644 --- a/demos/cosmology/GravityParticleContainer.hpp +++ b/demos/cosmology/GravityParticleContainer.hpp @@ -21,7 +21,7 @@ class ParticleContainer : public ippl::ParticleBase< /** * @brief Mass of the particle. */ - ippl::ParticleAttrib m; + ippl::ParticleAttrib m; /** * @brief Velocity of the particle. @@ -91,4 +91,4 @@ class ParticleContainer : public ippl::ParticleBase< void setBCAllPeriodic() { this->setParticleBC(ippl::BC::PERIODIC); } }; -#endif // IPPL_PARTICLE_CONTAINER_H \ No newline at end of file +#endif // IPPL_PARTICLE_CONTAINER_H diff --git a/demos/cosmology/StructureFormation.cpp b/demos/cosmology/StructureFormation.cpp index 07b71babf4..387dc59b41 100644 --- a/demos/cosmology/StructureFormation.cpp +++ b/demos/cosmology/StructureFormation.cpp @@ -1,34 +1,73 @@ -// StructureFormation Test -// +// StructureFormation.cpp // The StructureFormation simulation is designed to model the formation and evolution of cosmic // structures, such as galaxies and clusters of galaxies, in the universe. It uses initial // conditions, including particle positions and velocities, to simulate the gravitational // interactions and dynamics of a large number of particles over time. The goal is to understand how // initial density fluctuations grow and evolve under the influence of gravity, leading to the -// large-scale structure observed in the universe today. +// large-scale structure observed in the universe today. The old mc-4 initializer based in Zaria's +// old and new initializer is added. // // Usage: // srun ./StructureFormation -// [...] -// --overallocate --info 10 -// path = path to initial conditions folder containing the file Data.csv -// nx = No. cell-centered points in the x-direction -// ny... = No. cell-centered points in the y-, z-, ...-direction -// Np = Total no. of macro-particles in the simulation. Needs to match the IC file -// Data.csv Nt = Number of time steps stype = Field solver type (FFT and CG supported) -// lbthres = Load balancing threshold i.e., lbthres*100 is the maximum load imbalance -// percentage which can be tolerated and beyond which -// particle load balancing occurs. A value of 0.01 is good for many typical -// simulations. -// ovfactor = Over-allocation factor for the buffers used in the communication. Typical -// values are 1.0, 2.0. Value 1.0 means no over-allocation. +// +// --overallocate --info 10 +// inputfile = describes the simulation +// tfFn = transfer function +// outFn = base file name for snapshots +// outDir = base directory for output data +// fsType = Field solver type (FFT and CG supported) +// lbthres = Load balancing threshold i.e., lbthres*100 is the maximum load imbalance +// percentage which can be tolerated and beyond which +// particle load balancing occurs. A value of 0.01 is good for many typical +// simulations. +// integrator= LeapFrog +// ovfactor = Over-allocation factor for the buffers used in the communication. Typical +// values are 1.0, 2.0. Value 1.0 means no over-allocation. // Example: -// srun ./StructureFormation data/lsf_32/ 32 32 32 32768 10 FFT 1.0 LeapFrog --overallocate 1.0 +// srun StructureFormation input.par tf.dat out.data datadir FFT 1.0 LeapFrog --overallocate 1.0 // --info 5 +// +// more stuff and todos at: +// https://docs.google.com/document/d/1PhQ2Wo6QhdeS7U3lN1q2TEaJm2XjD4eXdjGCMKfPac4/edit?usp=sharing + +/* Example of an input file + +np=128 +nt=0 +box_size=64.0 // In Mpc/h +seed=9854373 +z_in=50.0 +z_fi=0.0 + +// Cosmology: +hubble=0.7 +Omega_m=0.25 // Total matter content; 1-this will be DE +Omega_nu=0.0 +Omega_bar=0.04 +Omega_r=0.0 + +Sigma_8=0.8 +n_s=0.9 +w_de=-1.0 +N_nu=3 +nu_pairs=4 +f_NL=0.0 +TFFlag=2 + +// Code stuff: +PrintFormat=0 + +// new +Omega_L=0.7 + +// Read from file (0) or create (1) +ReadInParticles=1 + + +*/ constexpr unsigned Dim = 3; using T = double; -const char* TestName = "StructureFormation"; #include "Ippl.h" @@ -42,12 +81,12 @@ const char* TestName = "StructureFormation"; #include #include -#include "Manager/datatypes.h" - #include "Utility/IpplTimings.h" #include "Manager/PicManager.h" #include "StructureFormationManager.h" +#include "mc-4-Initializer/DataBase.h" +#include "mc-4-Initializer/InputParser.h" using size_type = ippl::detail::size_type; @@ -57,25 +96,42 @@ using Vector_t = ippl::Vector; int main(int argc, char* argv[]) { ippl::initialize(argc, argv); { - Inform msg(TestName); - Inform msg2all(TestName, INFORM_ALL_NODES); + Inform msg("StructureFormation "); + Inform msg2all("StructureFormation ", INFORM_ALL_NODES); static IpplTimings::TimerRef mainTimer = IpplTimings::getTimer("total"); IpplTimings::startTimer(mainTimer); - // Read input parameters, assign them to the corresponding members of manager - int arg = 1; - // initial conditions folder - std::string ic_folder = argv[arg++]; + std::string indatName = argv[1]; + std::string tfName = argv[2]; + std::string outBase = argv[3]; + std::string ic_folder = argv[4]; + + int arg = 5; + initializer::InputParser par(indatName); + initializer::GlobalStuff::instance().GetParameters(par); + + size_t totalP; + if ((totalP = static_cast(initializer::GlobalStuff::instance().ngrid) + * static_cast(initializer::GlobalStuff::instance().ngrid) + * static_cast(initializer::GlobalStuff::instance().ngrid)) + > std::numeric_limits::max()) { + throw std::overflow_error("Index space exceeds size_t capacity"); + } + // Number of gridpoints in each dimension Vector_t nr; for (unsigned d = 0; d < Dim; d++) { - nr[d] = std::atoi(argv[arg++]); + nr[d] = initializer::GlobalStuff::instance().ngrid; } - // Total number of particles - size_type totalP = std::atoll(argv[arg++]); + // Number of time steps - int nt = std::atoi(argv[arg++]); + int nt = 0; + par.getByName("nt", nt); + + int readInParticles = 0; + par.getByName("ReadInParticles", readInParticles); + bool readICs = (readInParticles == 0); // Solver method std::string solver = argv[arg++]; @@ -89,7 +145,8 @@ int main(int argc, char* argv[]) { std::string step_method = argv[arg++]; // Create an instance of a manager for the considered application - StructureFormationManager manager(totalP, nt, nr, lbt, solver, step_method); + StructureFormationManager manager(totalP, nt, nr, lbt, solver, step_method, par, + tfName, readICs); // set initial conditions folder manager.setIC(ic_folder); @@ -97,7 +154,7 @@ int main(int argc, char* argv[]) { // Perform pre-run operations, including creating mesh, particles,... manager.pre_run(); - msg << "Starting iterations ..." << endl; + msg << "Starting iterations ... up to " << manager.getNt() << endl; manager.run(manager.getNt()); diff --git a/demos/cosmology/StructureFormationManager.h b/demos/cosmology/StructureFormationManager.h index 50ea332606..98ae5226f3 100644 --- a/demos/cosmology/StructureFormationManager.h +++ b/demos/cosmology/StructureFormationManager.h @@ -19,9 +19,35 @@ #include "Random/InverseTransformSampling.h" #include "Random/NormalDistribution.h" #include "Random/Randn.h" +#include "mc-4-Initializer/Cosmology.h" +#include "mc-4-Initializer/DataBase.h" +#include "mc-4-Initializer/InputParser.h" + +// #define DOUBLE_INIT_COND + +// #define KOKKOS_PRINT // Kokkos::printf of interesting quantities. Does not work multirank using view_type = typename ippl::detail::ViewType, 1>::view_type; +#ifdef DOUBLE_INIT_COND +typedef ippl::Field, Dim, Mesh_t, Mesh_t::DefaultCentering> + field_type; +typedef ippl::FFT CFFT_type; +typedef Field, Dim> CField_t; +typedef Field RField_t; +#else +typedef ippl::Field, Dim, Mesh_t, Mesh_t::DefaultCentering> + field_type; +typedef ippl::FFT CFFT_type; +typedef Field, Dim> CField_t; +typedef Field RField_t; +#endif + +struct HermitianPkg { + int kx, ky, kz; + double re, im; +}; + /** * @brief Construct a new StructureFormationManager object. * @@ -31,18 +57,38 @@ using view_type = typename ippl::detail::ViewType, 1>: * @param lbt_ Lookback time. * @param solver_ Solver method. * @param stepMethod_ Time stepping method. + * @param par_ the parser to read the input file + * @param tfname_ filename for transfer function + * @param readICs_ read or create initial conditions */ template class StructureFormationManager : public GravityManager { + /// all for the initializer + std::unique_ptr Cfft_m; + CField_t cfield_m; + RField_t Pk_m; + bool readICs_m; + + initializer::CosmoClass cosmo_m; + public: using ParticleContainer_t = ParticleContainer; using FieldContainer_t = FieldContainer; using FieldSolver_t = FieldSolver; using LoadBalancer_t = LoadBalancer; + using index_array_type = typename ippl::RangePolicy::index_array_type; + using index_type = typename ippl::RangePolicy::index_type; + StructureFormationManager(size_type totalP_, int nt_, Vector_t& nr_, double lbt_, - std::string& solver_, std::string& stepMethod_) - : GravityManager(totalP_, nt_, nr_, lbt_, solver_, stepMethod_) {} + std::string& solver_, std::string& stepMethod_, + initializer::InputParser par_, std::string tfname_, bool readICs_) + : GravityManager(totalP_, nt_, nr_, lbt_, solver_, stepMethod_, par_) + , readICs_m(readICs_) { + Inform m("StructureFormationManager "); + m << "totalP = " << totalP_ << endl; + cosmo_m.SetParameters(initializer::GlobalStuff::instance(), tfname_.c_str()); + } /** * @brief Destructor for StructureFormationManager. @@ -53,7 +99,7 @@ class StructureFormationManager : public GravityManager { * @brief Pre-run setup for the simulation. */ void pre_run() override { - Inform mes("Pre Run"); + Inform msg("Pre Run"); if (this->solver_m == "OPEN") { throw IpplException("StructureFormation", @@ -68,23 +114,32 @@ class StructureFormationManager : public GravityManager { this->decomp_m.fill(true); this->Hubble0 = 0.1; // h * km/sec/kpc (h = 0.7, H = 0.07) this->G = 4.30071e04; // kpc km^2 /s^2 / M_Sun e10 - this->z_m = 63; // initial z - this->z_f = 0; // final z + + float zm; + float zf; + this->parser_m.getByName("z_in", zm); // initial z + this->parser_m.getByName("z_fi", zf); // final z + this->z_m = zm; + this->z_f = zf; + this->InitialiseTime(); + float box_size; + this->parser_m.getByName("box_size", box_size); this->rmin_m = 0.0; - this->rmax_m = 50000.0; // kpc/h + this->rmax_m = box_size * 1000.0; // kpc/h + double Vol = std::reduce(this->rmax_m.begin(), this->rmax_m.end(), 1., std::multiplies()); this->M_m = this->rho_crit0 * Vol * this->O_m; // 1e10 M_Sun - mes << "total mass: " << this->M_m << endl; - mes << "mass of a single particle " << this->M_m / this->totalP_m << endl; + msg << "total mass: " << this->M_m << endl; + msg << "mass of a single particle " << this->M_m / this->totalP_m << endl; this->hr_m = this->rmax_m / this->nr_m; this->origin_m = this->rmin_m; this->it_m = 0; - mes << "Discretization:" << endl + msg << "Discretization:" << endl << "nt " << this->nt_m << ", Np = " << this->totalP_m << ", grid = " << this->nr_m << endl; @@ -103,12 +158,43 @@ class StructureFormationManager : public GravityManager { this->solver_m, &this->fcontainer_m->getRho(), &this->fcontainer_m->getF(), &this->fcontainer_m->getPhi())); + msg << "Set FieldSolver done ..." << endl; + this->fsolver_m->initSolver(); + msg << "Init FieldSolver done ..." << endl; + this->setLoadBalancer(std::make_shared( this->lbt_m, this->fcontainer_m, this->pcontainer_m, this->fsolver_m)); - readParticlesDomain(); // defines particle positions, velocities + if (readICs_m) { + msg << "Read in particles ..." << endl; + readParticlesDomain(); // defines particle positions, velocities + msg << "Read particles done" << endl; + } else { + msg << "Create Particles ....." << endl; + ippl::ParameterList fftParams; + fftParams.add("use_heffte_defaults", true); + + /** ADA + make string to annotate pushRegion + */ + + const ippl::NDIndex& lDom = this->fcontainer_m->getFL().getLocalNDIndex(); + std::ostringstream oss; + oss << "CFFT CField Pk nx= " << lDom[0].length(); + + Kokkos::Profiling::pushRegion(oss.str()); + Cfft_m = std::make_unique(this->fcontainer_m->getFL(), fftParams); + cfield_m.initialize(this->fcontainer_m->getMesh(), this->fcontainer_m->getFL()); + Pk_m.initialize(this->fcontainer_m->getMesh(), this->fcontainer_m->getFL()); + msg << "FFT and Pk structures initialized" << endl; + Kokkos::Profiling::popRegion(); + createParticles(); + msg << "Create particles done" << endl; + } + + return; static IpplTimings::TimerRef DummySolveTimer = IpplTimings::getTimer("solveWarmup"); IpplTimings::startTimer(DummySolveTimer); @@ -118,28 +204,419 @@ class StructureFormationManager : public GravityManager { this->fsolver_m->runSolver(); IpplTimings::stopTimer(DummySolveTimer); + this->par2grid(); static IpplTimings::TimerRef SolveTimer = IpplTimings::getTimer("solve"); - IpplTimings::startTimer(SolveTimer); this->fsolver_m->runSolver(); IpplTimings::stopTimer(SolveTimer); this->grid2par(); this->dump(); + msg << "Done creating initial conditions" << endl; + + this->savePositions(0); + } - mes << "Done"; + void initPwrSpec() { + Inform msg("initspec "); + + constexpr auto tpi = 2 * Kokkos::numbers::pi_v; + const double n_s = initializer::GlobalStuff::instance().n_s; + const double sigma8 = initializer::GlobalStuff::instance().Sigma_8; + const double f_NL = initializer::GlobalStuff::instance().f_NL; + const int ngrid = initializer::GlobalStuff::instance().ngrid; + const int nq = ngrid / 2; + const double tpiL = + tpi / initializer::GlobalStuff::instance().box_size; // k0, physical units + + /* TFFlag == 2) Hu-Sugiyama transfer function + const double Omega_m = initializer::GlobalStuff::instance().Omega_m; + const double Omega_bar = initializer::GlobalStuff::instance().Omega_bar; + const double h = initializer::GlobalStuff::instance().Hubble; + const double akh1=pow(46.9*Omega_m*h*h, 0.670)*(1.0+pow(32.1*Omega_m*h*h, -0.532)); + const double akh2=pow(12.0*Omega_m*h*h, 0.424)*(1.0+pow(45.0*Omega_m*h*h, -0.582)); + const double alpha=pow(akh1, -1.0*Omega_bar/Omega_m)*pow(akh2, + pow(-1.0*Omega_bar/Omega_m, 3.0)); const double kh_tmp = Omega_m*h*Kokkos::sqrt(alpha); + */ + // TFFlag == 4 BBKS transfer function + + const float Omega_m = initializer::GlobalStuff::instance().Omega_m; + const float h = initializer::GlobalStuff::instance().Hubble; + const float kh_tmp = Omega_m * h; + + /* Set P(k)=T^2*k^n array. + Linear array, taking care of the mirror symmetry + (reality condition for the density field). + */ + + msg << "Pulled all needed physics quantities" << endl; + msg << "kh_tmp= " << kh_tmp << " n_s= " << n_s << " tpiL= " << tpiL << endl; + + auto pkview = Pk_m.getView(); + auto* FL = &this->fcontainer_m->getFL(); + const ippl::NDIndex& ldom = + FL->getLocalNDIndex(); // local processor domain coordinates + const int nghost = Pk_m.getNghost(); + + ippl::parallel_for( + "Compute Pk_m", ippl::getRangePolicy(pkview, nghost), + KOKKOS_LAMBDA(const index_array_type& idx) { + int i = idx[0] - nghost + ldom[0].first(); + int j = idx[1] - nghost + ldom[1].first(); + int k = idx[2] - nghost + ldom[2].first(); + + int k_i = (i >= nq) ? -(nq - i % nq) : i; + int k_j = (j >= nq) ? -(nq - j % nq) : j; + int k_k = (k >= nq) ? -(nq - k % nq) : k; + + int k2 = k_i * k_i + k_j * k_j + k_k * k_k; + float kk = tpiL * Kokkos::sqrt(k2); + double trans_f = StructureFormationManager::bbks(kk, kh_tmp); + // double trans_f = StructureFormationManager::hu_sugiyama(kk, kh_tmp); + double val = trans_f * trans_f * Kokkos::pow(kk, n_s); + ippl::apply(pkview, idx) = val; + +#ifdef KOKKOS_PRINT + int index = + (k_i) + (k_j)*ldom[0].length() + (k_k)*ldom[0].length() * ldom[1].length(); + + if ((k_i == 0) && (k_j == 0) && (k_k == 0)) + Kokkos::printf("Pk: i \t j \t k \t index \t kk \t transf \t Pk \n"); + Kokkos::printf("Pk: %d \t %d \t %d \t %d \t %f \t %f \t %f \n", k_i, k_j, k_k, + index, kk, trans_f, val); +#endif + }); + + msg << "Pk created using BBKS TFFlag ==4" << endl; + + double s8 = cosmo_m.Sigma_r(8.0, 1.0); + const double norm = sigma8 * sigma8 / (s8 * s8); + s8 = cosmo_m.Sigma_r(8.0, norm); + msg << "sigma_8=" << s8 << ", target was " << sigma8 << endl; + msg << "norm= " << norm << endl; + + // For non-Gaussian initial conditions: P(k)=A*k^n, trasfer function will come later + if (f_NL != 0.0) { + msg << "Non-Gaussian initial conditions, f_NL=" << f_NL << endl; + ippl::parallel_for( + "Norm Pk_m (non-Gaussian)", ippl::getRangePolicy(pkview, nghost), + KOKKOS_LAMBDA(const index_array_type& args) { + int k_i = args[0] - nghost + ldom[0].first(); + int k_j = args[1] - nghost + ldom[1].first(); + int k_k = args[2] - nghost + ldom[2].first(); + if (k_k >= nq) { + k_k = -fmod(ngrid - k_k, ngrid); + } + double kk = tpiL * Kokkos::sqrt(k_i * k_i + k_j * k_j + k_k * k_k); + ippl::apply(pkview, args) *= norm * pow(kk, n_s); + }); + } else { + // For Gaussian initial conditions: P(k)=A*T^2*k^n + msg << "Gaussian initial conditions, f_NL=" << f_NL << endl; + ippl::parallel_for( + "Norm Pk_m (Gaussian)", ippl::getRangePolicy(pkview, nghost), + KOKKOS_LAMBDA(const index_array_type& args) { ippl::apply(pkview, args) *= norm; }); + } + } + + /** + @brief bbks hardcoded TFFLAG==4 + + */ + + KOKKOS_FUNCTION static float bbks(float k, float kh_tmp) { + float qkh, t_f; + + if (k == 0.0) + return (0.0); + qkh = k / kh_tmp; + t_f = Kokkos::log(1.0 + 2.34 * qkh) / (2.34 * qkh) + * Kokkos::pow(1.0 + 3.89 * qkh + Kokkos::pow(16.1 * qkh, 2.0) + + Kokkos::pow(5.46 * qkh, 3.0) + Kokkos::pow(6.71 * qkh, 4.0), + -0.25); + return (t_f); + } + + void LinearZeldoInitMP() { + // After creating the field layout (cfield_m) and determining global grid sizes Nx, Ny, Nz: + { + Inform msg("LinearZeldoInitMP "); + Inform m2a("LinearZeldoInitMP ", INFORM_ALL_NODES); + + typename CField_t::view_type& cfview = cfield_m.getView(); + typename RField_t::view_type& pkview = Pk_m.getView(); + + auto rView = this->pcontainer_m->R.getView(); + auto vView = this->pcontainer_m->V.getView(); + + // auto rView = this->pcontainer_m->R.getHostMirror(); + // auto vView = this->pcontainer_m->V.getHostMirror(); + + const int ngh = cfield_m.getNghost(); + const ippl::NDIndex& lDom = this->fcontainer_m->getFL().getLocalNDIndex(); + + index_type lgridsize = 1; + for (unsigned d = 0; d < Dim; d++) { + lgridsize *= lDom[d].length(); + } + const uint64_t global_seed = 12345ULL; // Shared global seed for reproducibility + + const int Nx = this->nr_m[0]; + const int Ny = this->nr_m[1]; + const int Nz = this->nr_m[2]; + const Vector_t hr = this->hr_m; + const double Lx = this->rmax_m[0]; + const double Ly = this->rmax_m[1]; + const double Lz = this->rmax_m[2]; + + msg << "Lx= " << Lx << " Ly= " << Ly << " Lz= " << Lz << endl; + msg << "Nx= " << Nx << " Ny= " << Ny << " Nz= " << Nz << endl; + msg << "h = " << this->hr_m << endl; + + float d_z; + float ddot; + float z_in = initializer::GlobalStuff::instance().z_in; + cosmo_m.GrowthFactor(z_in, &d_z, &ddot); + msg << "z_in= " << z_in << " d_z= " << d_z << " ddot= " << ddot << endl; + + Kokkos::Random_XorShift64_Pool<> rand_pool(global_seed); + + // Initialize the Fourier density field with Gaussian random modes (Hermitian symmetric) + ippl::parallel_for( + "InitDeltaField", ippl::getRangePolicy(cfview, ngh), + KOKKOS_LAMBDA(const index_array_type& idx) { + const double pi = Kokkos::numbers::pi_v; + // Compute global coordinates (i,j,k) for this local index + int i = idx[0] - ngh + lDom[0].first(); + int j = idx[1] - ngh + lDom[1].first(); + int k = idx[2] - ngh + lDom[2].first(); + + // DC mode (k=0 vector) set to 0 (no DC offset) + if (i == 0 && j == 0 && k == 0) { + ippl::apply(cfview, idx) = Kokkos::complex(0.0, 0.0); + } else { + // Compute the global “negative” indices for Hermitian pair + int i_neg = (i == 0 ? 0 : Nx - i); + int j_neg = (j == 0 ? 0 : Ny - j); + int k_neg = (k == 0 ? 0 : Nz - k); + + // Determine if this index is its own conjugate (self-Hermitian case) + bool self = (i_neg == i && j_neg == j && k_neg == k); + // Determine lexicographically which of (i,j,k) and its negative is smaller + bool is_conjugate = (!self + && (i_neg < i || (i_neg == i && j_neg < j) + || (i_neg == i && j_neg == j && k_neg < k))); + + auto rand_gen = rand_pool.get_state(); + // Generate a uniform random number in [0,1) + double u1 = rand_gen.drand(); + double u2 = rand_gen.drand(); + // Return the state to the pool + rand_pool.free_state(rand_gen); + + // Convert uniforms to Gaussian via Box-Muller + double R = Kokkos::sqrt(-2.0 * Kokkos::log(u1)); + double theta = 2.0 * pi * u2; + double gauss_re = R * Kokkos::cos(theta); // Gaussian(0,1) for real part + double gauss_im = + R * Kokkos::sin(theta); // Gaussian(0,1) for imaginary part + double Pk = ippl::apply(pkview, idx); // power spectrum P(k) + // Set amplitude: for self-conjugate modes use sqrt(Pk), otherwise + // sqrt(Pk/2) + double amp = self ? Kokkos::sqrt(Pk) : Kokkos::sqrt(Pk / 2.0); + double val_re = amp * gauss_re; + double val_im = amp * gauss_im; + if (self) { + // For self-conjugate (Nyquist) modes, enforce the mode is real + val_im = 0.0; + } else if (is_conjugate) { + // If this index is the "conjugate partner" (lexicographically larger), + // flip the imaginary sign + val_im = -val_im; + } +#ifdef KOKKOS_PRINT + Kokkos::printf("rho: %g %g \n", val_re, val_im); +#endif + // Assign the complex value to this local mode + ippl::apply(cfview, idx) = Kokkos::complex(val_re, val_im); + } + }); + + // Store delta(k) for reuse + auto tmpcfield = cfield_m.deepCopy(); + + /* + Now we can delete Pk and allocate the particles + */ + + // 2–4. Loop over displacement components x(0), y(1), z(2) + double xx, yy, zz; + int nq = Nx / 2; + + for (int dim = 0; dim < 3; ++dim) { + switch (dim) { + case 0: + xx = 1.0; + yy = 0.0; + zz = 0.0; + break; + case 1: + xx = 0.0; + yy = 1.0; + zz = 0.0; + break; + case 2: + xx = 0.0; + yy = 0.0; + zz = 1.0; + break; + } + + // Compute displacement component in k-space + ippl::parallel_for( + "force_in_k_space", ippl::getRangePolicy(cfview, ngh), + KOKKOS_LAMBDA(const index_array_type& idx) { + const double pi = Kokkos::numbers::pi_v; + int i = idx[0] - ngh + lDom[0].first(); + int j = idx[1] - ngh + lDom[1].first(); + int k = idx[2] - ngh + lDom[2].first(); + + int k_i = (i >= nq) ? -(nq - i % nq) : i; + int k_j = (j >= nq) ? -(nq - j % nq) : j; + int k_k = (k >= nq) ? -(nq - k % nq) : k; + + int k2 = k_i * k_i + k_j * k_j + k_k * k_k; + double Grad = -2 * pi / Nx * (k_i * xx + k_j * yy + k_k * zz); + double kk = 2 * pi / Nx * sqrt(k2); + double Green = (kk != 0.0) ? -1.0 / (kk * kk) : 0.0; + + Kokkos::complex tmp = ippl::apply(tmpcfield, idx); + Kokkos::complex res(-Grad * Green * tmp.imag(), + Grad * Green * tmp.real()); + ippl::apply(cfview, idx) = res; +#ifdef KOKKOS_PRINT + Kokkos::printf("feval: %i %g %g %g %g %g %g \n", dim, Grad, Green, + tmp.imag(), tmp.real(), res.real(), res.imag()); +#endif + }); + + // Inverse FFT to real space + Cfft_m->transform(ippl::BACKWARD, cfield_m); + + const double scale = Kokkos::pow(initializer::GlobalStuff::instance().ngrid, 1.5); + ippl::parallel_for( + "scale_rho_inv", ippl::getRangePolicy(cfview, ngh), + KOKKOS_LAMBDA(const index_array_type& idx) { + const Kokkos::complex tmp = ippl::apply(cfview, idx); + ippl::apply(cfview, idx) = + Kokkos::complex(tmp.real() / scale, tmp.imag() / scale); +#ifdef KOKKOS_PRINT + Kokkos::printf("rhoinvfft: %i %g %g \n", dim, tmp.real() / scale, + tmp.imag() / scale); +#endif + }); + + ippl::parallel_for( + "set_particle", ippl::getRangePolicy(cfview, ngh), + KOKKOS_LAMBDA(const index_array_type& idx) { + const unsigned int i = idx[0] - ngh + lDom[0].first(); + const unsigned int j = idx[1] - ngh + lDom[1].first(); + const unsigned int k = idx[2] - ngh + lDom[2].first(); + + const double x = ippl::apply(cfview, idx).real(); + const unsigned int d = (dim == 0) ? i : (dim == 1) ? j : k; + + const unsigned int il = idx[0] - ngh; + const unsigned int jl = idx[1] - ngh; + const unsigned int kl = idx[2] - ngh; + const index_type n = il + lDom[0].length() * (jl + kl * lDom[1].length()); + + rView(n)[dim] = (d * hr[dim]) - (d_z * x); + vView(n)[dim] = -ddot * x; +#ifdef KOKKOS_PRINT + if (dim == 0) + Kokkos::printf("setparticle: %i %g \n", d, x); + + if (n == 0) + Kokkos::printf("zeldo: dim, pos0, x, v, re: d_z= %g ddot= %g \n", d_z, + ddot); + Kokkos::printf("zeldo: %i %g %g %g %g \n", dim, d + 0.5, -d_z * x, + -ddot * x, x); +#endif + }); + cfield_m = tmpcfield; + } + Kokkos::fence(); + } + } + + /** + * @brief Create particles using Zarijas initializer method + */ + + void createParticles() { + Inform m2a("createParticles ", INFORM_ALL_NODES); + Inform msg("createParticles "); + + static IpplTimings::TimerRef creaIpplParts = IpplTimings::getTimer("creaIpplParts"); + IpplTimings::startTimer(creaIpplParts); + size_type nloc = this->totalP_m / ippl::Comm->size(); + std::shared_ptr pc = this->pcontainer_m; + + std::ostringstream oss; + oss << "pc->create nloc= " << nloc; + + Kokkos::Profiling::pushRegion(oss.str()); + pc->create(nloc); + pc->m = this->M_m / this->totalP_m; + IpplTimings::stopTimer(creaIpplParts); + msg << "pc->create(nloc) done" << endl; + Kokkos::Profiling::popRegion(); + + /** + construct Pk_m + */ + static IpplTimings::TimerRef iniPwrSpec = IpplTimings::getTimer("initPwrSpec"); + IpplTimings::startTimer(iniPwrSpec); + initPwrSpec(); + IpplTimings::stopTimer(iniPwrSpec); + + msg << "Construct Pk_m done" << endl; + + /** + the following code can be found + as standalone test in test/particles/zeldo-test-mp1.cpp + */ + + static IpplTimings::TimerRef linZeldo = IpplTimings::getTimer("LinearZeldoInitMP"); + IpplTimings::startTimer(linZeldo); + LinearZeldoInitMP(); + IpplTimings::stopTimer(linZeldo); + msg << "LinearZeldoInitMP() done" << endl; + // Load Balancer Initialisation + /* + auto* mesh = &this->fcontainer_m->getMesh(); + auto* FL = &this->fcontainer_m->getFL(); + if ((this->lbt_m != 1.0) && (ippl::Comm->size() > 1)) { + this->isFirstRepartition_m = true; + this->loadbalancer_m->initializeORB(FL, mesh); + this->loadbalancer_m->repartition(FL, mesh, this->isFirstRepartition_m); + } + */ + pc->update(); + m2a << "local number of galaxies after initializer " << pc->getLocalNum() << endl; } /** * @brief Read particle data from a file. */ void readParticles() { - Inform mes("Reading Particles"); + Inform msg("Reading Particles"); size_type nloc = this->totalP_m / ippl::Comm->size(); - mes << "Local number of particles: " << nloc << endl; + msg << "Local number of particles: " << nloc << endl; std::shared_ptr pc = this->pcontainer_m; pc->create(nloc); pc->m = this->M_m / this->totalP_m; @@ -150,7 +627,7 @@ class StructureFormationManager : public GravityManager { auto* mesh = &this->fcontainer_m->getMesh(); auto* FL = &this->fcontainer_m->getFL(); if ((this->lbt_m != 1.0) && (ippl::Comm->size() > 1)) { - mes << "Starting first repartition" << endl; + msg << "Starting first repartition" << endl; this->isFirstRepartition_m = true; this->loadbalancer_m->initializeORB(FL, mesh); this->loadbalancer_m->repartition(FL, mesh, this->isFirstRepartition_m); @@ -209,9 +686,9 @@ class StructureFormationManager : public GravityManager { } // Boundaries of Particle Positions - mes << "Minimum Position: " << MinPos << endl; - mes << "Maximum Position: " << MaxPos << endl; - mes << "Defined maximum: " << this->rmax_m << endl; + msg << "Minimum Position: " << MinPos << endl; + msg << "Maximum Position: " << MaxPos << endl; + msg << "Defined maximum: " << this->rmax_m << endl; // Number of Particles if (nloc != ParticlePositions.size()) { @@ -222,7 +699,7 @@ class StructureFormationManager : public GravityManager { // Particle positions and velocities, which are read in above from the initial // conditions file, are assigned to the particle attributes R and V in the particle // container. - mes << "successfully done." << endl; + msg << "successfully done." << endl; auto R_host = pc->R.getHostMirror(); auto V_host = pc->V.getHostMirror(); @@ -257,14 +734,14 @@ class StructureFormationManager : public GravityManager { printf("first repartition works \n"); } - mes << "Assignment of positions and velocities done." << endl; + msg << "Assignment of positions and velocities done." << endl; } /** * @brief Read particle data from a file and assign to the domain. */ void readParticlesDomain() { - Inform mes("Reading Particles"); + Inform msg("Reading Particles"); this->fcontainer_m->getRho() = 0.0; @@ -272,7 +749,7 @@ class StructureFormationManager : public GravityManager { auto* mesh = &this->fcontainer_m->getMesh(); auto* FL = &this->fcontainer_m->getFL(); if ((this->lbt_m != 1.0) && (ippl::Comm->size() > 1)) { - mes << "Starting first repartition" << endl; + msg << "Starting first repartition" << endl; this->isFirstRepartition_m = true; this->loadbalancer_m->initializeORB(FL, mesh); this->loadbalancer_m->repartition(FL, mesh, this->isFirstRepartition_m); @@ -320,12 +797,12 @@ class StructureFormationManager : public GravityManager { // positions. This avoids double-counting of particles and ensures the total // number of particles is conserved. if (Pos == Max[j]) { - mes << "Particle was on edge. Shift position from " << Pos << " to " + msg << "Particle was on edge. Shift position from " << Pos << " to " << Pos * 0.9999 << endl; Pos = 0.9999 * Pos; } if (Pos == 0) { - mes << "Particle was on edge. Shift position from " << Pos << " to " + msg << "Particle was on edge. Shift position from " << Pos << " to " << 0.0001 * Max[j] << endl; Pos = 0.0001 * Max[j]; } @@ -386,7 +863,7 @@ class StructureFormationManager : public GravityManager { printf("first repartition works \n"); } - mes << "Assignment of positions and velocities done." << endl; + msg << "Assignment of positions and velocities done." << endl; } /** @@ -396,7 +873,7 @@ class StructureFormationManager : public GravityManager { if (this->stepMethod_m == "LeapFrog") { LeapFrogStep(); } else { - throw IpplException(TestName, "Step method is not set/recognized!"); + throw IpplException("StructureFormation ", "Step method is not set/recognized!"); } } @@ -480,12 +957,12 @@ class StructureFormationManager : public GravityManager { * @param index Current time step number */ void savePositions(unsigned int index) { - Inform mes("Saving Particles"); + Inform msg("Saving Particles"); static IpplTimings::TimerRef SavingTimer = IpplTimings::getTimer("Save Data"); IpplTimings::startTimer(SavingTimer); - mes << "snapshot " << this->it_m << endl; + msg << "snapshot " << this->it_m << endl; std::stringstream ss; if (ippl::Comm->size() == 1) @@ -501,7 +978,9 @@ class StructureFormationManager : public GravityManager { if (!file.is_open()) { std::cerr << "Error opening saving file!" << std::endl; return; - } + } else + file << "# rx \t ry \t rz \n"; + std::shared_ptr pc = this->pcontainer_m; auto Rview = this->pcontainer_m->R.getView(); @@ -516,23 +995,26 @@ class StructureFormationManager : public GravityManager { Kokkos::deep_copy(V_host, Vview); Kokkos::deep_copy(F_host, Fview); - double a = this->a_m; + // double a = this->a_m; // Write data to the file for (unsigned int i = 0; i < pc->getLocalNum(); ++i) { for (unsigned int d = 0; d < Dim; ++d) - file << R_host(i)[d] << ","; + file << R_host(i)[d] << " \t"; for (unsigned int d = 0; d < Dim; ++d) - file << V_host(i)[d] << ","; + file << V_host(i)[d] << " \t"; + /* + for (unsigned int d = 0; d < Dim; ++d) file << -4 * M_PI * this->G / (a * a) * F_host(i)[d] << ","; + */ file << "\n"; } ippl::Comm->barrier(); // Close the file stream file.close(); - mes << "done." << endl; + msg << "done." << endl; IpplTimings::stopTimer(SavingTimer); } diff --git a/demos/cosmology/mc-4-Initializer/Cosmology.cpp b/demos/cosmology/mc-4-Initializer/Cosmology.cpp new file mode 100644 index 0000000000..66d31c9041 --- /dev/null +++ b/demos/cosmology/mc-4-Initializer/Cosmology.cpp @@ -0,0 +1,773 @@ +/* + Initializer: + Cosmology.cpp + + This file has implementation for routines: + + GetParameters(DataBase) + which gets cosmology parameters from the + DataBase (defined in DataBase.h), and stores them + with this class as well. + + Sigma_r(R, AA) + which returns mass variance on a scale R, of the + linear density field with normalization AA. + + TransferFunction(k) + which returns transfer function for the mode k in + Fourier space. + + GrowthFactor(z) + which returns the linear growth factor at redshift z. + + GrowthDeriv(z) + which returns the derivative of the linear growth factor + at redshift z. + + Zarija Lukic, February 2009 + zarija@lanl.gov +*/ +#include "Ippl.h" + +#include "Cosmology.h" + +#include +#include +#include +#include +#include +#include + +#include "TypesAndDefs.h" + +#include "DataBase.h" +#include "Random/Distribution.h" +#include "Random/InverseTransformSampling.h" +#include "Random/NormalDistribution.h" +#include "Random/Randn.h" + +#ifdef USENAMESPACE +namespace initializer { +#endif + + void CosmoClass::SetParameters(GlobalStuff& DataBase, const char* tfName) { + std::ifstream inputFile; + int i, ln; + real tmp, tfcdm, tfbar, norm; + real akh1, akh2, alpha; + // real file_kmax = M_PI/DataBase.box_size * DataBase.ngrid; + + Omega_m = DataBase.Omega_m; + Omega_bar = DataBase.Omega_bar; + Omega_nu = DataBase.Omega_nu; + Omega_r = DataBase.Omega_r; + + h = DataBase.Hubble; + n_s = DataBase.n_s; + TFFlag = DataBase.TFFlag; + w_de = DataBase.w_de; + N_nu = DataBase.N_nu; + z_in = DataBase.z_in; + + cobe_temp = 2.728; // COBE/FIRAS CMB temperature in K + tt = cobe_temp / 2.7 * cobe_temp / 2.7; + + if (TFFlag == 0) { // CMBFast transfer function + // inputFile.open("cmb.tf", std::ios::in); + inputFile.open(tfName, std::ios::in); + if (!inputFile) { + printf("Cannot open CMBFast file!"); + abort(); + } + ln = 0; + while (!inputFile.eof()) { + inputFile >> tmp >> tfcdm >> tfbar >> tmp >> tmp >> tmp >> tmp; + ++ln; + } + table_size = ln - 1; // It also reads newline... + table_kk = (real*)malloc(table_size * sizeof(real)); + table_tf = (real*)malloc(table_size * sizeof(real)); + inputFile.close(); + + inputFile.open(tfName, std::ios::in); + inputFile >> table_kk[0] >> tfcdm >> tfbar >> tmp >> tmp >> tmp >> tmp; + table_tf[0] = tfbar * Omega_bar / Omega_m + tfcdm * (Omega_m - Omega_bar) / Omega_m; + norm = table_tf[0]; + for (i = 1; i < table_size; ++i) { + inputFile >> table_kk[i] >> tfcdm >> tfbar >> tmp >> tmp >> tmp >> tmp; + table_tf[i] = tfbar * Omega_bar / Omega_m + tfcdm * (Omega_m - Omega_bar) / Omega_m; + table_tf[i] = + table_tf[i] + / norm; // Normalization + // std::cout << table_kk[0] << " " << table_tf[0] << std::endl; + } + last_k = table_kk[table_size - 1]; + + /* + + Brrrrrr adafixme: need to comment this out because of abort + + if (last_k < file_kmax) { + std::cout << "CMBFast file goes only to k = " << last_k << std::endl; + std::cout << "Aborting" << std::endl << std::flush; + int rc=801; + MPI_Abort(MPI_COMM_WORLD, rc); + } + */ + + inputFile.close(); + tan_f = &CosmoClass::cmbfast; + } else if (TFFlag == 1) { // Klypin-Holtzmann transfer function + akh1 = pow(46.9 * Omega_m * h * h, 0.670) * (1.0 + pow(32.1 * Omega_m * h * h, -0.532)); + akh2 = pow(12.0 * Omega_m * h * h, 0.424) * (1.0 + pow(45.0 * Omega_m * h * h, -0.582)); + alpha = pow(akh1, -1.0 * Omega_bar / Omega_m) + * pow(akh2, pow(-1.0 * Omega_bar / Omega_m, 3.0)); + kh_tmp = Omega_m * h * sqrt(alpha) * pow(1.0 - Omega_bar / Omega_m, 0.6); + tan_f = &CosmoClass::klypin_holtzmann; + last_k = 10.0; + } else if (TFFlag == 2) { // Hu-Sugiyama transfer function + akh1 = pow(46.9 * Omega_m * h * h, 0.670) * (1.0 + pow(32.1 * Omega_m * h * h, -0.532)); + akh2 = pow(12.0 * Omega_m * h * h, 0.424) * (1.0 + pow(45.0 * Omega_m * h * h, -0.582)); + alpha = pow(akh1, -1.0 * Omega_bar / Omega_m) + * pow(akh2, pow(-1.0 * Omega_bar / Omega_m, 3.0)); + kh_tmp = Omega_m * h * sqrt(alpha); + tan_f = &CosmoClass::hu_sugiyama; + last_k = 10.0; + } else if (TFFlag == 3) { // Peacock-Dodds transfer function + kh_tmp = Omega_m * h * exp(-2.0 * Omega_bar); + tan_f = &CosmoClass::peacock_dodds; + last_k = 10.0; + } else if (TFFlag == 4) { // BBKS transfer function + kh_tmp = Omega_m * h; + tan_f = &CosmoClass::bbks; + last_k = 10.0; + } else if (TFFlag == 5) { // Eisenstein-Hu transfer function + real f_bar, f_nub, f_c, f_cb, p_c, p_cb; + real omhh, obhh, z_drag, z_equality, y_d, R_drag, R_equality, k_equality; + f_nu = DataBase.Omega_nu / Omega_m; + f_bar = Omega_bar / Omega_m; + f_nub = f_nu + f_bar; + f_c = 1.0 - f_nu - f_bar; + f_cb = 1.0 - f_nu; + p_c = (5.0 - sqrt(1.0 + 24.0 * f_c)) / 4.0; + p_cb = (5.0 - sqrt(1.0 + 24.0 * f_cb)) / 4.0; + + omhh = Omega_m * h * h; + obhh = Omega_bar * h * h; + z_equality = 2.50e4 * omhh * pow(tt, -2.0); + z_drag = 0.313 * pow(omhh, -0.419) * (1. + 0.607 * pow(omhh, 0.674)); + z_drag = 1.0 + z_drag * pow(obhh, 0.238 * pow(omhh, 0.223)); + z_drag = 1291.0 * pow(omhh, 0.251) / (1.0 + 0.659 * pow(omhh, 0.828)) * z_drag; + y_d = (1.0 + z_equality) / (1. + z_drag); + + alpha_nu = (f_c / f_cb) * (5.0 - 2.0 * (p_c + p_cb)) / (5.0 - 4.0 * p_cb); + alpha_nu = alpha_nu * (1.0 - 0.553 * f_nub + 0.126 * pow(f_nub, 3.0)); + alpha_nu = alpha_nu / (1.0 - 0.193 * sqrt(f_nu * N_nu) + 0.169 * f_nu * pow(N_nu, 0.2)); + alpha_nu = alpha_nu * pow(1.0 + y_d, p_cb - p_c); + alpha_nu = alpha_nu + * (1.0 + + (p_c - p_cb) / 2. * (1. + 1. / (3. - 4. * p_c) / (7. - 4. * p_cb)) + / (1. + y_d)); + + k_equality = 0.0746 * omhh / tt; + R_drag = 31.5 * obhh * pow(tt, -2.) * 1000.0 / (1.0 + z_drag); + R_equality = 31.5 * obhh * pow(tt, -2.) * 1000.0 / (1.0 + z_equality); + sound_horizon = + 2. / 3. / k_equality * sqrt(6. / R_equality) + * log((sqrt(1. + R_drag) + sqrt(R_drag + R_equality)) / (1. + sqrt(R_equality))); + + beta_c = 1. / (1. - 0.949 * f_nub); + tan_f = &CosmoClass::eisenstein_hu; + last_k = 10.0; + } else { + printf("Cosmology::SetParameters: TFFlag has to be 0, 1, 2, 3, 4, or 5!\n"); + abort(); + } + + // neutrino stuff below + + /** + * v_0(z) = (1+z) c k_B T_\nu / (m_\nu c^2) + * m_\nu c^2 = 94 \Omega_\nu h^2 [eV] + * v_0(z) = (1+z) / (\Omega_nu h^2) \times c k T_\nu / 94 + * = (1+z) / (\Omega_nu h^2) \times 0.536 + * + * where T_\nu = 1.95 K + * for more details see: http://arxiv.org/abs/astro-ph/9305011 + * Note: we used 94 instead of 91.4 + */ + neut_vv0 = (1.0 + z_in) * 0.536 / (Omega_nu * h * h); + neut_p = (real*)(malloc(neut_nmax * sizeof(real))); + neut_c = (real*)(malloc(neut_nmax * sizeof(real))); + // generate cumulative array + genFDDist(); + + return; + } + + CosmoClass::~CosmoClass() { + if (TFFlag == 0) { + free(table_kk); + free(table_tf); + } + free(neut_p); + free(neut_c); + } + +#define pi M_PI + + /* Transfer functions */ + + real CosmoClass::cmbfast(real k) { + real t_f; + t_f = interpolate(table_kk, table_tf, table_size, k); + return (t_f); + } + + real CosmoClass::eisenstein_hu(real k) { + // t_f is the T_master(k) in Eisenstein & Hu 1999, eq. 24 + real q, q_eff, q_nu, gamma_eff, omhh, t_f; + k = k * h; + omhh = Omega_m * h * h; + q = k * tt / omhh; + gamma_eff = + omhh + * (sqrt(alpha_nu) + (1. - sqrt(alpha_nu)) / (1. + pow(0.43 * k * sound_horizon, 4.0))); + q_eff = k * tt / gamma_eff; + t_f = log(exp(1.0) + 1.84 * beta_c * sqrt(alpha_nu) * q_eff); + t_f = t_f / (t_f + q_eff * q_eff * (14.4 + 325. / (1. + 60.5 * pow(q_eff, 1.11)))); + q_nu = 3.92 * q * sqrt(N_nu) / f_nu; + t_f = t_f + * (1. + + (1.2 * pow(f_nu, 0.64) * pow(N_nu, 0.3 + 0.6 * f_nu)) + / (pow(q_nu, -1.6) + pow(q_nu, 0.8))); + return (t_f); + } + + real CosmoClass::klypin_holtzmann(real k) { + real qkh, t_f; + + if (k == 0.0) + return (0.0); + qkh = k * tt / kh_tmp; + /* NOTE: the following line has 0/0 for k=0. + This was taken care of at the beginning of the routine. */ + t_f = log(1.0 + 2.34 * qkh) / (2.34 * qkh) + * pow(1.0 + 13.0 * qkh + pow(10.5 * qkh, 2.0) + pow(10.4 * qkh, 3.0) + + pow(6.51 * qkh, 4.0), + -0.25); + return (t_f); + } + + real CosmoClass::hu_sugiyama(real k) { + real qkh, t_f; + + if (k == 0.0) + return (0.0); + qkh = k * tt / kh_tmp; + /* NOTE: the following line has 0/0 for k=0. + This was taken care of at the beginning of the routine. */ + t_f = log(1.0 + 2.34 * qkh) / (2.34 * qkh) + * pow(1.0 + 3.89 * qkh + pow(16.1 * qkh, 2.0) + pow(5.46 * qkh, 3.0) + + pow(6.71 * qkh, 4.0), + -0.25); + return (t_f); + } + + real CosmoClass::peacock_dodds(real k) { + real qkh, t_f; + + if (k == 0.0) + return (0.0); + qkh = k / kh_tmp; + /* NOTE: the following line has 0/0 for k=0. + This was taken care of at the beginning of the routine. */ + t_f = log(1.0 + 2.34 * qkh) / (2.34 * qkh) + * pow(1.0 + 3.89 * qkh + pow(16.1 * qkh, 2.0) + pow(5.46 * qkh, 3.0) + + pow(6.71 * qkh, 4.0), + -0.25); + return (t_f); + } + + real CosmoClass::bbks(real k) { + real qkh, t_f; + + if (k == 0.0) + return (0.0); + qkh = k / kh_tmp; + /* NOTE: the following line has 0/0 for k=0. + This was taken care of at the beginning of the routine. */ + t_f = log(1.0 + 2.34 * qkh) / (2.34 * qkh) + * pow(1.0 + 3.89 * qkh + pow(16.1 * qkh, 2.0) + pow(5.46 * qkh, 3.0) + + pow(6.71 * qkh, 4.0), + -0.25); + return (t_f); + } + + real CosmoClass::TransferFunction(real k) { + return ((this->*tan_f)(k)); + } + + /* Sigma(R) routines */ + + real CosmoClass::sigma2(real k) { + real prim_ps, w_f, t_f, s2; + + prim_ps = Pk_norm * pow(k, n_s); + w_f = 3.0 * (sin(k * R_M) - k * R_M * cos(k * R_M)) / pow(k * R_M, 3.0); + t_f = TransferFunction(k); + s2 = k * k * t_f * t_f * w_f * w_f * prim_ps; + return (s2); + } + + real CosmoClass::Sigma_r(real R, real AA) { + real sigma; + const real k_min = 0.0, k_max = last_k; + + R_M = R; + Pk_norm = AA; + + sigma = 1.0 / (2.0 * pi * pi) * integrate(&CosmoClass::sigma2, k_min, k_max); + sigma = sqrt(sigma); + return (sigma); + } + +#undef pi + + /* Thermal velocity for neutrinos, following Fermi-Dirac distribution */ + + // Generates a table of velocities (varray) and the + // probability that the velocity is lass than this + // velocity (neut_c) + + void CosmoClass::genFDDist() { + float neut_pmax = 3.0E5; + neut_c[0] = 0.0; + neut_p[0] = 0.0; + real norm = integrate(&CosmoClass::FD, 0, neut_pmax); + + for (int j = 1; j < neut_nmax; j++) { + neut_p[j] = (j)*neut_pmax / (1.0 * neut_nmax - 1); + neut_c[j] = integrate(&CosmoClass::FD, 0, neut_p[j]) / norm; + } + + return; + } + + real CosmoClass::FD(real vv) { + return vv * vv / (exp(vv / neut_vv0) + 1.0); + } + + real CosmoClass::Maxwell(real vv) { + real T = 1.95; //[K] + real k = 8.617343e-5; //[eV/K] + real c = 299792.458; //[km/s] + real m = 94 * Omega_nu * h * h / (c * c); + + return 4 * M_PI * pow(m / 2 / M_PI / k / T, 1.5) * vv * vv * exp(-m * vv * vv / 2 / k / T); + } + + void CosmoClass::FDVelocity(real& x, real& y, real& z) { + real aran, bran, pp, mu, phi; + Kokkos::Random_XorShift64_Pool<> rand_pool(12345); + auto rand_gen = rand_pool.get_state(); + aran = rand_gen.drand(); + bran = rand_gen.drand(); + rand_pool.free_state(rand_gen); + + pp = interpolate(neut_c, neut_p, neut_nmax, aran); + mu = 2.0 * (aran - 0.5); + phi = 8 * atan(1.0) * bran; + + x = pp * sqrt(1.0 - mu * mu) * cos(phi); + y = pp * sqrt(1.0 - mu * mu) * sin(phi); + z = pp * mu; + + return; + } + + /* Growth factor for flat wCDM cosmologies: */ + + void CosmoClass::GrowthFactor(real z, real* gf, real* g_dot) { + real x1, x2, dplus, ddot; + // const float redshift=200.0; + real* ystart; + + x1 = 1.0 / (1.0 + 100000.0); + x2 = 1.0 / (1.0 + z); + ystart = (real*)malloc(2 * sizeof(real)); + ystart[0] = x1; + ystart[1] = 0.0; + + odesolve(ystart, 2, x1, x2, 1.0e-6, 1.0e-6, &CosmoClass::growths, false); + // printf("Dplus = %f; Ddot = %f \n", ystart[0], ystart[1]); + + dplus = ystart[0]; + ddot = ystart[1]; + x1 = 1.0 / (1.0 + 100000.0); + x2 = 1.0; + ystart[0] = x1; + ystart[1] = 0.0; + + odesolve(ystart, 2, x1, x2, 1.0e-6, 1.0e-6, &CosmoClass::growths, false); + // printf("Dplus = %f; Ddot = %f \n", ystart[0], ystart[1]); + + *gf = dplus / ystart[0]; + *g_dot = ddot / ystart[0]; + // printf("\nGrowth factor = %f; Derivative = %f \n", dplus/ystart[0], ddot/ystart[0]); + free(ystart); + + return; + } + + void CosmoClass::growths(real a, real y[], real dydx[]) { + real H; + H = sqrt(Omega_r / (a * a * a * a) + Omega_m / (a * a * a) + + (1.0 - Omega_m - Omega_r) * pow(a, -3.0 * (1.0 + w_de))); + dydx[0] = y[1] / (a * H); + dydx[1] = -2.0 * y[1] / a + 1.5 * Omega_m * y[0] / (H * pow(a, 4.0)); + return; + } + +#define MAXSTP 10000 +#define TINY 1.0e-30 +#define SIGN(a, b) ((b) >= 0.0 ? fabs(a) : -fabs(a)) + void CosmoClass::odesolve(real ystart[], int nvar, real x1, real x2, real eps, real h1, + void (CosmoClass::*derivs)(real, real[], real[]), bool print_stat) { + int i, nstp, nok, nbad, feval; + real x, hnext, hdid, h; + real *yscal, *y, *dydx; + const real hmin = 0.0; + + feval = 0; + yscal = (real*)malloc(nvar * sizeof(real)); + y = (real*)malloc(nvar * sizeof(real)); + dydx = (real*)malloc(nvar * sizeof(real)); + x = x1; + h = SIGN(h1, x2 - x1); + nok = nbad = 0; + for (i = 0; i < nvar; ++i) { + y[i] = ystart[i]; + } + + for (nstp = 0; nstp < MAXSTP; ++nstp) { + (this->*derivs)(x, y, dydx); + ++feval; + for (i = 0; i < nvar; ++i) { + yscal[i] = fabs(y[i]) + fabs(dydx[i] * h) + TINY; + } + if ((x + h - x2) * (x + h - x1) > 0.0) + h = x2 - x; + rkqs(y, dydx, nvar, &x, h, eps, yscal, &hdid, &hnext, &feval, derivs); + if (hdid == h) + ++nok; + else + ++nbad; + if ((x - x2) * (x2 - x1) >= 0.0) { + for (i = 0; i < nvar; ++i) { + ystart[i] = y[i]; + } + free(dydx); + free(y); + free(yscal); + if (print_stat) { + printf("ODEsolve:\n"); + printf(" Evolved from x = %f to x = %f\n", x1, x2); + printf(" successful steps: %d\n", nok); + printf(" bad steps: %d\n", nbad); + printf(" function evaluations: %d\n", feval); + } + return; + } + if (fabs(hnext) <= hmin) { + printf("Step size too small in ODEsolve"); + exit(1); + } + h = hnext; + } + printf("Too many steps in ODEsolve"); + exit(1); + } +#undef MAXSTP +#undef TINY +#undef SIGN + +#define SAFETY 0.9 +#define PGROW -0.2 +#define PSHRNK -0.25 +#define ERRCON 1.89e-4 + static real maxarg1, maxarg2, minarg1, minarg2; +#define FMAX(a, b) (maxarg1 = (a), maxarg2 = (b), (maxarg1) > (maxarg2) ? (maxarg1) : (maxarg2)) +#define FMIN(a, b) (minarg1 = (a), minarg2 = (b), (minarg1) < (minarg2) ? (minarg1) : (minarg2)) + void CosmoClass::rkqs(real y[], real dydx[], int n, real* x, real htry, real eps, real yscal[], + real* hdid, real* hnext, int* feval, + void (CosmoClass::*derivs)(real, real[], real[])) { + int i; + real errmax, h, htemp, xnew, *yerr, *ytemp; + + yerr = (real*)malloc(n * sizeof(real)); + ytemp = (real*)malloc(n * sizeof(real)); + h = htry; + + for (;;) { + rkck(y, dydx, n, *x, h, ytemp, yerr, derivs); + *feval += 5; + errmax = 0.0; + for (i = 0; i < n; ++i) { + errmax = FMAX(errmax, fabs(yerr[i] / yscal[i])); + } + errmax /= eps; + if (errmax <= 1.0) + break; + htemp = SAFETY * h * pow(errmax, PSHRNK); + h = (h >= 0.0 ? FMAX(htemp, 0.1 * h) : FMIN(htemp, 0.1 * h)); + xnew = (*x) + h; + if (xnew == *x) { + printf("Stepsize underflow in ODEsolve rkqs"); + exit(1); + } + } + if (errmax > ERRCON) + *hnext = SAFETY * h * pow(errmax, PGROW); + else + *hnext = 5.0 * h; + *x += (*hdid = h); + for (i = 0; i < n; ++i) { + y[i] = ytemp[i]; + } + free(ytemp); + free(yerr); + } +#undef SAFETY +#undef PGROW +#undef PSHRNK +#undef ERRCON +#undef FMAX +#undef FMIN + + /* Cash-Karp Runge-Kutta Step, based on the + work of Fehlberg, who figured out that six function evaluations could + be used to determine two Runge-Kutta steps, one fourth-order and one + fifth-order. The difference between these estimates can be used as a + truncation error for adjusting the stepsize. */ + void CosmoClass::rkck(real y[], real dydx[], int n, real x, real h, real yout[], real yerr[], + void (CosmoClass::*derivs)(real, real[], real[])) { + int i; + static real a2 = 0.2, a3 = 0.3, a4 = 0.6, a5 = 1.0, a6 = 0.875, b21 = 0.2, b31 = 3.0 / 40.0, + b32 = 9.0 / 40.0, b41 = 0.3, b42 = -0.9, b43 = 1.2, b51 = -11.0 / 54.0, + b52 = 2.5, b53 = -70.0 / 27.0, b54 = 35.0 / 27.0, b61 = 1631.0 / 55296.0, + b62 = 175.0 / 512.0, b63 = 575.0 / 13824.0, b64 = 44275.0 / 110592.0, + b65 = 253.0 / 4096.0, c1 = 37.0 / 378.0, c3 = 250.0 / 621.0, c4 = 125.0 / 594.0, + c6 = 512.0 / 1771.0, dc5 = -277.00 / 14336.0; + real dc1 = c1 - 2825.0 / 27648.0, dc3 = c3 - 18575.0 / 48384.0, + dc4 = c4 - 13525.0 / 55296.0, dc6 = c6 - 0.25; + real *ak2, *ak3, *ak4, *ak5, *ak6, *ytemp; + + ak2 = (real*)malloc(n * sizeof(real)); + ak3 = (real*)malloc(n * sizeof(real)); + ak4 = (real*)malloc(n * sizeof(real)); + ak5 = (real*)malloc(n * sizeof(real)); + ak6 = (real*)malloc(n * sizeof(real)); + ytemp = (real*)malloc(n * sizeof(real)); + + for (i = 0; i < n; ++i) + ytemp[i] = y[i] + b21 * h * dydx[i]; + (this->*derivs)(x + a2 * h, ytemp, ak2); + for (i = 0; i < n; ++i) + ytemp[i] = y[i] + h * (b31 * dydx[i] + b32 * ak2[i]); + (this->*derivs)(x + a3 * h, ytemp, ak3); + for (i = 0; i < n; ++i) + ytemp[i] = y[i] + h * (b41 * dydx[i] + b42 * ak2[i] + b43 * ak3[i]); + (this->*derivs)(x + a4 * h, ytemp, ak4); + for (i = 0; i < n; ++i) + ytemp[i] = y[i] + h * (b51 * dydx[i] + b52 * ak2[i] + b53 * ak3[i] + b54 * ak4[i]); + (this->*derivs)(x + a5 * h, ytemp, ak5); + for (i = 0; i < n; ++i) + ytemp[i] = + y[i] + + h * (b61 * dydx[i] + b62 * ak2[i] + b63 * ak3[i] + b64 * ak4[i] + b65 * ak5[i]); + (this->*derivs)(x + a6 * h, ytemp, ak6); + for (i = 0; i < n; ++i) + yout[i] = y[i] + h * (c1 * dydx[i] + c3 * ak3[i] + c4 * ak4[i] + c6 * ak6[i]); + for (i = 0; i < n; ++i) + yerr[i] = + h * (dc1 * dydx[i] + dc3 * ak3[i] + dc4 * ak4[i] + dc5 * ak5[i] + dc6 * ak6[i]); + + free(ytemp); + free(ak6); + free(ak5); + free(ak4); + free(ak3); + free(ak2); + } + + /* Numerical integration routines from Numerical Recipes + calling: + x = integrate(&CosmoClass::func, a, b) + + where + func -- is the (real) function whose integral is calculated, + and has to be member of the CosmoClass + a -- is the (real) lower boundary for integration + b -- is the (real) upper boundary for integration + */ + +#define FUNC(x) ((this->*func)(x)) + real CosmoClass::midpoint(real (CosmoClass::*func)(real), real a, real b, int n) { + real x, tnm, sum, del, ddel; + static real s; + int it, j; + + if (n == 1) { + return (s = (b - a) * FUNC(0.5 * (a + b))); + } else { + it = 1; + for (j = 1; j < n - 1; ++j) + it *= 3; + tnm = it; + del = (b - a) / (3.0 * tnm); + ddel = del + del; + x = a + 0.5 * del; + sum = 0.0; + for (j = 1; j <= it; ++j) { + sum += FUNC(x); + x += ddel; + sum += FUNC(x); + x += del; + } + s = (s + (b - a) * sum / tnm) / 3.0; + return (s); + } + } +#undef FUNC + +#define EPS 1.0e-4 +#define JMAX 20 +#define JMIN 5 + real CosmoClass::integrate(real (CosmoClass::*func)(real), real a, real b) { + real sol, old; + int j; + + sol = 0.0; + old = -1.0e26; + for (j = 1; j <= JMAX; ++j) { + sol = midpoint(func, a, b, j); + if (j > JMIN) { + if (fabs(sol - old) < EPS * fabs(old) || (sol == 0.0 && old == 0.0)) + return (sol); + } + old = sol; + } + printf("integrate: no convergence!"); + abort(); + } +#undef EPS +#undef JMAX +#undef JMIN + + /* Linear interpolation routine + calling: + y = interpolate(xx, yy, n, x) + + where + xx -- is the (real) x array of ordered data + yy -- is the (real) y array + n -- size of above arrays + x -- is the (real) point whose y value should be interpolated + */ + + real CosmoClass::interpolate(real xx[], real yy[], unsigned long n, real x) { + real y, dx, dy; + unsigned long jlo; + + locate(xx, n, x, &jlo); + // Linear interpolation: + dx = xx[jlo] - xx[jlo + 1]; + dy = yy[jlo] - yy[jlo + 1]; + y = dy / dx * (x - xx[jlo]) + yy[jlo]; + + return (y); + } + + /* Routines for locating a value in an ordered table + from Numerical Recipes */ + + void CosmoClass::locate(real xx[], unsigned long n, real x, unsigned long* j) { + unsigned long ju, jm, jl; + int ascnd; + + jl = 0; + ju = n - 1; + ascnd = (xx[n - 1] >= xx[0]); + while (ju - jl > 1) { + jm = (ju + jl) / 2; + if ((x >= xx[jm]) == ascnd) + jl = jm; + else + ju = jm; + } + if (x == xx[0]) + *j = 0; + else if (x == xx[n - 1]) + *j = n - 2; + else + *j = jl; + return; + } + + void CosmoClass::hunt(real xx[], unsigned long n, real x, unsigned long* jlo) { + unsigned long jm, jhi, inc; + int ascnd; + + ascnd = (xx[n - 1] >= xx[0]); + if (*jlo <= 0 || *jlo > n) { + *jlo = 0; + jhi = n + 1; + } else { + inc = 1; + if (x >= xx[*jlo] == ascnd) { + if (*jlo == n) + return; + jhi = (*jlo) + 1; + while (x >= xx[jhi] == ascnd) { + *jlo = jhi; + inc += inc; + jhi = (*jlo) + inc; + if (jhi > n) { + jhi = n + 1; + break; + } + } + } else { + if (*jlo == 1) { + *jlo = 0; + return; + } + jhi = (*jlo)--; + while (x < xx[*jlo] == ascnd) { + jhi = (*jlo); + inc <<= 1; + if (inc >= jhi) { + *jlo = 0; + break; + } else + *jlo = jhi - inc; + } + } + } + while (jhi - (*jlo) != 1) { + jm = (jhi + (*jlo)) >> 1; + if (x >= xx[jm] == ascnd) + *jlo = jm; + else + jhi = jm; + } + if (x == xx[n]) + *jlo = n - 1; + if (x == xx[1]) + *jlo = 1; + } + +#ifdef USENAMESPACE +} +#endif diff --git a/demos/cosmology/mc-4-Initializer/Cosmology.h b/demos/cosmology/mc-4-Initializer/Cosmology.h new file mode 100644 index 0000000000..c161759a67 --- /dev/null +++ b/demos/cosmology/mc-4-Initializer/Cosmology.h @@ -0,0 +1,85 @@ +/* + Initializer: + Cosmology.h + + Defines CosmoStuff class. This will be used + for calculating many quantities -- transfer function, + mass variance on certain scale, growth factor. + See Cosmology.cpp for more explanations and + actual implementation. + + Zarija Lukic, February 2009 + zarija@lanl.gov +*/ + +#ifndef Cosmology_Header_Included +#define Cosmology_Header_Included + +#include "TypesAndDefs.h" + +#ifdef USENAMESPACE +namespace initializer { +#endif + class GlobalStuff; + + class CosmoClass { + public: + CosmoClass() {}; + void SetParameters(GlobalStuff& DataBase, const char* tfName); + real Sigma_r(real R, real norm); + void GrowthFactor(real z, real* gf, real* g_dot); + real TransferFunction(real k); + void FDVelocity(real& x, real& y, real& z); + ~CosmoClass(); + + private: + real cobe_temp, tt; + real Omega_m, Omega_bar, Omega_nu, Omega_r, h, n_s, w_de, z_in; + real sound_horizon, alpha_nu, beta_c, f_nu, N_nu; + int TFFlag; + real *table_kk, *table_tf; + unsigned long table_size; + + real kh_tmp, R_M, Pk_norm, last_k; + + // neutrino stuff + real *neut_p, *neut_c; + real neut_vv0; + static const unsigned int neut_nmax = 1000; + // FIXME: check upper limit for integral + // static const unsigned int neut_pmax = 3.0e5; + + real (CosmoClass::*tan_f)(real k); + real integrate(real (CosmoClass::*func)(real), real a, real b); + real midpoint(real (CosmoClass::*func)(real), real a, real b, int n); + real interpolate(real xx[], real yy[], unsigned long n, real x); + void locate(real xx[], unsigned long n, real x, unsigned long* j); + void hunt(real xx[], unsigned long n, real x, unsigned long* jlo); + real eisenstein_hu(real k); + real klypin_holtzmann(real k); + real hu_sugiyama(real k); + real peacock_dodds(real k); + real bbks(real k); + real cmbfast(real k); + real sigma2(real k); + + void growths(real a, real y[], real dydx[]); + void odesolve(real ystart[], int nvar, real x1, real x2, real eps, real h1, + void (CosmoClass::*derivs)(real, real[], real[]), bool print_stat); + void rkqs(real y[], real dydx[], int n, real* x, real htry, real eps, real yscal[], + real* hdid, real* hnext, int* feval, + void (CosmoClass::*derivs)(real, real[], real[])); + void rkck(real y[], real dydx[], int n, real x, real h, real yout[], real yerr[], + void (CosmoClass::*derivs)(real, real[], real[])); + + // neutrino stuff + void genFDDist(); + real FD(real vv); + real Maxwell(real vv); + }; + +#ifdef USENAMESPACE +} +#endif + +#endif diff --git a/demos/cosmology/mc-4-Initializer/DataBase.cpp b/demos/cosmology/mc-4-Initializer/DataBase.cpp new file mode 100644 index 0000000000..006ab56ead --- /dev/null +++ b/demos/cosmology/mc-4-Initializer/DataBase.cpp @@ -0,0 +1,131 @@ +/* + Initializer: + DataBase.cpp + + Has only + GetParameters(Basedata& bdata, ParallelTools& Parallel) + routine which reads in simulation and cosmology + parameters from the main code's Basedata and stors them into + DataBase class (defined in DataBase.h). + + Zarija Lukic, February 2009 + zarija@lanl.gov +*/ + +#include "Ippl.h" + +#include "DataBase.h" + +#include +#include + +#include "TypesAndDefs.h" + +initializer::GlobalStuff& initializer::GlobalStuff::instance() { + static initializer::GlobalStuff instance; + return instance; +} + +#include "InputParser.h" + +#ifdef USENAMESPACE +namespace initializer { +#endif + + void GlobalStuff::GetParameters(InputParser& par) { + Inform msg("DataBase "); + + int np; + if (!par.getByName("np", np)) // Grid size is equal to np^3 + msg << "Error: np not found!" << endl; + ngrid = np; + + if (!par.getByName("box_size", box_size)) + msg << "Error: box_size not found!" << endl; + + int decomposition_type = 3; + + dim = decomposition_type; + + int sseed; + if (!par.getByName("seed", sseed)) + msg << "Error: seed not found!" << endl; + seed = sseed; + + if (!par.getByName("z_in", z_in)) + msg << "Error: z_in not found!" << endl; + + // Get cosmology parameters: + if (!par.getByName("hubble", Hubble)) + msg << "Error: hubble not found!" << endl; + + if (!par.getByName("Omega_m", Omega_m)) + msg << "Error: Omega_m not found!" << endl; + + if (!par.getByName("Omega_bar", Omega_bar)) + msg << "Error: Omega_bar not found!" << endl; + + if (!par.getByName("Omega_nu", Omega_nu)) + msg << "Error: Omega_nu not found!" << endl; + + if (!par.getByName("Omega_r", Omega_r)) + msg << "Error: Omega_r not found set Omega_r to zero!" << endl; + else + Omega_r = 0.0; + + if (!par.getByName("Sigma_8", Sigma_8)) + msg << "Error: Sigma_8 not found!" << endl; + + if (!par.getByName("n_s", n_s)) + msg << "Error: n_s not found!" << endl; + + if (!par.getByName("w_de", w_de)) + msg << "Error: w_de not found!" << endl; + + if (!par.getByName("TFFlag", TFFlag)) + msg << "Error: TFFlag not found!" << endl; + + // Get neutrino parameters: + if (!par.getByName("N_nu", N_nu)) + msg << "Error: N_nu not found!" << endl; + + if (!par.getByName("nu_pairs", nu_pairs)) + msg << "Error: nu_pairs not found!" << endl; + + // Get non-gaussianity parameters: + if (!par.getByName("f_NL", f_NL)) { + msg << "Warning: f_NL not found, assuming it is zero!" << endl; + f_NL = 0.0; + } + + // Get code parameters: + if (!par.getByName("PrintFormat", PrintFormat)) + msg << "Error: PrintFormat not found!" << endl; + + // Make some sanity checks: + // Errors: + if (Omega_m > 1.0) + throw IpplException("DataBase", "Cannot initialize closed universe!"); + + if (Omega_bar > Omega_m) + throw IpplException("DataBase", "More baryons than total matter content!"); + if (dim < 1 || dim > 3) + throw IpplException("DataBase", "1, 2 or 3D decomposition is possible!"); + if (Omega_nu > 0 && nu_pairs < 1) + throw IpplException("DataBase", "Omega_nu > 0, yet nu_pairs < 1!"); + + // Warnings: + if (z_in > 1000.0) + msg << "Starting before CMB epoch. Make sure you want that." << endl; + if (Sigma_8 > 1.0) + msg << "Sigma(8) > 1. Make sure you want that." << endl; + if (Hubble > 1.0) + msg << "H > 100 km/s/Mpc. Make sure you want that." << endl; + + msg << "Done loading parameters." << endl; + return; + } + +#ifdef USENAMESPACE +} +#endif diff --git a/demos/cosmology/mc-4-Initializer/DataBase.h b/demos/cosmology/mc-4-Initializer/DataBase.h new file mode 100644 index 0000000000..6309d1e64d --- /dev/null +++ b/demos/cosmology/mc-4-Initializer/DataBase.h @@ -0,0 +1,77 @@ +/* + Initializer: + DataBase.h + + Defines GlobalStuff class. (main makes a globaly visible + instance of that class -- DataBase.) + This class will hold simulation + and cosmology parameters needed by most of routines. + + Zarija Lukic, February 2009 + zarija@lanl.gov + */ + +#ifndef DataBase_Header_Included +#define DataBase_Header_Included + +#include + +#include "TypesAndDefs.h" + +#include "InputParser.h" + +#ifdef USENAMESPACE +namespace initializer { +#endif + + // class ParallelTools; + + class GlobalStuff { + public: + static GlobalStuff& instance(); + GlobalStuff(const GlobalStuff&) = delete; + GlobalStuff& operator=(const GlobalStuff&) = delete; + GlobalStuff(GlobalStuff&&) = delete; + GlobalStuff& operator=(GlobalStuff&&) = delete; + + private: + GlobalStuff() = default; + + public: + void GetParameters(InputParser& par); + + // Simulation parameters: + int ngrid; // Number of mesh points in any direction + real box_size; // Lenght of the box size in Mpc/h + int dim; // dimensionality of domain decomposition + unsigned long seed; // Master seed for random numbers + real z_in; // Starting redshift + + // Cosmological parameters: + real Omega_m; // Total matter fraction (today) + real Omega_bar; // Baryon fraction (today) + real Omega_nu; // neutrino fraction + real Omega_r; // radiation fraction + + real Hubble; // Hubble constant + real Sigma_8; // Mass RMS in 8Mpc/h spheres + real n_s; // Spectral index + real w_de; // Dark energy EOS parameter + real f_NL; // f_NL for non-gaussianity + int TFFlag; /* Transfer function used: 0=CMBFast, 1=KH, + 2=HS, 3=PD, 4=BBKS */ + + // Neutrino parameters: + int N_nu, nu_pairs; + + // Code stuff: + int PrintFormat; /* 0 - serial ASCII, 1 - serial binary, + 2 - parallel binary, 3 - HDF5, + 4 - Gadget */ + }; + +#ifdef USENAMESPACE +} +#endif + +#endif diff --git a/demos/cosmology/mc-4-Initializer/InputParser.cpp b/demos/cosmology/mc-4-Initializer/InputParser.cpp new file mode 100644 index 0000000000..e727eb72a9 --- /dev/null +++ b/demos/cosmology/mc-4-Initializer/InputParser.cpp @@ -0,0 +1,73 @@ +#include "InputParser.h" + +#ifdef USENAMESPACE +namespace initializer { +#endif + + InputParser::InputParser(std::string filename) { + filename_m = filename; + parseFile(); + } + + bool InputParser::getByName(std::string name, int& param) { + std::map::iterator it; + it = IntDict_m.find(name); + if (it == IntDict_m.end()) + return false; + else { + param = (*it).second; + return true; + } + } + + bool InputParser::getByName(std::string name, real& param) { + std::map::iterator it; + it = RealDict_m.find(name); + if (it == RealDict_m.end()) + return false; + else { + param = (*it).second; + return true; + } + } + + void InputParser::parseFile() { + std::ifstream fin; + std::string line = ""; + + fin.open(filename_m.c_str()); + while (std::getline(fin, line)) { + if (line.find("//") == 0 || line.empty()) + continue; + else { + int pos = line.find("="); + std::string name = line.substr(0, pos); + std::string svalue = line.substr(pos + 1); + if (isInt(svalue)) { + std::istringstream instr(svalue); + int value; + instr >> value; + IntDict_m.insert(std::pair(name, value)); + } else { + std::istringstream instr(svalue); + real value; + instr >> value; + RealDict_m.insert(std::pair(name, value)); + } + } + } + } + +#ifdef USENAMESPACE +} +#endif + +#ifdef TEST_PARSER +int main(int argc, char** argv) { + initializer::InputParser* par = new initializer::InputParser("ExampleMCInputFile.in"); + real zin = 0.0; + bool found = par->getByName("zin", zin); + std::cout << "found= " << found << " value= " << zin << std::endl; + return 0; +} +#endif diff --git a/demos/cosmology/mc-4-Initializer/InputParser.h b/demos/cosmology/mc-4-Initializer/InputParser.h new file mode 100644 index 0000000000..4d674a4cbd --- /dev/null +++ b/demos/cosmology/mc-4-Initializer/InputParser.h @@ -0,0 +1,35 @@ +#ifndef INPUT_PARSER_H +#define INPUT_PARSER_H + +#include +#include +#include +#include +#include + +#include "TypesAndDefs.h" + +#ifdef USENAMESPACE +namespace initializer { +#endif + + class InputParser { + public: + InputParser(std::string filename); + bool getByName(std::string name, int& param); + bool getByName(std::string name, real& param); + + private: + std::map RealDict_m; + std::map IntDict_m; + std::string filename_m; + + void parseFile(); + bool isInt(std::string str) { return str.find(".") == std::string::npos; } + }; + +#ifdef USENAMESPACE +} +#endif + +#endif diff --git a/demos/cosmology/mc-4-Initializer/Output.cpp b/demos/cosmology/mc-4-Initializer/Output.cpp new file mode 100644 index 0000000000..e82a8c01a3 --- /dev/null +++ b/demos/cosmology/mc-4-Initializer/Output.cpp @@ -0,0 +1,333 @@ +/* + * Output.cpp + * Initializer + * + * Created by Zarija on 1/27/10. + * + */ + +#include "Output.h" + +#include +#include +#include +#include +#include +#include +#include + +#include "TypesAndDefs.h" + +#include "InputParser.h" + +#ifdef USENAMESPACE +namespace initializer { +#endif + using namespace std; + + void Output::grid2phys(real* pos_x, real* pos_y, real* pos_z, real* vel_x, real* vel_y, + real* vel_z, int Npart, int np, float rL) { + long i; + + float grid2phys_pos = 1.0 * rL / np; + float grid2phys_vel = 100.0 * rL / np; + + for (i = 0; i < Npart; i++) { + pos_x[i] *= grid2phys_pos; + pos_y[i] *= grid2phys_pos; + pos_z[i] *= grid2phys_pos; + vel_x[i] *= grid2phys_vel; + vel_y[i] *= grid2phys_vel; + vel_z[i] *= grid2phys_vel; + } + + return; + } + + void Output::write_hcosmo_serial(real* pos_x, real* pos_y, real* pos_z, real* vel_x, + real* vel_y, real* vel_z, int* id, int Npart, string outBase) { + FILE* outFile; + ostringstream outName; + long i; + int MyPE, NumPEs, MasterPE, proc; + MPI_Status status; + + MPI_Comm_size(MPI_COMM_WORLD, &NumPEs); + MPI_Comm_rank(MPI_COMM_WORLD, &MyPE); + MasterPE = 0; + + outName << outBase << ".bin"; + if (MyPE == MasterPE) { + outFile = fopen(outName.str().c_str(), "wb"); + for (i = 0; i < Npart; i++) { + fwrite(&pos_x[i], sizeof(real), 1, outFile); + fwrite(&vel_x[i], sizeof(real), 1, outFile); + fwrite(&pos_y[i], sizeof(real), 1, outFile); + fwrite(&vel_y[i], sizeof(real), 1, outFile); + fwrite(&pos_z[i], sizeof(real), 1, outFile); + fwrite(&vel_z[i], sizeof(real), 1, outFile); + fwrite(&id[i], sizeof(int), 1, outFile); + } + for (proc = 0; proc < NumPEs; ++proc) { // Get data from other processors: + if (proc == MasterPE) + continue; + MPI_Recv(pos_x, Npart, MY_MPI_REAL, proc, 101, MPI_COMM_WORLD, &status); + MPI_Recv(pos_y, Npart, MY_MPI_REAL, proc, 102, MPI_COMM_WORLD, &status); + MPI_Recv(pos_z, Npart, MY_MPI_REAL, proc, 103, MPI_COMM_WORLD, &status); + MPI_Recv(vel_x, Npart, MY_MPI_REAL, proc, 104, MPI_COMM_WORLD, &status); + MPI_Recv(vel_y, Npart, MY_MPI_REAL, proc, 105, MPI_COMM_WORLD, &status); + MPI_Recv(vel_z, Npart, MY_MPI_REAL, proc, 106, MPI_COMM_WORLD, &status); + MPI_Recv(id, Npart, MY_MPI_INT, proc, 107, MPI_COMM_WORLD, &status); + for (i = 0; i < Npart; ++i) { + fwrite(&pos_x[i], sizeof(real), 1, outFile); + fwrite(&vel_x[i], sizeof(real), 1, outFile); + fwrite(&pos_y[i], sizeof(real), 1, outFile); + fwrite(&vel_y[i], sizeof(real), 1, outFile); + fwrite(&pos_z[i], sizeof(real), 1, outFile); + fwrite(&vel_z[i], sizeof(real), 1, outFile); + fwrite(&id[i], sizeof(int), 1, outFile); + } + } + fclose(outFile); + } else { // Send data to master processor: + MPI_Send(pos_x, Npart, MY_MPI_REAL, MasterPE, 101, MPI_COMM_WORLD); + MPI_Send(pos_y, Npart, MY_MPI_REAL, MasterPE, 102, MPI_COMM_WORLD); + MPI_Send(pos_z, Npart, MY_MPI_REAL, MasterPE, 103, MPI_COMM_WORLD); + MPI_Send(vel_x, Npart, MY_MPI_REAL, MasterPE, 104, MPI_COMM_WORLD); + MPI_Send(vel_y, Npart, MY_MPI_REAL, MasterPE, 105, MPI_COMM_WORLD); + MPI_Send(vel_z, Npart, MY_MPI_REAL, MasterPE, 106, MPI_COMM_WORLD); + MPI_Send(id, Npart, MY_MPI_INT, MasterPE, 107, MPI_COMM_WORLD); + } + + return; + } + + void Output::write_hcosmo_parallel(real* pos_x, real* pos_y, real* pos_z, real* vel_x, + real* vel_y, real* vel_z, int* id, int Npart, + string outBase) { + FILE* outFile; + ostringstream outName; + long i; + int MyPE; + + MPI_Comm_rank(MPI_COMM_WORLD, &MyPE); + + outName << outBase << ".bin." << MyPE; + outFile = fopen(outName.str().c_str(), "wb"); + for (i = 0; i < Npart; i++) { + fwrite(&pos_x[i], sizeof(real), 1, outFile); + fwrite(&vel_x[i], sizeof(real), 1, outFile); + fwrite(&pos_y[i], sizeof(real), 1, outFile); + fwrite(&vel_y[i], sizeof(real), 1, outFile); + fwrite(&pos_z[i], sizeof(real), 1, outFile); + fwrite(&vel_z[i], sizeof(real), 1, outFile); + fwrite(&id[i], sizeof(int), 1, outFile); + } + fclose(outFile); + + return; + } + + void Output::write_hcosmo_ascii(real* pos_x, real* pos_y, real* pos_z, real* vel_x, real* vel_y, + real* vel_z, int Npart, string outBase) { + long i; + int MyPE, NumPEs, MasterPE, proc; + MPI_Status status; + ostringstream outName; + ofstream of; + + MPI_Comm_size(MPI_COMM_WORLD, &NumPEs); + MPI_Comm_rank(MPI_COMM_WORLD, &MyPE); + MasterPE = 0; + + outName << outBase << ".txt"; + + if (MyPE == MasterPE) { + of.open(outName.str().c_str(), ios::out); + for (i = 0; i < Npart; ++i) { + of << pos_x[i] << ' '; + of << vel_x[i] << ' '; + of << pos_y[i] << ' '; + of << vel_y[i] << ' '; + of << pos_z[i] << ' '; + of << vel_z[i] << endl; + } + for (proc = 0; proc < NumPEs; ++proc) { // Get data from other processors: + if (proc == MasterPE) + continue; + MPI_Recv(pos_x, Npart, MY_MPI_REAL, proc, 101, MPI_COMM_WORLD, &status); + MPI_Recv(pos_y, Npart, MY_MPI_REAL, proc, 102, MPI_COMM_WORLD, &status); + MPI_Recv(pos_z, Npart, MY_MPI_REAL, proc, 103, MPI_COMM_WORLD, &status); + MPI_Recv(vel_x, Npart, MY_MPI_REAL, proc, 104, MPI_COMM_WORLD, &status); + MPI_Recv(vel_y, Npart, MY_MPI_REAL, proc, 105, MPI_COMM_WORLD, &status); + MPI_Recv(vel_z, Npart, MY_MPI_REAL, proc, 106, MPI_COMM_WORLD, &status); + for (i = 0; i < Npart; ++i) { + of << pos_x[i] << ' '; + of << vel_x[i] << ' '; + of << pos_y[i] << ' '; + of << vel_y[i] << ' '; + of << pos_z[i] << ' '; + of << vel_z[i] << endl; + } + } + of.close(); + } else { // Send data to master processor: + MPI_Send(pos_x, Npart, MY_MPI_REAL, MasterPE, 101, MPI_COMM_WORLD); + MPI_Send(pos_y, Npart, MY_MPI_REAL, MasterPE, 102, MPI_COMM_WORLD); + MPI_Send(pos_z, Npart, MY_MPI_REAL, MasterPE, 103, MPI_COMM_WORLD); + MPI_Send(vel_x, Npart, MY_MPI_REAL, MasterPE, 104, MPI_COMM_WORLD); + MPI_Send(vel_y, Npart, MY_MPI_REAL, MasterPE, 105, MPI_COMM_WORLD); + MPI_Send(vel_z, Npart, MY_MPI_REAL, MasterPE, 106, MPI_COMM_WORLD); + } + + return; + } + + void Output::write_gadget(InputParser& par, real* pos_x, real* pos_y, real* pos_z, real* vel_x, + real* vel_y, real* vel_z, int* id, int Npart, string outBase) { + FILE* outFile; + ostringstream outName; + long i; + int np; + real rL, Om, h, redshift; + gadget_header ghead; + int MyPE, NumPEs, MasterPE, proc; + MPI_Status status; + float x, y, z, vx, vy, vz; + unsigned int gID; + + MPI_Comm_size(MPI_COMM_WORLD, &NumPEs); + MPI_Comm_rank(MPI_COMM_WORLD, &MyPE); + MasterPE = 0; + outName << outBase << ".gadget"; + + // Gadget header -- dark matter only: + ghead.flag_sfr = 0; + ghead.flag_feedback = 0; + ghead.flag_cooling = 0; + ghead.num_files = 1; + ghead.flag_stellarage = 0; + ghead.flag_metals = 0; + for (i = 0; i < 6; ++i) { + ghead.npartTotalHighWord[i] = 0; + } + ghead.flag_entropy_instead_u = 0; + + par.getByName("box_size", rL); + par.getByName("Omega_m", Om); + par.getByName("hubble", h); + par.getByName("z_in", redshift); + ghead.BoxSize = (double)rL; // In Mpc/h + ghead.Omega0 = (double)Om; + ghead.OmegaLambda = (double)(1.0 - Om); // Flat universe + ghead.HubbleParam = (double)h; + ghead.redshift = (double)redshift; + ghead.time = (double)(1.0 / (1.0 + redshift)); + + for (i = 0; i < 6; ++i) { + ghead.npart[i] = 0; + ghead.npartTotal[i] = 0; + ghead.mass[i] = 0.0; + } + par.getByName("np", np); + ghead.npart[1] = np * np * np; + ghead.npartTotal[1] = np * np * np; + ghead.mass[1] = (double)27.75197 * Om * pow(rL / np, 3); // In 10^10 M_sun/h + + // Write it: + int blksize; +#define SKIP \ + { \ + fwrite(&blksize, sizeof(int), 1, outFile); \ + } + blksize = sizeof(ghead); + if (MyPE == MasterPE) { + outFile = fopen(outName.str().c_str(), "wb"); + SKIP fwrite(&ghead, sizeof(ghead), 1, outFile); + SKIP blksize = np * np * np * 3 * sizeof(float); + SKIP for (i = 0; i < Npart; i++) { + x = (float)pos_x[i]; // In Mpc/h + y = (float)pos_y[i]; + z = (float)pos_z[i]; + fwrite(&x, sizeof(float), 1, outFile); + fwrite(&y, sizeof(float), 1, outFile); + fwrite(&z, sizeof(float), 1, outFile); + } + for (proc = 0; proc < NumPEs; ++proc) { // Get data from other processors: + if (proc == MasterPE) + continue; + MPI_Recv(pos_x, Npart, MY_MPI_REAL, proc, 101, MPI_COMM_WORLD, &status); + MPI_Recv(pos_y, Npart, MY_MPI_REAL, proc, 102, MPI_COMM_WORLD, &status); + MPI_Recv(pos_z, Npart, MY_MPI_REAL, proc, 103, MPI_COMM_WORLD, &status); + for (i = 0; i < Npart; ++i) { + x = (float)pos_x[i]; // In Mpc/h + y = (float)pos_y[i]; + z = (float)pos_z[i]; + fwrite(&x, sizeof(float), 1, outFile); + fwrite(&y, sizeof(float), 1, outFile); + fwrite(&z, sizeof(float), 1, outFile); + } + } + SKIP + } else { // Send data to master processor: + MPI_Send(pos_x, Npart, MY_MPI_REAL, MasterPE, 101, MPI_COMM_WORLD); + MPI_Send(pos_y, Npart, MY_MPI_REAL, MasterPE, 102, MPI_COMM_WORLD); + MPI_Send(pos_z, Npart, MY_MPI_REAL, MasterPE, 103, MPI_COMM_WORLD); + } + if (MyPE == MasterPE) { + blksize = np * np * np * 3 * sizeof(float); + SKIP for (i = 0; i < Npart; i++) { + vx = (float)vel_x[i]; // In km/s + vy = (float)vel_y[i]; + vz = (float)vel_z[i]; + fwrite(&vx, sizeof(float), 1, outFile); + fwrite(&vy, sizeof(float), 1, outFile); + fwrite(&vz, sizeof(float), 1, outFile); + } + for (proc = 0; proc < NumPEs; ++proc) { // Get data from other processors: + if (proc == MasterPE) + continue; + MPI_Recv(vel_x, Npart, MY_MPI_REAL, proc, 104, MPI_COMM_WORLD, &status); + MPI_Recv(vel_y, Npart, MY_MPI_REAL, proc, 105, MPI_COMM_WORLD, &status); + MPI_Recv(vel_z, Npart, MY_MPI_REAL, proc, 106, MPI_COMM_WORLD, &status); + for (i = 0; i < Npart; ++i) { + vx = (float)vel_x[i]; // In km/s + vy = (float)vel_y[i]; + vz = (float)vel_z[i]; + fwrite(&vx, sizeof(float), 1, outFile); + fwrite(&vy, sizeof(float), 1, outFile); + fwrite(&vz, sizeof(float), 1, outFile); + } + } + SKIP + } else { // Send data to master processor: + MPI_Send(vel_x, Npart, MY_MPI_REAL, MasterPE, 104, MPI_COMM_WORLD); + MPI_Send(vel_y, Npart, MY_MPI_REAL, MasterPE, 105, MPI_COMM_WORLD); + MPI_Send(vel_z, Npart, MY_MPI_REAL, MasterPE, 106, MPI_COMM_WORLD); + } + if (MyPE == MasterPE) { + blksize = np * np * np * sizeof(int); + SKIP for (i = 0; i < Npart; i++) { + gID = (int)id[i]; + fwrite(&gID, sizeof(int), 1, outFile); + } + for (proc = 0; proc < NumPEs; ++proc) { // Get data from other processors: + if (proc == MasterPE) + continue; + MPI_Recv(id, Npart, MY_MPI_INT, proc, 107, MPI_COMM_WORLD, &status); + for (i = 0; i < Npart; ++i) { + gID = (int)id[i]; + fwrite(&gID, sizeof(int), 1, outFile); + } + } + SKIP fclose(outFile); + } else { // Send data to master processor: + MPI_Send(id, Npart, MY_MPI_INT, MasterPE, 107, MPI_COMM_WORLD); + } + + return; + } + +#ifdef USENAMESPACE +} +#endif diff --git a/demos/cosmology/mc-4-Initializer/Output.h b/demos/cosmology/mc-4-Initializer/Output.h new file mode 100644 index 0000000000..adad6f31d5 --- /dev/null +++ b/demos/cosmology/mc-4-Initializer/Output.h @@ -0,0 +1,61 @@ +/* + * Output.h + * Initializer + * + * Created by Zarija on 1/27/10. + * + */ + +#ifndef Output_Header_Included +#define Output_Header_Included + +#include + +#include "InputParser.h" + +#ifdef USENAMESPACE +namespace initializer { +#endif + using namespace std; + + class Output { + public: + void grid2phys(real* pos_x, real* pos_y, real* pos_z, real* vel_x, real* vel_y, real* vel_z, + int Npart, int np, float rL); + void write_hcosmo_ascii(real* pos_x, real* pos_y, real* pos_z, real* vel_x, real* vel_y, + real* vel_z, int Npart, string outBase); + void write_hcosmo_serial(real* pos_x, real* pos_y, real* pos_z, real* vel_x, real* vel_y, + real* vel_z, int* id, int Npart, string outBase); + void write_hcosmo_parallel(real* pos_x, real* pos_y, real* pos_z, real* vel_x, real* vel_y, + real* vel_z, int* id, int Npart, string outBase); + void write_gadget(InputParser& par, real* pos_x, real* pos_y, real* pos_z, real* vel_x, + real* vel_y, real* vel_z, int* id, int Npart, string outBase); + + private: + struct gadget_header { + int npart[6]; + double mass[6]; + double time; + double redshift; + int flag_sfr; + int flag_feedback; + unsigned int npartTotal[6]; + int flag_cooling; + int num_files; + double BoxSize; + double Omega0; + double OmegaLambda; + double HubbleParam; + int flag_stellarage; + int flag_metals; + unsigned int npartTotalHighWord[6]; + int flag_entropy_instead_u; + char fill[60]; /* fills to 256 Bytes */ + }; + }; + +#ifdef USENAMESPACE +} +#endif + +#endif diff --git a/demos/cosmology/mc-4-Initializer/TypesAndDefs.h b/demos/cosmology/mc-4-Initializer/TypesAndDefs.h new file mode 100644 index 0000000000..7d8a05b6cb --- /dev/null +++ b/demos/cosmology/mc-4-Initializer/TypesAndDefs.h @@ -0,0 +1,46 @@ +/* Definitions of real and integer variables: + if the code should be in double precision, here's + the only place where change should be made. + Also, some F90 intrinsics are defined here. + + Zarija Lukic, February 2009 + zarija@lanl.gov +*/ + +#ifndef TypesDefs_Header_Included +#define TypesDefs_Header_Included + +#ifdef USENAMESPACE +namespace initializer { +#endif + +#ifdef DOUBLE_REAL + typedef double real; +#define MY_MPI_REAL MPI_DOUBLE +#else +typedef float real; +#define MY_MPI_REAL MPI_FLOAT +#endif + +#ifdef LONG_INTEGER + typedef long integer; +#define MY_MPI_INTEGER MPI_DOUBLE +#else +// typedef int integer; +// #define MY_MPI_INTEGER MPI_INT +#endif + + typedef struct { + double re; + double im; + } my_fftw_complex; + + inline int MOD(int x, int y) { + return (x - y * (int)(x / y)); + } + +#ifdef USENAMESPACE +} +#endif + +#endif diff --git a/src/Decomposition/OrthogonalRecursiveBisection.hpp b/src/Decomposition/OrthogonalRecursiveBisection.hpp index 0ada6cfca7..2e2ad60689 100644 --- a/src/Decomposition/OrthogonalRecursiveBisection.hpp +++ b/src/Decomposition/OrthogonalRecursiveBisection.hpp @@ -1,3 +1,5 @@ +#include + #include "Utility/IpplTimings.h" namespace ippl { @@ -294,7 +296,7 @@ namespace ippl { Vector args = index - lDom.first() + nghost; // Scatter - scatterToField(std::make_index_sequence<1 << Dim>{}, view, wlo, whi, args); + scatterToField(std::make_index_sequence<1 << Dim>{}, view, wlo, whi, args, 1); }); bf_m.accumulateHalo(); diff --git a/src/Interpolation/CIC.h b/src/Interpolation/CIC.h index 5adac510d3..570a170651 100644 --- a/src/Interpolation/CIC.h +++ b/src/Interpolation/CIC.h @@ -68,6 +68,7 @@ namespace ippl { * @tparam View the field view type * @tparam T the field data type * @tparam IndexType the index type for accessing the field (default size_t) + $ @tparam Val the type of the value to be interpolated to the grid * @param view the field view on which to scatter * @param wlo lower weights for interpolation * @param whi upper weights for interpolation @@ -75,11 +76,11 @@ namespace ippl { * @param val the value to interpolate */ template + typename IndexType = size_t, typename Val = T> KOKKOS_INLINE_FUNCTION constexpr void scatterToField( const std::index_sequence&, const View& view, const Vector& wlo, const Vector& whi, - const Vector& args, T val = 1); + const Vector& args, Val val); /*! * Gathers from a field at a single point diff --git a/src/Interpolation/CIC.hpp b/src/Interpolation/CIC.hpp index 0f8d1732c6..568915c271 100644 --- a/src/Interpolation/CIC.hpp +++ b/src/Interpolation/CIC.hpp @@ -1,6 +1,5 @@ #include - namespace ippl { namespace detail { template @@ -33,14 +32,22 @@ namespace ippl { val * (interpolationWeight(wlo, whi) * ...)); } - template + template KOKKOS_INLINE_FUNCTION constexpr void scatterToField( const std::index_sequence&, const View& view, const Vector& wlo, const Vector& whi, - const Vector& args, T val) { + const Vector& args, Val val) { // The number of indices is equal to the view rank + + /** If scatterToPoint ultimately writes into `view`, convert once to its element type. + Now Val can c.f be double or floar + */ + using out_t = + std::remove_cv_t>; + const out_t v_out = static_cast(val); (scatterToPoint(std::make_index_sequence{}, view, wlo, whi, - args, val), + args, v_out), ...); } diff --git a/src/Particle/ParticleAttrib.h b/src/Particle/ParticleAttrib.h index 80b9dad20a..56e1b85353 100644 --- a/src/Particle/ParticleAttrib.h +++ b/src/Particle/ParticleAttrib.h @@ -137,14 +137,14 @@ namespace ippl { * particle range [0, size()). */ view_type getView() { - return Kokkos::subview(dview_m, - Kokkos::make_pair(size_type(0), - static_cast(*(this->localNum_mp)))); + return Kokkos::subview( + dview_m, + Kokkos::make_pair(size_type(0), static_cast(*(this->localNum_mp)))); } const view_type getView() const { - return Kokkos::subview(dview_m, - Kokkos::make_pair(size_type(0), - static_cast(*(this->localNum_mp)))); + return Kokkos::subview( + dview_m, + Kokkos::make_pair(size_type(0), static_cast(*(this->localNum_mp)))); } host_mirror_type getHostMirror() const { return Kokkos::create_mirror(getView()); } diff --git a/src/Particle/ParticleAttrib.hpp b/src/Particle/ParticleAttrib.hpp index 15cb155d9c..ca477f21d0 100644 --- a/src/Particle/ParticleAttrib.hpp +++ b/src/Particle/ParticleAttrib.hpp @@ -16,7 +16,6 @@ #include "Ippl.h" #include - #include #include "Communicate/DataTypes.h" @@ -31,6 +30,72 @@ #include "Particle/SortBuffer.h" namespace ippl { + namespace detail { + template + KOKKOS_INLINE_FUNCTION size_t scatterMappedIndex(const size_t idx, + const HashView& hashView) { + if constexpr (UseHashView) { + return static_cast(hashView(idx)); + } else { + (void)hashView; + return idx; + } + } + + template + void particleAttribScatterImpl(Field& f, const ValuesView& dview, const PositionAttrib& pp, + policy_type iteration_policy, HashView hash_array) { + constexpr unsigned Dim = Field::dim; + using PositionType = typename Field::Mesh_t::value_type; + + static IpplTimings::TimerRef scatterTimer = IpplTimings::getTimer("scatter"); + IpplTimings::startTimer(scatterTimer); + using view_type = typename Field::view_type; + view_type view = f.getView(); + + using mesh_type = typename Field::Mesh_t; + const mesh_type& mesh = f.get_mesh(); + + using vector_type = typename mesh_type::vector_type; + using value_type = typename ValuesView::non_const_value_type; + + const vector_type& dx = mesh.getMeshSpacing(); + const vector_type& origin = mesh.getOrigin(); + const vector_type invdx = 1.0 / dx; + + const FieldLayout& layout = f.getLayout(); + const NDIndex& lDom = layout.getLocalNDIndex(); + const int nghost = f.getNghost(); + + auto ppview = pp.getView(); + auto hashView = hash_array; + Kokkos::parallel_for( + "ParticleAttrib::scatter", iteration_policy, KOKKOS_LAMBDA(const size_t idx) { + // map index to possible hash_map + const size_t mapped_idx = + detail::scatterMappedIndex(idx, hashView); + + vector_type l = (ppview(mapped_idx) - origin) * invdx + 0.5; + Vector index = l; + Vector whi = l - index; + Vector wlo = 1.0 - whi; + + Vector args = index - lDom.first() + nghost; + + const value_type& val = dview(mapped_idx); + detail::scatterToField(std::make_index_sequence<1 << Field::dim>{}, view, wlo, + whi, args, val); + }); + IpplTimings::stopTimer(scatterTimer); + + static IpplTimings::TimerRef accumulateHaloTimer = + IpplTimings::getTimer("accumulateHalo"); + IpplTimings::startTimer(accumulateHaloTimer); + f.accumulateHalo(); + IpplTimings::stopTimer(accumulateHaloTimer); + } + } // namespace detail template void ParticleAttrib::create(size_type n, bool non_destructive) { @@ -133,60 +198,20 @@ namespace ippl { void ParticleAttrib::scatter( Field& f, const ParticleAttrib, Properties...>& pp, policy_type iteration_policy, hash_type hash_array) const { - constexpr unsigned Dim = Field::dim; - using PositionType = typename Field::Mesh_t::value_type; - - static IpplTimings::TimerRef scatterTimer = IpplTimings::getTimer("scatter"); - IpplTimings::startTimer(scatterTimer); - using view_type = typename Field::view_type; - view_type view = f.getView(); - - using mesh_type = typename Field::Mesh_t; - const mesh_type& mesh = f.get_mesh(); - - using vector_type = typename mesh_type::vector_type; - using value_type = typename ParticleAttrib::value_type; - - const vector_type& dx = mesh.getMeshSpacing(); - const vector_type& origin = mesh.getOrigin(); - const vector_type invdx = 1.0 / dx; - - const FieldLayout& layout = f.getLayout(); - const NDIndex& lDom = layout.getLocalNDIndex(); - const int nghost = f.getNghost(); - - // using policy_type = Kokkos::RangePolicy; - const bool useHashView = hash_array.extent(0) > 0; - if (useHashView && (iteration_policy.end() > hash_array.extent(0))) { + const auto hashExtent = static_cast(hash_array.extent(0)); + const bool useHashView = hashExtent > 0; + if (useHashView && (iteration_policy.end() > hashExtent)) { Inform m("scatter"); m << "Hash array was passed to scatter, but size does not match iteration policy." << endl; ippl::Comm->abort(); } - auto dview = dview_m; - auto ppview = pp.getView(); - Kokkos::parallel_for( - "ParticleAttrib::scatter", iteration_policy, KOKKOS_LAMBDA(const size_t idx) { - // map index to possible hash_map - size_t mapped_idx = useHashView ? hash_array(idx) : idx; - - vector_type l = (ppview(mapped_idx) - origin) * invdx + 0.5; - Vector index = l; - Vector whi = l - index; - Vector wlo = 1.0 - whi; - Vector args = index - lDom.first() + nghost; - - const value_type& val = dview(mapped_idx); - detail::scatterToField(std::make_index_sequence<1 << Field::dim>{}, view, wlo, whi, - args, val); - }); - IpplTimings::stopTimer(scatterTimer); - - static IpplTimings::TimerRef accumulateHaloTimer = IpplTimings::getTimer("accumulateHalo"); - IpplTimings::startTimer(accumulateHaloTimer); - f.accumulateHalo(); - IpplTimings::stopTimer(accumulateHaloTimer); + if (useHashView) { + detail::particleAttribScatterImpl(f, dview_m, pp, iteration_policy, hash_array); + } else { + detail::particleAttribScatterImpl(f, dview_m, pp, iteration_policy, hash_array); + } } template @@ -410,8 +435,7 @@ namespace ippl { Field& f, Field& Sk, const ParticleAttrib, Properties...>& pp, FFT>* nufft, ParticleAttrib& q) { - static IpplTimings::TimerRef gatherPIFNUFFTTimer = - IpplTimings::getTimer("GatherPIFNUFFT"); + static IpplTimings::TimerRef gatherPIFNUFFTTimer = IpplTimings::getTimer("GatherPIFNUFFT"); IpplTimings::startTimer(gatherPIFNUFFTTimer); typename Field::uniform_type tempField; @@ -452,8 +476,7 @@ namespace ippl { "Gather NUFFT", mdrange_type( {nghost, nghost, nghost}, - {fview.extent(0) - nghost, fview.extent(1) - nghost, - fview.extent(2) - nghost}), + {fview.extent(0) - nghost, fview.extent(1) - nghost, fview.extent(2) - nghost}), KOKKOS_LAMBDA(const int i, const int j, const int k) { Vector iVec = {i, j, k}; for (unsigned d = 0; d < Dim; ++d) { diff --git a/test/particle/TestScatter.cpp b/test/particle/TestScatter.cpp index bfa11e9daf..8ca2403a13 100644 --- a/test/particle/TestScatter.cpp +++ b/test/particle/TestScatter.cpp @@ -6,13 +6,17 @@ template struct Bunch : public ippl::ParticleBase { Bunch(PLayout& playout) : ippl::ParticleBase(playout) { - this->addAttribute(Q); + this->addAttribute(Q1); + this->addAttribute(Q2); } ~Bunch() {} - typedef ippl::ParticleAttrib charge_container_type; - charge_container_type Q; + typedef ippl::ParticleAttrib charge_container_typeF; + charge_container_typeF Q1; + + typedef ippl::ParticleAttrib charge_container_typeD; + charge_container_typeD Q2; }; int main(int argc, char* argv[]) { @@ -85,26 +89,46 @@ int main(int argc, char* argv[]) { std::cout << "Sum coord: " << global_sum_coord << std::endl; } - bunch.Q = 1.0; + bunch.Q1 = 1.0; + + bunch.update(); + + field = 0.0; + + scatter(bunch.Q1, field, bunch.R); + + // Check charge conservation + try { + double Total_charge_field = field.sum(); + + std::cout << "Float:: Total charge in the field:" << Total_charge_field << std::endl; + std::cout << "Float:: Total charge of the particles:" << bunch.Q1.sum() << std::endl; + std::cout << "Float:: Error:" << std::fabs(bunch.Q1.sum() - Total_charge_field) + << std::endl; + } catch (const std::exception& e) { + std::cout << e.what() << std::endl; + } + + bunch.Q2 = 1.0; bunch.update(); field = 0.0; - scatter(bunch.Q, field, bunch.R); + scatter(bunch.Q2, field, bunch.R); // Check charge conservation try { double Total_charge_field = field.sum(); - std::cout << "Total charge in the field:" << Total_charge_field << std::endl; - std::cout << "Total charge of the particles:" << bunch.Q.sum() << std::endl; - std::cout << "Error:" << std::fabs(bunch.Q.sum() - Total_charge_field) << std::endl; + std::cout << "Double:: Total charge in the field:" << Total_charge_field << std::endl; + std::cout << "Double:: Total charge of the particles:" << bunch.Q2.sum() << std::endl; + std::cout << "Double:: Error:" << std::fabs(bunch.Q2.sum() - Total_charge_field) + << std::endl; } catch (const std::exception& e) { std::cout << e.what() << std::endl; } } ippl::finalize(); - return 0; } diff --git a/unit_tests/Particle/GatherScatterTest.cpp b/unit_tests/Particle/GatherScatterTest.cpp index 5a05f16189..387ce03ec0 100644 --- a/unit_tests/Particle/GatherScatterTest.cpp +++ b/unit_tests/Particle/GatherScatterTest.cpp @@ -9,41 +9,60 @@ // #include "Ippl.h" -#include "TestUtils.h" -#include "gtest/gtest.h" -#include -#include -#include #include +#include #include +#include +#include +#include +#include + +#include "TestUtils.h" +#include "gtest/gtest.h" // A helper needed to reduce over a hash_type. // This is needed, since Kokkos kernels apparently // cannot be called inside a TYPED_TEST on device. +template struct ComputeTotalChargeLambda { - Kokkos::View viewQ; - Kokkos::View hash; + ViewType viewQ; + HashType hash; - ComputeTotalChargeLambda(Kokkos::View viewQ_, Kokkos::View hash_) - : viewQ(viewQ_), hash(hash_) {} + ComputeTotalChargeLambda(ViewType viewQ_, HashType hash_) + : viewQ(viewQ_) + , hash(hash_) {} KOKKOS_INLINE_FUNCTION void operator()(const size_t i, double& val) const { - val += viewQ(hash(i)); + val += static_cast(viewQ(hash(i))); } }; -// A simple bunch_type holding a charge attribute +template +double scatterConservationTolerance(const double reference) { + using tolerance_type = std::conditional_t<(std::numeric_limits::digits + < std::numeric_limits::digits), + ScatterType, FieldType>; + return 100.0 * static_cast(std::numeric_limits::epsilon()) + * std::max(1.0, std::abs(reference)); +} + +// A simple bunch_type holding a charge attribute template struct Bunch : public ippl::ParticleBase { Bunch(PLayout& playout) : ippl::ParticleBase(playout) { this->addAttribute(Q); + this->addAttribute(QFloat); } ~Bunch() = default; - typedef ippl::ParticleAttrib charge_container_type; + typedef ippl::ParticleAttrib + charge_container_type; + typedef ippl::ParticleAttrib + mixed_charge_container_type; charge_container_type Q; + mixed_charge_container_type QFloat; }; template @@ -52,30 +71,30 @@ class GatherScatterTest; template class GatherScatterTest>> : public ::testing::Test { public: - using scalar_type = T; - using exec_space = ExecSpace; + using scalar_type = T; + using exec_space = ExecSpace; static const unsigned dim = Dim; - using flayout_type = ippl::FieldLayout; - using mesh_type = ippl::UniformCartesian; - using playout_type = ippl::ParticleSpatialLayout; - using bunch_type = Bunch; + using flayout_type = ippl::FieldLayout; + using mesh_type = ippl::UniformCartesian; + using playout_type = ippl::ParticleSpatialLayout; + using bunch_type = Bunch; // Domain parameters: use a high resolution grid so that cells are small. std::array nPoints; std::array domain; flayout_type layout; mesh_type mesh; - std::shared_ptr playout; + std::shared_ptr playout; std::shared_ptr bunch; // Particle counts for the tests. - size_t nGather = 10; // for gather test: local particles per rank + size_t nGather = 10; // for gather test: local particles per rank size_t nScatter = static_cast(std::pow(64, Dim)); // for scatter tests // Store cell sizes (hx) for use in generating positions. T hx[Dim]; - GatherScatterTest() { } + GatherScatterTest() {} void SetUp() override { // Use a high-resolution grid (e.g. 512 cells per dimension) @@ -91,13 +110,13 @@ class GatherScatterTest>> : public ::testing: std::array isParallel; isParallel.fill(true); auto owned_tu = std::make_from_tuple>(owned); - layout = flayout_type(MPI_COMM_WORLD, owned_tu, isParallel); + layout = flayout_type(MPI_COMM_WORLD, owned_tu, isParallel); ippl::Vector hx_vec; ippl::Vector origin; for (size_t d = 0; d < Dim; d++) { hx_vec[d] = domain[d] / nPoints[d]; - hx[d] = hx_vec[d]; // store cell size for distribution + hx[d] = hx_vec[d]; // store cell size for distribution origin[d] = 0; } mesh = mesh_type(owned_tu, hx_vec, origin); @@ -136,32 +155,42 @@ class GatherScatterTest>> : public ::testing: Kokkos::deep_copy(bunch->Q.getView(), Q_host); ippl::Comm->barrier(); } + + void fillAttributeQFloat(float value) { + auto Q_host = bunch->QFloat.getHostMirror(); + for (size_t i = 0; i < Q_host.size(); ++i) { + Q_host(i) = value; + } + Kokkos::deep_copy(bunch->QFloat.getView(), Q_host); + ippl::Comm->barrier(); + } }; -using TestTypes = ::testing::Types< - Parameters>, - Parameters>, - Parameters>//, - //Parameters>, - //Parameters>, - //Parameters> ->; +using TestTypes = ::testing::Types>, + Parameters>, + Parameters> //, + // Parameters>, + // Parameters>, + // Parameters> + >; TYPED_TEST_SUITE(GatherScatterTest, TestTypes); // -// GatherTest: +// GatherTest: // First, set each local Q to 10.0. -// Then, call gather with addToAttribute = false so that Q becomes 1.0 (should replace value with 1.0). -// If Q != 0, then the values were 1. not replaced and 2. not correctly gathered from the field. -// Note: for a constant field, there should not be an error during linear interpolation. +// Then, call gather with addToAttribute = false so that Q becomes 1.0 (should replace value +// with 1.0). If Q != 0, then the values were 1. not replaced and 2. not correctly gathered from the +// field. Note: for a constant field, there should not be an error during linear interpolation. // TYPED_TEST(GatherScatterTest, GatherTestReplace) { const size_t n = this->nGather; this->fillRandomPositions(n); this->fillAttributeQ(10.0); - using Mesh_t = typename TestFixture::mesh_type; - using FieldType = ippl::Field; + using Mesh_t = typename TestFixture::mesh_type; + using FieldType = + ippl::Field; FieldType field; field.initialize(this->mesh, this->layout); field = 1.0; @@ -178,7 +207,7 @@ TYPED_TEST(GatherScatterTest, GatherTestReplace) { } // -// GatherTest: +// GatherTest: // First, set each local Q to 1.0. // Then, call gather with addToAttribute = true so that Q becomes 2.0 (should add 1.0 per particle). // @@ -187,8 +216,10 @@ TYPED_TEST(GatherScatterTest, GatherTestIncrement) { this->fillRandomPositions(n); this->fillAttributeQ(1.0); - using Mesh_t = typename TestFixture::mesh_type; - using FieldType = ippl::Field; + using Mesh_t = typename TestFixture::mesh_type; + using FieldType = + ippl::Field; FieldType field; field.initialize(this->mesh, this->layout); field = 1.0; @@ -216,46 +247,80 @@ TYPED_TEST(GatherScatterTest, ScatterSimpleTest) { this->fillAttributeQ(1.0); // Create and initialize a field/mesh. - using Mesh_t = typename TestFixture::mesh_type; - using FieldType = ippl::Field; + using Mesh_t = typename TestFixture::mesh_type; + using FieldType = + ippl::Field; FieldType field; field.initialize(this->mesh, this->layout); - + field = 0.0; // Perform the simple scatter operation (extended functionality is tested below). scatter(this->bunch->Q, field, this->bunch->R); // Compute the total charge in the field and from the particles. - double total_field = field.sum(); + double total_field = field.sum(); double total_particles = this->bunch->Q.sum(); // Check that the scattered field conserves charge. ASSERT_NEAR(total_field, total_particles, 1e-6); } +// +// ScatterMixedValueTypeTest: +// Scatter a lower-precision particle attribute into a higher-precision field. +// This covers ParticleAttrib -> Field, which is the mixed type +// use case enabled by ParticleAttrib::scatter accepting a field value type that +// differs from the attribute value type. +// +TYPED_TEST(GatherScatterTest, ScatterMixedValueTypeTest) { + using ScatterType = float; + using FieldType = double; + + const unsigned int n = this->nScatter; + this->fillRandomPositions(n); + this->fillAttributeQFloat(ScatterType(1.0)); + + using Mesh_t = typename TestFixture::mesh_type; + using Field = ippl::Field; + Field field; + field.initialize(this->mesh, this->layout); + field = FieldType(0.0); + + scatter(this->bunch->QFloat, field, this->bunch->R); + + const double total_field = field.sum(); + const double total_particles = static_cast(this->bunch->QFloat.sum()); + const double tolerance = scatterConservationTolerance(total_particles); + + ASSERT_NEAR(total_field, total_particles, tolerance); +} // -// ScatterCustomRangeTest: +// ScatterCustomRangeTest: // Set Q = 1.0 for all particles and scatter only a subset defined by a custom range policy. // Then compare the total charge in the field to the expected value. // TYPED_TEST(GatherScatterTest, ScatterCustomRangeTest) { const size_t n = this->nScatter; - if(n % ippl::Comm->size() != 0) { + if (n % ippl::Comm->size() != 0) { GTEST_SKIP() << "nScatter not divisible by number of ranks."; } this->fillRandomPositions(n); this->fillAttributeQ(1.0); - using Mesh_t = typename TestFixture::mesh_type; - using FieldType = ippl::Field; + using Mesh_t = typename TestFixture::mesh_type; + using FieldType = + ippl::Field; FieldType field; field.initialize(this->mesh, this->layout); field = 0.0; - size_t rank = ippl::Comm->rank(); - size_t nLoc = this->bunch->getLocalNum(); + size_t rank = ippl::Comm->rank(); + size_t nLoc = this->bunch->getLocalNum(); size_t NScattered = nLoc / 2 + rank; double Q_total = 1.0 * NScattered; @@ -269,7 +334,7 @@ TYPED_TEST(GatherScatterTest, ScatterCustomRangeTest) { } // -// ScatterCustomHashTest: +// ScatterCustomHashTest: // Assign random charges (in [0.5, 1.5]), create and shuffle an index array, // use it as a custom hash, scatter the first NScattered particles accordingly, // and compare the field’s total charge to the expected total. @@ -280,10 +345,10 @@ TYPED_TEST(GatherScatterTest, ScatterCustomHashTest) { GTEST_SKIP() << "nScatter not divisible by number of ranks."; } this->fillRandomPositions(n); - + size_t rank = ippl::Comm->rank(); - size_t nLoc = this->bunch->getLocalNum(); // since update() might change number of particles - size_t NScattered = nLoc / 2 + rank; // can be anything + size_t nLoc = this->bunch->getLocalNum(); // since update() might change number of particles + size_t NScattered = nLoc / 2 + rank; // can be anything // Assign random charges to particles std::mt19937_64 eng(42); @@ -295,8 +360,10 @@ TYPED_TEST(GatherScatterTest, ScatterCustomHashTest) { Kokkos::deep_copy(this->bunch->Q.getView(), Q_host); // Create and initialize a field/mesh. - using Mesh_t = typename TestFixture::mesh_type; - using FieldType = ippl::Field; + using Mesh_t = typename TestFixture::mesh_type; + using FieldType = + ippl::Field; FieldType field; field.initialize(this->mesh, this->layout); field = 0.0; @@ -315,15 +382,16 @@ TYPED_TEST(GatherScatterTest, ScatterCustomHashTest) { } Kokkos::deep_copy(hash, hash_host); - // First compute the total charge of the first NScattered particles as determined by the hash map + // First compute the total charge of the first NScattered particles as determined by the hash + // map double Q_total = 0.0; - auto viewQ = this->bunch->Q.getView(); + auto viewQ = this->bunch->Q.getView(); - ComputeTotalChargeLambda lambda(viewQ, hash); - Kokkos::parallel_reduce("computeTotalCharge", - Kokkos::RangePolicy(0, NScattered), - lambda, Q_total); - /*Kokkos::parallel_reduce("computeTotalCharge", + ComputeTotalChargeLambda lambda(viewQ, hash); + Kokkos::parallel_reduce("computeTotalCharge", + Kokkos::RangePolicy(0, NScattered), + lambda, Q_total); + /*Kokkos::parallel_reduce("computeTotalCharge", Kokkos::RangePolicy(0, NScattered), KOKKOS_LAMBDA(const size_t i, double& val) { val += viewQ(hash(i)); @@ -339,6 +407,69 @@ TYPED_TEST(GatherScatterTest, ScatterCustomHashTest) { ASSERT_NEAR(Q_total, Total_charge_field, 1e-6); } +// +// ScatterMixedValueTypeCustomHashTest: +// Scatter ParticleAttrib into Field through a custom hash array. +// This verifies that both the mixed value type and the hash-dispatched scatter +// path conserve the charge represented by the lower-precision attribute. +// +TYPED_TEST(GatherScatterTest, ScatterMixedValueTypeCustomHashTest) { + using ScatterType = float; + using FieldType = double; + + const size_t n = this->nScatter / ippl::Comm->size(); + if (this->nScatter % ippl::Comm->size() > 0) { + GTEST_SKIP() << "nScatter not divisible by number of ranks."; + } + this->fillRandomPositions(n); + + const size_t rank = ippl::Comm->rank(); + const size_t nLoc = this->bunch->getLocalNum(); + const size_t NScattered = nLoc / 2 + rank; + + std::mt19937_64 eng(42); + std::uniform_real_distribution unif_charge(0.5f, 1.5f); + auto Q_host = this->bunch->QFloat.getHostMirror(); + for (size_t i = 0; i < nLoc; ++i) { + Q_host(i) = unif_charge(eng); + } + Kokkos::deep_copy(this->bunch->QFloat.getView(), Q_host); + + using Mesh_t = typename TestFixture::mesh_type; + using Field = ippl::Field; + Field field; + field.initialize(this->mesh, this->layout); + field = FieldType(0.0); + + using hash_type = typename TestFixture::bunch_type::mixed_charge_container_type::hash_type; + hash_type hash("mixedIndexArray", nLoc); + std::vector host_indices(nLoc); + std::iota(host_indices.begin(), host_indices.end(), 0); + std::shuffle(host_indices.begin(), host_indices.end(), eng); + + auto hash_host = Kokkos::create_mirror_view(hash); + for (size_t i = 0; i < nLoc; ++i) { + hash_host(i) = host_indices[i]; + } + Kokkos::deep_copy(hash, hash_host); + + double Q_total = 0.0; + auto viewQ = this->bunch->QFloat.getView(); + ComputeTotalChargeLambda lambda(viewQ, hash); + Kokkos::parallel_reduce("computeMixedTotalCharge", + Kokkos::RangePolicy(0, NScattered), + lambda, Q_total); + ippl::Comm->allreduce(Q_total, 1, std::plus()); + + Kokkos::RangePolicy policy(0, NScattered); + scatter(this->bunch->QFloat, field, this->bunch->R, policy, hash); + + const double Total_charge_field = field.sum(); + const double tolerance = scatterConservationTolerance(Q_total); + ASSERT_NEAR(Q_total, Total_charge_field, tolerance); +} + int main(int argc, char* argv[]) { ippl::initialize(argc, argv); int result = 1;