lesgourg / class_public

Public repository of the Cosmic Linear Anisotropy Solving System (master for the most recent version of the standard code; GW_CLASS to include Cosmic Gravitational Wave Background anisotropies; classnet branch for acceleration with neutral networks; ExoCLASS branch for exotic energy injection; class_matter branch for FFTlog)
223 stars 291 forks source link

How to add delta_fld to the perturbations dictionary? #285

Closed rodrigovonmarttens closed 5 years ago

rodrigovonmarttens commented 5 years ago

I'm trying to print DE perturbation for a dynamical DE EoS (fld approach), so, I tried to add two columns in the perturbations output related to delta_fld and theta_fld. To do that, I added two new titles for the fld quantities:

class_store_columntitle(ppt->scalar_titles, "delta_fld", pba->has_fld); class_store_columntitle(ppt->scalar_titles, "theta_fld", pba->has_fld);

and I stored the solutions

class_store_double(dataptr, delta_fld, pba->has_fld, storeidx); class_store_double(dataptr, theta_fld, pba->has_fld, storeidx);

where delta_fld and theta_fld are defined with their respective solutions y[ppw->pv->...].

However, when I run the code I got "Segmentation fault (core dumped)". Am I missing something?

Stefan-Heimersheim commented 5 years ago

Hi,

I'm not sure where changed the code but I tried to reproduce your problem and get no segmentation fault. Did you check that the index you use for y is valid?

delta_fld = y[ppw->pv->index_pt_delta_fld];
theta_fld = y[ppw->pv->index_pt_theta_fld];

I just implemented everything analogous to delta_b and theta_b, here is a git diff:

@@ -2517,6 +2517,9 @@ int perturb_prepare_output(struct background * pba,
       class_store_columntitle(ppt->scalar_titles,"pol2_g",_TRUE_);
       class_store_columntitle(ppt->scalar_titles,"delta_b",_TRUE_);
       class_store_columntitle(ppt->scalar_titles,"theta_b",_TRUE_);
+      /* fld */
+      class_store_columntitle(ppt->scalar_titles,"delta_fld",pba->has_fld);
+      class_store_columntitle(ppt->scalar_titles,"theta_fld",pba->has_fld);
       class_store_columntitle(ppt->scalar_titles,"psi",_TRUE_);
       class_store_columntitle(ppt->scalar_titles,"phi",_TRUE_);
       /* Perturbed recombination */
@@ -6357,6 +6360,7 @@ int perturb_print_variables(double tau,

   double delta_g,theta_g,shear_g,l4_g,pol0_g,pol1_g,pol2_g,pol4_g;
   double delta_b,theta_b;
+  double theta_fld, delta_fld;
   double delta_cdm=0.,theta_cdm=0.;
   double delta_dcdm=0.,theta_dcdm=0.;
   double delta_dr=0.,theta_dr=0.,shear_dr=0., f_dr=1.0;
@@ -6503,6 +6507,9 @@ int perturb_print_variables(double tau,
     delta_b = y[ppw->pv->index_pt_delta_b];
     theta_b = y[ppw->pv->index_pt_theta_b];

+    delta_fld = y[ppw->pv->index_pt_delta_fld];
+    theta_fld = y[ppw->pv->index_pt_theta_fld];
+
     if (pba->has_cdm == _TRUE_) {

       delta_cdm = y[ppw->pv->index_pt_delta_cdm];
@@ -6698,6 +6705,8 @@ int perturb_print_variables(double tau,
     class_store_double(dataptr, pol2_g, _TRUE_, storeidx);
     class_store_double(dataptr, delta_b, _TRUE_, storeidx);
     class_store_double(dataptr, theta_b, _TRUE_, storeidx);
+    class_store_double(dataptr, delta_fld, _TRUE_, storeidx);
+    class_store_double(dataptr, theta_fld, _TRUE_, storeidx);
     class_store_double(dataptr, psi, _TRUE_, storeidx);
     class_store_double(dataptr, phi, _TRUE_, storeidx);
     /* perturbed recombination */

Cheers, Stefan

rodrigovonmarttens commented 5 years ago

Hi Stefan!

Thanks for your response!

I did exactly what you showed in your comment, but I still get a segmentation fault error. The fact I'm using a mac (clang compiler) would be an issue?

best wishes,

Rodrigo

Stefan-Heimersheim commented 5 years ago

Hi Rodrigo,

I just compiled my code with clang and it still works without problems (with explanatory.ini). I don't have a Mac so I cannot test your case but I think it is unlikely that this is the cause.

What ini file are you using for your test? I just tested the default options.

Best, Stefan

rodrigovonmarttens commented 5 years ago

Hi Stefan,

I have used the explonatory.ini, but I set Omega_Lambda to zero and commented Omega_fld. I also set w0_fld to -0.9 and wa_fld to -0.1. Have you used the fld description for DE? Indeed I only have segmentation fault if I use fld description for DE, otherwise, delta_fld and theta_fld are identically zero in the output.

best wishes,

Rodrigo

Stefan-Heimersheim commented 5 years ago

Hi Rodrigo,

I tried the same changes to the ini file:

@@ -165,8 +165,8 @@ Omega_k = 0.
     (default: 'Omega_fld' and 'Omega_scf' set to 0 and 'Omega_Lambda' inferred
     by code)

-# Omega_Lambda = 0.7
-Omega_fld = 0
+Omega_Lambda = 0
+# Omega_fld = 0
 Omega_scf = 0

 8b) The flag 'use_ppf' is 'yes' by default, to use the PPF approximation
@@ -196,7 +196,7 @@ fluid_equation_of_state = CLP
     'wa_fld' to 0, 'cs2_fls' to 1)

 w0_fld = -0.9
-wa_fld = 0.
+wa_fld = -0.1
 cs2_fld = 1

 8c2) equation of state of the fluid in 'EDE' case and squared sound speed 'cs2_fld' of the fluid
@@ -207,9 +207,9 @@ cs2_fld = 1
     the function background_w_fld().  (default: 'w0_fld' set to -1,
     'Omega_EDE' to 0, 'cs2_fls' to 1)

-w0_fld = -0.9
-Omega_EDE = 0.
-cs2_fld = 1
+# w0_fld = -0.9
+# Omega_EDE = 0.
+# cs2_fld = 1

 8d) Scalar field (scf) initial conditions from attractor solution (assuming
     pure exponential potential). (default: yes)

and the code does still run without segfault, did I forget an option there or is this the same config you are using?

rodrigovonmarttens commented 5 years ago

Hi Stefan,

Indeed there is one more option I forgot to mention, in order to compute the perturbations output, you must set the option k_output to some value. In my case, I set k_output to 1. Could you please try this?

best wishes,

Rodrigo

Stefan-Heimersheim commented 5 years ago

Ah, I forgot that option, now I can reproduce the error. You have to add use_ppf = no to your ini file to fix it.

The reason is the following: The delta_fld and theta_fld indices are only defined if the ppf approximation is not used:

/* fluid */

if (pba->use_ppf == _FALSE_) {
  class_define_index(ppv->index_pt_delta_fld,pba->has_fld,index_pt,1); /* fluid density */
  class_define_index(ppv->index_pt_theta_fld,pba->has_fld,index_pt,1); /* fluid velocity */
}
else {
  class_define_index(ppv->index_pt_Gamma_fld,pba->has_fld,index_pt,1); /* Gamma variable of PPF scheme */
}

Therefore you have to add use_ppf = no to the ini file. Then the indices are allocated and you can use them. Otherwise accessing y[ppw->pv->index_pt_delta_fld] produces the segfault.

Best, Stefan

rodrigovonmarttens commented 5 years ago

Thanks once again Stefan!

It would be possible to compute the delta_fld for a case where PPF is necessary? For example, some choices of w0_fld and wa_fld for CPL parameterization.

best wishes,

Rodrigo

Stefan-Heimersheim commented 5 years ago

I'm not sure how delta and theta are defined in the PPF description and if the fld equations implemented in CLASS are valid in these cases.

The PPF description of dark energy replaces the den- sity and momentum components with a single joint dy- namical variable Γ [...]

To check this you could search for pba->use_ppf in the code and change those places, defining the indices and setting delta in a meaningful way. I have no experience with this though.

rodrigovonmarttens commented 5 years ago

I understand. Thanks a lot for your help Stefan!