Fixissue/#415

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})
 
 

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<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;
-        }
+	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;
     }
 
 private:
@@ -180,4 +177,4 @@ class LoadBalancer {
     ORB<T, Dim> orb;                                  ///< ORB for domain partitioning
 };
 
-#endif  // IPPL_LOAD_BALANCER_H
+#endif  // IPPL_LOAD_BALANCER_H

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<ippl::Vector<double, Dim>, 1>::view_type;
 
 /**
@@ -41,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_)
@@ -64,6 +68,12 @@ 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.
      *
@@ -210,8 +220,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));
@@ -240,8 +252,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;
     }
 
     /**
@@ -261,8 +273,8 @@ 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
+        mes << "Step: " << this->it_m;
+        mes << " comological time: " << this->time_m << ", dt: " << this->dt_m << ", z: " << this->z_m
             << ", a: " << this->a_m << endl;
     }
 
@@ -290,10 +302,10 @@ class GravityManager : public ippl::PicManager<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;
@@ -302,15 +314,14 @@ 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;
@@ -350,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]
@@ -367,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

cosmology/GravityParticleContainer.hpp

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