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

wip. fixed some bugs. appears to work now for arbitrary margins. needs more testing

Browse source
This commit is contained in:
Oliver Hahn 2023-02-05 12:27:55 -08:00
parent 1e8eeef76a
commit 6f05ca2e7b
3 changed files with 202 additions and 25 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

25
src/densities.cc
View file

@ -173,6 +173,8 @@ void fft_interpolate(m1 &V, m2 &v, bool from_basegrid = false)
size_t nxf = v.size(0), nyf = v.size(1), nzf = v.size(2), nzfp = nzf + 2;
size_t nxF = V.size(0), nyF = V.size(1), nzF = V.size(2);
size_t mxf = v.margin(0), myf = v.margin(1), mzf = v.margin(2);
if (!from_basegrid)
{
#ifdef NO_COARSE_OVERLAP
@ -180,19 +182,20 @@ void fft_interpolate(m1 &V, m2 &v, bool from_basegrid = false)
oyf += nyF / 4;
ozf += nzF / 4;
#else
oxf += nxF / 4 - nxf / 8;
oyf += nyF / 4 - nyf / 8;
ozf += nzF / 4 - nzf / 8;
oxf += nxF / 4 - nxf / 8; //mxf / 2;
oyf += nyF / 4 - nyf / 8; //myf / 2;
ozf += nzF / 4 - nzf / 8; //mzf / 2;
}
else
{
oxf -= nxf / 8;
oyf -= nyf / 8;
ozf -= nzf / 8;
oxf -= mxf/2; //nxf / 8;
oyf -= myf/2; //nyf / 8;
ozf -= mzf/2; //nzf / 8;
#endif
}
LOGUSER("FFT interpolate: offset=%d,%d,%d size=%d,%d,%d", oxf, oyf, ozf, nxf, nyf, nzf);
// LOGUSER("FFT interpolate: offset=%d,%d,%d size=%d,%d,%d", oxf, oyf, ozf, nxf, nyf, nzf);
LOGINFO("FFT interpolate: offset=%d,%d,%d size=%d,%d,%d", oxf, oyf, ozf, nxf, nyf, nzf);
// cut out piece of coarse grid that overlaps the fine:
assert(nxf % 2 == 0 && nyf % 2 == 0 && nzf % 2 == 0);
@ -632,7 +635,7 @@ void GenerateDensityHierarchy(config_file &cf, transfer_function *ptf, tf_type t
LOGUSER(" size =(%5d,%5d,%5d)", refh.size(levelmin + i, 0),
refh.size(levelmin + i, 1), refh.size(levelmin + i, 2));
if( refh.get_margin() > 0 )
if( refh.get_margin() > 0 ){
fine = new PaddedDensitySubGrid<real_t>(refh.offset(levelmin + i, 0),
refh.offset(levelmin + i, 1),
refh.offset(levelmin + i, 2),
@ -640,14 +643,16 @@ void GenerateDensityHierarchy(config_file &cf, transfer_function *ptf, tf_type t
refh.size(levelmin + i, 1),
refh.size(levelmin + i, 2),
refh.get_margin(), refh.get_margin(), refh.get_margin() );
else
LOGUSER(" margin = %d",refh.get_margin());
}else{
fine = new PaddedDensitySubGrid<real_t>(refh.offset(levelmin + i, 0),
refh.offset(levelmin + i, 1),
refh.offset(levelmin + i, 2),
refh.size(levelmin + i, 0),
refh.size(levelmin + i, 1),
refh.size(levelmin + i, 2));
LOGUSER(" margin = %d",refh.size(levelmin + i, 0)/2);
}
/////////////////////////////////////////////////////////////////////////
// load white noise for patch

35
src/densities.hh
View file

@ -245,7 +245,7 @@ public:
using DensityGrid<real_t>::oz_;
using DensityGrid<real_t>::data_;
size_t padx_, pady_, padz_;
std::array<size_t,3> pad_;
using DensityGrid<real_t>::fill_rand;
using DensityGrid<real_t>::get_data_ptr;
@ -253,26 +253,31 @@ public:
public:
PaddedDensitySubGrid(int ox, int oy, int oz, unsigned nx, unsigned ny, unsigned nz )
: DensityGrid<real_t>(nx*2, ny*2, nz*2, ox, oy, oz),
padx_(nx / 2), pady_(ny / 2), padz_(nz / 2)
pad_{{ nx / 2, ny / 2, nz / 2 }}
{ }
PaddedDensitySubGrid(int ox, int oy, int oz, unsigned nx, unsigned ny, unsigned nz, unsigned padx, unsigned pady, unsigned padz )
: DensityGrid<real_t>(nx + 2 * padx, ny + 2 * pady, nz + 2 * padz, ox, oy, oz),
padx_(padx), pady_(pady), padz_(padz)
pad_{{ padx, pady, padz }}
{ }
PaddedDensitySubGrid(const PaddedDensitySubGrid<real_t> &o)
: DensityGrid<real_t>(o)
{ }
size_t margin(int i) const
{
return pad_[i];
}
template <class array3>
void copy_unpad(array3 &v)
{
#pragma omp parallel for
for (size_t ix = padx_; ix < nx_-padx_; ++ix){
const size_t ixu = ix - padx_;
for (size_t iy = pady_, iyu = 0; iy < ny_-pady_; ++iy, ++iyu)
for (size_t iz = padz_, izu = 0; iz < nz_-padz_; ++iz, ++izu)
for (size_t ix = pad_[0]; ix < nx_-pad_[0]; ++ix){
const size_t ixu = ix - pad_[0];
for (size_t iy = pad_[1], iyu = 0; iy < ny_-pad_[1]; ++iy, ++iyu)
for (size_t iz = pad_[2], izu = 0; iz < nz_-pad_[2]; ++iz, ++izu)
v(ixu, iyu, izu) = (*this)(ix, iy, iz);
}
}
@ -281,10 +286,10 @@ public:
void copy_add_unpad(array3 &v)
{
#pragma omp parallel for
for (size_t ix = padx_; ix < nx_-padx_; ++ix){
const size_t ixu = ix - padx_;
for (size_t iy = pady_, iyu = 0; iy < ny_-pady_; ++iy, ++iyu)
for (size_t iz = padz_, izu = 0; iz < nz_-padz_; ++iz, ++izu)
for (size_t ix = pad_[0]; ix < nx_-pad_[0]; ++ix){
const size_t ixu = ix - pad_[0];
for (size_t iy = pad_[1], iyu = 0; iy < ny_-pad_[1]; ++iy, ++iyu)
for (size_t iz = pad_[2], izu = 0; iz < nz_-pad_[2]; ++iz, ++izu)
v(ixu, iyu, izu) += (*this)(ix, iy, iz);
}
}
@ -293,10 +298,10 @@ public:
void copy_subtract_unpad(array3 &v)
{
#pragma omp parallel for
for (size_t ix = padx_; ix < nx_-padx_; ++ix){
const size_t ixu = ix - padx_;
for (size_t iy = pady_, iyu = 0; iy < ny_-pady_; ++iy, ++iyu)
for (size_t iz = padz_, izu = 0; iz < nz_-padz_; ++iz, ++izu)
for (size_t ix = pad_[0]; ix < nx_-pad_[0]; ++ix){
const size_t ixu = ix - pad_[0];
for (size_t iy = pad_[1], iyu = 0; iy < ny_-pad_[1]; ++iy, ++iyu)
for (size_t iz = pad_[2], izu = 0; iz < nz_-pad_[2]; ++iz, ++izu)
v(ixu, iyu, izu) -= (*this)(ix, iy, iz);
}
}

167
tools/check_output.ipynb Normal file
View file

File diff suppressed because one or more lines are too long