1
0
Fork 0
mirror of https://github.com/cosmo-sims/monofonIC.git synced 2024-09-19 17:03:45 +02:00

on the fly output of power spectra

This commit is contained in:
Oliver Hahn 2019-05-19 12:05:04 +02:00
parent 7e9bbc6dc6
commit 06ed21cae0
3 changed files with 46 additions and 6 deletions

View file

@ -430,9 +430,11 @@ public:
void Write_to_HDF5(std::string fname, std::string datasetname);
void ComputePowerSpectrum(std::vector<double> &bin_k, std::vector<double> &bin_P, std::vector<double> &bin_eP, int nbins);
void Write_PowerSpectrum( std::string ofname );
void Compute_PDF(std::string ofname, int nbins = 1000, double scale = 1.0, double rhomin = 1e-3, double rhomax = 1e3);
void Compute_PowerSpectrum(std::vector<double> &bin_k, std::vector<double> &bin_P, std::vector<double> &bin_eP, std::vector<size_t> &bin_count, int nbins);
void Write_PDF(std::string ofname, int nbins = 1000, double scale = 1.0, double rhomin = 1e-3, double rhomax = 1e3);
void zero_DC_mode(void)
{

View file

@ -501,7 +501,7 @@ void Grid_FFT<data_t>::Write_to_HDF5(std::string fname, std::string datasetname)
#include <iomanip>
template <typename data_t>
void Grid_FFT<data_t>::Compute_PDF( std::string ofname, int nbins, double scale, double vmin, double vmax )
void Grid_FFT<data_t>::Write_PDF( std::string ofname, int nbins, double scale, double vmin, double vmax )
{
double logvmin = std::log10(vmin);
double logvmax = std::log10(vmax);
@ -551,15 +551,44 @@ void Grid_FFT<data_t>::Compute_PDF( std::string ofname, int nbins, double scale,
}
template <typename data_t>
void Grid_FFT<data_t>::ComputePowerSpectrum(std::vector<double> &bin_k, std::vector<double> &bin_P, std::vector<double> &bin_eP, int nbins)
void Grid_FFT<data_t>::Write_PowerSpectrum( std::string ofname )
{
std::vector<double> bin_k, bin_P, bin_eP;
std::vector<size_t> bin_count;
int nbins = 4*std::max( nhalf_[0], std::max(nhalf_[1], nhalf_[2]) );
this->Compute_PowerSpectrum( bin_k, bin_P, bin_eP, bin_count, nbins );
#if defined(USE_MPI)
if (CONFIG::MPI_task_rank == 0)
{
#endif
std::ofstream ofs( ofname.c_str());
std::size_t numcells = size(0) * size(1) * size(2);
//ofs << "# a = " << aexp << std::endl;
ofs << "# " << std::setw(14) << "k" << std::setw(16) << "P(k)" << std::setw(16) << "err. P(k)"
<< std::setw(16) <<"#modes" << "\n";
for( int ibin=0; ibin<nbins; ++ibin ){
if( bin_count[ibin] > 0 )
ofs << std::setw(16) << bin_k[ibin]
<< std::setw(16) << bin_P[ibin]
<< std::setw(16) << bin_eP[ibin]
<< std::setw(16) << bin_count[ibin]
<< std::endl;
}
#if defined(USE_MPI)
}
#endif
}
template <typename data_t>
void Grid_FFT<data_t>::Compute_PowerSpectrum(std::vector<double> &bin_k, std::vector<double> &bin_P, std::vector<double> &bin_eP, std::vector<size_t> &bin_count, int nbins)
{
real_t kmax = std::max(std::max(kfac_[0] * nhalf_[0], kfac_[1] * nhalf_[1]),
kfac_[2] * nhalf_[2]),
kmin = std::min(std::min(kfac_[0], kfac_[1]), kfac_[2]),
dklog = log10(kmax / kmin) / nbins;
std::vector<size_t> bin_count;
bin_count.assign(nbins,0);
bin_k.assign(nbins, 0);
bin_P.assign(nbins, 0);

View file

@ -15,6 +15,7 @@
#include <output_plugin.hh>
#include <cosmology_calculator.hh>
// initialise with "default" values
namespace CONFIG{
int MPI_thread_support = -1;
int MPI_task_rank = 0;
@ -94,6 +95,7 @@ int main( int argc, char** argv )
//...
const std::string fname_hdf5 = the_config.GetValueSafe<std::string>("output", "fname_hdf5", "output.hdf5");
const std::string fname_analysis = the_config.GetValueSafe<std::string>("output", "fbase_analysis", "output");
//////////////////////////////////////////////////////////////////////////////////////////////
std::unique_ptr<CosmologyCalculator> the_cosmo_calc;
@ -276,6 +278,13 @@ int main( int argc, char** argv )
}
}
}
delta.Write_PowerSpectrum(fname_analysis+"_"+"power_delta1.txt");
delta2.Write_PowerSpectrum(fname_analysis+"_"+"power_delta2.txt");
delta3a.Write_PowerSpectrum(fname_analysis+"_"+"power_delta3a.txt");
delta3b.Write_PowerSpectrum(fname_analysis+"_"+"power_delta3b.txt");
delta3.Write_PowerSpectrum(fname_analysis+"_"+"power_delta3.txt");
phi.FourierTransformBackward();
phi2.FourierTransformBackward();
phi3a.FourierTransformBackward();