perturbations.h

Go to the documentation of this file.

 
#ifndef __PERTURBATIONS__
#define __PERTURBATIONS__
 
#include <memory>
#include <vector>
 
#include "../species/base_species.h"
#include "evolver_ndf15.h"
#include "evolver_rkck.h"
#include "evolver_rkdp45.h"
#include "thermodynamics.h"
 
#define _scalars_ (ppt->has_scalars && (index_md == index_md_scalars_))
#define _vectors_ (ppt->has_vectors && (index_md == index_md_vectors_))
#define _tensors_ (ppt->has_tensors && (index_md == index_md_tensors_))
 
//TODO: Remove those when possible!
#define _scalarsEXT_ (ppt->has_scalars && (index_md == perturbations_module_->index_md_scalars_))
#define _vectorsEXT_ (ppt->has_vectors && (index_md == perturbations_module_->index_md_vectors_))
#define _tensorsEXT_ (ppt->has_tensors && (index_md == perturbations_module_->index_md_tensors_))
 
#define _set_source_(index) \
  sources_[index_md][index_ic * tp_size_[index_md] + index][index_tau * k_size_[index_md] + index_k]
 
 
enum tca_flags { tca_on, tca_off };
enum rsa_flags { rsa_off, rsa_on };
enum tca_idm_dr_flags { tca_idm_dr_on, tca_idm_dr_off };
enum tca_idm_drmd_flags { tca_idm_drmd_on, tca_idm_drmd_off };
enum rsa_idr_flags { rsa_idr_off, rsa_idr_on };
enum ufa_flags { ufa_off, ufa_on };
enum ncdmfa_flags { ncdmfa_off, ncdmfa_on };
 
 
 
enum idr_method { idr_free_streaming, idr_fluid }; /* for the idm-idr case */
enum tensor_methods { tm_photons_only, tm_massless_approximation, tm_exact };
 
 
 
enum class possible_gauges {
  newtonian,  
  synchronous 
};
 
 
 
#define _SELECTION_NUM_MAX_ 100
enum selection_type { gaussian, tophat, dirac };
 
 
 
#define _MAX_NUMBER_OF_K_FILES_ 30
 
 
struct perturbs {
 
  bool has_perturbations; 
  bool
      has_cls; 
  bool has_scalars = true;  
  bool has_vectors = false; 
  bool has_tensors = false; 
  bool has_ad  = true;  
  bool has_bi  = false; 
  bool has_cdi = false; 
  bool has_nid = false; 
  bool has_niv = false; 
  /* perturbed recombination */
  bool has_perturbed_recombination = false;
  enum tensor_methods tensor_method =
      tm_massless_approximation; 
  bool has_cl_cmb_temperature  = false; 
  bool has_cl_cmb_polarization = false; 
  bool has_cl_cmb_lensing_potential =
      false; 
  bool want_lcmb_full_limber =
      true; 
  bool has_cl_lensing_potential =
      false;                        
  bool has_cl_number_count = false; 
  bool has_pk_matter       = false; 
  bool has_density_transfers =
      false; 
  bool has_velocity_transfers =
      false; 
  bool has_metricpotential_transfers =
      false; 
  bool has_Nbody_gauge_transfers =
      false; 
  bool has_nl_corrections_based_on_delta_m =
      false; 
  bool has_nc_density = false; 
  bool has_nc_rsd     = false; 
  bool has_nc_lens    = false; 
  bool has_nc_gr      = false; 
  int l_scalar_max = 2500; 
  int l_vector_max = 500;  
  int l_tensor_max = 500;  
  int l_lss_max =
      300; 
  double k_max_for_pk =
      1.; 
  int selection_num                          = 1;        
  enum selection_type selection              = gaussian; 
  double selection_mean[_SELECTION_NUM_MAX_] = {
      1.}; 
  double selection_width[_SELECTION_NUM_MAX_] = {
      0.1}; 
  int switch_sw =
      1; 
  int switch_eisw =
      1; 
  int switch_lisw =
      1; 
  int switch_dop = 1; 
  int switch_pol =
      1; 
  double eisw_lisw_split_z =
      120; 
  int store_perturbations = false; 
  int k_output_values_num = 0;     
  double k_output_values
      [_MAX_NUMBER_OF_K_FILES_]; 
  double three_ceff2_ur =
      1.; 
  double three_cvis2_ur =
      1.; 
  double G_eff_ur = 0.;
 
  double z_max_pk =
      0.; 
 
 
  possible_gauges gauge =
      possible_gauges::synchronous; 
 
 
  double
      selection_min_of_tau_min; 
  double
      selection_max_of_tau_max; 
  double
      selection_delta_tau; 
  double*
      selection_tau_min; 
  double*
      selection_tau_max; 
  double* selection_tau; 
  double*
      selection_function; 
 
 
  short perturbations_verbose =
      0; 
};
 
struct perturb_vector {
  // Per-species perturbation layouts, parallel to all_species_ (lex-key order).
  // Each thread allocates its own pv, so this storage is per-thread.
  std::vector<std::unique_ptr<BaseSpecies::PerturbLayout>> species_layouts;
 
  // Non-owning view into species_layouts: the species that registered ≥1
  // perturbation variable in THIS pv's mode (no-ops like Lambda excluded).
  // Scalar-mode construction keeps scalar fields last so their Newtonian
  // stress-energy evaluation sees the complete shear accumulated from other
  // species. Consumed by the scalar/vector/tensor RHS dispatch loops.
  struct ActiveSpecies {
    const BaseSpecies* species;
    const BaseSpecies::PerturbLayout* layout;
  };
  std::vector<ActiveSpecies> active_species;
 
  // Non-owning view into species_layouts: the species that emit a per-species
  // transfer source (delta_g, delta_cdm, delta_ncdm, ...), in lex-key order.
  // Source-column membership is k-independent (fixed in
  // perturb_indices_of_perturbs); this list materialises it per pv so
  // perturb_sources_member dispatches FillSources by pointer. Empty in runs
  // without dTk/mTk/k_output, so that per-sample loop is skipped entirely.
  struct SourceSpecies {
    const BaseSpecies* species;
    const BaseSpecies::PerturbLayout* layout;
  };
  std::vector<SourceSpecies> source_species;
 
  // Always-present species, resolved once per pv (non-owning, into
  // species_layouts). Base-typed because perturbations.h cannot see the nested
  // PhotonsSpecies/BaryonsSpecies::PerturbLayout types (photons.h/baryons.h
  // include perturbations.h). Read sites static_cast to the concrete type.
  BaseSpecies::PerturbLayout* photon_layout = nullptr;
  BaseSpecies::PerturbLayout* baryon_layout = nullptr;
  BaseSpecies::PerturbLayout* cdm_layout    = nullptr;
 
  /* Photons bare fields removed — use pv->species_layouts[g_i] (PhotonsSpecies::PerturbLayout). */
  /* Baryons bare fields removed — use pv->species_layouts[b_i] (BaryonsSpecies::PerturbLayout). */
  /* CDM bare fields removed — use pv->species_layouts[cdm_i] (CDMSpecies::PerturbLayout). */
  /* IDM_DR bare fields removed — use
     IDM_DR_IDR_Species::idm_dr_layout(*pv->species_layouts[idm_dr_idr_i]). */
  /* IDM_DRMD bare fields removed — use
     IDM_DRMD_IDR_DRMD_Species::idm_drmd_layout(*pv->species_layouts[idm_drmd_idr_drmd_i]). */
  /* DCDM bare fields removed — use DCDM_DR_Species::dcdm_layout(*pv->species_layouts[dcdm_dr_i]). */
  /* UR bare fields removed — use pv->species_layouts[ur_i] (UltraRelativisticSpecies::PerturbLayout). */
  /* IDR bare fields removed — use
     IDM_DR_IDR_Species::idr_layout(*pv->species_layouts[idm_dr_idr_i]). */
  /* IDR_DRMD bare fields removed — use
     IDM_DRMD_IDR_DRMD_Species::idr_drmd_layout(*pv->species_layouts[idm_drmd_idr_drmd_i]). */
 
  /* perturbed recombination */
  int index_pt_perturbed_recombination_delta_temp; 
  int index_pt_perturbed_recombination_delta_chi;  
  int index_pt_eta;      
  int index_pt_phi;      
  int index_pt_hv_prime; 
  int index_pt_V;        
  int index_pt_gw;    
  int index_pt_gwdot; 
  struct TensorRelativisticLayout {
    int idx_delta = -1; 
    int idx_theta = -1; 
    int idx_shear = -1; 
    int idx_l3    = -1; 
    int l_max     = -1; 
  } tensor_ur_layout;
 
  int pt_size; 
  std::vector<double> y;  
  std::vector<double> dy; 
  std::vector<int> used_in_sources; 
};
 
struct perturb_workspace {
 
  int index_mt_psi;           
  int index_mt_phi_prime;     
  int index_mt_h_prime;       
  int index_mt_h_prime_prime; 
  int index_mt_eta_prime;     
  int index_mt_alpha;         
  int index_mt_alpha_prime;   
  int index_mt_gw_prime_prime; 
  int index_mt_V_prime; 
  int index_mt_hv_prime_prime; 
  int mt_size;                 
 
 
  std::vector<double> pvecback;   
  std::vector<double> pvecthermo; 
  std::vector<double> pvecmetric; 
  /* Single owner of the perturb_vector; species methods receive a
     non-owning raw pointer via pv.get() (used transiently while the
     workspace is alive, never stored or shared). */
  std::unique_ptr<perturb_vector> pv; 
  double delta_rho;        
  double rho_plus_p_theta; 
  double rho_plus_p_shear; 
  double delta_p;          
  double rho_plus_p_tot; 
  double
      gw_source; 
  double vector_source_pi; 
  double vector_source_v;  
  /* Module-owned approximation closures. Computed by the perturbations module
     (perturb_tca_slip_and_shear and perturb_rsa_* / perturb_rsa_idr_*, the
     latter run inside perturb_einstein after the metric solve because they
     depend on h'/phi') and then read — read-only — by the relevant species
     (photons, ultra-relativistic, IDM_DR). These are not per-species private
     state despite the species-suffixed names: they live here because the module
     owns the approximation machinery and several consumers cross species. */
  double tca_shear_g;  
  double tca_slip;     
  double rsa_delta_g;  
  double rsa_theta_g;  
  double rsa_delta_ur; 
  double rsa_theta_ur; 
  double
      rsa_delta_idr; 
  double
      rsa_theta_idr; 
  double delta_m; 
  double theta_m; 
  double delta_cb; 
  double theta_cb; 
  /* Module-owned PPF dark-energy closure. Written by FluidSpecies::ComputePpf
     during stress-energy assembly (it depends on the totals of every other
     species), then read back by the module totals and by the fluid's own
     derivs/sources. PPF fluid is defined relative to the global stress-energy,
     so this is genuine module-level coupling state, not fluid-private scratch.
     (The S quantity of eq. 15 in 0808.3125 is a ComputePpf-internal temporary,
     no longer stored here.) */
  double delta_rho_fld; 
  double delta_p_fld;   
  double rho_plus_p_theta_fld; 
  double Gamma_prime_fld;      
  FILE* perturb_output_file; 
  int index_ikout;           
 
  PerturbScalarContext scalar_ctx;
 
 
  short
      inter_mode; 
  int last_index_back; 
  int last_index_thermo; 
 
 
  int index_ap_tca;          
  int index_ap_rsa;          
  int index_ap_tca_idm_dr;   
  int index_ap_tca_idm_drmd; 
  int index_ap_rsa_idr;      
  int index_ap_ufa;          
  int index_ap_ncdmfa;       
  int ap_size;               
  std::vector<int>
      approx; 
 
 
  int max_l_max; 
  std::vector<double>
      s_l; 
  double
      cotKgen; 
};
 
#endif