mirror of
https://github.com/cosmo-sims/MUSIC.git
synced 2024-09-19 17:03:46 +02:00
2383f8b68d
Fixed a bug with top grid kernel normalization, this reduces errors a bit. A non-cubic grid is now extended to a cubic grid in the density calculation. This should reduce artefacts due to TF truncation. Added a new (undocumented and internal) parameter that forces grid centering on the specified refinement region also in unigrid mode. This should in the future also obey that the shift comes in multiples of the top grid cells (which are however unknown). Fixed a bug in the generic output plugin that caused non-cubic grids to be cut to cubic grids.
831 lines
23 KiB
C++
831 lines
23 KiB
C++
/*
|
|
|
|
densities.hh - This file is part of MUSIC -
|
|
a code to generate multi-scale initial conditions
|
|
for cosmological simulations
|
|
|
|
Copyright (C) 2010 Oliver Hahn
|
|
|
|
This program is free software: you can redistribute it and/or modify
|
|
it under the terms of the GNU General Public License as published by
|
|
the Free Software Foundation, either version 3 of the License, or
|
|
(at your option) any later version.
|
|
|
|
This program is distributed in the hope that it will be useful,
|
|
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|
GNU General Public License for more details.
|
|
|
|
You should have received a copy of the GNU General Public License
|
|
along with this program. If not, see <http://www.gnu.org/licenses/>.
|
|
|
|
*/
|
|
|
|
#ifndef __DENSITIES_HH
|
|
#define __DENSITIES_HH
|
|
|
|
|
|
#ifdef SINGLE_PRECISION
|
|
#ifdef SINGLETHREAD_FFTW
|
|
#include <srfftw.h>
|
|
#else
|
|
#include <srfftw_threads.h>
|
|
#endif
|
|
#else
|
|
#ifdef SINGLETHREAD_FFTW
|
|
#include <drfftw.h>
|
|
#else
|
|
#include <drfftw_threads.h>
|
|
#endif
|
|
#endif
|
|
|
|
#include "config_file.hh"
|
|
#include "random.hh"
|
|
#include "cosmology.hh"
|
|
#include "transfer_function.hh"
|
|
|
|
|
|
void GenerateDensityHierarchy( config_file& cf, transfer_function *ptf, tf_type type,
|
|
refinement_hierarchy& refh, grid_hierarchy& delta );
|
|
|
|
void GenerateDensityUnigrid( config_file& cf, transfer_function *ptf, tf_type type,
|
|
refinement_hierarchy& refh, grid_hierarchy& delta, bool kspace=false );
|
|
|
|
void normalize_density( grid_hierarchy& delta );
|
|
|
|
|
|
/*!
|
|
* @class DensityGrid
|
|
* @brief provides infrastructure for computing the initial density field
|
|
*
|
|
* This class provides access and data management member functions that
|
|
* are used when computing the initial density field by convolution with
|
|
* transfer functions.
|
|
*/
|
|
template< typename real_t >
|
|
class DensityGrid
|
|
{
|
|
public:
|
|
|
|
int nx_; //!< number of grid cells in x-direction
|
|
int ny_; //!< number of grid cells in y-direction
|
|
int nz_; //!< number of grid cells in z-direction
|
|
int nzp_; //!< number of cells in memory (z-dir), used for Nyquist padding
|
|
|
|
//! the actual data container in the form of a 1D array
|
|
std::vector< real_t > data_;
|
|
|
|
//! constructor
|
|
/*! constructs an instance given the dimensions of the density field
|
|
* @params nx the number of cells in x
|
|
* @params ny the number of cells in y
|
|
* @params nz the number of cells in z
|
|
*/
|
|
DensityGrid( unsigned nx, unsigned ny, unsigned nz )
|
|
: nx_(nx), ny_(ny), nz_(nz)
|
|
{
|
|
data_.assign(nx_*ny_*2*(nz_/2+1),0.0);
|
|
nzp_ = 2*(nz_/2+1);
|
|
}
|
|
|
|
//! copy constructor
|
|
explicit DensityGrid( const DensityGrid<real_t> & g )
|
|
: nx_(g.nx_), ny_(g.ny_), nz_(g.nz_), nzp_(g.nzp_)
|
|
{
|
|
data_ = g.data_;
|
|
}
|
|
|
|
//!destructor
|
|
~DensityGrid()
|
|
{ }
|
|
|
|
//! clears the density object
|
|
/*! sets all dimensions to zero and frees the memory
|
|
*/
|
|
void clear( void )
|
|
{
|
|
nx_ = ny_ = nz_ = nzp_ = 0;
|
|
data_.clear();
|
|
std::vector<real_t>().swap(data_);
|
|
}
|
|
|
|
//! assigns the contents of another DensityGrid to this
|
|
DensityGrid& operator=( const DensityGrid<real_t>& g )
|
|
{
|
|
nx_ = g.nx_;
|
|
ny_ = g.ny_;
|
|
nz_ = g.nz_;
|
|
nzp_= g.nzp_;
|
|
data_ = g.data_;
|
|
|
|
return *this;
|
|
}
|
|
|
|
//! 3D index based data access operator
|
|
inline real_t& operator()( int i, int j, int k )
|
|
{ return data_[(i*ny_+j)*nzp_+k]; }
|
|
|
|
//! 3D index based const data access operator
|
|
inline const real_t& operator()( int i, int j, int k ) const
|
|
{ return data_[(i*ny_+j)*nzp_+k]; }
|
|
|
|
//! recover the pointer to the 1D data array
|
|
inline real_t * get_data_ptr( void )
|
|
{ return &data_[0]; }
|
|
|
|
|
|
//! fills the density field with random number values
|
|
/*! given a pointer to a random_numbers object, fills the field with random values
|
|
* @params prc pointer to a random_numbers object
|
|
* @params variance the variance of the random numbers (the values returned by prc are multiplied by this)
|
|
* @params i0 x-offset (shift) in cells of the density field with respect to the random number field
|
|
* @params j0 y-offset (shift) in cells of the density field with respect to the random number field
|
|
* @params k0 z-offset (shift) in cells of the density field with respect to the random number field
|
|
* @params zero boolean, if true, the global mean will be subtracted
|
|
*/
|
|
void fill_rand( /*const*/ random_numbers<real_t>* prc, real_t variance, int i0, int j0, int k0, bool zero=false )
|
|
{
|
|
double sum = 0.0;
|
|
|
|
#pragma omp parallel for reduction(+:sum)
|
|
for( int i=0; i<nx_; ++i )
|
|
for( int j=0; j<ny_; ++j )
|
|
for( int k=0; k<nz_; ++k )
|
|
{
|
|
(*this)(i,j,k) = (*prc)(i0+i,j0+j,k0+k) * variance;
|
|
sum += (*this)(i,j,k);
|
|
}
|
|
|
|
sum /= nx_*ny_*nz_;
|
|
|
|
if( zero )
|
|
{
|
|
#pragma omp parallel for
|
|
for( int i=0; i<nx_; ++i )
|
|
for( int j=0; j<ny_; ++j )
|
|
for( int k=0; k<nz_; ++k )
|
|
(*this)(i,j,k) -= sum;
|
|
|
|
}
|
|
|
|
}
|
|
|
|
//! copies the data from another field with 3D index-based access operator
|
|
template< class array3 >
|
|
void copy( array3& v )
|
|
{
|
|
for( int ix=0; ix<(int)nx_; ++ix )
|
|
for( int iy=0; iy<(int)ny_; ++iy )
|
|
for( int iz=0; iz<(int)nz_; ++iz )
|
|
v(ix,iy,iz) = (*this)(ix,iy,iz);
|
|
|
|
}
|
|
|
|
//! adds the data from another field with 3D index-based access operator
|
|
template< class array3 >
|
|
void copy_add( array3& v )
|
|
{
|
|
for( int ix=0; ix<(int)nx_; ++ix )
|
|
for( int iy=0; iy<(int)ny_; ++iy )
|
|
for( int iz=0; iz<(int)nz_; ++iz )
|
|
v(ix,iy,iz) += (*this)(ix,iy,iz);
|
|
}
|
|
|
|
//! subtract the mean of each oct of the field
|
|
/*! subtracting the mean of each oct implies that a restriction of the field to the coarser level
|
|
* is zero everywhere (needed for Hoffman-Ribak constraints)
|
|
*/
|
|
void subtract_oct_mean( void )
|
|
{
|
|
for( int i=0; i<nx_; i+=2 )
|
|
for( int j=0; j<ny_; j+=2 )
|
|
for( int k=0; k<nz_; k+=2 )
|
|
{
|
|
real_t avg =
|
|
0.125*((*this)(i,j,k)+(*this)(i+1,j,k)+(*this)(i,j+1,k)+(*this)(i,j,k+1)+
|
|
(*this)(i+1,j+1,k)+(*this)(i+1,j,k+1)+(*this)(i,j+1,k+1)+(*this)(i+1,j+1,k+1));
|
|
|
|
(*this)(i,j,k) -= avg;
|
|
(*this)(i+1,j,k) -= avg;
|
|
(*this)(i,j+1,k) -= avg;
|
|
(*this)(i,j,k+1) -= avg;
|
|
(*this)(i+1,j+1,k) -= avg;
|
|
(*this)(i+1,j,k+1) -= avg;
|
|
(*this)(i,j+1,k+1) -= avg;
|
|
(*this)(i+1,j+1,k+1) -= avg;
|
|
|
|
}
|
|
|
|
}
|
|
|
|
//! replaces the value of all cells in an oct with the average value of the oct
|
|
void set_to_oct_mean( void )
|
|
{
|
|
for( int i=0; i<nx_; i+=2 )
|
|
for( int j=0; j<ny_; j+=2 )
|
|
for( int k=0; k<nz_; k+=2 )
|
|
{
|
|
real_t avg =
|
|
0.125*((*this)(i,j,k)+(*this)(i+1,j,k)+(*this)(i,j+1,k)+(*this)(i,j,k+1)+
|
|
(*this)(i+1,j+1,k)+(*this)(i+1,j,k+1)+(*this)(i,j+1,k+1)+(*this)(i+1,j+1,k+1));
|
|
|
|
(*this)(i,j,k) = avg;
|
|
(*this)(i+1,j,k) = avg;
|
|
(*this)(i,j+1,k) = avg;
|
|
(*this)(i,j,k+1) = avg;
|
|
(*this)(i+1,j+1,k) = avg;
|
|
(*this)(i+1,j,k+1) = avg;
|
|
(*this)(i,j+1,k+1) = avg;
|
|
(*this)(i+1,j+1,k+1) = avg;
|
|
|
|
}
|
|
|
|
}
|
|
|
|
//! sets the field to zero inside of a rectangular box
|
|
void zero_subgrid( unsigned oxsub, unsigned oysub, unsigned ozsub, unsigned lxsub, unsigned lysub, unsigned lzsub )
|
|
{
|
|
//... correct offsets for padding (not needed for top grid)
|
|
for( int ix=oxsub; ix<(int)(oxsub+lxsub); ++ix )
|
|
for( int iy=oysub; iy<(int)(oysub+lysub); ++iy )
|
|
for( int iz=ozsub; iz<(int)(ozsub+lzsub); ++iz )
|
|
(*this)(ix,iy,iz) = 0.0;
|
|
|
|
}
|
|
|
|
|
|
//! sets the field to zero inside of a rectangular box including the boundary of the box (used for isolated FFT)
|
|
void zero_padded_subgrid( unsigned oxsub, unsigned oysub, unsigned ozsub, unsigned lxsub, unsigned lysub, unsigned lzsub )
|
|
{
|
|
//... correct offsets for padding (not needed for top grid)
|
|
oxsub -= lxsub/2;
|
|
oysub -= lysub/2;
|
|
ozsub -= lzsub/2;
|
|
lxsub *= 2;
|
|
lysub *= 2;
|
|
lzsub *= 2;
|
|
for( int ix=oxsub; ix<(int)(oxsub+lxsub); ++ix )
|
|
for( int iy=oysub; iy<(int)(oysub+lysub); ++iy )
|
|
for( int iz=ozsub; iz<(int)(ozsub+lzsub); ++iz )
|
|
(*this)(ix,iy,iz) = 0.0;
|
|
|
|
}
|
|
|
|
//! sets the field to zero outside of a rectangular box including the boundary of the box (used for isolated FFT)
|
|
void zero_but_padded_subgrid( int oxsub, int oysub, int ozsub, int lxsub, int lysub, int lzsub )
|
|
{
|
|
//... correct offsets for padding (not needed for top grid)
|
|
oxsub -= lxsub/2;
|
|
oysub -= lysub/2;
|
|
ozsub -= lzsub/2;
|
|
lxsub *= 2;
|
|
lysub *= 2;
|
|
lzsub *= 2;
|
|
|
|
for( int ix=0; ix<nx_; ++ix )
|
|
for( int iy=0; iy<ny_; ++iy )
|
|
for( int iz=0; iz<nz_; ++iz )
|
|
{
|
|
if( (ix<oxsub||ix>=oxsub+lxsub)
|
|
|| (iy<oysub||iy>=oysub+lysub)
|
|
|| (iz<ozsub||iz>=ozsub+lzsub) )
|
|
(*this)(ix,iy,iz) = 0.0;
|
|
|
|
}
|
|
}
|
|
|
|
//! sets the field to zero outside of a rectangular box
|
|
void zero_but_subgrid( int oxsub, int oysub, int ozsub, int lxsub, int lysub, int lzsub )
|
|
{
|
|
for( int ix=0; ix<nx_; ++ix )
|
|
for( int iy=0; iy<ny_; ++iy )
|
|
for( int iz=0; iz<nz_; ++iz )
|
|
{
|
|
if( (ix<oxsub||ix>=oxsub+lxsub)
|
|
|| (iy<oysub||iy>=oysub+lysub)
|
|
|| (iz<ozsub||iz>=ozsub+lzsub) )
|
|
(*this)(ix,iy,iz) = 0.0;
|
|
|
|
}
|
|
}
|
|
|
|
//! sets the field to zero except on the boundary of a rectangular box
|
|
void zero_but_subgrid_bnd( unsigned oxsub, unsigned oysub, unsigned ozsub, unsigned lxsub, unsigned lysub, unsigned lzsub )
|
|
{
|
|
//... correct offsets for padding (not needed for top grid)
|
|
|
|
// zero the subgrid
|
|
for( int ix=oxsub; ix<(int)(oxsub+lxsub); ++ix )
|
|
for( int iy=oysub; iy<(int)(oysub+lysub); ++iy )
|
|
for( int iz=ozsub; iz<(int)(ozsub+lzsub); ++iz )
|
|
(*this)(ix,iy,iz) = 0.0;
|
|
|
|
oxsub -= lxsub/2;
|
|
oysub -= lysub/2;
|
|
ozsub -= lzsub/2;
|
|
lxsub *= 2;
|
|
lysub *= 2;
|
|
lzsub *= 2;
|
|
|
|
// zero the outside, except the boundary
|
|
//#pragma parallel nowait
|
|
{
|
|
|
|
for( int ix=0; ix<oxsub; ++ix )
|
|
for( int iy=0; iy<ny_; ++iy )
|
|
for( int iz=0; iz<nz_; ++iz )
|
|
(*this)(ix,iy,iz) = 0.0;
|
|
|
|
for( int ix=oxsub+lxsub; ix<nx_; ++ix )
|
|
for( int iy=0; iy<ny_; ++iy )
|
|
for( int iz=0; iz<nz_; ++iz )
|
|
(*this)(ix,iy,iz) = 0.0;
|
|
|
|
for( int ix=0; ix<nx_; ++ix )
|
|
for( int iy=0; iy<oysub; ++iy )
|
|
for( int iz=0; iz<nz_; ++iz )
|
|
(*this)(ix,iy,iz) = 0.0;
|
|
for( int ix=0; ix<nx_; ++ix )
|
|
for( int iy=oysub+lysub; iy<ny_; ++iy )
|
|
for( int iz=0; iz<nz_; ++iz )
|
|
(*this)(ix,iy,iz) = 0.0;
|
|
|
|
for( int ix=0; ix<nx_; ++ix )
|
|
for( int iy=0; iy<ny_; ++iy )
|
|
for( int iz=0; iz<ozsub; ++iz )
|
|
(*this)(ix,iy,iz) = 0.0;
|
|
for( int ix=0; ix<nx_; ++ix )
|
|
for( int iy=0; iy<ny_; ++iy )
|
|
for( int iz=ozsub+lzsub; iz<nz_; ++iz )
|
|
(*this)(ix,iy,iz) = 0.0;
|
|
|
|
}
|
|
|
|
}
|
|
|
|
};
|
|
|
|
template< typename real_t >
|
|
class PaddedDensitySubGrid : public DensityGrid<real_t>
|
|
{
|
|
public:
|
|
using DensityGrid<real_t>::nx_;
|
|
using DensityGrid<real_t>::ny_;
|
|
using DensityGrid<real_t>::nz_;
|
|
using DensityGrid<real_t>::data_;
|
|
|
|
using DensityGrid<real_t>::fill_rand;
|
|
using DensityGrid<real_t>::get_data_ptr;
|
|
|
|
int ox_, oy_, oz_;
|
|
public:
|
|
|
|
PaddedDensitySubGrid( int ox, int oy, int oz, unsigned nx, unsigned ny, unsigned nz)
|
|
: DensityGrid<real_t>(2*nx,2*ny,2*nz), ox_(ox), oy_(oy), oz_(oz)
|
|
{
|
|
if( ox-nx/4 < 0 || oy-ny/4 < 0 || oz-nz/4 < 0 )
|
|
throw std::runtime_error("subgrid extends across top grid");
|
|
|
|
//.. the size in top grid cells is nx/2,ny/2,nz/2, so padding starts at ox-nx/4...
|
|
//.. loop over relevant part of top grid and copy down grid values
|
|
/*for( int ix=ox-nx/4,ixs=0; ix<ox+3*nx/4; ix++,ixs+=2 )
|
|
for( int iy=oy-ny/4,iys=0; iy<oy+3*ny/4; iy++,iys+=2 )
|
|
for( int iz=oz-nz/4,izs=0; iz<oz+3*nz/4; iz++,izs+=2 )
|
|
{
|
|
(*this)(ixs,iys,izs) =
|
|
(*this)(ixs+1,iys,izs) =
|
|
(*this)(ixs,iys+1,izs) =
|
|
(*this)(ixs,iys,izs+1) =
|
|
(*this)(ixs+1,iys+1,izs) =
|
|
(*this)(ixs+1,iys,izs+1) =
|
|
(*this)(ixs,iys+1,izs+1) =
|
|
(*this)(ixs+1,iys+1,izs+1) = top(ix,iy,iz);
|
|
}*/
|
|
}
|
|
|
|
PaddedDensitySubGrid( const PaddedDensitySubGrid<real_t>& o )
|
|
: DensityGrid<real_t>( o ), ox_(o.ox_), oy_(o.oy_), oz_(o.oz_)
|
|
{ }
|
|
|
|
void zero_padded_subgrid( unsigned oxsub, unsigned oysub, unsigned ozsub, unsigned lxsub, unsigned lysub, unsigned lzsub )
|
|
{
|
|
//... correct offsets for padding (not needed for top grid)
|
|
|
|
oxsub += nx_/4;
|
|
oysub += ny_/4;
|
|
ozsub += nz_/4;
|
|
|
|
//... correct offsets for padding (not needed for top grid)
|
|
for( int ix=oxsub; ix<(int)(oxsub+lxsub); ++ix )
|
|
for( int iy=oysub; iy<(int)(oysub+lysub); ++iy )
|
|
for( int iz=ozsub; iz<(int)(ozsub+lzsub); ++iz )
|
|
(*this)(ix,iy,iz) = 0.0;
|
|
|
|
}
|
|
|
|
void zero_subgrid( unsigned oxsub, unsigned oysub, unsigned ozsub, unsigned lxsub, unsigned lysub, unsigned lzsub )
|
|
{
|
|
//... correct offsets for padding (not needed for top grid)
|
|
|
|
oxsub += nx_/4;
|
|
oysub += ny_/4;
|
|
ozsub += nz_/4;
|
|
|
|
for( int ix=oxsub; ix<(int)(oxsub+lxsub); ++ix )
|
|
for( int iy=oysub; iy<(int)(oysub+lysub); ++iy )
|
|
for( int iz=ozsub; iz<(int)(ozsub+lzsub); ++iz )
|
|
(*this)(ix,iy,iz) = 0.0;
|
|
|
|
}
|
|
|
|
void zero_but_subgrid_bnd( unsigned oxsub, unsigned oysub, unsigned ozsub, unsigned lxsub, unsigned lysub, unsigned lzsub )
|
|
{
|
|
//... correct offsets for padding (not needed for top grid)
|
|
|
|
oxsub += nx_/4;
|
|
oysub += ny_/4;
|
|
ozsub += nz_/4;
|
|
|
|
//lxsub += nx_/4;
|
|
//lysub += nx_/4;
|
|
//lzsub += nx_/4;
|
|
|
|
// zero the subgrid
|
|
for( int ix=oxsub; ix<(int)(oxsub+lxsub); ++ix )
|
|
for( int iy=oysub; iy<(int)(oysub+lysub); ++iy )
|
|
for( int iz=ozsub; iz<(int)(ozsub+lzsub); ++iz )
|
|
(*this)(ix,iy,iz) = 0.0;
|
|
|
|
oxsub -= lxsub/2;
|
|
oysub -= lysub/2;
|
|
ozsub -= lzsub/2;
|
|
lxsub *= 2;
|
|
lysub *= 2;
|
|
lzsub *= 2;
|
|
|
|
// zero the outside, except the boundary
|
|
//#pragma parallel block nowait
|
|
{
|
|
|
|
for( int ix=0; ix<oxsub; ++ix )
|
|
for( int iy=0; iy<ny_; ++iy )
|
|
for( int iz=0; iz<nz_; ++iz )
|
|
(*this)(ix,iy,iz) = 0.0;
|
|
|
|
std::cerr << oxsub+lxsub << " -> " << nx_ << std::endl;
|
|
for( int ix=oxsub+lxsub; ix<nx_; ++ix )
|
|
for( int iy=0; iy<ny_; ++iy )
|
|
for( int iz=0; iz<nz_; ++iz )
|
|
(*this)(ix,iy,iz) = 0.0;
|
|
|
|
for( int ix=0; ix<nx_; ++ix )
|
|
for( int iy=0; iy<oysub; ++iy )
|
|
for( int iz=0; iz<nz_; ++iz )
|
|
(*this)(ix,iy,iz) = 0.0;
|
|
for( int ix=0; ix<nx_; ++ix )
|
|
for( int iy=oysub+lysub; iy<ny_; ++iy )
|
|
for( int iz=0; iz<nz_; ++iz )
|
|
(*this)(ix,iy,iz) = 0.0;
|
|
|
|
for( int ix=0; ix<nx_; ++ix )
|
|
for( int iy=0; iy<ny_; ++iy )
|
|
for( int iz=0; iz<ozsub; ++iz )
|
|
(*this)(ix,iy,iz) = 0.0;
|
|
for( int ix=0; ix<nx_; ++ix )
|
|
for( int iy=0; iy<ny_; ++iy )
|
|
for( int iz=ozsub+lzsub; iz<nz_; ++iz )
|
|
(*this)(ix,iy,iz) = 0.0;
|
|
|
|
}
|
|
}
|
|
|
|
|
|
template< class array3 >
|
|
void upload_bnd_add( array3& top )
|
|
{
|
|
int ox=ox_-nx_/8, oy=oy_-ny_/8, oz=oz_-nz_/8;
|
|
int lx=ox+nx_/2, ly=oy+ny_/2, lz=oz+nz_/2;
|
|
ox = std::max(ox,0);
|
|
oy = std::max(oy,0);
|
|
oz = std::max(oz,0);
|
|
lx = std::min(lx,(int)top.size(0));
|
|
ly = std::min(ly,(int)top.size(1));
|
|
lz = std::min(lz,(int)top.size(2));
|
|
|
|
|
|
if( ox < 0 || oy < 0 || oz < 0 ){
|
|
|
|
std::cerr << "offset = " << ox << ", " << oy << ", " << oz << std::endl;
|
|
std::cerr << "nx = " << nx_ << ", " << ny_ << ", " << nz_ << std::endl;
|
|
|
|
throw std::runtime_error("Subgrid extends beyond top grid. Increase levelmin or reduce size of refinement region!");
|
|
}
|
|
for( int ix=0,ixu=ox;ix<nx_&&ixu<lx;ix+=2,++ixu )
|
|
for( int iy=0,iyu=oy;iy<ny_&&iyu<ly;iy+=2,++iyu )
|
|
for( int iz=0,izu=oz;iz<nz_&&izu<lz;iz+=2,++izu )
|
|
{
|
|
if( (ix<nx_/4||ix>=3*nx_/4)
|
|
|| (iy<ny_/4||iy>=3*ny_/4)
|
|
|| (iz<nz_/4||iz>=3*nz_/4) )
|
|
{
|
|
top(ixu,iyu,izu) += 0.125* ((*this)(ix,iy,iz)+(*this)(ix+1,iy,iz)+(*this)(ix,iy+1,iz)+(*this)(ix,iy,iz+1)
|
|
+(*this)(ix+1,iy+1,iz)+(*this)(ix+1,iy,iz+1)+(*this)(ix,iy+1,iz+1)+(*this)(ix+1,iy+1,iz+1) );
|
|
|
|
|
|
}
|
|
}
|
|
}
|
|
|
|
void constrain( const DensityGrid<real_t>& top )
|
|
{
|
|
int ox=ox_-nx_/8, oy=oy_-ny_/8, oz=oz_-nz_/8;
|
|
|
|
if( ox < 0 || oy < 0 || oz < 0 ){
|
|
|
|
std::cerr << "offset = " << ox << ", " << oy << ", " << oz << std::endl;
|
|
std::cerr << "nx = " << nx_ << ", " << ny_ << ", " << nz_ << std::endl;
|
|
|
|
throw std::runtime_error("Subgrid extends beyond top grid. Increase levelmin or reduce size of refinement region!");
|
|
}
|
|
|
|
//...top grid is not padded
|
|
|
|
//.. fine grid will be multiplied by sqrt(8) later to account for increase in variance
|
|
//.. so, we need to divide the coarse grid value by it
|
|
|
|
real_t fac = 1.0/sqrt(8);
|
|
for( int ixs=0,ix=ox;ixs<nx_;ixs+=2,++ix )
|
|
for( int iys=0,iy=oy;iys<ny_;iys+=2,++iy )
|
|
for( int izs=0,iz=oz;izs<nz_;izs+=2,++iz )
|
|
{
|
|
real_t mean = 0.125 *
|
|
((*this)(ixs,iys,izs)+
|
|
(*this)(ixs+1,iys,izs)+
|
|
(*this)(ixs,iys+1,izs)+
|
|
(*this)(ixs,iys,izs+1)+
|
|
(*this)(ixs+1,iys+1,izs)+
|
|
(*this)(ixs+1,iys,izs+1)+
|
|
(*this)(ixs,iys+1,izs+1)+
|
|
(*this)(ixs+1,iys+1,izs+1));
|
|
|
|
|
|
double aa = top(ix,iy,iz)*fac - mean;
|
|
|
|
(*this)(ixs,iys,izs) += aa;
|
|
(*this)(ixs+1,iys,izs) += aa;
|
|
(*this)(ixs,iys+1,izs) += aa;
|
|
(*this)(ixs,iys,izs+1) += aa;
|
|
(*this)(ixs+1,iys+1,izs) += aa;
|
|
(*this)(ixs+1,iys,izs+1) += aa;
|
|
(*this)(ixs,iys+1,izs+1) += aa;
|
|
(*this)(ixs+1,iys+1,izs+1) += aa;
|
|
}
|
|
|
|
#if 0
|
|
unsigned nx=nx_/2, ny=ny_/2, nz=nz_/2;
|
|
|
|
int ox=ox_, oy=oy_, oz=oz_;
|
|
|
|
if( ox-nx/4 < 0 || oy-ny/4 < 0 || oz-nz/4 < 0 )
|
|
throw("subgrid extends across top grid");
|
|
|
|
//.. the size in top grid cells is nx/2,ny/2,nz/2, so padding starts at ox-nx/4...
|
|
//.. loop over relevant part of top grid and copy down grid values
|
|
for( int ix=ox-nx/4,ixs=0; ix<ox+3*nx/4; ix++,ixs+=2 )
|
|
for( int iy=oy-ny/4,iys=0; iy<oy+3*ny/4; iy++,iys+=2 )
|
|
for( int iz=oz-nz/4,izs=0; iz<oz+3*nz/4; iz++,izs+=2 )
|
|
{
|
|
//... removed -r, as the mean will be taken out right at the beginning
|
|
(*this)(ixs,iys,izs) += top(ix,iy,iz);
|
|
(*this)(ixs+1,iys,izs) += top(ix,iy,iz);
|
|
(*this)(ixs,iys+1,izs) += top(ix,iy,iz);
|
|
(*this)(ixs,iys,izs+1) += top(ix,iy,iz);
|
|
(*this)(ixs+1,iys+1,izs) += top(ix,iy,iz);
|
|
(*this)(ixs+1,iys,izs+1) += top(ix,iy,iz);
|
|
(*this)(ixs,iys+1,izs+1) += top(ix,iy,iz);
|
|
(*this)(ixs+1,iys+1,izs+1) += top(ix,iy,iz);
|
|
}
|
|
#endif
|
|
}
|
|
|
|
void constrain( const PaddedDensitySubGrid<real_t>& top )
|
|
{
|
|
///... top grid is padded too...
|
|
//.. ox is shifted by nx_/4, padded overlap starts at -nx_/8
|
|
int ox=ox_+top.nx_/4-nx_/8, oy=oy_+top.ny_/4-ny_/8, oz=oz_+top.nz_/4-nz_/8;
|
|
if( ox < 0 || oy < 0 || oz < 0 )
|
|
throw std::runtime_error("subgrid extends beyond top grid");
|
|
double fac = 1.0/sqrt(8.0);
|
|
for( int ixs=0,ix=ox;ixs<nx_;ixs+=2,++ix )
|
|
for( int iys=0,iy=oy;iys<ny_;iys+=2,++iy )
|
|
for( int izs=0,iz=oz;izs<nz_;izs+=2,++iz )
|
|
{
|
|
|
|
double mean = 0.125 *
|
|
((*this)(ixs,iys,izs)+
|
|
(*this)(ixs+1,iys,izs)+
|
|
(*this)(ixs,iys+1,izs)+
|
|
(*this)(ixs,iys,izs+1)+
|
|
(*this)(ixs+1,iys+1,izs)+
|
|
(*this)(ixs+1,iys,izs+1)+
|
|
(*this)(ixs,iys+1,izs+1)+
|
|
(*this)(ixs+1,iys+1,izs+1));
|
|
|
|
|
|
double aa = top(ix,iy,iz)*fac - mean;
|
|
|
|
(*this)(ixs,iys,izs) += aa;
|
|
(*this)(ixs+1,iys,izs) += aa;
|
|
(*this)(ixs,iys+1,izs) += aa;
|
|
(*this)(ixs,iys,izs+1) += aa;
|
|
(*this)(ixs+1,iys+1,izs) += aa;
|
|
(*this)(ixs+1,iys,izs+1) += aa;
|
|
(*this)(ixs,iys+1,izs+1) += aa;
|
|
(*this)(ixs+1,iys+1,izs+1) += aa;
|
|
}
|
|
|
|
}
|
|
|
|
|
|
template< class array3 >
|
|
void upload_bnd( unsigned nbnd, array3& v )
|
|
{
|
|
for( int ix=nx_/4-2*nbnd,ixu=-nbnd; ix<3*nx_/4+2*nbnd; ix+=2,++ixu)
|
|
for( int iy=ny_/4-2*nbnd,iyu=-nbnd; iy<3*ny_/4+2*nbnd; iy+=2,++iyu )
|
|
for( int iz=nz_/4-2*nbnd,izu=-nbnd; iz<3*nz_/4+2*nbnd; iz+=2,++izu )
|
|
{
|
|
if( (ix<nx_/4||ix>=3*nx_/4)
|
|
|| (iy<ny_/4||iy>=3*ny_/4)
|
|
|| (iz<nz_/4||iz>=3*nz_/4) )
|
|
{
|
|
v(ixu+ox_,iyu+oy_,izu+oz_)
|
|
= 1.0;
|
|
}
|
|
}
|
|
}
|
|
|
|
template< class array3 >
|
|
void copy_unpad( array3& v )
|
|
{
|
|
for( int ix=nx_/4,ixu=0; ix<3*nx_/4; ++ix,++ixu )
|
|
for( int iy=ny_/4,iyu=0; iy<3*ny_/4; ++iy,++iyu )
|
|
for( int iz=nz_/4,izu=0; iz<3*nz_/4; ++iz,++izu )
|
|
v(ixu,iyu,izu) = (*this)(ix,iy,iz);
|
|
}
|
|
|
|
template< class array3 >
|
|
void copy_add_unpad( array3& v )
|
|
{
|
|
for( int ix=nx_/4,ixu=0; ix<3*nx_/4; ++ix,++ixu )
|
|
for( int iy=ny_/4,iyu=0; iy<3*ny_/4; ++iy,++iyu )
|
|
for( int iz=nz_/4,izu=0; iz<3*nz_/4; ++iz,++izu )
|
|
v(ixu,iyu,izu) += (*this)(ix,iy,iz);
|
|
}
|
|
|
|
template< class array3 >
|
|
void copy_subtract_unpad( array3& v )
|
|
{
|
|
for( int ix=nx_/4,ixu=0; ix<3*nx_/4; ++ix,++ixu )
|
|
for( int iy=ny_/4,iyu=0; iy<3*ny_/4; ++iy,++iyu )
|
|
for( int iz=nz_/4,izu=0; iz<3*nz_/4; ++iz,++izu )
|
|
v(ixu,iyu,izu) -= (*this)(ix,iy,iz);
|
|
}
|
|
|
|
void zero_boundary( void )
|
|
{
|
|
for( int ix=0; ix<nx_/4; ++ix )
|
|
for( int iy=0; iy<ny_; ++iy )
|
|
for( int iz=0; iz<nz_; ++iz )
|
|
{
|
|
(*this)(ix,iy,iz) = 0.0;
|
|
(*this)(ix+3*nx_/4,iy,iz) = 0.0;
|
|
}
|
|
|
|
for( int ix=0; ix<nx_; ++ix )
|
|
for( int iy=0; iy<ny_/4; ++iy )
|
|
for( int iz=0; iz<nz_; ++iz )
|
|
{
|
|
(*this)(ix,iy,iz) = 0.0;
|
|
(*this)(ix,iy+3*ny_/4,iz) = 0.0;
|
|
}
|
|
|
|
|
|
for( int ix=0; ix<nx_; ++ix )
|
|
for( int iy=0; iy<ny_; ++iy )
|
|
for( int iz=0; iz<nz_/4; ++iz )
|
|
{
|
|
(*this)(ix,iy,iz) = 0.0;
|
|
(*this)(ix,iy,iz+3*nz_/4) = 0.0;
|
|
}
|
|
|
|
}
|
|
|
|
|
|
};
|
|
|
|
template< typename M >
|
|
inline void enforce_coarse_mean( M& v, M& V )
|
|
{
|
|
double outersum =0.0, innersum = 0.0;
|
|
int
|
|
nx = v.size(0)/2,
|
|
ny = v.size(1)/2,
|
|
nz = v.size(2)/2,
|
|
ox = v.offset(0),
|
|
oy = v.offset(1),
|
|
oz = v.offset(2);
|
|
|
|
unsigned innercount = 0, outercount = 0;
|
|
for( int i=0; i<V.size(0); ++i )
|
|
for( int j=0; j<V.size(1); ++j )
|
|
for( int k=0; k<V.size(2); ++k )
|
|
{
|
|
if( (i > ox && i < ox+nx) &&
|
|
(j > oy && j < oy+ny) &&
|
|
(k > oz && k < oz+nz ))
|
|
{
|
|
innersum += V(i,j,k);
|
|
++innercount;
|
|
}
|
|
else
|
|
{
|
|
outersum += V(i,j,k);
|
|
++outercount;
|
|
}
|
|
|
|
}
|
|
innersum /= innercount;
|
|
outersum /= outercount;
|
|
|
|
double dcorr = innersum * innercount / outercount;
|
|
|
|
for( int i=0; i<V.size(0); ++i )
|
|
for( int j=0; j<V.size(1); ++j )
|
|
for( int k=0; k<V.size(2); ++k )
|
|
{
|
|
if( !((i > ox && i < ox+nx) &&
|
|
(j > oy && j < oy+ny) &&
|
|
(k > oz && k < oz+nz )))
|
|
V(i,j,k) -= outersum + dcorr;//innersum;
|
|
}
|
|
|
|
}
|
|
|
|
template< typename M >
|
|
inline void enforce_mean( M& v, M& V )
|
|
{
|
|
|
|
int
|
|
nx = v.size(0)/2,
|
|
ny = v.size(1)/2,
|
|
nz = v.size(2)/2,
|
|
ox = v.offset(0),
|
|
oy = v.offset(1),
|
|
oz = v.offset(2);
|
|
|
|
|
|
double coarsemean = 0.0, finemean = 0.0;
|
|
unsigned count = 0;
|
|
#pragma omp parallel for reduction(+:coarsemean,finemean,count)
|
|
for( int i=0; i<nx; ++i )
|
|
{
|
|
int i2 = 2*i;
|
|
for( int j=0,j2=0; j<ny; ++j,j2+=2 )
|
|
for( int k=0,k2=0; k<nz; ++k,k2+=2 )
|
|
{
|
|
|
|
finemean += 0.125 * ( v(i2+1,j2,k2) + v(i2,j2+1,k2) + v(i2,j2,k2+1) + v(i2+1,j2+1,k2) +
|
|
v(i2+1,j2,k2+1) + v(i2+1,j2+1,k2+1) + v(i2,j2+1,k2+1) + v(i2,j2,k2) );
|
|
|
|
coarsemean += V(i+ox,j+oy,k+oz);
|
|
++count;
|
|
}
|
|
}
|
|
|
|
coarsemean /= count;
|
|
finemean /= count;
|
|
|
|
double dmean = coarsemean-finemean;
|
|
|
|
#pragma omp parallel for reduction(+:coarsemean,finemean)
|
|
for( int i=0; i<nx; ++i )
|
|
{
|
|
int i2 = 2*i;
|
|
for( int j=0,j2=0; j<ny; ++j,j2+=2 )
|
|
for( int k=0,k2=0; k<nz; ++k,k2+=2 )
|
|
{
|
|
|
|
v(i2+1,j2,k2) += dmean;
|
|
v(i2,j2+1,k2) += dmean;
|
|
v(i2,j2,k2+1) += dmean;
|
|
v(i2+1,j2+1,k2) += dmean;
|
|
v(i2+1,j2,k2+1) += dmean;
|
|
v(i2+1,j2+1,k2+1) += dmean;
|
|
v(i2,j2+1,k2+1) += dmean;
|
|
v(i2,j2,k2) += dmean;
|
|
}
|
|
}
|
|
|
|
}
|
|
#endif
|
|
|