mirror of
https://github.com/cosmo-sims/monofonIC.git
synced 2024-09-19 17:03:45 +02:00
Added theoretical convergence test.
This commit is contained in:
parent
e9b11b1bae
commit
01b22b76a3
3 changed files with 99 additions and 9 deletions
|
@ -4,6 +4,7 @@
|
|||
#include <general.hh>
|
||||
#include <config_file.hh>
|
||||
#include <grid_fft.hh>
|
||||
#include <cosmology_calculator.hh>
|
||||
|
||||
namespace testing{
|
||||
void output_potentials_and_densities(
|
||||
|
@ -27,6 +28,7 @@ namespace testing{
|
|||
|
||||
void output_convergence(
|
||||
ConfigFile &the_config,
|
||||
CosmologyCalculator* the_cosmo_calc,
|
||||
std::size_t ngrid, real_t boxlen, real_t vfac, real_t dplus,
|
||||
Grid_FFT<real_t> &phi,
|
||||
Grid_FFT<real_t> &phi2,
|
||||
|
|
|
@ -335,7 +335,7 @@ int Run( ConfigFile& the_config )
|
|||
} else if(testing == "velocity_displacement_symmetries") {
|
||||
testing::output_velocity_displacement_symmetries(the_config, ngrid, boxlen, vfac, Dplus0, phi, phi2, phi3a, phi3b, A3);
|
||||
} else if(testing == "convergence") {
|
||||
testing::output_convergence(the_config, ngrid, boxlen, vfac, Dplus0, phi, phi2, phi3a, phi3b, A3);
|
||||
testing::output_convergence(the_config, the_cosmo_calc.get(), ngrid, boxlen, vfac, Dplus0, phi, phi2, phi3a, phi3b, A3);
|
||||
} else {
|
||||
csoca::flog << "unknown test '" << testing << "'" << std::endl;
|
||||
std::abort();
|
||||
|
|
104
src/testing.cc
104
src/testing.cc
|
@ -242,6 +242,7 @@ void output_velocity_displacement_symmetries(
|
|||
|
||||
void output_convergence(
|
||||
ConfigFile &the_config,
|
||||
CosmologyCalculator* the_cosmo_calc,
|
||||
std::size_t ngrid, real_t boxlen, real_t vfac, real_t dplus,
|
||||
Grid_FFT<real_t> &phi,
|
||||
Grid_FFT<real_t> &phi2,
|
||||
|
@ -249,7 +250,6 @@ void output_convergence(
|
|||
Grid_FFT<real_t> &phi3b,
|
||||
std::array<Grid_FFT<real_t> *, 3> &A3)
|
||||
{
|
||||
|
||||
// scale all potentials to remove dplus0
|
||||
phi /= dplus;
|
||||
phi2 /= dplus * dplus;
|
||||
|
@ -259,6 +259,90 @@ void output_convergence(
|
|||
(*A3[1]) /= dplus * dplus * dplus;
|
||||
(*A3[2]) /= dplus * dplus * dplus;
|
||||
|
||||
////////////////////// theoretical convergence radius //////////////////////
|
||||
|
||||
// compute phi_code
|
||||
Grid_FFT<real_t> phi_code({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen});
|
||||
phi_code.FourierTransformForward(false);
|
||||
#pragma omp parallel for collapse(3)
|
||||
for (std::size_t i = 0; i < phi_code.size(0); ++i) {
|
||||
for (std::size_t j = 0; j < phi_code.size(1); ++j) {
|
||||
for (std::size_t k = 0; k < phi_code.size(2); ++k) {
|
||||
std::size_t idx = phi_code.get_idx(i, j, k);
|
||||
phi_code.kelem(idx) = -phi.kelem(idx);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// initialize norm to 0
|
||||
Grid_FFT<real_t> nabla_vini_norm({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen});
|
||||
#pragma omp parallel for collapse(3)
|
||||
for (std::size_t i = 0; i < nabla_vini_norm.size(0); ++i) {
|
||||
for (std::size_t j = 0; j < nabla_vini_norm.size(1); ++j) {
|
||||
for (std::size_t k = 0; k < nabla_vini_norm.size(2); ++k) {
|
||||
std::size_t idx = nabla_vini_norm.get_idx(i, j, k);
|
||||
nabla_vini_norm.relem(idx) = 0.0;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
Grid_FFT<real_t> nabla_vini_mn({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen});
|
||||
for(std::size_t m = 0; m < 3; m++) {
|
||||
for(std::size_t n = m; n < 3; n++) {
|
||||
nabla_vini_mn.FourierTransformForward(false);
|
||||
#pragma omp parallel for collapse(3)
|
||||
for (std::size_t i = 0; i < phi_code.size(0); ++i) {
|
||||
for (std::size_t j = 0; j < phi_code.size(1); ++j) {
|
||||
for (std::size_t k = 0; k < phi_code.size(2); ++k) {
|
||||
std::size_t idx = phi_code.get_idx(i, j, k);
|
||||
auto kk = phi_code.get_k<real_t>(i, j, k);
|
||||
nabla_vini_mn.kelem(idx) = phi_code.kelem(idx) * (kk[m] * kk[n]);
|
||||
}
|
||||
}
|
||||
}
|
||||
nabla_vini_mn.FourierTransformBackward();
|
||||
nabla_vini_mn *= (3.2144004915 / the_cosmo_calc->CalcGrowthFactor(1.0));
|
||||
// sum of squares
|
||||
#pragma omp parallel for collapse(3)
|
||||
for (std::size_t i = 0; i < nabla_vini_norm.size(0); ++i) {
|
||||
for (std::size_t j = 0; j < nabla_vini_norm.size(1); ++j) {
|
||||
for (std::size_t k = 0; k < nabla_vini_norm.size(2); ++k) {
|
||||
std::size_t idx = nabla_vini_norm.get_idx(i, j, k);
|
||||
if(m != n) {
|
||||
nabla_vini_norm.relem(idx) += (2.0 * nabla_vini_mn.relem(idx) * nabla_vini_mn.relem(idx));
|
||||
} else {
|
||||
nabla_vini_norm.relem(idx) += (nabla_vini_mn.relem(idx) * nabla_vini_mn.relem(idx));
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
// square root
|
||||
#pragma omp parallel for collapse(3)
|
||||
for (std::size_t i = 0; i < nabla_vini_norm.size(0); ++i) {
|
||||
for (std::size_t j = 0; j < nabla_vini_norm.size(1); ++j) {
|
||||
for (std::size_t k = 0; k < nabla_vini_norm.size(2); ++k) {
|
||||
std::size_t idx = nabla_vini_norm.get_idx(i, j, k);
|
||||
nabla_vini_norm.relem(idx) = std::sqrt(nabla_vini_norm.relem(idx));
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// get t_eds
|
||||
Grid_FFT<real_t> t_eds({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen});
|
||||
#pragma omp parallel for collapse(3)
|
||||
for (std::size_t i = 0; i < t_eds.size(0); ++i) {
|
||||
for (std::size_t j = 0; j < t_eds.size(1); ++j) {
|
||||
for (std::size_t k = 0; k < t_eds.size(2); ++k) {
|
||||
std::size_t idx = t_eds.get_idx(i, j, k);
|
||||
t_eds.relem(idx) = 0.0204 / nabla_vini_norm.relem(idx);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
////////////////////////// 3lpt convergence test ///////////////////////////
|
||||
|
||||
// initialize grids to 0
|
||||
Grid_FFT<real_t> psi_1({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen});
|
||||
Grid_FFT<real_t> psi_2({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen});
|
||||
|
@ -351,13 +435,17 @@ void output_convergence(
|
|||
}
|
||||
}
|
||||
|
||||
// write results
|
||||
unlink("convergence_test.hdf5");
|
||||
inv_convergence_radius.Write_to_HDF5("convergence_test.hdf5", "inv_convergence_radius");
|
||||
psi_1.Write_to_HDF5("convergence_test.hdf5", "psi_1_norm");
|
||||
psi_2.Write_to_HDF5("convergence_test.hdf5", "psi_2_norm");
|
||||
psi_3.Write_to_HDF5("convergence_test.hdf5", "psi_3_norm");
|
||||
|
||||
////////////////////////////// write results ///////////////////////////////
|
||||
std::string convergence_test_filename("convergence_test.hdf5");
|
||||
unlink(convergence_test_filename.c_str());
|
||||
#if defined(USE_MPI)
|
||||
MPI_Barrier(MPI_COMM_WORLD);
|
||||
#endif
|
||||
t_eds.Write_to_HDF5(convergence_test_filename, "t_eds");
|
||||
inv_convergence_radius.Write_to_HDF5(convergence_test_filename, "inv_convergence_radius");
|
||||
// psi_1.Write_to_HDF5(convergence_test_filename, "psi_1_norm");
|
||||
// psi_2.Write_to_HDF5(convergence_test_filename, "psi_2_norm");
|
||||
// psi_3.Write_to_HDF5(convergence_test_filename, "psi_3_norm");
|
||||
}
|
||||
|
||||
} // namespace testing
|
||||
|
|
Loading…
Reference in a new issue