diff --git a/cosmology/CMakeLists.txt b/cosmology/CMakeLists.txt index df5df0b1e..66c101726 100644 --- a/cosmology/CMakeLists.txt +++ b/cosmology/CMakeLists.txt @@ -1,22 +1,24 @@ file (RELATIVE_PATH _relPath "${CMAKE_SOURCE_DIR}" "${CMAKE_CURRENT_SOURCE_DIR}") -message (STATUS "Adding index test found in ${_relPath}") +message (STATUS "Adding StructureFormation ${_relPath}") include_directories ( ${CMAKE_SOURCE_DIR}/src + ${CMAKE_SOURCE_DIR}/cosmology/mc-4-Initializer ) - - link_directories ( ${CMAKE_CURRENT_SOURCE_DIR} ${Kokkos_DIR}/.. ) - set (IPPL_LIBS ippl ${MPI_CXX_LIBRARIES}) -set (COMPILE_FLAGS ${OPAL_CXX_FLAGS}) +set (COMPILE_FLAGS ${OPAL_CXX_FLAGS}) + +add_definitions(-DUSENAMESPACE) + +# Add the main executable and link the object file +add_executable(StructureFormation StructureFormation.cpp mc-4-Initializer/Cosmology.cpp mc-4-Initializer/DataBase.cpp mc-4-Initializer/InputParser.cpp) -add_executable (StructureFormation StructureFormation.cpp) target_link_libraries (StructureFormation ${IPPL_LIBS}) diff --git a/cosmology/GravityLoadBalancer.hpp b/cosmology/GravityLoadBalancer.hpp index 2bcc7a6f8..dbfff961c 100644 --- a/cosmology/GravityLoadBalancer.hpp +++ b/cosmology/GravityLoadBalancer.hpp @@ -147,26 +147,23 @@ 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; - } - } - return false; - } + 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; } private: @@ -180,4 +177,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/cosmology/GravityManager.h b/cosmology/GravityManager.h index 80ca4fd87..222b25eca 100644 --- a/cosmology/GravityManager.h +++ b/cosmology/GravityManager.h @@ -13,6 +13,8 @@ #include "Random/NormalDistribution.h" #include "Random/Randn.h" +#include "mc-4-Initializer/InputParser.h" + using view_type = typename ippl::detail::ViewType, 1>::view_type; /** @@ -41,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_) @@ -64,6 +68,12 @@ 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. * @@ -210,8 +220,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)); @@ -240,8 +252,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; } /** @@ -261,8 +273,8 @@ 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 + mes << "Step: " << this->it_m; + mes << " comological time: " << this->time_m << ", dt: " << this->dt_m << ", z: " << this->z_m << ", a: " << this->a_m << endl; } @@ -290,10 +302,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; @@ -302,7 +314,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(); @@ -310,7 +321,7 @@ 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; @@ -350,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] @@ -367,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/cosmology/GravityParticleContainer.hpp b/cosmology/GravityParticleContainer.hpp index ab000c6e4..dcb3a0f24 100644 --- a/cosmology/GravityParticleContainer.hpp +++ b/cosmology/GravityParticleContainer.hpp @@ -19,7 +19,7 @@ class ParticleContainer : public ippl::ParticleBase m; + ippl::ParticleAttrib m; /** * @brief Velocity of the particle. @@ -89,4 +89,4 @@ class ParticleContainer : public ippl::ParticleBasesetParticleBC(ippl::BC::PERIODIC); } }; -#endif // IPPL_PARTICLE_CONTAINER_H \ No newline at end of file +#endif // IPPL_PARTICLE_CONTAINER_H diff --git a/cosmology/StructureFormation.cpp b/cosmology/StructureFormation.cpp index 07b71babf..51c8f2e4d 100644 --- a/cosmology/StructureFormation.cpp +++ b/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 -// --info 5 +// 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,13 +81,14 @@ 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; template @@ -57,26 +97,44 @@ 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++]; - // Number of gridpoints in each dimension + 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++]; // Check if the solver type is valid @@ -89,7 +147,7 @@ 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 +155,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/cosmology/StructureFormationManager.h b/cosmology/StructureFormationManager.h index 50ea33260..d761e370c 100644 --- a/cosmology/StructureFormationManager.h +++ b/cosmology/StructureFormationManager.h @@ -20,8 +20,36 @@ #include "Random/NormalDistribution.h" #include "Random/Randn.h" +#include "mc-4-Initializer/InputParser.h" +#include "mc-4-Initializer/DataBase.h" +#include "mc-4-Initializer/Cosmology.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 +59,40 @@ 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 { -public: - using ParticleContainer_t = ParticleContainer; - using FieldContainer_t = FieldContainer; - using FieldSolver_t = FieldSolver; - using LoadBalancer_t = LoadBalancer; - StructureFormationManager(size_type totalP_, int nt_, Vector_t& nr_, double lbt_, - std::string& solver_, std::string& stepMethod_) - : GravityManager(totalP_, nt_, nr_, lbt_, solver_, stepMethod_) {} + /// 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_, + 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 +103,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 +118,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(); - this->rmin_m = 0.0; - this->rmax_m = 50000.0; // kpc/h + float box_size; + this->parser_m.getByName("box_size", box_size); + this->rmin_m = 0.0; + 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,43 +162,452 @@ 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 - - static IpplTimings::TimerRef DummySolveTimer = IpplTimings::getTimer("solveWarmup"); - IpplTimings::startTimer(DummySolveTimer); - - this->fcontainer_m->getRho() = 0.0; - - this->fsolver_m->runSolver(); + 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); + + this->fcontainer_m->getRho() = 0.0; + + 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); + } - IpplTimings::stopTimer(DummySolveTimer); - this->par2grid(); - static IpplTimings::TimerRef SolveTimer = IpplTimings::getTimer("solve"); + 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(); + } + } - IpplTimings::startTimer(SolveTimer); - this->fsolver_m->runSolver(); - IpplTimings::stopTimer(SolveTimer); + /** + * @brief Create particles using Zarijas initializer method + */ - this->grid2par(); - this->dump(); + 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 + */ - mes << "Done"; + 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 +618,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 +677,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 +690,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 +725,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 +740,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 +788,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 +854,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 +864,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 +948,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) @@ -502,6 +970,9 @@ class StructureFormationManager : public GravityManager { 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 +987,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/cosmology/mc-4-Initializer/Cosmology.cpp b/cosmology/mc-4-Initializer/Cosmology.cpp new file mode 100644 index 000000000..e41772107 --- /dev/null +++ b/cosmology/mc-4-Initializer/Cosmology.cpp @@ -0,0 +1,742 @@ +/* + 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 "TypesAndDefs.h" +#include +#include +#include +#include +#include +#include +#include "TypesAndDefs.h" +#include "DataBase.h" +#include "Cosmology.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_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*derivs)(x, y, dydx); + ++feval; + for (i=0; i 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 (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= 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*derivs)(x+a2*h,ytemp,ak2); + for (i=0; i*derivs)(x+a3*h,ytemp,ak3); + for (i=0; i*derivs)(x+a4*h,ytemp,ak4); + for (i=0; i*derivs)(x+a5*h,ytemp,ak5); + for (i=0; i*derivs)(x+a6*h,ytemp,ak6); + for (i=0; i*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 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/cosmology/mc-4-Initializer/Cosmology.h b/cosmology/mc-4-Initializer/Cosmology.h new file mode 100644 index 000000000..20ff34997 --- /dev/null +++ b/cosmology/mc-4-Initializer/Cosmology.h @@ -0,0 +1,86 @@ +/* + 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/cosmology/mc-4-Initializer/DataBase.cpp b/cosmology/mc-4-Initializer/DataBase.cpp new file mode 100644 index 000000000..12ad12503 --- /dev/null +++ b/cosmology/mc-4-Initializer/DataBase.cpp @@ -0,0 +1,130 @@ +/* + 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 +#include +#include "TypesAndDefs.h" +#include "Ippl.h" +#include "DataBase.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/cosmology/mc-4-Initializer/DataBase.h b/cosmology/mc-4-Initializer/DataBase.h new file mode 100644 index 000000000..5cc588c93 --- /dev/null +++ b/cosmology/mc-4-Initializer/DataBase.h @@ -0,0 +1,74 @@ +/* + 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 "TypesAndDefs.h" +#include "InputParser.h" +#include + +#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/cosmology/mc-4-Initializer/InputParser.cpp b/cosmology/mc-4-Initializer/InputParser.cpp new file mode 100644 index 000000000..681c80351 --- /dev/null +++ b/cosmology/mc-4-Initializer/InputParser.cpp @@ -0,0 +1,82 @@ +#include "InputParser.h" + + +#ifdef USENAMESPACE +namespace initializer { +#endif + +InputParser::InputParser(std::string filename) +{ + filename_m = filename; + parseFile(); +} + + +bool InputParser::getByName(std::string name, int ¶m) +{ + 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 ¶m) +{ + 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/cosmology/mc-4-Initializer/InputParser.h b/cosmology/mc-4-Initializer/InputParser.h new file mode 100644 index 000000000..f88e805d5 --- /dev/null +++ b/cosmology/mc-4-Initializer/InputParser.h @@ -0,0 +1,39 @@ +#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 ¶m); + bool getByName(std::string name, real ¶m); + +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/cosmology/mc-4-Initializer/Output.cpp b/cosmology/mc-4-Initializer/Output.cpp new file mode 100644 index 000000000..f5868d107 --- /dev/null +++ b/cosmology/mc-4-Initializer/Output.cpp @@ -0,0 +1,340 @@ +/* + * Output.cpp + * Initializer + * + * Created by Zarija on 1/27/10. + * + */ + +#include +#include +#include +#include +#include +#include +#include +#include "TypesAndDefs.h" +#include "Output.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 +#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/cosmology/mc-4-Initializer/TypesAndDefs.h b/cosmology/mc-4-Initializer/TypesAndDefs.h new file mode 100644 index 000000000..5abd66b6c --- /dev/null +++ b/cosmology/mc-4-Initializer/TypesAndDefs.h @@ -0,0 +1,44 @@ +/* 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 c00aa696b..069c70a31 100644 --- a/src/Decomposition/OrthogonalRecursiveBisection.hpp +++ b/src/Decomposition/OrthogonalRecursiveBisection.hpp @@ -1,4 +1,5 @@ #include "Utility/IpplTimings.h" +#include namespace ippl { @@ -288,9 +289,9 @@ namespace ippl { Vector wlo = 1.0 - whi; 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 5adac510d..4ad5117ea 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 cbe8064f1..6080128fb 100644 --- a/src/Interpolation/CIC.hpp +++ b/src/Interpolation/CIC.hpp @@ -43,14 +43,21 @@ 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) { - // The number of indices is equal to the view rank + 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.hpp b/src/Particle/ParticleAttrib.hpp index 78df54058..a38c34210 100644 --- a/src/Particle/ParticleAttrib.hpp +++ b/src/Particle/ParticleAttrib.hpp @@ -151,7 +151,7 @@ namespace ippl { // scatter const value_type& val = dview_m(mapped_idx); - detail::scatterToField(std::make_index_sequence<1 << Field::dim>{}, view, wlo, whi, + detail::scatterToField(std::make_index_sequence<1 << Field::dim>{}, view, wlo, whi, args, val); }); IpplTimings::stopTimer(scatterTimer); diff --git a/test/particle/TestScatter.cpp b/test/particle/TestScatter.cpp index d92f6a9b1..3c78c73da 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[]) { @@ -84,26 +88,44 @@ 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.Q, field, bunch.R); + scatter(bunch.Q1, 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 << "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.Q2, field, bunch.R); + + // Check charge conservation + try { + double Total_charge_field = field.sum(); + + 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; }