From 5beed1b498b8a8f4b723a9473031c9e5508ff566 Mon Sep 17 00:00:00 2001 From: Oliver Hahn Date: Wed, 7 Aug 2019 18:43:00 +0200 Subject: [PATCH] some cleanup, moved testing output routines outside of main loop --- example.conf | 4 +- include/testing.hh | 11 +++++ src/ic_generator.cc | 109 ++++++++------------------------------------ src/testing.cc | 89 ++++++++++++++++++++++++++++++++++++ 4 files changed, 121 insertions(+), 92 deletions(-) create mode 100644 include/testing.hh create mode 100644 src/testing.cc diff --git a/example.conf b/example.conf index bdbd76e..37e7375 100644 --- a/example.conf +++ b/example.conf @@ -37,5 +37,5 @@ sigma_8 = 0.811 nspec = 0.961 [sch] -hbar = 100.0 -dt = 5.0 +hbar = 0.5 +dt = 1.0 diff --git a/include/testing.hh b/include/testing.hh new file mode 100644 index 0000000..7796a21 --- /dev/null +++ b/include/testing.hh @@ -0,0 +1,11 @@ +#pragma once + +namespace testing{ + void output_potentials_and_densities( + size_t ngrid, real_t boxlen, + const Grid_FFT& phi, + const Grid_FFT& phi2, + const Grid_FFT& phi3a, + const Grid_FFT& phi3b, + const std::array< Grid_FFT*,3 >& A3 ); +} diff --git a/src/ic_generator.cc b/src/ic_generator.cc index ddc875d..500572b 100644 --- a/src/ic_generator.cc +++ b/src/ic_generator.cc @@ -2,6 +2,7 @@ #include #include #include +#include #include @@ -74,8 +75,6 @@ int Run( ConfigFile& the_config ) Grid_FFT A3x({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen}); Grid_FFT A3y({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen}); Grid_FFT A3z({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen}); - Grid_FFT psi({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen}); - Grid_FFT rho({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen}); //... array [.] access to components of A3: std::array< Grid_FFT*,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("output", "fname_hdf5", "output.hdf5"); - const std::string fname_analysis = the_config.GetValueSafe("output", "fbase_analysis", "output"); - - Grid_FFT delta({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen}); - Grid_FFT delta2({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen}); - Grid_FFT delta3a({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen}); - Grid_FFT delta3b({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen}); - Grid_FFT 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(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 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 psi({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen}); + Grid_FFT rho({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen}); + //====================================================================== // initialise psi = exp(i Phi(1)/hbar) //====================================================================== - real_t hbar= the_config.GetValueSafe("sch", "hbar", 0.01); + real_t hbar= the_config.GetValueSafe("sch", "hbar", 0.000001); real_t dt = the_config.GetValueSafe("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 ); diff --git a/src/testing.cc b/src/testing.cc new file mode 100644 index 0000000..57be3a3 --- /dev/null +++ b/src/testing.cc @@ -0,0 +1,89 @@ +#include + +void output_potentials_and_densities( + size_t ngrid, real_t boxlen, + const Grid_FFT &phi, + const Grid_FFT &phi2, + const Grid_FFT &phi3a, + const Grid_FFT &phi3b, + const std::array *, 3> &A3) +{ + const std::string fname_hdf5 = the_config.GetValueSafe("output", "fname_hdf5", "output.hdf5"); + const std::string fname_analysis = the_config.GetValueSafe("output", "fbase_analysis", "output"); + + Grid_FFT delta({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen}); + Grid_FFT delta2({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen}); + Grid_FFT delta3a({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen}); + Grid_FFT delta3b({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen}); + Grid_FFT 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(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"); +} \ No newline at end of file