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

fixed mixup between time coeffs and 3a 3b terms

Browse source
This commit is contained in:
Oliver Hahn 2019-05-20 17:23:21 +02:00
parent 2b704c3a1f
commit 28e072dea3
2 changed files with 27 additions and 46 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

2
src/grid_fft.cc
View file

@ -564,9 +564,7 @@ void Grid_FFT<data_t>::Write_PowerSpectrum( std::string ofname )
{
#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";

71
src/main.cc
View file

@ -147,6 +147,15 @@ int main( int argc, char** argv )
OrszagConvolver<real_t> Conv({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen});
//--------------------------------------------------------------------
// Some operators to add or subtract terms
auto assign_op = []( ccomplex_t res, ccomplex_t val ) -> ccomplex_t{ return res; };
auto add_op = []( ccomplex_t res, ccomplex_t val ) -> ccomplex_t{ return val+res; };
auto add2_op = []( ccomplex_t res, ccomplex_t val ) -> ccomplex_t{ return val+2.0*res; };
auto sub_op = []( ccomplex_t res, ccomplex_t val ) -> ccomplex_t{ return val-res; };
auto sub2_op = []( ccomplex_t res, ccomplex_t val ) -> ccomplex_t{ return val-2.0*res; };
//--------------------------------------------------------------------
//phi.FillRandomReal(6519);
the_random_number_generator->Fill_Grid( phi );
@ -164,19 +173,15 @@ int main( int argc, char** argv )
ccomplex_t delta = x * the_cosmo_calc->GetAmplitude(kmod, total);
return -delta / (kmod * kmod) * phifac / volfac;
});
phi.zero_DC_mode();
csoca::ilog << " took " << get_wtime()-wtime << "s" << std::endl;
//////////////////////////////////////////////////////////////////////////////////////////////////
auto assign_op = []( ccomplex_t res, ccomplex_t val ) -> ccomplex_t{ return res; };
auto add_op = []( ccomplex_t res, ccomplex_t val ) -> ccomplex_t{ return val+res; };
auto add2_op = []( ccomplex_t res, ccomplex_t val ) -> ccomplex_t{ return val+2.0*res; };
auto sub_op = []( ccomplex_t res, ccomplex_t val ) -> ccomplex_t{ return val-res; };
auto sub2_op = []( ccomplex_t res, ccomplex_t val ) -> ccomplex_t{ return val-2.0*res; };
//======================================================================
//... compute 2LPT displacement potential ....
wtime = get_wtime();
csoca::ilog << "Computing phi(2) term..." << std::flush;
// Compute the source term for phi(2)
Conv.convolve_SumHessians( phi, {0,0}, phi, {1,1}, {2,2}, phi2, assign_op );
Conv.convolve_Hessians( phi, {1,1}, phi, {2,2}, phi2, add_op );
Conv.convolve_Hessians( phi, {0,1}, phi, {0,1}, phi2, sub_op );
@ -186,57 +191,35 @@ int main( int argc, char** argv )
phi2 /= phifac;
csoca::ilog << " took " << get_wtime()-wtime << "s" << std::endl;
// phi2.FourierTransformForward();
// phi2.apply_function_k_dep([&](auto x, auto k) {
// real_t kmod2 = k.norm_squared();
// return x * (-1.0 / kmod2) * phifac / phifac / phifac;
// });
// phi2.zero_DC_mode();
//======================================================================
//... compute 3LPT displacement potential
wtime = get_wtime();
csoca::ilog << "Computing phi(3a) term..." << std::flush;
Conv.convolve_SumHessians( phi, {0,0}, phi2, {1,1}, {2,2}, phi3a, assign_op );
Conv.convolve_SumHessians( phi, {1,1}, phi2, {2,2}, {0,0}, phi3a, add_op );
Conv.convolve_SumHessians( phi, {2,2}, phi2, {0,0}, {1,1}, phi3a, add_op );
Conv.convolve_Hessians( phi, {0,1}, phi2, {0,1}, phi3a, sub2_op );
Conv.convolve_Hessians( phi, {0,2}, phi2, {0,2}, phi3a, sub2_op );
Conv.convolve_Hessians( phi, {1,2}, phi2, {1,2}, phi3a, sub2_op );
Conv.convolve_Hessians( phi, {0,0}, phi, {1,1}, phi, {2,2}, phi3a, assign_op );
Conv.convolve_Hessians( phi, {0,1}, phi, {0,2}, phi, {1,2}, phi3a, add2_op );
Conv.convolve_Hessians( phi, {1,2}, phi, {1,2}, phi, {0,0}, phi3a, sub_op );
Conv.convolve_Hessians( phi, {0,2}, phi, {0,2}, phi, {1,1}, phi3a, sub_op );
Conv.convolve_Hessians( phi, {0,1}, phi, {0,1}, phi, {2,2}, phi3a, sub_op );
phi3a.apply_InverseLaplacian();
phi3a *= 0.5/phifac; // factor 1/2 from definition of phi(3a)!
// phi3a.apply_function_k_dep([&](auto x, auto k) {
// return 0.5 * x;
// });
phi3a /= phifac*phifac;
csoca::ilog << " took " << get_wtime()-wtime << "s" << std::endl;
// phi3a.FourierTransformForward();
// phi3a.apply_function_k_dep([&](auto x, auto k) {
// real_t kmod2 = k.norm_squared();
// return x * (-1.0 / kmod2) * phifac / phifac / phifac;
// });
// phi3a.zero_DC_mode();
//...
wtime = get_wtime();
csoca::ilog << "Computing phi(3b) term..." << std::flush;
Conv.convolve_Hessians( phi, {0,0}, phi, {1,1}, phi, {2,2}, phi3b, assign_op );
Conv.convolve_Hessians( phi, {0,1}, phi, {0,2}, phi, {1,2}, phi3b, add2_op );
Conv.convolve_Hessians( phi, {1,2}, phi, {1,2}, phi, {0,0}, phi3b, sub_op );
Conv.convolve_Hessians( phi, {0,2}, phi, {0,2}, phi, {1,1}, phi3b, sub_op );
Conv.convolve_Hessians( phi, {0,1}, phi, {0,1}, phi, {2,2}, phi3b, sub_op );
// phi3b.FourierTransformForward();
// phi3b.apply_function_k_dep([&](auto x, auto k) {
// real_t kmod2 = k.norm_squared();
// return x * (-1.0 / kmod2) * phifac / phifac / phifac/phifac;
// });
// phi3b.zero_DC_mode();
Conv.convolve_SumHessians( phi, {0,0}, phi2, {1,1}, {2,2}, phi3b, assign_op );
Conv.convolve_SumHessians( phi, {1,1}, phi2, {2,2}, {0,0}, phi3b, add_op );
Conv.convolve_SumHessians( phi, {2,2}, phi2, {0,0}, {1,1}, phi3b, add_op );
Conv.convolve_Hessians( phi, {0,1}, phi2, {0,1}, phi3b, sub2_op );
Conv.convolve_Hessians( phi, {0,2}, phi2, {0,2}, phi3b, sub2_op );
Conv.convolve_Hessians( phi, {1,2}, phi2, {1,2}, phi3b, sub2_op );
phi3b.apply_InverseLaplacian();
phi3b /= phifac*phifac;
phi3b *= 0.5/phifac; // factor 1/2 from definition of phi(3a)!
csoca::ilog << " took " << get_wtime()-wtime << "s" << std::endl;
///////////////////////////////////////////////////////////////////////
// we store the densities here if we compute them
const bool compute_densities = true;