mirror of
https://github.com/cosmo-sims/monofonIC.git
synced 2024-09-19 17:03:45 +02:00
Added convergence test and made testings activable from the config file.
This commit is contained in:
Michael Michaux
parent
8bc6e5acac
commit
b064928763
4 changed files with 148 additions and 7 deletions
|
|
@ -8,6 +8,9 @@ DoBaryons = no
|
|||
DoFixing = no
|
||||
BCClattice = no
|
||||
|
||||
[testing]
|
||||
test = convergence
|
||||
|
||||
[execution]
|
||||
NumThreads = 4
|
||||
|
||||
|
|
|
|||
|
|
@ -24,4 +24,13 @@ namespace testing{
|
|||
Grid_FFT<real_t> &phi3b,
|
||||
std::array<Grid_FFT<real_t> *, 3> &A3,
|
||||
bool bwrite_out_fields=false);
|
||||
|
||||
void output_convergence(
|
||||
ConfigFile &the_config,
|
||||
std::size_t ngrid, real_t boxlen, real_t vfac, real_t dplus,
|
||||
Grid_FFT<real_t> &phi,
|
||||
Grid_FFT<real_t> &phi2,
|
||||
Grid_FFT<real_t> &phi3a,
|
||||
Grid_FFT<real_t> &phi3b,
|
||||
std::array<Grid_FFT<real_t> *, 3> &A3);
|
||||
}
|
||||
|
|
|
|||
|
|
@ -323,15 +323,23 @@ int Run( ConfigFile& the_config )
|
|||
///////////////////////////////////////////////////////////////////////
|
||||
// we store the densities here if we compute them
|
||||
//======================================================================
|
||||
const bool testing_compute_densities = false;
|
||||
if( testing_compute_densities )
|
||||
{
|
||||
|
||||
// Testing
|
||||
const std::string testing = the_config.GetValueSafe<std::string>("testing", "test", "none");
|
||||
|
||||
if(testing != "none") {
|
||||
csoca::wlog << "you are running in testing mode. No ICs, only diagnostic output will be written out!" << std::endl;
|
||||
// testing::output_potentials_and_densities( the_config, ngrid, boxlen, phi, phi2, phi3a, phi3b, A3 );
|
||||
testing::output_velocity_displacement_symmetries( the_config, ngrid, boxlen, vfac, Dplus0, phi, phi2, phi3a, phi3b, A3 );
|
||||
if(testing == "potentials_and_densities") {
|
||||
testing::output_potentials_and_densities(the_config, ngrid, boxlen, phi, phi2, phi3a, phi3b, A3);
|
||||
} 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);
|
||||
} else {
|
||||
csoca::flog << "unknown test '" << testing << "'" << std::endl;
|
||||
std::abort();
|
||||
}
|
||||
else
|
||||
{
|
||||
} else {
|
||||
// temporary storage of data
|
||||
Grid_FFT<real_t> tmp({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen});
|
||||
|
||||
|
|
|
|||
121
src/testing.cc
121
src/testing.cc
|
|
@ -1,5 +1,6 @@
|
|||
#include <testing.hh>
|
||||
#include <unistd.h> // for unlink
|
||||
#include <memory>
|
||||
|
||||
#include <operators.hh>
|
||||
#include <convolution.hh>
|
||||
|
|
@ -239,4 +240,124 @@ void output_velocity_displacement_symmetries(
|
|||
|
||||
}
|
||||
|
||||
void output_convergence(
|
||||
ConfigFile &the_config,
|
||||
std::size_t ngrid, real_t boxlen, real_t vfac, real_t dplus,
|
||||
Grid_FFT<real_t> &phi,
|
||||
Grid_FFT<real_t> &phi2,
|
||||
Grid_FFT<real_t> &phi3a,
|
||||
Grid_FFT<real_t> &phi3b,
|
||||
std::array<Grid_FFT<real_t> *, 3> &A3)
|
||||
{
|
||||
|
||||
// scale all potentials to remove dplus0
|
||||
phi /= dplus;
|
||||
phi2 /= dplus * dplus;
|
||||
phi3a /= dplus * dplus * dplus;
|
||||
phi3b /= dplus * dplus * dplus;
|
||||
(*A3[0]) /= dplus * dplus * dplus;
|
||||
(*A3[1]) /= dplus * dplus * dplus;
|
||||
(*A3[2]) /= dplus * dplus * dplus;
|
||||
|
||||
// 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});
|
||||
Grid_FFT<real_t> psi_3({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen});
|
||||
#pragma omp parallel for collapse(3)
|
||||
for (std::size_t i = 0; i < psi_1.size(0); ++i) {
|
||||
for (std::size_t j = 0; j < psi_1.size(1); ++j) {
|
||||
for (std::size_t k = 0; k < psi_1.size(2); ++k) {
|
||||
std::size_t idx = psi_1.get_idx(i, j, k);
|
||||
psi_1.relem(idx) = 0.0;
|
||||
psi_2.relem(idx) = 0.0;
|
||||
psi_3.relem(idx) = 0.0;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// temporaries
|
||||
Grid_FFT<real_t> psi_1_tmp({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen});
|
||||
Grid_FFT<real_t> psi_2_tmp({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen});
|
||||
Grid_FFT<real_t> psi_3_tmp({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen});
|
||||
|
||||
// compute psi 1, 2 and 3
|
||||
for (int idim = 0; idim < 3; ++idim) {
|
||||
// cyclic rotations of indices
|
||||
int idimp = (idim + 1) % 3, idimpp = (idim + 2) % 3;
|
||||
|
||||
psi_1_tmp.FourierTransformForward(false);
|
||||
psi_2_tmp.FourierTransformForward(false);
|
||||
psi_3_tmp.FourierTransformForward(false);
|
||||
|
||||
#pragma omp parallel for collapse(3)
|
||||
for (std::size_t i = 0; i < phi.size(0); ++i) {
|
||||
for (std::size_t j = 0; j < phi.size(1); ++j) {
|
||||
for (std::size_t k = 0; k < phi.size(2); ++k) {
|
||||
auto kk = phi.get_k<real_t>(i, j, k);
|
||||
std::size_t idx = phi.get_idx(i, j, k);
|
||||
psi_1_tmp.kelem(idx) = ccomplex_t(0.0, 1.0) * (kk[idim] * phi.kelem(idx));
|
||||
psi_2_tmp.kelem(idx) = ccomplex_t(0.0, 1.0) * (kk[idim] * phi2.kelem(idx));
|
||||
psi_3_tmp.kelem(idx) = ccomplex_t(0.0, 1.0) * (
|
||||
kk[idim] * (phi3a.kelem(idx) + phi3b.kelem(idx)) +
|
||||
kk[idimp] * A3[idimpp]->kelem(idx) -
|
||||
kk[idimpp] * A3[idimp]->kelem(idx)
|
||||
);
|
||||
}
|
||||
}
|
||||
}
|
||||
psi_1_tmp.FourierTransformBackward();
|
||||
psi_2_tmp.FourierTransformBackward();
|
||||
psi_3_tmp.FourierTransformBackward();
|
||||
|
||||
// sum of squares
|
||||
#pragma omp parallel for collapse(3)
|
||||
for (std::size_t i = 0; i < psi_1.size(0); ++i) {
|
||||
for (std::size_t j = 0; j < psi_1.size(1); ++j) {
|
||||
for (std::size_t k = 0; k < psi_1.size(2); ++k) {
|
||||
std::size_t idx = psi_1.get_idx(i, j, k);
|
||||
psi_1.relem(idx) += psi_1_tmp.relem(idx) * psi_1_tmp.relem(idx);
|
||||
psi_2.relem(idx) += psi_2_tmp.relem(idx) * psi_2_tmp.relem(idx);
|
||||
psi_3.relem(idx) += psi_3_tmp.relem(idx) * psi_3_tmp.relem(idx);
|
||||
}
|
||||
}
|
||||
}
|
||||
} // loop on dimensions
|
||||
|
||||
// apply square root for the L2 norm
|
||||
#pragma omp parallel for collapse(3)
|
||||
for (std::size_t i = 0; i < psi_1.size(0); ++i) {
|
||||
for (std::size_t j = 0; j < psi_1.size(1); ++j) {
|
||||
for (std::size_t k = 0; k < psi_1.size(2); ++k) {
|
||||
std::size_t idx = psi_1.get_idx(i, j, k);
|
||||
psi_1.relem(idx) = std::sqrt(psi_1.relem(idx));
|
||||
psi_2.relem(idx) = std::sqrt(psi_2.relem(idx));
|
||||
psi_3.relem(idx) = std::sqrt(psi_3.relem(idx));
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// convergence radius
|
||||
Grid_FFT<real_t> inv_convergence_radius({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen});
|
||||
#pragma omp parallel for collapse(3)
|
||||
for (std::size_t i = 0; i < psi_1.size(0); ++i) {
|
||||
for (std::size_t j = 0; j < psi_1.size(1); ++j) {
|
||||
for (std::size_t k = 0; k < psi_1.size(2); ++k) {
|
||||
std::size_t idx = psi_1.get_idx(i, j, k);
|
||||
inv_convergence_radius.relem(idx) =
|
||||
3.0 * (std::abs(psi_3.relem(idx)) / std::abs(psi_2.relem(idx))) -
|
||||
2.0 * (std::abs(psi_2.relem(idx)) / std::abs(psi_1.relem(idx)));
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// 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");
|
||||
|
||||
}
|
||||
|
||||
} // namespace testing
|
||||
|
|
|
|||
Loading…
Reference in a new issue