From bb466f8590b2765c59ebb12a4f4151a0ee2b49ca Mon Sep 17 00:00:00 2001 From: Andreas Adelmann Date: Fri, 12 Sep 2025 23:56:13 +0300 Subject: [PATCH 01/18] step 1: augment the original code with mc-4 compatible inputfile, change command line --- cosmology/CMakeLists.txt | 14 +- cosmology/GravityLoadBalancer.hpp | 37 +- cosmology/StructureFormation.cpp | 113 +++- cosmology/StructureFormationManager.h | 22 +- cosmology/mc-4-Initializer/Cosmology.cpp | 742 +++++++++++++++++++++ cosmology/mc-4-Initializer/Cosmology.h | 86 +++ cosmology/mc-4-Initializer/DataBase.cpp | 130 ++++ cosmology/mc-4-Initializer/DataBase.h | 74 ++ cosmology/mc-4-Initializer/InputParser.cpp | 82 +++ cosmology/mc-4-Initializer/InputParser.h | 39 ++ cosmology/mc-4-Initializer/Output.cpp | 340 ++++++++++ cosmology/mc-4-Initializer/Output.h | 68 ++ cosmology/mc-4-Initializer/TypesAndDefs.h | 44 ++ 13 files changed, 1732 insertions(+), 59 deletions(-) create mode 100644 cosmology/mc-4-Initializer/Cosmology.cpp create mode 100644 cosmology/mc-4-Initializer/Cosmology.h create mode 100644 cosmology/mc-4-Initializer/DataBase.cpp create mode 100644 cosmology/mc-4-Initializer/DataBase.h create mode 100644 cosmology/mc-4-Initializer/InputParser.cpp create mode 100644 cosmology/mc-4-Initializer/InputParser.h create mode 100644 cosmology/mc-4-Initializer/Output.cpp create mode 100644 cosmology/mc-4-Initializer/Output.h create mode 100644 cosmology/mc-4-Initializer/TypesAndDefs.h diff --git a/cosmology/CMakeLists.txt b/cosmology/CMakeLists.txt index df5df0b1e3..66c1017267 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 2bcc7a6f8b..dbfff961c3 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/StructureFormation.cpp b/cosmology/StructureFormation.cpp index 07b71babf4..ca03e85ee6 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" @@ -49,6 +88,9 @@ const char* TestName = "StructureFormation"; #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 +99,37 @@ 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(argv[0]); + Inform msg2all(argv[0], 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); + // 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++]); + size_type totalP = nr[0]*nr[1]*nr[2]; + // 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 +142,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 +150,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 50ea332606..bc2d20053a 100644 --- a/cosmology/StructureFormationManager.h +++ b/cosmology/StructureFormationManager.h @@ -20,6 +20,13 @@ #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 KOKKOS_PRINT // Kokkos::printf of interesting quantities. Does not work multirank + using view_type = typename ippl::detail::ViewType, 1>::view_type; /** @@ -31,9 +38,15 @@ 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 { + + bool readICs_m; + public: using ParticleContainer_t = ParticleContainer; using FieldContainer_t = FieldContainer; @@ -41,8 +54,11 @@ class StructureFormationManager : public GravityManager { 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_) {} + std::string& solver_, std::string& stepMethod_, + initializer::InputParser par_, std::string tfname_, bool readICs_) + : GravityManager(totalP_, nt_, nr_, lbt_, solver_, stepMethod_), + readICs_m(readICs_) + {} /** * @brief Destructor for StructureFormationManager. @@ -396,7 +412,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!"); } } diff --git a/cosmology/mc-4-Initializer/Cosmology.cpp b/cosmology/mc-4-Initializer/Cosmology.cpp new file mode 100644 index 0000000000..e417721072 --- /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 0000000000..20ff349976 --- /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 0000000000..12ad12503b --- /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 0000000000..5cc588c936 --- /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 0000000000..681c80351b --- /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 0000000000..f88e805d58 --- /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 0000000000..f5868d107b --- /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 0000000000..5abd66b6c0 --- /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 From de2c42609bfaf4c805879775f514811b1f6771b0 Mon Sep 17 00:00:00 2001 From: Andreas Adelmann Date: Sat, 13 Sep 2025 00:29:14 +0300 Subject: [PATCH 02/18] make output more to the point --- cosmology/GravityManager.h | 15 +++++++-------- 1 file changed, 7 insertions(+), 8 deletions(-) diff --git a/cosmology/GravityManager.h b/cosmology/GravityManager.h index 80ca4fd87a..776a6bb63d 100644 --- a/cosmology/GravityManager.h +++ b/cosmology/GravityManager.h @@ -240,8 +240,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 +261,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,7 +290,7 @@ 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; @@ -302,7 +302,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 +309,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; @@ -367,4 +366,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 From b53bf914b643d94315dafa5c94c0adfb7a85c33e Mon Sep 17 00:00:00 2001 From: Andreas Adelmann Date: Sat, 13 Sep 2025 00:52:16 +0300 Subject: [PATCH 03/18] cleanup and read all parameters from inutfile --- cosmology/GravityManager.h | 22 +++++++++++++++++----- cosmology/StructureFormation.cpp | 2 -- cosmology/StructureFormationManager.h | 2 +- 3 files changed, 18 insertions(+), 8 deletions(-) diff --git a/cosmology/GravityManager.h b/cosmology/GravityManager.h index 776a6bb63d..939be52663 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)); @@ -349,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] diff --git a/cosmology/StructureFormation.cpp b/cosmology/StructureFormation.cpp index ca03e85ee6..88b51292aa 100644 --- a/cosmology/StructureFormation.cpp +++ b/cosmology/StructureFormation.cpp @@ -81,8 +81,6 @@ using T = double; #include #include -#include "Manager/datatypes.h" - #include "Utility/IpplTimings.h" #include "Manager/PicManager.h" diff --git a/cosmology/StructureFormationManager.h b/cosmology/StructureFormationManager.h index bc2d20053a..2ca0654e2c 100644 --- a/cosmology/StructureFormationManager.h +++ b/cosmology/StructureFormationManager.h @@ -56,7 +56,7 @@ class StructureFormationManager : public GravityManager { 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_), + : GravityManager(totalP_, nt_, nr_, lbt_, solver_, stepMethod_, par_), readICs_m(readICs_) {} From 73669bd853f8492d7426e64e1b340869ca6832c3 Mon Sep 17 00:00:00 2001 From: Andreas Adelmann Date: Sat, 13 Sep 2025 13:23:45 +0300 Subject: [PATCH 04/18] read all constants from inputfile --- cosmology/StructureFormationManager.h | 88 ++++++++++++++++----------- 1 file changed, 52 insertions(+), 36 deletions(-) diff --git a/cosmology/StructureFormationManager.h b/cosmology/StructureFormationManager.h index 2ca0654e2c..cac46d9c48 100644 --- a/cosmology/StructureFormationManager.h +++ b/cosmology/StructureFormationManager.h @@ -47,18 +47,25 @@ class StructureFormationManager : public GravityManager { bool readICs_m; + initializer::CosmoClass cosmo_m; + 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_, - initializer::InputParser par_, std::string tfname_, bool readICs_) - : GravityManager(totalP_, nt_, nr_, lbt_, solver_, stepMethod_, par_), - readICs_m(readICs_) - {} + 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_) + { + cosmo_m.SetParameters(initializer::GlobalStuff::instance(), tfname_.c_str()); + } /** * @brief Destructor for StructureFormationManager. @@ -69,7 +76,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", @@ -84,23 +91,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; @@ -145,17 +161,17 @@ class StructureFormationManager : public GravityManager { this->grid2par(); this->dump(); - mes << "Done"; + msg << "Done"; } /** * @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; @@ -166,7 +182,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); @@ -225,9 +241,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()) { @@ -238,7 +254,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(); @@ -273,14 +289,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; @@ -288,7 +304,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); @@ -336,12 +352,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]; } @@ -402,7 +418,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; } /** @@ -496,12 +512,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) @@ -548,7 +564,7 @@ class StructureFormationManager : public GravityManager { // Close the file stream file.close(); - mes << "done." << endl; + msg << "done." << endl; IpplTimings::stopTimer(SavingTimer); } From 2a2e3729732e36e21914cc8fbcf14a52d0a058c0 Mon Sep 17 00:00:00 2001 From: Andreas Adelmann Date: Sun, 14 Sep 2025 22:10:56 +0300 Subject: [PATCH 05/18] Add structures for creating particles --- cosmology/StructureFormationManager.h | 64 ++++++++++++++++++--------- 1 file changed, 43 insertions(+), 21 deletions(-) diff --git a/cosmology/StructureFormationManager.h b/cosmology/StructureFormationManager.h index cac46d9c48..b83711cd15 100644 --- a/cosmology/StructureFormationManager.h +++ b/cosmology/StructureFormationManager.h @@ -29,6 +29,16 @@ using view_type = typename ippl::detail::ViewType, 1>::view_type; +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; + +struct HermitianPkg { + int kx, ky, kz; + double re, im; +}; + /** * @brief Construct a new StructureFormationManager object. * @@ -45,6 +55,10 @@ using view_type = typename ippl::detail::ViewType, 1>: 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; @@ -140,28 +154,36 @@ class StructureFormationManager : public GravityManager { 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(); - - 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"; + msg << "About to generate ... readICs_m= " << readICs_m << endl; + + if (readICs_m) { + msg << "Read in particles ..." << endl; + readParticlesDomain(); // defines particle positions, velocities + msg << "Read particles done" << endl; + 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 reading initial conditions"; + } else { + msg << "Create Particles" << endl; + msg << "Done creating initial conditions"; + } } /** From 18161c66a5d9b78542f42def74e9ad1a2b876106 Mon Sep 17 00:00:00 2001 From: Andreas Adelmann Date: Wed, 17 Sep 2025 15:08:13 +0300 Subject: [PATCH 06/18] make sure we do not overflow --- cosmology/StructureFormation.cpp | 23 ++++++++++++++++------- 1 file changed, 16 insertions(+), 7 deletions(-) diff --git a/cosmology/StructureFormation.cpp b/cosmology/StructureFormation.cpp index 88b51292aa..0b03c94498 100644 --- a/cosmology/StructureFormation.cpp +++ b/cosmology/StructureFormation.cpp @@ -97,8 +97,8 @@ using Vector_t = ippl::Vector; int main(int argc, char* argv[]) { ippl::initialize(argc, argv); { - Inform msg(argv[0]); - Inform msg2all(argv[0], INFORM_ALL_NODES); + Inform msg("StructureFormation "); + Inform msg2all("StructureFormation ", INFORM_ALL_NODES); static IpplTimings::TimerRef mainTimer = IpplTimings::getTimer("total"); IpplTimings::startTimer(mainTimer); @@ -111,15 +111,24 @@ int main(int argc, char* argv[]) { int arg = 5; initializer::InputParser par(indatName); initializer::GlobalStuff::instance().GetParameters(par); - - // Number of gridpoints in each dimension + + 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] = initializer::GlobalStuff::instance().ngrid; } - // Total number of particles - size_type totalP = nr[0]*nr[1]*nr[2]; + msg << "Np = " << totalP << endl; + // Number of time steps int nt = 0; par.getByName("nt", nt); @@ -150,7 +159,7 @@ int main(int argc, char* argv[]) { msg << "Starting iterations ... up to " << manager.getNt() << endl; - manager.run(manager.getNt()); + // manager.run(manager.getNt()); msg << "End." << endl; From 322bed5ec4fcd92bbdddcb93dc5448493a35d45d Mon Sep 17 00:00:00 2001 From: Andreas Adelmann Date: Wed, 17 Sep 2025 15:18:24 +0300 Subject: [PATCH 07/18] diverse fixes, runs now up to 1024**3 --- cosmology/StructureFormationManager.h | 386 +++++++++++++++++++++++++- 1 file changed, 380 insertions(+), 6 deletions(-) diff --git a/cosmology/StructureFormationManager.h b/cosmology/StructureFormationManager.h index b83711cd15..09afd7d259 100644 --- a/cosmology/StructureFormationManager.h +++ b/cosmology/StructureFormationManager.h @@ -25,7 +25,7 @@ #include "mc-4-Initializer/Cosmology.h" -#define KOKKOS_PRINT // Kokkos::printf of interesting quantities. Does not work multirank +//#define KOKKOS_PRINT // Kokkos::printf of interesting quantities. Does not work multirank using view_type = typename ippl::detail::ViewType, 1>::view_type; @@ -78,6 +78,8 @@ class StructureFormationManager : public GravityManager { : 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()); } @@ -154,9 +156,6 @@ class StructureFormationManager : public GravityManager { this->setLoadBalancer(std::make_shared( this->lbt_m, this->fcontainer_m, this->pcontainer_m, this->fsolver_m)); - - msg << "About to generate ... readICs_m= " << readICs_m << endl; - if (readICs_m) { msg << "Read in particles ..." << endl; readParticlesDomain(); // defines particle positions, velocities @@ -181,11 +180,386 @@ class StructureFormationManager : public GravityManager { this->dump(); msg << "Done reading initial conditions"; } else { - msg << "Create Particles" << endl; - msg << "Done creating initial conditions"; + msg << "Create Particles" << endl; + ippl::ParameterList fftParams; + fftParams.add("use_heffte_defaults", true); + 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"; + createParticles(); } } + + 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(); + + 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 + 0.5) - d_z*x); + vView(n)[dim] = -ddot*x; + + //Kokkos::printf("set_particle: rank= %i \t dim= %i \t d=%i \t n=%i \n", myrank, dim, d, n); + +#ifdef KOKKOS_PRINT + 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; + } + } + + /** + * @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; + pc->create(nloc); + pc->m = this->M_m / this->totalP_m; + IpplTimings::stopTimer(creaIpplParts); + + msg << "pc->create(nloc) done" << endl; + + /** + construct Pk_m + */ + static IpplTimings::TimerRef iniPwrSpec = IpplTimings::getTimer("initPwrSpec"); + IpplTimings::startTimer(iniPwrSpec); + initPwrSpec(); + IpplTimings::stopTimer(iniPwrSpec); + + /** + 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); + + // 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; + this->dumpParticles(); + */ +} + + + + + /** * @brief Read particle data from a file. */ From fc882b674253092a35146b46fdb797a69acafe9c Mon Sep 17 00:00:00 2001 From: Andreas Adelmann Date: Wed, 17 Sep 2025 22:26:01 +0300 Subject: [PATCH 08/18] positions are on the same order of magnitude --- cosmology/StructureFormationManager.h | 50 +++++++++++++++++++++------ 1 file changed, 39 insertions(+), 11 deletions(-) diff --git a/cosmology/StructureFormationManager.h b/cosmology/StructureFormationManager.h index 09afd7d259..c7ead80215 100644 --- a/cosmology/StructureFormationManager.h +++ b/cosmology/StructureFormationManager.h @@ -180,15 +180,37 @@ class StructureFormationManager : public GravityManager { this->dump(); msg << "Done reading initial conditions"; } else { - msg << "Create Particles" << endl; + msg << "Create Particles ....." << endl; ippl::ParameterList fftParams; fftParams.add("use_heffte_defaults", true); 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"; + msg << "FFT and Pk structures initialized" << endl;; createParticles(); + msg << "Create particles done" << endl; + + 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); } @@ -319,6 +341,9 @@ class StructureFormationManager : public GravityManager { 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(); @@ -434,6 +459,7 @@ class StructureFormationManager : public GravityManager { zz = 1.0; break; } + // Compute displacement component in k-space ippl::parallel_for("force_in_k_space", ippl::getRangePolicy(cfview, ngh), @@ -487,10 +513,11 @@ class StructureFormationManager : public GravityManager { const unsigned int kl = idx[2] - ngh; const index_type n = il + lDom[0].length()*(jl + kl*lDom[1].length()); - rView(n)[dim] = ((d + 0.5) - d_z*x); + rView(n)[dim] = (d*hr[dim]) - (d_z*x); vView(n)[dim] = -ddot*x; - - //Kokkos::printf("set_particle: rank= %i \t dim= %i \t d=%i \t n=%i \n", myrank, dim, d, n); + + if (dim==0) + Kokkos::printf("setparticle: %i %g \n", d, x); #ifdef KOKKOS_PRINT if (n==0) @@ -498,7 +525,6 @@ class StructureFormationManager : public GravityManager { Kokkos::printf("zeldo: %i %g %g %g %g \n", dim, d+0.5,-d_z*x, -ddot*x, x); #endif }); - cfield_m = tmpcfield; } } @@ -541,7 +567,7 @@ class StructureFormationManager : public GravityManager { IpplTimings::stopTimer(linZeldo); // 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)) { @@ -549,11 +575,8 @@ class StructureFormationManager : public GravityManager { 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; - this->dumpParticles(); - */ } @@ -930,6 +953,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(); @@ -949,11 +975,13 @@ class StructureFormationManager : public GravityManager { // 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] << ","; 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(); From ca57e1bcd6f13a86175e3b08530567ee48ffae9b Mon Sep 17 00:00:00 2001 From: Andreas Adelmann Date: Wed, 17 Sep 2025 23:00:51 +0300 Subject: [PATCH 09/18] add velocities to csv file, they also match in the picture norm --- cosmology/StructureFormationManager.h | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/cosmology/StructureFormationManager.h b/cosmology/StructureFormationManager.h index c7ead80215..eafd972957 100644 --- a/cosmology/StructureFormationManager.h +++ b/cosmology/StructureFormationManager.h @@ -970,15 +970,16 @@ 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] << " \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] << ","; */ From ef51d57ee2cd6299278163193769c48e1d78e1ac Mon Sep 17 00:00:00 2001 From: Andreas Adelmann Date: Wed, 17 Sep 2025 23:13:49 +0300 Subject: [PATCH 10/18] remove duplicated code, enable integration --- cosmology/StructureFormation.cpp | 2 +- cosmology/StructureFormationManager.h | 61 +++++++++------------------ 2 files changed, 22 insertions(+), 41 deletions(-) diff --git a/cosmology/StructureFormation.cpp b/cosmology/StructureFormation.cpp index 0b03c94498..2646ff4009 100644 --- a/cosmology/StructureFormation.cpp +++ b/cosmology/StructureFormation.cpp @@ -159,7 +159,7 @@ int main(int argc, char* argv[]) { msg << "Starting iterations ... up to " << manager.getNt() << endl; - // manager.run(manager.getNt()); + manager.run(manager.getNt()); msg << "End." << endl; diff --git a/cosmology/StructureFormationManager.h b/cosmology/StructureFormationManager.h index eafd972957..88af15435a 100644 --- a/cosmology/StructureFormationManager.h +++ b/cosmology/StructureFormationManager.h @@ -160,25 +160,6 @@ class StructureFormationManager : public GravityManager { msg << "Read in particles ..." << endl; readParticlesDomain(); // defines particle positions, velocities msg << "Read particles done" << endl; - 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 reading initial conditions"; } else { msg << "Create Particles ....." << endl; ippl::ParameterList fftParams; @@ -188,28 +169,29 @@ class StructureFormationManager : public GravityManager { Pk_m.initialize(this->fcontainer_m->getMesh(), this->fcontainer_m->getFL()); msg << "FFT and Pk structures initialized" << endl;; createParticles(); - msg << "Create particles done" << endl; + msg << "Create particles done" << endl; + } - static IpplTimings::TimerRef DummySolveTimer = IpplTimings::getTimer("solveWarmup"); - IpplTimings::startTimer(DummySolveTimer); - - this->fcontainer_m->getRho() = 0.0; - - this->fsolver_m->runSolver(); + 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::stopTimer(DummySolveTimer); + + this->par2grid(); - IpplTimings::startTimer(SolveTimer); - this->fsolver_m->runSolver(); - IpplTimings::stopTimer(SolveTimer); + 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->grid2par(); + this->dump(); + msg << "Done creating initial conditions" << endl; + this->savePositions(0); } @@ -515,11 +497,10 @@ class StructureFormationManager : public GravityManager { 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); - -#ifdef KOKKOS_PRINT + 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); From 35b89f338dbc6d64ad749a9848814d2b157c3d13 Mon Sep 17 00:00:00 2001 From: Andreas Adelmann Date: Wed, 24 Sep 2025 10:10:29 +0300 Subject: [PATCH 11/18] add float for mass and greate Pk in single precision --- cosmology/GravityManager.h | 4 +- cosmology/GravityParticleContainer.hpp | 4 +- cosmology/StructureFormation.cpp | 2 - cosmology/StructureFormationManager.h | 56 +++++++++++++++++++++----- 4 files changed, 50 insertions(+), 16 deletions(-) diff --git a/cosmology/GravityManager.h b/cosmology/GravityManager.h index 939be52663..222b25ecad 100644 --- a/cosmology/GravityManager.h +++ b/cosmology/GravityManager.h @@ -305,7 +305,7 @@ class GravityManager : public ippl::PicManager 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; @@ -321,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 > 10.*Kokkos::Experimental::epsilon_v ) { + 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; diff --git a/cosmology/GravityParticleContainer.hpp b/cosmology/GravityParticleContainer.hpp index ab000c6e46..dcb3a0f249 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 2646ff4009..51c8f2e4d5 100644 --- a/cosmology/StructureFormation.cpp +++ b/cosmology/StructureFormation.cpp @@ -126,8 +126,6 @@ int main(int argc, char* argv[]) { for (unsigned d = 0; d < Dim; d++) { nr[d] = initializer::GlobalStuff::instance().ngrid; } - - msg << "Np = " << totalP << endl; // Number of time steps int nt = 0; diff --git a/cosmology/StructureFormationManager.h b/cosmology/StructureFormationManager.h index 88af15435a..d761e370c3 100644 --- a/cosmology/StructureFormationManager.h +++ b/cosmology/StructureFormationManager.h @@ -24,15 +24,26 @@ #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; @@ -151,8 +162,12 @@ 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)); @@ -164,14 +179,27 @@ class StructureFormationManager : public GravityManager { 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;; + 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); @@ -314,6 +342,7 @@ class StructureFormationManager : public GravityManager { 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); @@ -441,7 +470,6 @@ class StructureFormationManager : public GravityManager { zz = 1.0; break; } - // Compute displacement component in k-space ippl::parallel_for("force_in_k_space", ippl::getRangePolicy(cfview, ngh), @@ -508,6 +536,8 @@ class StructureFormationManager : public GravityManager { }); cfield_m = tmpcfield; } + Kokkos::fence(); + } } /** @@ -519,24 +549,32 @@ class StructureFormationManager : public GravityManager { 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 @@ -546,9 +584,9 @@ class StructureFormationManager : public GravityManager { 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)) { @@ -556,13 +594,11 @@ class StructureFormationManager : public GravityManager { 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. From 5ed6fa6194e9cf8eaf59e931c2597fc0018c8ebb Mon Sep 17 00:00:00 2001 From: Andreas Adelmann Date: Wed, 24 Sep 2025 10:17:41 +0300 Subject: [PATCH 12/18] allow type of attribute be different than field in scatter --- .../OrthogonalRecursiveBisection.hpp | 5 ++- src/Interpolation/CIC.h | 5 ++- src/Interpolation/CIC.hpp | 15 +++++-- src/Particle/ParticleAttrib.hpp | 2 +- test/particle/TestScatter.cpp | 40 ++++++++++++++----- 5 files changed, 49 insertions(+), 18 deletions(-) diff --git a/src/Decomposition/OrthogonalRecursiveBisection.hpp b/src/Decomposition/OrthogonalRecursiveBisection.hpp index c00aa696b9..069c70a317 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 5adac510d3..4ad5117ea1 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 cbe8064f13..6080128fb0 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 78df540589..a38c342108 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 d92f6a9b13..3c78c73da5 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; } From 05c6898968d37fef41639dd528b8f7988c67a8bb Mon Sep 17 00:00:00 2001 From: Andreas Adelmann Date: Thu, 6 Aug 2026 16:30:43 +0200 Subject: [PATCH 13/18] Make ParticleAttrib scatter hash dispatch CUDA-safe Split ParticleAttrib::scatter into a host-side hash/no-hash dispatcher and a templated scatterImpl so the non-hashed kernel does not instantiate hash-array indexing. Route the mapped-index selection through a small KOKKOS_INLINE_FUNCTION helper. This avoids the CUDA extended-lambda restriction where a variable is first captured inside an if constexpr context, while still keeping the hash access compiled out for the non-hashed instantiation. Also normalize the hash extent to the range-policy end type before checking the policy range, removing the signed/unsigned comparison warning. --- src/Particle/ParticleAttrib.h | 4 +++ src/Particle/ParticleAttrib.hpp | 49 +++++++++++++++++++++++++-------- 2 files changed, 42 insertions(+), 11 deletions(-) diff --git a/src/Particle/ParticleAttrib.h b/src/Particle/ParticleAttrib.h index 80b9dad20a..d2725bea2c 100644 --- a/src/Particle/ParticleAttrib.h +++ b/src/Particle/ParticleAttrib.h @@ -272,6 +272,10 @@ namespace ippl { void internalCopy(const hash_type& indices) override; private: + template + void scatterImpl(Field& f, const ParticleAttrib, Properties...>& pp, + policy_type iteration_policy, hash_type hash_array) const; + view_type dview_m{"ParticleAttrib::dview", 0}; view_type buf_m{"ParticleAttrib::buf", 0}; }; diff --git a/src/Particle/ParticleAttrib.hpp b/src/Particle/ParticleAttrib.hpp index 15cb155d9c..c400bc3e57 100644 --- a/src/Particle/ParticleAttrib.hpp +++ b/src/Particle/ParticleAttrib.hpp @@ -31,6 +31,18 @@ #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; + } + } + } // namespace detail template void ParticleAttrib::create(size_type n, bool non_destructive) { @@ -131,6 +143,27 @@ namespace ippl { template template void ParticleAttrib::scatter( + Field& f, const ParticleAttrib, Properties...>& pp, + policy_type iteration_policy, hash_type hash_array) const { + 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(); + } + + if (useHashView) { + scatterImpl(f, pp, iteration_policy, hash_array); + } else { + scatterImpl(f, pp, iteration_policy, hash_array); + } + } + + template + template + void ParticleAttrib::scatterImpl( Field& f, const ParticleAttrib, Properties...>& pp, policy_type iteration_policy, hash_type hash_array) const { constexpr unsigned Dim = Field::dim; @@ -155,20 +188,14 @@ namespace ippl { 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))) { - 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(); + auto dview = dview_m; + 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 - size_t mapped_idx = useHashView ? hash_array(idx) : idx; + const size_t mapped_idx = + detail::scatterMappedIndex(idx, hashView); vector_type l = (ppview(mapped_idx) - origin) * invdx + 0.5; Vector index = l; From ba26b1cdf0c7f66df174224d90ac6ad932e6f28c Mon Sep 17 00:00:00 2001 From: Andreas Adelmann Date: Thu, 6 Aug 2026 17:06:40 +0200 Subject: [PATCH 14/18] Make ParticleAttrib scatter dispatch CUDA-safe Move the hash-aware ParticleAttrib scatter kernel body out of the ParticleAttrib private member scope and into a namespace-scope detail helper. CUDA extended host-device lambdas cannot be enclosed by private or protected class member functions. The previous scatterImpl helper fixed the hash-array constexpr dispatch but still placed the KOKKOS_LAMBDA inside a private member, which NVCC rejects for CUDA builds. Keep ParticleAttrib::scatter as the public host-side dispatcher. It validates the optional hash view and selects either the hashed or non-hashed compile-time implementation: particleAttribScatterImpl particleAttribScatterImpl The mapped-index selection remains in a KOKKOS_INLINE_FUNCTION helper so the hash access is compiled only for the hashed instantiation, avoiding CUDA's first-capture-in-constexpr-if extended-lambda restriction. Tested locally with a SERIAL debug build: cmake --build build-fixissue-415-merge -j 8 ctest --test-dir build-fixissue-415-merge --output-on-failure -R '^(TestScatter|GatherScatterTest)$' Both GatherScatterTest and TestScatter passed. --- src/Particle/ParticleAttrib.h | 4 -- src/Particle/ParticleAttrib.hpp | 114 ++++++++++++++++---------------- 2 files changed, 57 insertions(+), 61 deletions(-) diff --git a/src/Particle/ParticleAttrib.h b/src/Particle/ParticleAttrib.h index d2725bea2c..80b9dad20a 100644 --- a/src/Particle/ParticleAttrib.h +++ b/src/Particle/ParticleAttrib.h @@ -272,10 +272,6 @@ namespace ippl { void internalCopy(const hash_type& indices) override; private: - template - void scatterImpl(Field& f, const ParticleAttrib, Properties...>& pp, - policy_type iteration_policy, hash_type hash_array) const; - view_type dview_m{"ParticleAttrib::dview", 0}; view_type buf_m{"ParticleAttrib::buf", 0}; }; diff --git a/src/Particle/ParticleAttrib.hpp b/src/Particle/ParticleAttrib.hpp index c400bc3e57..92e54daedd 100644 --- a/src/Particle/ParticleAttrib.hpp +++ b/src/Particle/ParticleAttrib.hpp @@ -42,6 +42,61 @@ namespace ippl { 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 @@ -155,67 +210,12 @@ namespace ippl { } if (useHashView) { - scatterImpl(f, pp, iteration_policy, hash_array); + detail::particleAttribScatterImpl(f, dview_m, pp, iteration_policy, hash_array); } else { - scatterImpl(f, pp, iteration_policy, hash_array); + detail::particleAttribScatterImpl(f, dview_m, pp, iteration_policy, hash_array); } } - template - template - void ParticleAttrib::scatterImpl( - 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(); - - auto dview = dview_m; - 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); - } - template template void ParticleAttrib::gather( From b28412aba727ea3ada382062c27ffad8a9121a74 Mon Sep 17 00:00:00 2001 From: Andreas Adelmann Date: Thu, 6 Aug 2026 17:16:18 +0200 Subject: [PATCH 15/18] clang formatting --- src/Particle/ParticleAttrib.h | 12 ++++++------ src/Particle/ParticleAttrib.hpp | 16 ++++++---------- 2 files changed, 12 insertions(+), 16 deletions(-) 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 92e54daedd..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" @@ -45,9 +44,8 @@ namespace ippl { template - void particleAttribScatterImpl(Field& f, const ValuesView& dview, - const PositionAttrib& pp, policy_type iteration_policy, - HashView hash_array) { + 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; @@ -78,8 +76,8 @@ namespace ippl { const size_t mapped_idx = detail::scatterMappedIndex(idx, hashView); - vector_type l = (ppview(mapped_idx) - origin) * invdx + 0.5; - Vector index = l; + vector_type l = (ppview(mapped_idx) - origin) * invdx + 0.5; + Vector index = l; Vector whi = l - index; Vector wlo = 1.0 - whi; @@ -437,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; @@ -479,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) { From 537c539edb1a47f77359b7fafa36cb461b10f0d5 Mon Sep 17 00:00:00 2001 From: Andreas Adelmann Date: Thu, 6 Aug 2026 19:32:11 +0200 Subject: [PATCH 16/18] Add mixed-value ParticleAttrib scatter regression tests Add explicit GatherScatterTest coverage for scattering a lower-precision particle attribute into a higher-precision field: ParticleAttrib -> Field This exercises the scatter API change that allows the particle attribute value type to differ from the field value type. The new coverage includes both the plain scatter path and the custom hash-array scatter path, so the host-side hash/no-hash dispatch and the CUDA-safe scatter implementation are both covered for mixed value types. Extend the test bunch with an additional float charge attribute and add a generic hash-reduction helper so expected charges can be accumulated in double while reading either float or double particle attributes. Add a scatterConservationTolerance helper that chooses the tolerance from the weaker precision of the scatter attribute and field value type. This avoids testing float -> double scatter with an unrealistically strict double tolerance while still scaling the allowed error with the conserved charge norm. Tested locally with: cmake --build build-fixissue-415-merge --target GatherScatterTest -j 8 ctest --test-dir build-fixissue-415-merge --output-on-failure -R '^GatherScatterTest$' mpirun -np 2 build-fixissue-415-merge/unit_tests/Particle/GatherScatterTest --gtest_filter='*MixedValueType*' mpirun -np 4 build-fixissue-415-merge/unit_tests/Particle/GatherScatterTest --gtest_filter='*MixedValueType*' All tests passed. --- unit_tests/Particle/GatherScatterTest.cpp | 257 ++++++++++++++++------ 1 file changed, 194 insertions(+), 63 deletions(-) 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; From f2d409ad996085e2ffcaee6aa01cc1cd7ac5357a Mon Sep 17 00:00:00 2001 From: Andreas Adelmann Date: Thu, 6 Aug 2026 19:54:42 +0200 Subject: [PATCH 17/18] clang formating --- demos/cosmology/GravityLoadBalancer.hpp | 31 +- demos/cosmology/GravityManager.h | 27 +- demos/cosmology/GravityParticleContainer.hpp | 2 +- demos/cosmology/StructureFormation.cpp | 67 +- demos/cosmology/StructureFormationManager.h | 938 ++++++------ .../cosmology/mc-4-Initializer/Cosmology.cpp | 1297 +++++++++-------- demos/cosmology/mc-4-Initializer/Cosmology.h | 119 +- demos/cosmology/mc-4-Initializer/DataBase.cpp | 203 +-- demos/cosmology/mc-4-Initializer/DataBase.h | 105 +- .../mc-4-Initializer/InputParser.cpp | 105 +- .../cosmology/mc-4-Initializer/InputParser.h | 36 +- demos/cosmology/mc-4-Initializer/Output.cpp | 617 ++++---- demos/cosmology/mc-4-Initializer/Output.h | 83 +- .../cosmology/mc-4-Initializer/TypesAndDefs.h | 28 +- .../OrthogonalRecursiveBisection.hpp | 5 +- src/Interpolation/CIC.h | 2 +- src/Interpolation/CIC.hpp | 18 +- 17 files changed, 1849 insertions(+), 1834 deletions(-) diff --git a/demos/cosmology/GravityLoadBalancer.hpp b/demos/cosmology/GravityLoadBalancer.hpp index 362d7128c1..fa51740753 100644 --- a/demos/cosmology/GravityLoadBalancer.hpp +++ b/demos/cosmology/GravityLoadBalancer.hpp @@ -150,22 +150,21 @@ class LoadBalancer { 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; + 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: diff --git a/demos/cosmology/GravityManager.h b/demos/cosmology/GravityManager.h index 20dadf12f8..3b9aa2defa 100644 --- a/demos/cosmology/GravityManager.h +++ b/demos/cosmology/GravityManager.h @@ -12,7 +12,6 @@ #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; @@ -44,13 +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 + * @param par_ Inputfile parser */ GravityManager(size_type totalP_, int nt_, Vector_t& nr_, double lbt_, std::string& solver_, std::string& stepMethod_, initializer::InputParser par_) : ippl::PicManager, FieldContainer, LoadBalancer>() - , parser_m(par_) + , parser_m(par_) , totalP_m(totalP_) , nt_m(nt_) , nr_m(nr_) @@ -74,7 +73,6 @@ class GravityManager : public ippl::PicManager */ initializer::InputParser parser_m; - /** * @brief Get the total number of particles. * @@ -221,9 +219,9 @@ class GravityManager : public ippl::PicManager */ void InitialiseTime() { Inform mes("Inititalise: "); - 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 + 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); @@ -253,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; } /** @@ -275,8 +273,8 @@ class GravityManager : public ippl::PicManager this->dt_m = this->Dloga / this->Hubble_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; + 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 @@ -322,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 > 10.*Kokkos::Experimental::epsilon_v ) { + 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; @@ -362,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] - float O_m; ///< Matter density parameter. [1] - float 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] diff --git a/demos/cosmology/GravityParticleContainer.hpp b/demos/cosmology/GravityParticleContainer.hpp index c7b738b223..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. diff --git a/demos/cosmology/StructureFormation.cpp b/demos/cosmology/StructureFormation.cpp index 51c8f2e4d5..387dc59b41 100644 --- a/demos/cosmology/StructureFormation.cpp +++ b/demos/cosmology/StructureFormation.cpp @@ -5,11 +5,11 @@ // 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. The old mc-4 initializer based in Zaria's -// old and new initializer is added. +// old and new initializer is added. // // Usage: // srun ./StructureFormation -// +// // --overallocate --info 10 // inputfile = describes the simulation // tfFn = transfer function @@ -24,9 +24,11 @@ // 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 input.par tf.dat out.data datadir 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 +// more stuff and todos at: +// https://docs.google.com/document/d/1PhQ2Wo6QhdeS7U3lN1q2TEaJm2XjD4eXdjGCMKfPac4/edit?usp=sharing /* Example of an input file @@ -56,7 +58,7 @@ TFFlag=2 PrintFormat=0 // new -Omega_L=0.7 +Omega_L=0.7 // Read from file (0) or create (1) ReadInParticles=1 @@ -64,8 +66,6 @@ ReadInParticles=1 */ - - constexpr unsigned Dim = 3; using T = double; @@ -85,7 +85,6 @@ using T = double; #include "Manager/PicManager.h" #include "StructureFormationManager.h" - #include "mc-4-Initializer/DataBase.h" #include "mc-4-Initializer/InputParser.h" @@ -103,38 +102,37 @@ int main(int argc, char* argv[]) { static IpplTimings::TimerRef mainTimer = IpplTimings::getTimer("total"); IpplTimings::startTimer(mainTimer); - std::string indatName = argv[1]; - std::string tfName = argv[2]; - std::string outBase = argv[3]; + 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 + + 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] = initializer::GlobalStuff::instance().ngrid; + nr[d] = initializer::GlobalStuff::instance().ngrid; } - + // Number of time steps int nt = 0; - par.getByName("nt", nt); + par.getByName("nt", nt); + + int readInParticles = 0; + par.getByName("ReadInParticles", readInParticles); + bool readICs = (readInParticles == 0); - 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 @@ -147,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, par, tfName, readICs); + StructureFormationManager manager(totalP, nt, nr, lbt, solver, step_method, par, + tfName, readICs); // set initial conditions folder manager.setIC(ic_folder); diff --git a/demos/cosmology/StructureFormationManager.h b/demos/cosmology/StructureFormationManager.h index d761e370c3..98ae5226f3 100644 --- a/demos/cosmology/StructureFormationManager.h +++ b/demos/cosmology/StructureFormationManager.h @@ -19,35 +19,33 @@ #include "Random/InverseTransformSampling.h" #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" +#include "mc-4-Initializer/DataBase.h" +#include "mc-4-Initializer/InputParser.h" -//#define DOUBLE_INIT_COND +// #define DOUBLE_INIT_COND -//#define KOKKOS_PRINT // Kokkos::printf of interesting quantities. Does not work multirank +// #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::Field, Dim, Mesh_t, Mesh_t::DefaultCentering> + field_type; typedef ippl::FFT CFFT_type; typedef Field, Dim> CField_t; -typedef Field RField_t; +typedef Field RField_t; #else -typedef ippl::Field, Dim, Mesh_t, Mesh_t::DefaultCentering> field_type; +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; + int kx, ky, kz; + double re, im; }; /** @@ -65,34 +63,32 @@ struct HermitianPkg { */ 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; - /// all for the initializer - std::unique_ptr Cfft_m; - CField_t cfield_m; - RField_t Pk_m; - bool readICs_m; + initializer::CosmoClass cosmo_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()); - } + 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. @@ -119,19 +115,19 @@ class StructureFormationManager : public GravityManager { 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 - 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; - + 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 = box_size*1000.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()); @@ -162,444 +158,457 @@ class StructureFormationManager : public GravityManager { this->solver_m, &this->fcontainer_m->getRho(), &this->fcontainer_m->getF(), &this->fcontainer_m->getPhi())); - msg << "Set FieldSolver done ..." << endl; - + msg << "Set FieldSolver done ..." << endl; + this->fsolver_m->initSolver(); - msg << "Init FieldSolver done ..." << endl; - + msg << "Init FieldSolver done ..." << endl; + this->setLoadBalancer(std::make_shared( this->lbt_m, this->fcontainer_m, this->pcontainer_m, this->fsolver_m)); - 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); + 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); } + 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; - 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; - }); + 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(); + + /** + @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); } - 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; + + 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("feval: %i %g %g %g %g %g %g \n", dim, Grad, Green, tmp.imag(), tmp.real(), res.real(), res.imag()); + Kokkos::printf("rho: %g %g \n", val_re, val_im); #endif - }); - - // Inverse FFT to real space - Cfft_m->transform(ippl::BACKWARD, cfield_m); + // 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; + } - 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); + // 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("rhoinvfft: %i %g %g \n", dim, tmp.real()/scale, tmp.imag()/scale); + 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(); + 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 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()); + 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; + 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; + 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(); + } } - 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 - */ + void createParticles() { + Inform m2a("createParticles ", INFORM_ALL_NODES); + Inform msg("createParticles "); - 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); + 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; } - */ - pc->update(); - m2a << "local number of galaxies after initializer " << pc->getLocalNum() << endl; -} - /** * @brief Read particle data from a file. */ @@ -969,10 +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"; - + } else + file << "# rx \t ry \t rz \n"; + std::shared_ptr pc = this->pcontainer_m; auto Rview = this->pcontainer_m->R.getView(); @@ -987,7 +995,7 @@ 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) { @@ -995,11 +1003,11 @@ class StructureFormationManager : public GravityManager { file << R_host(i)[d] << " \t"; for (unsigned int d = 0; d < Dim; ++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(); diff --git a/demos/cosmology/mc-4-Initializer/Cosmology.cpp b/demos/cosmology/mc-4-Initializer/Cosmology.cpp index e417721072..66d31c9041 100644 --- a/demos/cosmology/mc-4-Initializer/Cosmology.cpp +++ b/demos/cosmology/mc-4-Initializer/Cosmology.cpp @@ -6,38 +6,41 @@ GetParameters(DataBase) which gets cosmology parameters from the - DataBase (defined in DataBase.h), and stores them - with this class as well. + 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 + 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 + 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 + 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 "Cosmology.h" + +#include +#include #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" @@ -47,470 +50,494 @@ 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(); + 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]; - return; -} + /* + + Brrrrrr adafixme: need to comment this out because of abort -CosmoClass::~CosmoClass(){ - if (TFFlag == 0) { - free(table_kk); - free(table_tf); + 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); } - free(neut_p); - free(neut_c); -} #define pi M_PI -/* Transfer functions */ + /* 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::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::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::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::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::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::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)); -} + real CosmoClass::TransferFunction(real k) { + return ((this->*tan_f)(k)); + } -/* Sigma(R) routines */ + /* 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::sigma2(real k) { + real prim_ps, w_f, t_f, 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); -} + 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 */ + /* 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) -//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); - 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; + } - 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; } - return; -} + real CosmoClass::FD(real vv) { + return vv * vv / (exp(vv / neut_vv0) + 1.0); + } -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); + } -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); + 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); - return 4*M_PI*pow(m/2/M_PI/k/T, 1.5)*vv*vv*exp(-m*vv*vv/2/k/T); -} + 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; - 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; -} + return; + } + /* Growth factor for flat wCDM cosmologies: */ -/* 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::GrowthFactor(real z, real* gf, real* g_dot) { + real x1, x2, dplus, ddot; + // const float redshift=200.0; + real* ystart; - 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; -} + 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; -#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= 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 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= 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 (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 @@ -518,224 +545,228 @@ static real maxarg1,maxarg2, minarg1, minarg2; #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*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*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) + /* 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 -*/ + 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 JMIN){ - if (fabs(sol-old) < EPS*fabs(old) || - (sol==0.0 && old==0.0)) return(sol); - } - old = sol; - } - printf("integrate: no convergence!"); - abort(); -} + 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) + /* 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 -*/ + 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; + 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]; + 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); -} + 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; -} + /* 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; + 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; + } } - 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; -} + } + 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 } diff --git a/demos/cosmology/mc-4-Initializer/Cosmology.h b/demos/cosmology/mc-4-Initializer/Cosmology.h index 20ff349976..c161759a67 100644 --- a/demos/cosmology/mc-4-Initializer/Cosmology.h +++ b/demos/cosmology/mc-4-Initializer/Cosmology.h @@ -2,11 +2,11 @@ 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. + 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 @@ -20,64 +20,63 @@ #ifdef USENAMESPACE namespace initializer { #endif -class GlobalStuff; + 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; + 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(); - 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 [])); + 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; - //neutrino stuff - void genFDDist(); - real FD(real vv); - real Maxwell(real vv); -}; + 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 } diff --git a/demos/cosmology/mc-4-Initializer/DataBase.cpp b/demos/cosmology/mc-4-Initializer/DataBase.cpp index 12ad12503b..006ab56ead 100644 --- a/demos/cosmology/mc-4-Initializer/DataBase.cpp +++ b/demos/cosmology/mc-4-Initializer/DataBase.cpp @@ -2,22 +2,25 @@ Initializer: DataBase.cpp - Has only + 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 + 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" +#include +#include + +#include "TypesAndDefs.h" + initializer::GlobalStuff& initializer::GlobalStuff::instance() { static initializer::GlobalStuff instance; return instance; @@ -29,101 +32,99 @@ initializer::GlobalStuff& initializer::GlobalStuff::instance() { 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; -} + 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 } diff --git a/demos/cosmology/mc-4-Initializer/DataBase.h b/demos/cosmology/mc-4-Initializer/DataBase.h index 5cc588c936..6309d1e64d 100644 --- a/demos/cosmology/mc-4-Initializer/DataBase.h +++ b/demos/cosmology/mc-4-Initializer/DataBase.h @@ -1,12 +1,12 @@ /* 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. - + + 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 */ @@ -14,58 +14,61 @@ #ifndef DataBase_Header_Included #define DataBase_Header_Included +#include + #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 */ -}; + // 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 } diff --git a/demos/cosmology/mc-4-Initializer/InputParser.cpp b/demos/cosmology/mc-4-Initializer/InputParser.cpp index 681c80351b..e727eb72a9 100644 --- a/demos/cosmology/mc-4-Initializer/InputParser.cpp +++ b/demos/cosmology/mc-4-Initializer/InputParser.cpp @@ -1,81 +1,72 @@ #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; + InputParser::InputParser(std::string filename) { + filename_m = filename; + parseFile(); } -} - -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; + 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 = ""; + 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)); + 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); +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; } diff --git a/demos/cosmology/mc-4-Initializer/InputParser.h b/demos/cosmology/mc-4-Initializer/InputParser.h index f88e805d58..4d674a4cbd 100644 --- a/demos/cosmology/mc-4-Initializer/InputParser.h +++ b/demos/cosmology/mc-4-Initializer/InputParser.h @@ -1,11 +1,11 @@ #ifndef INPUT_PARSER_H #define INPUT_PARSER_H -#include -#include +#include #include +#include #include -#include +#include #include "TypesAndDefs.h" @@ -13,27 +13,23 @@ namespace initializer { #endif -class InputParser { + class InputParser { + public: + InputParser(std::string filename); + bool getByName(std::string name, int& param); + bool getByName(std::string name, real& param); -public: + private: + std::map RealDict_m; + std::map IntDict_m; + std::string filename_m; - 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; } - -}; + 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 index f5868d107b..e82a8c01a3 100644 --- a/demos/cosmology/mc-4-Initializer/Output.cpp +++ b/demos/cosmology/mc-4-Initializer/Output.cpp @@ -6,335 +6,328 @@ * */ -#include -#include +#include "Output.h" + #include +#include #include +#include #include -#include +#include #include + #include "TypesAndDefs.h" -#include "Output.h" -#include "InputParser.h" +#include "InputParser.h" #ifdef USENAMESPACE namespace initializer { #endif -using namespace std; + 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; -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 { + using namespace std; -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); + 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 */ - }; -}; + 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 index 5abd66b6c0..7d8a05b6cb 100644 --- a/demos/cosmology/mc-4-Initializer/TypesAndDefs.h +++ b/demos/cosmology/mc-4-Initializer/TypesAndDefs.h @@ -1,8 +1,8 @@ -/* Definitions of real and integer variables: - if the code should be in double precision, here's - the only place where change should be made. +/* 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 */ @@ -15,7 +15,7 @@ namespace initializer { #endif #ifdef DOUBLE_REAL -typedef double real; + typedef double real; #define MY_MPI_REAL MPI_DOUBLE #else typedef float real; @@ -23,19 +23,21 @@ typedef float real; #endif #ifdef LONG_INTEGER -typedef long integer; + typedef long integer; #define MY_MPI_INTEGER MPI_DOUBLE #else - // typedef int integer; - // #define MY_MPI_INTEGER MPI_INT +// typedef int integer; +// #define MY_MPI_INTEGER MPI_INT #endif -typedef struct{ - double re; - double im; -} my_fftw_complex; + typedef struct { + double re; + double im; + } my_fftw_complex; -inline int MOD(int x, int y) { return (x - y*(int)(x/y));} + inline int MOD(int x, int y) { + return (x - y * (int)(x / y)); + } #ifdef USENAMESPACE } diff --git a/src/Decomposition/OrthogonalRecursiveBisection.hpp b/src/Decomposition/OrthogonalRecursiveBisection.hpp index 3d67887444..2e2ad60689 100644 --- a/src/Decomposition/OrthogonalRecursiveBisection.hpp +++ b/src/Decomposition/OrthogonalRecursiveBisection.hpp @@ -1,6 +1,7 @@ -#include "Utility/IpplTimings.h" #include +#include "Utility/IpplTimings.h" + namespace ippl { template @@ -293,7 +294,7 @@ 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, 1); }); diff --git a/src/Interpolation/CIC.h b/src/Interpolation/CIC.h index 4ad5117ea1..570a170651 100644 --- a/src/Interpolation/CIC.h +++ b/src/Interpolation/CIC.h @@ -68,7 +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 + $ @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 diff --git a/src/Interpolation/CIC.hpp b/src/Interpolation/CIC.hpp index 288165afd2..568915c271 100644 --- a/src/Interpolation/CIC.hpp +++ b/src/Interpolation/CIC.hpp @@ -1,6 +1,5 @@ #include - namespace ippl { namespace detail { template @@ -33,19 +32,20 @@ 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, Val val) { - // The number of indices is equal to the view rank + // 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); + /** 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, v_out), ...); From 8e7f2a94f0c75c823c3cf42510830f51d5932bad Mon Sep 17 00:00:00 2001 From: Andreas Adelmann Date: Thu, 6 Aug 2026 20:45:01 +0200 Subject: [PATCH 18/18] add test with float/double scatter --- test/particle/TestScatter.cpp | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/test/particle/TestScatter.cpp b/test/particle/TestScatter.cpp index 49d4606c56..8ca2403a13 100644 --- a/test/particle/TestScatter.cpp +++ b/test/particle/TestScatter.cpp @@ -7,14 +7,14 @@ struct Bunch : public ippl::ParticleBase { Bunch(PLayout& playout) : ippl::ParticleBase(playout) { this->addAttribute(Q1); - this->addAttribute(Q2); + this->addAttribute(Q2); } ~Bunch() {} typedef ippl::ParticleAttrib charge_container_typeF; charge_container_typeF Q1; - + typedef ippl::ParticleAttrib charge_container_typeD; charge_container_typeD Q2; }; @@ -103,12 +103,13 @@ int main(int argc, char* argv[]) { 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; + 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.Q2 = 1.0; bunch.update(); @@ -122,7 +123,8 @@ int main(int argc, char* argv[]) { 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; + 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; }