Simplified the structure of main_program.py by removing redundant steps.

Also updated:
 - model files
 - linearized evolution re-initialization
 - README

README

@@ -9,7 +9,12 @@ Please cite arXiv:
 if you use this code in your research.
 See also http://www.physics.utu.fi/theory/particlecosmology/pycool/
 
-Please submit any errors at https://github.com/jtksai/PyCOOL/issues
+I have done my best to add comments and info texts into different
+functions in the Python and CUDA codes in order to make the program
+more understandable. If you still have questions please email me or
+post them at https://github.com/jtksai/PyCOOL/ .
+
+Please submit any errors at https://github.com/jtksai/PyCOOL/issues .
 
 ------------------------------------------------------
 
@@ -146,25 +151,26 @@ also in other operating systems.
 
 Typing 'python main_program.py' runs the code for a specific model that is
 imported from the models directory. The current program prints to the screen
-information of current step number, scale factor, canonical momentum p and
-the accumulated numerical error.
+information of current step number, scale factor, canonical momentum p,
+physical time and numerical error.
+
+The initial field perturbations are created with the method that was
+also used in Defrost.
 
-The main program is divided into four parts. Thse are:
+The main program is divided into
  - Homogeneous system solver
- - Homogeneous + linearized perturbation sovler
- - Non-linear solver
- - Non-linear solver for multiple simulations e.g. for non-Gaussianity simulations
-Which of these to use have to be specified in the imported model object.
+ - Non-linear solver that includes a linearized perturbation solver
 
 The program creates a new folder by time stamp in the data folder.
-For large lattices and frequent data saving generated data can be
-several hundred gigabytes. The folder also includes a simple info file
-that tells basic information about the model used to generate the data.
-These include the initial values and the potential functions.
+For large lattices and frequent data saving the generated data can be
+several hundred gigabytes if the fields are also stores. The folder also
+includes a simple info file that tells basic information about the model
+used to generate the data. These include the initial values and
+the potential functions.
 
-After the simulation is done PyCOOL can calculate spectrums used in
+During the simulation PyCOOL can also calculate various spectrums used in
 LatticeEasy and Defrost that are written into the generated SILO files.
-In VisIt these are available under Add -> Curve.
+In VisIt these are available under Add -> Curve. See section 5 for details.
 
 ------------------------------------------------------
 
@@ -200,26 +206,28 @@ Possible variables to change/consider include:
 5 Output and Post-processing functions
 
 PyCOOL writes variety of variables into the Silo files during run time. Which of these and how often
-write are determined in the scalar field model object.
+to write are determined in the scalar field model object. Note that the different spectra are
+calculated on the CPU and therefore the data must be transferred from the device to the host memory.
+This can significantly hinder the performance of the code.
 The variables include:
  - scale factor a
  - Hubble parameter H
  - Comoving horizon 1/(a*H)
  - field f
- - canonical momemtum of field pi
+ - canonical momentum of field pi
  - energy density rho
  - fractional energy densities of fields rho_field(n)/rho_total (omega_field(n) in Silo output)
  - fractional energy density of the interaction term between fields (omega_int in Silo output)
  - Absolute and relative numerical errors
  - spatial correlation length of the total energy density (l_p in Silo output) (See Defrost paper
-   for a definition.)
+   for the definition.)
 
 When solving homogeneous equations for the fields the program writes
  - scale factor a
  - Hubble parameter H
  - Comoving horizon 1/(a*H)
  - homogeneous field f
- - homogeneous canonical momemtum of field pi
+ - homogeneous canonical momentum of field pi
  - energy density rho
    These are labeled by adding '_hom' at the end of the variable name in Silo output.
  In addition the program writes
@@ -241,7 +249,8 @@ Most of these functions are defined/used in 'Dmitry I. Podolsky et al. : Equatio
 Beginning of Thermalization After Preheating' http://arxiv.org/abs/hep-ph/0507096v1 .
 
 N.B. These functions have been tested but not thoroughly. If you notice that the output
-is clearly wrong please submit this to the issues page at github.
+is clearly wrong please submit this to the issues page at github and/or create
+a functioning fork.
 
 ------------------------------------------------------
 

integrator/symp_integrator.py

@@ -402,7 +402,7 @@ class Simulation:
         - scale parameter a
         - canonical momentum p
         - physical time t
-        - number of steps taken i0
+        - number of steps taken in linearized kernel i0
 
         - initial radiation energy density rho_r0
         - initial matter energy density rho_m0
@@ -676,7 +676,7 @@ def calc_k2_bins(self, lat):
            that the last of the k2_bins are set to -1 and lead to clearly
            unphysical solutions and might become infinite. They are
            not however used when evolving the actual scalar field
-           perturbations andcause no problems."""
+           perturbations and cause no problems."""
 
         if new_len > len(k2_bins):
             for i in xrange(new_len-len(k2_bins)):
@@ -1143,13 +1143,21 @@ def reinit(self, model, lat, a, init_m = 'defrost_cpu'):
         self.p = -6.*lat.VL_reduced*self.H*self.a**2.
         self.p_list = [self.p]
 
-        self.t_gpu = gpuarray.to_gpu(np.array(self.t, dtype = lat.prec_real))
-        self.a_gpu = gpuarray.to_gpu(np.array(self.a, dtype = lat.prec_real))
-        self.p_gpu = gpuarray.to_gpu(np.array(self.p, dtype = lat.prec_real))
+        #self.t_gpu = gpuarray.to_gpu(np.array(self.t, dtype = lat.prec_real))
+        #self.a_gpu = gpuarray.to_gpu(np.array(self.a, dtype = lat.prec_real))
+        #self.p_gpu = gpuarray.to_gpu(np.array(self.p, dtype = lat.prec_real))
 
         for field in self.fields:
             field.sample_field(lat, a, init_m)
 
+        "Reset variables used in linearized evo:"
+        if self.lin_evo:
+            cuda.memcpy_htod(self.a_gpu.gpudata, np.array(a,lat.prec_real))
+            cuda.memcpy_htod(self.p_gpu.gpudata, np.array(self.p,lat.prec_real))
+            cuda.memcpy_htod(self.t_gpu.gpudata,
+                             np.array(model.t_in,lat.prec_real))
+
+
         "Clear different list:"
         self.omega_rad_list = []
         self.omega_mat_list = []
@@ -1190,6 +1198,11 @@ def reinit(self, model, lat, a, init_m = 'defrost_cpu'):
             field.kurt_list = []
             field.w_list = []
             field.omega_list = []
+            if self.lin_evo:
+                field.f0_field[0] = field.f0_in
+                field.pi0_field[0] = field.pi0_in
+                cuda.memcpy_htod(field.f0_gpu.gpudata, field.f0_field)
+                cuda.memcpy_htod(field.pi0_gpu.gpudata, field.pi0_field)
 
     def set_non_lin(self):
         """Set the values of a and p equal to the values
@@ -1273,6 +1286,12 @@ def __init__(self, lat, V, f0, pi0, field_i, m2_eff, a_in, init_m, hom_m):
         self.f0_gpu = gpuarray.to_gpu(self.f0_field)
         self.pi0_gpu = gpuarray.to_gpu(self.pi0_field)
 
+        """These arrays are used in the linearized kernel when solving
+           perturbation equations with intial values (0,1) and (1,0)
+           for f_i^(1) and pi_i^(1) respectively.
+           Note that these are set to zeros and ones
+           automatically in the CUDA kernel before
+           evolving the system:"""
         self.f_lin_01_gpu = gpuarray.to_gpu(self.zeros)
         self.pi_lin_01_gpu = gpuarray.to_gpu(self.zeros)
 
@@ -1711,6 +1730,7 @@ def k_to_x_space(self, lat, sim, unperturb = False):
             for field in sim.fields:
                 field.unperturb_field()
 
+        "Update memory ids:"
         self.update(lat, sim)
 
     def lin_evo_step(self, lat, V, sim):
@@ -1846,6 +1866,7 @@ def x_to_k_space(self, lat, sim, perturb = False):
         for field in sim.fields:
             field.fft()
 
+        "Update memory ids:"
         self.update(lat, sim)