mirror of
https://github.com/cosmo-sims/monofonIC.git
synced 2024-09-19 17:03:45 +02:00
some cleanup, moved testing output routines outside of main loop
This commit is contained in:
Oliver Hahn
parent
b29bb1ab18
commit
5beed1b498
4 changed files with 121 additions and 92 deletions
|
|
@ -37,5 +37,5 @@ sigma_8 = 0.811
|
|||
nspec = 0.961
|
||||
|
||||
[sch]
|
||||
hbar = 100.0
|
||||
dt = 5.0
|
||||
hbar = 0.5
|
||||
dt = 1.0
|
||||
|
|
|
|||
11
include/testing.hh
Normal file
11
include/testing.hh
Normal file
|
|
@ -0,0 +1,11 @@
|
|||
#pragma once
|
||||
|
||||
namespace testing{
|
||||
void output_potentials_and_densities(
|
||||
size_t ngrid, real_t boxlen,
|
||||
const Grid_FFT<real_t>& phi,
|
||||
const Grid_FFT<real_t>& phi2,
|
||||
const Grid_FFT<real_t>& phi3a,
|
||||
const Grid_FFT<real_t>& phi3b,
|
||||
const std::array< Grid_FFT<real_t>*,3 >& A3 );
|
||||
}
|
||||
|
|
@ -2,6 +2,7 @@
|
|||
#include <general.hh>
|
||||
#include <grid_fft.hh>
|
||||
#include <convolution.hh>
|
||||
#include <testing.hh>
|
||||
|
||||
#include <ic_generator.hh>
|
||||
|
||||
|
|
@ -74,8 +75,6 @@ int Run( ConfigFile& the_config )
|
|||
Grid_FFT<real_t> A3x({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen});
|
||||
Grid_FFT<real_t> A3y({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen});
|
||||
Grid_FFT<real_t> A3z({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen});
|
||||
Grid_FFT<ccomplex_t> psi({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen});
|
||||
Grid_FFT<real_t> rho({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen});
|
||||
//... array [.] access to components of A3:
|
||||
std::array< Grid_FFT<real_t>*,3 > A3({&A3x,&A3y,&A3z});
|
||||
|
||||
|
|
@ -232,105 +231,35 @@ int Run( ConfigFile& the_config )
|
|||
|
||||
csoca::ilog << "-----------------------------------------------------------------------------" << std::endl;
|
||||
|
||||
|
||||
///////////////////////////////////////////////////////////////////////
|
||||
// we store the densities here if we compute them
|
||||
//======================================================================
|
||||
const bool compute_densities = false;
|
||||
if( compute_densities ){
|
||||
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");
|
||||
|
||||
Grid_FFT<real_t> delta({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen});
|
||||
Grid_FFT<real_t> delta2({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen});
|
||||
Grid_FFT<real_t> delta3a({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen});
|
||||
Grid_FFT<real_t> delta3b({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen});
|
||||
Grid_FFT<real_t> delta3({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen});
|
||||
delta.FourierTransformForward(false);
|
||||
delta2.FourierTransformForward(false);
|
||||
delta3a.FourierTransformForward(false);
|
||||
delta3b.FourierTransformForward(false);
|
||||
delta3.FourierTransformForward(false);
|
||||
|
||||
#pragma omp parallel for
|
||||
for (size_t i = 0; i < phi.size(0); ++i)
|
||||
{
|
||||
for (size_t j = 0; j < phi.size(1); ++j)
|
||||
{
|
||||
for (size_t k = 0; k < phi.size(2); ++k)
|
||||
{
|
||||
auto kk = phi.get_k<real_t>(i,j,k);
|
||||
size_t idx = phi.get_idx(i,j,k);
|
||||
auto laplace = -kk.norm_squared();
|
||||
|
||||
// compute densities associated to respective potentials as well
|
||||
delta.kelem(idx) = laplace * phi.kelem(idx);
|
||||
delta2.kelem(idx) = laplace * phi2.kelem(idx);
|
||||
delta3a.kelem(idx) = laplace * phi3a.kelem(idx);
|
||||
delta3b.kelem(idx) = laplace * phi3b.kelem(idx);
|
||||
delta3.kelem(idx) = delta3a.kelem(idx) + delta3b.kelem(idx);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
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();
|
||||
phi3b.FourierTransformBackward();
|
||||
|
||||
delta.FourierTransformBackward();
|
||||
delta2.FourierTransformBackward();
|
||||
delta3a.FourierTransformBackward();
|
||||
delta3b.FourierTransformBackward();
|
||||
delta3.FourierTransformBackward();
|
||||
|
||||
A3[0]->FourierTransformBackward();
|
||||
A3[1]->FourierTransformBackward();
|
||||
A3[2]->FourierTransformBackward();
|
||||
|
||||
#if defined(USE_MPI)
|
||||
if( CONFIG::MPI_task_rank == 0 )
|
||||
unlink(fname_hdf5.c_str());
|
||||
MPI_Barrier( MPI_COMM_WORLD );
|
||||
#else
|
||||
unlink(fname_hdf5.c_str());
|
||||
#endif
|
||||
|
||||
phi.Write_to_HDF5(fname_hdf5, "phi");
|
||||
phi2.Write_to_HDF5(fname_hdf5, "phi2");
|
||||
phi3a.Write_to_HDF5(fname_hdf5, "phi3a");
|
||||
phi3b.Write_to_HDF5(fname_hdf5, "phi3b");
|
||||
|
||||
delta.Write_to_HDF5(fname_hdf5, "delta");
|
||||
delta2.Write_to_HDF5(fname_hdf5, "delta2");
|
||||
delta3a.Write_to_HDF5(fname_hdf5, "delta3a");
|
||||
delta3b.Write_to_HDF5(fname_hdf5, "delta3b");
|
||||
delta3.Write_to_HDF5(fname_hdf5, "delta3");
|
||||
|
||||
A3[0]->Write_to_HDF5(fname_hdf5, "A3x");
|
||||
A3[1]->Write_to_HDF5(fname_hdf5, "A3y");
|
||||
A3[2]->Write_to_HDF5(fname_hdf5, "A3z");
|
||||
|
||||
}else{
|
||||
const bool testing_compute_densities = false;
|
||||
if( testing_compute_densities )
|
||||
{
|
||||
testing::output_potentials_and_densities( ngrid, boxlen, phi, phi2, phi3a, phi3b, A3 );
|
||||
}
|
||||
else
|
||||
{
|
||||
// temporary storage of data
|
||||
Grid_FFT<real_t> tmp({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen});
|
||||
|
||||
|
||||
if( the_output_plugin->write_species_as( cosmo_species::dm ) == output_type::field_eulerian ){
|
||||
//======================================================================
|
||||
// use QPT to get density and velocity fields
|
||||
//======================================================================
|
||||
Grid_FFT<ccomplex_t> psi({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen});
|
||||
Grid_FFT<real_t> rho({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen});
|
||||
|
||||
//======================================================================
|
||||
// initialise psi = exp(i Phi(1)/hbar)
|
||||
//======================================================================
|
||||
real_t hbar= the_config.GetValueSafe<real_t>("sch", "hbar", 0.01);
|
||||
real_t hbar= the_config.GetValueSafe<real_t>("sch", "hbar", 0.000001);
|
||||
real_t dt = the_config.GetValueSafe<real_t>("sch", "dt", 1.0);
|
||||
phi.FourierTransformBackward();
|
||||
psi.assign_function_of_grids_r([&]( real_t p ){return std::exp(ccomplex_t(0.0,1.0)/hbar*p);}, phi );
|
||||
|
||||
psi.assign_function_of_grids_r([&]( real_t p ){return std::exp(-ccomplex_t(0.0,1.0)/hbar*p);}, phi );
|
||||
phi.FourierTransformForward();
|
||||
//======================================================================
|
||||
// evolve wave-function (one drift step) psi = psi *exp(-i hbar *k^2 dt / 2)
|
||||
//======================================================================
|
||||
|
|
@ -345,7 +274,7 @@ int Run( ConfigFile& the_config )
|
|||
// compute rho
|
||||
//======================================================================
|
||||
rho.assign_function_of_grids_r([&]( auto p ){
|
||||
auto pp = std::real(p)*std::real(p) + std::imag(p)*std::imag(p);
|
||||
auto pp = std::real(p)*std::real(p) + std::imag(p)*std::imag(p) - 1.0;
|
||||
return pp;
|
||||
}, psi);
|
||||
|
||||
|
|
@ -365,7 +294,7 @@ int Run( ConfigFile& the_config )
|
|||
grad_psi.FourierTransformBackward();
|
||||
|
||||
tmp.assign_function_of_grids_r([&](auto ppsi, auto pgrad_psi, auto prho) {
|
||||
return std::real((std::conj(ppsi) * pgrad_psi - ppsi * std::conj(pgrad_psi)) / ccomplex_t(0.0, 2.0 / hbar)/prho);
|
||||
return std::real((std::conj(ppsi) * pgrad_psi - ppsi * std::conj(pgrad_psi)) / ccomplex_t(0.0, 2.0 / hbar)/(1.0+prho));
|
||||
}, psi, grad_psi, rho);
|
||||
|
||||
fluid_component fc = (idim==0)? fluid_component::vx : ((idim==1)? fluid_component::vy : fluid_component::vz );
|
||||
|
|
|
|||
89
src/testing.cc
Normal file
89
src/testing.cc
Normal file
|
|
@ -0,0 +1,89 @@
|
|||
#include <testing.hh>
|
||||
|
||||
void output_potentials_and_densities(
|
||||
size_t ngrid, real_t boxlen,
|
||||
const Grid_FFT<real_t> &phi,
|
||||
const Grid_FFT<real_t> &phi2,
|
||||
const Grid_FFT<real_t> &phi3a,
|
||||
const Grid_FFT<real_t> &phi3b,
|
||||
const std::array<Grid_FFT<real_t> *, 3> &A3)
|
||||
{
|
||||
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");
|
||||
|
||||
Grid_FFT<real_t> delta({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen});
|
||||
Grid_FFT<real_t> delta2({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen});
|
||||
Grid_FFT<real_t> delta3a({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen});
|
||||
Grid_FFT<real_t> delta3b({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen});
|
||||
Grid_FFT<real_t> delta3({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen});
|
||||
delta.FourierTransformForward(false);
|
||||
delta2.FourierTransformForward(false);
|
||||
delta3a.FourierTransformForward(false);
|
||||
delta3b.FourierTransformForward(false);
|
||||
delta3.FourierTransformForward(false);
|
||||
|
||||
#pragma omp parallel for
|
||||
for (size_t i = 0; i < phi.size(0); ++i)
|
||||
{
|
||||
for (size_t j = 0; j < phi.size(1); ++j)
|
||||
{
|
||||
for (size_t k = 0; k < phi.size(2); ++k)
|
||||
{
|
||||
auto kk = phi.get_k<real_t>(i, j, k);
|
||||
size_t idx = phi.get_idx(i, j, k);
|
||||
auto laplace = -kk.norm_squared();
|
||||
|
||||
// compute densities associated to respective potentials as well
|
||||
delta.kelem(idx) = laplace * phi.kelem(idx);
|
||||
delta2.kelem(idx) = laplace * phi2.kelem(idx);
|
||||
delta3a.kelem(idx) = laplace * phi3a.kelem(idx);
|
||||
delta3b.kelem(idx) = laplace * phi3b.kelem(idx);
|
||||
delta3.kelem(idx) = delta3a.kelem(idx) + delta3b.kelem(idx);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
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();
|
||||
phi3b.FourierTransformBackward();
|
||||
|
||||
delta.FourierTransformBackward();
|
||||
delta2.FourierTransformBackward();
|
||||
delta3a.FourierTransformBackward();
|
||||
delta3b.FourierTransformBackward();
|
||||
delta3.FourierTransformBackward();
|
||||
|
||||
A3[0]->FourierTransformBackward();
|
||||
A3[1]->FourierTransformBackward();
|
||||
A3[2]->FourierTransformBackward();
|
||||
|
||||
#if defined(USE_MPI)
|
||||
if (CONFIG::MPI_task_rank == 0)
|
||||
unlink(fname_hdf5.c_str());
|
||||
MPI_Barrier(MPI_COMM_WORLD);
|
||||
#else
|
||||
unlink(fname_hdf5.c_str());
|
||||
#endif
|
||||
|
||||
phi.Write_to_HDF5(fname_hdf5, "phi");
|
||||
phi2.Write_to_HDF5(fname_hdf5, "phi2");
|
||||
phi3a.Write_to_HDF5(fname_hdf5, "phi3a");
|
||||
phi3b.Write_to_HDF5(fname_hdf5, "phi3b");
|
||||
|
||||
delta.Write_to_HDF5(fname_hdf5, "delta");
|
||||
delta2.Write_to_HDF5(fname_hdf5, "delta2");
|
||||
delta3a.Write_to_HDF5(fname_hdf5, "delta3a");
|
||||
delta3b.Write_to_HDF5(fname_hdf5, "delta3b");
|
||||
delta3.Write_to_HDF5(fname_hdf5, "delta3");
|
||||
|
||||
A3[0]->Write_to_HDF5(fname_hdf5, "A3x");
|
||||
A3[1]->Write_to_HDF5(fname_hdf5, "A3y");
|
||||
A3[2]->Write_to_HDF5(fname_hdf5, "A3z");
|
||||
}
|
||||
Loading…
Reference in a new issue