Fixissue/#415

demos/cosmology/CMakeLists.txt

@@ -5,6 +5,13 @@
 # -----------------------------------------------------------------------------
 message(STATUS "Configuring demos/cosmology/")
 
-add_executable(StructureFormation StructureFormation.cpp)
+add_executable(StructureFormation
+    StructureFormation.cpp
+    mc-4-Initializer/Cosmology.cpp
+    mc-4-Initializer/DataBase.cpp
+    mc-4-Initializer/InputParser.cpp
+)
+
+target_compile_definitions(StructureFormation PRIVATE USENAMESPACE)
 
 target_link_libraries(StructureFormation PRIVATE IPPL::ippl)

demos/cosmology/GravityLoadBalancer.hpp

@@ -149,26 +149,22 @@ class LoadBalancer {
         if (ippl::Comm->size() < 2) {
             return false;
         }
-        if (std::strcmp(TestName, "UniformPlasmaTest") == 0) {
-            return (nstep % loadbalancefreq_m == 0);
-        } else {
-            int local = 0;
-            std::vector<int> res(ippl::Comm->size());
-            double equalPart = static_cast<double>(totalP) / ippl::Comm->size();
-            double dev = std::abs(static_cast<double>(pc_m->getLocalNum()) - equalPart) / totalP;
-            if (dev > loadbalancethreshold_m) {
-                local = 1;
-            }
-            MPI_Allgather(&local, 1, MPI_INT, res.data(), 1, MPI_INT,
-                          ippl::Comm->getCommunicator());
 
-            for (unsigned int i = 0; i < res.size(); i++) {
-                if (res[i] == 1) {
-                    return true;
-                }
+        int local = 0;
+        std::vector<int> res(ippl::Comm->size());
+        double equalPart = static_cast<double>(totalP) / ippl::Comm->size();
+        double dev       = std::abs(static_cast<double>(pc_m->getLocalNum()) - equalPart) / totalP;
+        if (dev > loadbalancethreshold_m) {
+            local = 1;
+        }
+        MPI_Allgather(&local, 1, MPI_INT, res.data(), 1, MPI_INT, ippl::Comm->getCommunicator());
+
+        for (unsigned int i = 0; i < res.size(); i++) {
+            if (res[i] == 1) {
+                return true;
             }
-            return false;
         }
+        return false;
     }
 
 private:
@@ -182,4 +178,4 @@ class LoadBalancer {
     ORB<T, Dim> orb;                                  ///< ORB for domain partitioning
 };
 
-#endif  // IPPL_LOAD_BALANCER_H
+#endif  // IPPL_LOAD_BALANCER_H

demos/cosmology/GravityManager.h

@@ -12,6 +12,7 @@
 #include "Random/InverseTransformSampling.h"
 #include "Random/NormalDistribution.h"
 #include "Random/Randn.h"
+#include "mc-4-Initializer/InputParser.h"
 
 using view_type = typename ippl::detail::ViewType<ippl::Vector<double, Dim>, 1>::view_type;
 
@@ -42,11 +43,13 @@ class GravityManager : public ippl::PicManager<T, Dim, ParticleContainer<T, Dim>
      * @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<int, Dim>& nr_, double lbt_,
-                   std::string& solver_, std::string& stepMethod_)
+                   std::string& solver_, std::string& stepMethod_, initializer::InputParser par_)
         : ippl::PicManager<T, Dim, ParticleContainer<T, Dim>, FieldContainer<T, Dim>,
                            LoadBalancer<T, Dim>>()
+        , parser_m(par_)
         , totalP_m(totalP_)
         , nt_m(nt_)
         , nr_m(nr_)
@@ -65,6 +68,11 @@ class GravityManager : public ippl::PicManager<T, Dim, ParticleContainer<T, Dim>
      */
     std::string folder;
 
+    /**
+     * @brief Access to the input file with constants and simulation parameters.
+     */
+    initializer::InputParser parser_m;
+
     /**
      * @brief Get the total number of particles.
      *
@@ -211,8 +219,10 @@ class GravityManager : public ippl::PicManager<T, Dim, ParticleContainer<T, Dim>
      */
     void InitialiseTime() {
         Inform mes("Inititalise: ");
-        this->O_m      = 0.3;
-        this->O_L      = 0.7;
+        parser_m.getByName("Omega_m", this->O_m);
+        parser_m.getByName("Omega_L", this->O_L);
+        // this->O_m      = 0.3;       // \todo need to from input file
+        // this->O_L      = 0.7;       // \todo need to from input file
         this->t_L      = 2 / (3 * this->Hubble0 * sqrt(this->O_L));
         this->a_m      = 1 / (1 + this->z_m);
         this->Dloga    = 1. / (this->nt_m) * log((1 + this->z_m) / (1 + this->z_f));
@@ -241,8 +251,8 @@ class GravityManager : public ippl::PicManager<T, Dim, ParticleContainer<T, Dim>
      * @brief Pre-step method called before each simulation step.
      */
     void pre_step() override {
-        Inform mes("Pre-step");
-        mes << "Done" << endl;
+        //        Inform mes("Pre-step");
+        //  mes << "Done" << endl;
     }
 
     /**
@@ -262,9 +272,9 @@ class GravityManager : public ippl::PicManager<T, Dim, ParticleContainer<T, Dim>
         // dynamic time step
         this->dt_m = this->Dloga / this->Hubble_m;
 
-        mes << "Finished time step: " << this->it_m << endl;
-        mes << " time: " << this->time_m << ", timestep: " << this->dt_m << ", z: " << this->z_m
-            << ", a: " << this->a_m << endl;
+        mes << "Step: " << this->it_m;
+        mes << " comological time: " << this->time_m << ", dt: " << this->dt_m
+            << ", z: " << this->z_m << ", a: " << this->a_m << endl;
     }
 
     // Grid to particle and particle to grid methods
@@ -291,10 +301,10 @@ class GravityManager : public ippl::PicManager<T, Dim, ParticleContainer<T, Dim>
      */
     void scatterCIC() {
         Inform mes("scatter ");
-        mes << "starting ..." << endl;
+
         this->fcontainer_m->getRho() = 0.0;
 
-        ippl::ParticleAttrib<double>* m          = &this->pcontainer_m->m;
+        ippl::ParticleAttrib<float>* m           = &this->pcontainer_m->m;
         typename Base::particle_position_type* R = &this->pcontainer_m->R;
         Field_t<Dim>* rho                        = &this->fcontainer_m->getRho();
         Vector_t<double, Dim> rmin               = rmin_m;
@@ -303,15 +313,15 @@ class GravityManager : public ippl::PicManager<T, Dim, ParticleContainer<T, Dim>
 
         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();
 
         ippl::Comm->reduce(localParticles, TotalParticles, 1, std::plus<size_type>());
 
         if (ippl::Comm->rank() == 0) {
-            if (TotalParticles != totalP_m || relError > 1e-10) {
+            if (TotalParticles != totalP_m
+                || relError > 10. * Kokkos::Experimental::epsilon_v<float>) {
                 mes << "Time step: " << it_m << endl;
                 mes << "Total particles in the sim. " << totalP_m << " "
                     << "after update: " << TotalParticles << endl;
@@ -351,8 +361,8 @@ class GravityManager : public ippl::PicManager<T, Dim, ParticleContainer<T, Dim>
     double Hubble0;                  ///< Hubble constant today (73.8 km/sec/Mpc).
     double G;                        ///< Gravitational constant. [kpc^3/(Msun s^2)]
     double rho_crit0;                ///< Critical density today. [Msun/kpc^3]
-    double O_m;                      ///< Matter density parameter. [1]
-    double O_L;                      ///< Dark energy density parameter. [1]
+    float O_m;                       ///< Matter density parameter. [1]
+    float O_L;                       ///< Dark energy density parameter. [1]
     double t_L;                      ///< Characteristic time scale. [s]
     double z_m;                      ///< Initial redshift. [1]
     double z_f;                      ///< Final redshift.   [1]
@@ -368,4 +378,4 @@ class GravityManager : public ippl::PicManager<T, Dim, ParticleContainer<T, Dim>
     std::array<bool, Dim> decomp_m;  ///< Decomposition flags for each dimension. [bool]
     double rhoNorm_m;                ///< Normalized density. [1]
 };
-#endif
+#endif

demos/cosmology/GravityParticleContainer.hpp

@@ -21,7 +21,7 @@ class ParticleContainer : public ippl::ParticleBase<
     /**
      * @brief Mass of the particle.
      */
-    ippl::ParticleAttrib<double> m;
+    ippl::ParticleAttrib<float> m;
 
     /**
      * @brief Velocity of the particle.
@@ -91,4 +91,4 @@ class ParticleContainer : public ippl::ParticleBase<
     void setBCAllPeriodic() { this->setParticleBC(ippl::BC::PERIODIC); }
 };
 
-#endif  // IPPL_PARTICLE_CONTAINER_H
+#endif  // IPPL_PARTICLE_CONTAINER_H

demos/cosmology/StructureFormation.cpp

@@ -1,34 +1,73 @@
-// StructureFormation Test
-//
+// StructureFormation.cpp
 // The StructureFormation simulation is designed to model the formation and evolution of cosmic
 // structures, such as galaxies and clusters of galaxies, in the universe. It uses initial
 // conditions, including particle positions and velocities, to simulate the gravitational
 // interactions and dynamics of a large number of particles over time. The goal is to understand how
 // initial density fluctuations grow and evolve under the influence of gravity, leading to the
-// large-scale structure observed in the universe today.
+// large-scale structure observed in the universe today. The old mc-4 initializer based in Zaria's
+// old and new initializer is added.
 //
 //   Usage:
 //     srun ./StructureFormation
-//                  <path> <nx> [<ny>...] <Np> <Nt> <stype>
-//                  <lbthres> --overallocate <ovfactor> --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.
+//                  <inputfile> <tfFn> <outFn> <outDir> <fsType>
+//                  <lbthres> <integrator> --overallocate <ovfactor> --info 10
+//     inputfile = describes the simulation
+//     tfFn      = transfer function
+//     outFn     = base file name for snapshots
+//     outDir    = base directory for output data
+//     fsType    = Field solver type (FFT and CG supported)
+//     lbthres   = Load balancing threshold i.e., lbthres*100 is the maximum load imbalance
+//                 percentage which can be tolerated and beyond which
+//                 particle load balancing occurs. A value of 0.01 is good for many typical
+//                 simulations.
+//     integrator= LeapFrog
+//     ovfactor  = Over-allocation factor for the buffers used in the communication. Typical
+//                 values are 1.0, 2.0. Value 1.0 means no over-allocation.
 //     Example:
-//     srun ./StructureFormation data/lsf_32/ 32 32 32 32768 10 FFT 1.0 LeapFrog --overallocate 1.0
+//     srun StructureFormation input.par tf.dat out.data datadir FFT 1.0 LeapFrog --overallocate 1.0
 //     --info 5
+//
+//     more stuff and todos at:
+//     https://docs.google.com/document/d/1PhQ2Wo6QhdeS7U3lN1q2TEaJm2XjD4eXdjGCMKfPac4/edit?usp=sharing
+
+/*  Example of an input file
+
+np=128
+nt=0
+box_size=64.0   // In Mpc/h
+seed=9854373
+z_in=50.0
+z_fi=0.0
+
+// Cosmology:
+hubble=0.7
+Omega_m=0.25   // Total matter content; 1-this will be DE
+Omega_nu=0.0
+Omega_bar=0.04
+Omega_r=0.0
+
+Sigma_8=0.8
+n_s=0.9
+w_de=-1.0
+N_nu=3
+nu_pairs=4
+f_NL=0.0
+TFFlag=2
+
+// Code stuff:
+PrintFormat=0
+
+// new
+Omega_L=0.7
+
+// Read from file (0) or create (1)
+ReadInParticles=1
+
+
+*/
 
 constexpr unsigned Dim = 3;
 using T                = double;
-const char* TestName   = "StructureFormation";
 
 #include "Ippl.h"
 
@@ -42,12 +81,12 @@ const char* TestName   = "StructureFormation";
 #include <string>
 #include <vector>
 
-#include "Manager/datatypes.h"
-
 #include "Utility/IpplTimings.h"
 
 #include "Manager/PicManager.h"
 #include "StructureFormationManager.h"
+#include "mc-4-Initializer/DataBase.h"
+#include "mc-4-Initializer/InputParser.h"
 
 using size_type = ippl::detail::size_type;
 
@@ -57,25 +96,42 @@ using Vector_t = ippl::Vector<T, Dim>;
 int main(int argc, char* argv[]) {
     ippl::initialize(argc, argv);
     {
-        Inform msg(TestName);
-        Inform msg2all(TestName, INFORM_ALL_NODES);
+        Inform msg("StructureFormation ");
+        Inform msg2all("StructureFormation ", INFORM_ALL_NODES);
 
         static IpplTimings::TimerRef mainTimer = IpplTimings::getTimer("total");
         IpplTimings::startTimer(mainTimer);
 
-        // Read input parameters, assign them to the corresponding members of manager
-        int arg = 1;
-        // initial conditions folder
-        std::string ic_folder = argv[arg++];
+        std::string indatName = argv[1];
+        std::string tfName    = argv[2];
+        std::string outBase   = argv[3];
+        std::string ic_folder = argv[4];
+
+        int arg = 5;
+        initializer::InputParser par(indatName);
+        initializer::GlobalStuff::instance().GetParameters(par);
+
+        size_t totalP;
+        if ((totalP = static_cast<std::size_t>(initializer::GlobalStuff::instance().ngrid)
+                      * static_cast<std::size_t>(initializer::GlobalStuff::instance().ngrid)
+                      * static_cast<std::size_t>(initializer::GlobalStuff::instance().ngrid))
+            > std::numeric_limits<std::size_t>::max()) {
+            throw std::overflow_error("Index space exceeds size_t capacity");
+        }
+
         // Number of gridpoints in each dimension
         Vector_t<int, Dim> nr;
         for (unsigned d = 0; d < Dim; d++) {
-            nr[d] = std::atoi(argv[arg++]);
+            nr[d] = initializer::GlobalStuff::instance().ngrid;
         }
-        // Total number of particles
-        size_type totalP = std::atoll(argv[arg++]);
+
         // Number of time steps
-        int nt = std::atoi(argv[arg++]);
+        int nt = 0;
+        par.getByName("nt", nt);
+
+        int readInParticles = 0;
+        par.getByName("ReadInParticles", readInParticles);
+        bool readICs = (readInParticles == 0);
 
         // Solver method
         std::string solver = argv[arg++];
@@ -89,15 +145,16 @@ int main(int argc, char* argv[]) {
         std::string step_method = argv[arg++];
 
         // Create an instance of a manager for the considered application
-        StructureFormationManager<T, Dim> manager(totalP, nt, nr, lbt, solver, step_method);
+        StructureFormationManager<T, Dim> manager(totalP, nt, nr, lbt, solver, step_method, par,
+                                                  tfName, readICs);
 
         // set initial conditions folder
         manager.setIC(ic_folder);
 
         // 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());