From fb90897699757597576df5fbbdd4c146e8fe91ae Mon Sep 17 00:00:00 2001 From: Oliver Hahn Date: Thu, 16 Feb 2023 21:10:19 -0800 Subject: [PATCH] WIP some fixes to noise generators --- src/densities.cc | 52 ++-- src/main.cc | 2 +- src/mesh.hh | 266 ++++++++++--------- src/mg_solver.hh | 17 +- src/plugins/random_music.cc | 14 +- src/plugins/random_music_wnoise_generator.cc | 22 +- src/plugins/random_panphasia.cc | 14 +- 7 files changed, 199 insertions(+), 188 deletions(-) diff --git a/src/densities.cc b/src/densities.cc index f285e58..10b4ecc 100644 --- a/src/densities.cc +++ b/src/densities.cc @@ -23,7 +23,7 @@ double Meyer_scaling_function( double k, double kmax ) constexpr double fourpithirds{4.0*M_PI/3.0}; auto nu = []( double x ){ return x<0.0?0.0:(x<1.0?x:1.0); }; - k = k/kmax * fourpithirds; + k = std::abs(k)/kmax * fourpithirds; if( k < twopithirds ) return 1.0; else if( k< fourpithirds ){ @@ -126,15 +126,15 @@ void fft_interpolate(m1 &V, m2 &v, bool from_basegrid = false) if (!from_basegrid) { - oxf += mxf - mxf/2; - oyf += myf - myf/2; - ozf += mzf - mzf/2; + oxf += mxf/2; + oyf += myf/2; + ozf += mzf/2; } else { - oxf -= mxf/2; //nxf / 8; - oyf -= myf/2; //nyf / 8; - ozf -= mzf/2; //nzf / 8; + oxf -= mxf/2; + oyf -= myf/2; + ozf -= mzf/2; } music::ulog.Print("FFT interpolate: offset=%d,%d,%d size=%d,%d,%d", oxf, oyf, ozf, nxf, nyf, nzf); @@ -185,10 +185,9 @@ void fft_interpolate(m1 &V, m2 &v, bool from_basegrid = false) /*************************************************/ //.. perform actual interpolation double fftnorm = 1.0 / ((double)nxf * (double)nyf * (double)nzf); - double sqrt8 = 8.0; //sqrt(8.0); - double phasefac = -0.5; // this enables filtered splicing of coarse and fine modes + #pragma omp parallel for for (int i = 0; i < (int)nxc; i++) for (int j = 0; j < (int)nyc; j++) for (int k = 0; k < (int)nzc / 2 + 1; k++) @@ -210,20 +209,28 @@ void fft_interpolate(m1 &V, m2 &v, bool from_basegrid = false) double ky = (j <= (int)nyc / 2) ? (double)j : (double)(j - (int)nyc); double kz = (k <= (int)nzc / 2) ? (double)k : (double)(k - (int)nzc); - double phase = phasefac * (kx / nxc + ky / nyc + kz / nzc) * M_PI; + double phase = -0.5 * M_PI * (kx / nxc + ky / nyc + kz / nzc); std::complex val_phas(cos(phase), sin(phase)); std::complex val(RE(ccoarse[qc]), IM(ccoarse[qc])); - val *= sqrt8 * val_phas; + val *= val_phas * 8.0; + + if (false){ //(i != (int)nxc / 2 && j != (int)nyc / 2 && k != (int)nzc / 2){ + double blend_coarse_x = Meyer_scaling_function(kx, nxc / 2); + double blend_coarse_y = Meyer_scaling_function(ky, nyc / 2); + double blend_coarse_z = Meyer_scaling_function(kz, nzc / 2); + double blend_coarse = blend_coarse_x*blend_coarse_y*blend_coarse_z; + double blend_fine = 1.0-blend_coarse; - if (i != (int)nxc / 2 && j != (int)nyc / 2 && k != (int)nzc / 2){ - double blend_coarse = Meyer_scaling_function(kx, nxc / 2) * Meyer_scaling_function(ky, nyc / 2) * Meyer_scaling_function(kz, nzc / 2); - double blend_fine = 1.0 - blend_coarse; - RE(cfine[qf]) = blend_fine * RE(cfine[qf]) + blend_coarse * val.real(); IM(cfine[qf]) = blend_fine * IM(cfine[qf]) + blend_coarse * val.imag(); } + + if (i != (int)nxc / 2 && j != (int)nyc / 2 && k != (int)nzc / 2){ + RE(cfine[qf]) = val.real(); + IM(cfine[qf]) = val.imag(); + } } delete[] rcoarse; @@ -451,8 +458,7 @@ void normalize_density(grid_hierarchy &delta) ny = delta.get_grid(levelmin)->size(1); nz = delta.get_grid(levelmin)->size(2); -#pragma omp parallel for reduction(+ \ - : sum) +#pragma omp parallel for reduction(+ : sum) for (int ix = 0; ix < (int)nx; ++ix) for (size_t iy = 0; iy < ny; ++iy) for (size_t iz = 0; iz < nz; ++iz) @@ -481,9 +487,7 @@ void normalize_density(grid_hierarchy &delta) void coarsen_density(const refinement_hierarchy &rh, GridHierarchy &u, bool kspace) { unsigned levelmin_TF = u.levelmin(); - - /*for( int i=rh.levelmax(); i>0; --i ) - mg_straight().restrict( *(u.get_grid(i)), *(u.get_grid(i-1)) );*/ + bool benforce_coarse = true;//!kspace; if (kspace) { @@ -498,13 +502,11 @@ void coarsen_density(const refinement_hierarchy &rh, GridHierarchy &u, b for (unsigned i = 1; i <= rh.levelmax(); ++i) { - if (rh.offset(i, 0) != u.get_grid(i)->offset(0) || rh.offset(i, 1) != u.get_grid(i)->offset(1) || rh.offset(i, 2) != u.get_grid(i)->offset(2) || rh.size(i, 0) != u.get_grid(i)->size(0) || rh.size(i, 1) != u.get_grid(i)->size(1) || rh.size(i, 2) != u.get_grid(i)->size(2)) + if (rh.offset(i, 0) != u.get_grid(i)->offset(0) || rh.offset(i, 1) != u.get_grid(i)->offset(1) || rh.offset(i, 2) != u.get_grid(i)->offset(2) + || rh.size(i, 0) != u.get_grid(i)->size(0) || rh.size(i, 1) != u.get_grid(i)->size(1) || rh.size(i, 2) != u.get_grid(i)->size(2)) { u.cut_patch(i, rh.offset_abs(i, 0), rh.offset_abs(i, 1), rh.offset_abs(i, 2), - rh.size(i, 0), rh.size(i, 1), rh.size(i, 2)); + rh.size(i, 0), rh.size(i, 1), rh.size(i, 2), benforce_coarse); } } - - for (int i = rh.levelmax(); i > 0; --i) - mg_straight().restrict(*(u.get_grid(i)), *(u.get_grid(i - 1))); } diff --git a/src/main.cc b/src/main.cc index 64d341e..629b000 100644 --- a/src/main.cc +++ b/src/main.cc @@ -456,7 +456,7 @@ int main(int argc, const char *argv[]) if (!cf.contains_key("setup", "kspace_TF")) cf.insert_value("setup", "kspace_TF", "yes"); - bool bspectral_sampling = cf.get_value("setup", "kspace_TF"); + bool bspectral_sampling = true;//cf.get_value("setup", "kspace_TF"); if (bspectral_sampling) music::ilog.Print("Using k-space sampled transfer functions..."); diff --git a/src/mesh.hh b/src/mesh.hh index 18cfa3f..af43f90 100644 --- a/src/mesh.hh +++ b/src/mesh.hh @@ -8,8 +8,7 @@ */ -#ifndef __MESH_HH -#define __MESH_HH +#pragma once #include #include @@ -27,7 +26,6 @@ using index_t = ptrdiff_t; using index3_t = std::array; using vec3_t = std::array; - class refinement_mask { protected: @@ -147,9 +145,11 @@ public: m_pdata = new real_t[m_nx * m_ny * m_nz]; - if (copy_over) + if (copy_over){ + #pragma omp parallel for for (size_t i = 0; i < m_nx * m_ny * m_nz; ++i) m_pdata[i] = m.m_pdata[i]; + } } //! standard copy constructor @@ -165,6 +165,7 @@ public: m_pdata = new real_t[m_nx * m_ny * m_nz]; + #pragma omp parallel for for (size_t i = 0; i < m_nx * m_ny * m_nz; ++i) m_pdata[i] = m.m_pdata[i]; } @@ -217,6 +218,7 @@ public: //! set all the data to zero values void zero(void) { + #pragma omp parallel for for (size_t i = 0; i < m_nx * m_ny * m_nz; ++i) m_pdata[i] = 0.0; } @@ -258,6 +260,7 @@ public: //! direct multiplication of the whole data block with a number Meshvar &operator*=(real_t x) { + #pragma omp parallel for for (size_t i = 0; i < m_nx * m_ny * m_nz; ++i) m_pdata[i] *= x; return *this; @@ -266,6 +269,7 @@ public: //! direct addition of a number to the whole data block Meshvar &operator+=(real_t x) { + #pragma omp parallel for for (size_t i = 0; i < m_nx * m_ny * m_nz; ++i) m_pdata[i] += x; return *this; @@ -274,6 +278,7 @@ public: //! direct element-wise division of the whole data block by a number Meshvar &operator/=(real_t x) { + #pragma omp parallel for for (size_t i = 0; i < m_nx * m_ny * m_nz; ++i) m_pdata[i] /= x; return *this; @@ -282,6 +287,7 @@ public: //! direct subtraction of a number from the whole data block Meshvar &operator-=(real_t x) { + #pragma omp parallel for for (size_t i = 0; i < m_nx * m_ny * m_nz; ++i) m_pdata[i] -= x; return *this; @@ -295,6 +301,7 @@ public: music::elog.Print("Meshvar::operator*= : attempt to operate on incompatible data"); throw std::runtime_error("Meshvar::operator*= : attempt to operate on incompatible data"); } + #pragma omp parallel for for (size_t i = 0; i < m_nx * m_ny * m_nz; ++i) m_pdata[i] *= v.m_pdata[i]; @@ -310,6 +317,7 @@ public: throw std::runtime_error("Meshvar::operator/= : attempt to operate on incompatible data"); } + #pragma omp parallel for for (size_t i = 0; i < m_nx * m_ny * m_nz; ++i) m_pdata[i] /= v.m_pdata[i]; @@ -324,6 +332,7 @@ public: music::elog.Print("Meshvar::operator+= : attempt to operate on incompatible data"); throw std::runtime_error("Meshvar::operator+= : attempt to operate on incompatible data"); } + #pragma omp parallel for for (size_t i = 0; i < m_nx * m_ny * m_nz; ++i) m_pdata[i] += v.m_pdata[i]; @@ -338,6 +347,7 @@ public: music::elog.Print("Meshvar::operator-= : attempt to operate on incompatible data"); throw std::runtime_error("Meshvar::operator-= : attempt to operate on incompatible data"); } + #pragma omp parallel for for (size_t i = 0; i < m_nx * m_ny * m_nz; ++i) m_pdata[i] -= v.m_pdata[i]; @@ -360,6 +370,7 @@ public: m_pdata = new real_t[m_nx * m_ny * m_nz]; + #pragma omp parallel for for (size_t i = 0; i < m_nx * m_ny * m_nz; ++i) m_pdata[i] = m.m_pdata[i]; @@ -462,52 +473,13 @@ public: m_pdata = new real_t[m_nx * m_ny * m_nz]; } + #pragma omp parallel for for (size_t i = 0; i < m_nx * m_ny * m_nz; ++i) this->m_pdata[i] = m.m_pdata[i]; return *this; } - //! sets the value of all ghost zones to zero - void zero_bnd(void) - { - - int nx, ny, nz; - nx = this->size(0); - ny = this->size(1); - nz = this->size(2); - - for (int j = -m_nbnd; j < ny + m_nbnd; ++j) - for (int k = -m_nbnd; k < nz + m_nbnd; ++k) - { - for (int i = -m_nbnd; i < 0; ++i) - { - (*this)(i, j, k) = 0.0; - (*this)(nx - 1 - i, j, k) = 0.0; - } - } - - for (int i = -m_nbnd; i < nx + m_nbnd; ++i) - for (int k = -m_nbnd; k < nz + m_nbnd; ++k) - { - for (int j = -m_nbnd; j < 0; ++j) - { - (*this)(i, j, k) = 0.0; - (*this)(i, ny - j - 1, k) = 0.0; - } - } - - for (int i = -m_nbnd; i < nx + m_nbnd; ++i) - for (int j = -m_nbnd; j < ny + m_nbnd; ++j) - { - for (int k = -m_nbnd; k < 0; ++k) - { - (*this)(i, j, k) = 0.0; - (*this)(i, j, nz - k - 1) = 0.0; - } - } - } - //! outputs the data, for debugging only, not practical for large datasets void print(void) const { @@ -1121,8 +1093,9 @@ public: * @param nx new x-extent in fine grid cells * @param ny new y-extent in fine grid cells * @param nz new z-extent in fine grid cells + * @param enforce_coarse_mean enforces the average of 8 cells on a fine level to equal the coarse */ - void cut_patch(unsigned ilevel, unsigned xoff, unsigned yoff, unsigned zoff, unsigned nx, unsigned ny, unsigned nz) + void cut_patch(unsigned ilevel, unsigned xoff, unsigned yoff, unsigned zoff, unsigned nx, unsigned ny, unsigned nz, bool enforce_coarse_mean) { unsigned dx, dy, dz, dxtop, dytop, dztop; @@ -1139,10 +1112,19 @@ public: MeshvarBnd *mnew = new MeshvarBnd(m_nbnd, nx, ny, nz, dxtop, dytop, dztop); //... copy data - for (unsigned i = 0; i < nx; ++i) - for (unsigned j = 0; j < ny; ++j) - for (unsigned k = 0; k < nz; ++k) - (*mnew)(i, j, k) = (*m_pgrids[ilevel])(i + dx, j + dy, k + dz); + [[maybe_unused]] double coarsesum = 0.0, finesum = 0.0; + [[maybe_unused]] size_t coarsecount = 0, finecount = 0; + + //... copy data + #pragma omp parallel for reduction(+:finesum,finecount) collapse(3) + for( unsigned i=0; ioffset(2) -= dz; } + if( enforce_coarse_mean ) + { + //... enforce top mean density over same patch + if (ilevel > levelmin()) + { + int ox = m_pgrids[ilevel]->offset(0); + int oy = m_pgrids[ilevel]->offset(1); + int oz = m_pgrids[ilevel]->offset(2); + + #pragma omp parallel for reduction(+:coarsesum,coarsecount) collapse(3) + for (unsigned i = 0; i < nx / 2; ++i) + for (unsigned j = 0; j < ny / 2; ++j) + for (unsigned k = 0; k < nz / 2; ++k) + { + coarsesum += (*m_pgrids[ilevel - 1])(i + ox, j + oy, k + oz); + coarsecount++; + } + + coarsesum /= (double)coarsecount; + finesum /= (double)finecount; + + #pragma omp parallel for collapse(3) + for (unsigned i = 0; i < nx; ++i) + for (unsigned j = 0; j < ny; ++j) + for (unsigned k = 0; k < nz; ++k) + (*mnew)(i, j, k) += (coarsesum - finesum); + + music::ilog.Print(" .level %d : corrected patch overlap mean value by %f", ilevel, coarsesum - finesum); + } + } + find_new_levelmin(); } @@ -1207,14 +1220,14 @@ class refinement_hierarchy std::vector len_; unsigned - levelmin_, //!< minimum grid level for Poisson solver - levelmax_, //!< maximum grid level for all operations - levelmin_tf_, //!< minimum grid level for density calculation - padding_, //!< padding in number of coarse cells between refinement levels - blocking_factor_, //!< blocking factor of grids, necessary fo BoxLib codes such as NyX - gridding_unit_; //!< internal blocking factor of grids, necessary for Panphasia + levelmin_, //!< minimum grid level for Poisson solver + levelmax_, //!< maximum grid level for all operations + levelmin_tf_, //!< minimum grid level for density calculation + padding_, //!< padding in number of coarse cells between refinement levels + blocking_factor_, //!< blocking factor of grids, necessary fo BoxLib codes such as NyX + gridding_unit_; //!< internal blocking factor of grids, necessary for Panphasia - int margin_; //!< number of cells used for additional padding for convolutions with isolated boundaries (-1 = double padding) + int margin_; //!< number of cells used for additional padding for convolutions with isolated boundaries (-1 = double padding) config_file &cf_; //!< reference to config_file @@ -1232,28 +1245,29 @@ class refinement_hierarchy double rshift_[3]; //! calculates the greatest common divisor - int gcd(int a, int b) const { - return b == 0 ? a : gcd(b, a % b); - } + int gcd(int a, int b) const + { + return b == 0 ? a : gcd(b, a % b); + } - // calculates the cell shift in units of levelmin grid cells if there is an additional constraint to be - // congruent with another grid that partitions the same space in multiples of "base_unit" - int get_shift_unit( int base_unit, int levelmin ) const { - /*int Lp = 0; - while( base_unit * (1<( 1, (1<<(levelmin+1)) / (2*gcd(U,1<<(levelmin+1) )) );*/ + return std::max( 1, (1<<(levelmin+1)) / (2*gcd(U,1<<(levelmin+1) )) );*/ - int level_m = 0; - while( base_unit * (1<( 1, (1<(1, (1 << levelmin) / gcd(base_unit * (1 << level_m), (1 << levelmin))); + } public: //! copy constructor @@ -1275,21 +1289,23 @@ public: preserve_dims_ = cf_.get_value_safe("setup", "preserve_dims", false); equal_extent_ = cf_.get_value_safe("setup", "force_equal_extent", false); blocking_factor_ = cf.get_value_safe("setup", "blocking_factor", 0); - margin_ = cf.get_value_safe("setup","convolution_margin",32); + margin_ = cf.get_value_safe("setup", "convolution_margin", 32); bool bnoshift = cf_.get_value_safe("setup", "no_shift", false); bool force_shift = cf_.get_value_safe("setup", "force_shift", false); gridding_unit_ = cf.get_value_safe("setup", "gridding_unit", 2); - if (gridding_unit_ != 2 && blocking_factor_==0) { - blocking_factor_ = gridding_unit_; // THIS WILL LIKELY CAUSE PROBLEMS WITH NYX - }else if (gridding_unit_ != 2 && blocking_factor_!=0 && gridding_unit_!=blocking_factor_ ) { - music::elog.Print("incompatible gridding unit %d and blocking factor specified", gridding_unit_, blocking_factor_ ); + if (gridding_unit_ != 2 && blocking_factor_ == 0) + { + blocking_factor_ = gridding_unit_; // THIS WILL LIKELY CAUSE PROBLEMS WITH NYX + } + else if (gridding_unit_ != 2 && blocking_factor_ != 0 && gridding_unit_ != blocking_factor_) + { + music::elog.Print("incompatible gridding unit %d and blocking factor specified", gridding_unit_, blocking_factor_); throw std::runtime_error("Incompatible gridding unit and blocking factor!"); } - //... call the region generator if (levelmin_ != levelmax_) { @@ -1312,8 +1328,8 @@ public: // if not doing any refinement levels, set extent to full box if (levelmin_ == levelmax_) { - x0ref_ = { 0.0, 0.0, 0.0 }; - lxref_ = { 1.0, 1.0, 1.0 }; + x0ref_ = {0.0, 0.0, 0.0}; + lxref_ = {1.0, 1.0, 1.0}; } unsigned ncoarse = 1 << levelmin_; @@ -1327,18 +1343,19 @@ public: if ((levelmin_ != levelmax_) && (!bnoshift || force_shift)) { - int random_base_grid_unit = cf.get_value_safe("random","base_unit",1); - int shift_unit = get_shift_unit( random_base_grid_unit, levelmin_ ); - if( shift_unit != 1 ){ - music::ilog.Print("volume can only be shifted by multiples of %d coarse cells.",shift_unit); - } - xshift_[0] = (int)((0.5-xc[0]) * (double)ncoarse / shift_unit + 0.5) * shift_unit;//ARJ(int)((0.5 - xc[0]) * ncoarse); - xshift_[1] = (int)((0.5-xc[1]) * (double)ncoarse / shift_unit + 0.5) * shift_unit;//ARJ(int)((0.5 - xc[1]) * ncoarse); - xshift_[2] = (int)((0.5-xc[2]) * (double)ncoarse / shift_unit + 0.5) * shift_unit;//ARJ(int)((0.5 - xc[2]) * ncoarse); + int random_base_grid_unit = cf.get_value_safe("random", "base_unit", 1); + int shift_unit = get_shift_unit(random_base_grid_unit, levelmin_); + if (shift_unit != 1) + { + music::ilog.Print("volume can only be shifted by multiples of %d coarse cells.", shift_unit); + } + xshift_[0] = (int)((0.5 - xc[0]) * (double)ncoarse / shift_unit + 0.5) * shift_unit; // ARJ(int)((0.5 - xc[0]) * ncoarse); + xshift_[1] = (int)((0.5 - xc[1]) * (double)ncoarse / shift_unit + 0.5) * shift_unit; // ARJ(int)((0.5 - xc[1]) * ncoarse); + xshift_[2] = (int)((0.5 - xc[2]) * (double)ncoarse / shift_unit + 0.5) * shift_unit; // ARJ(int)((0.5 - xc[2]) * ncoarse); // xshift_[0] = (int)((0.5 - xc[0]) * ncoarse); - // xshift_[1] = (int)((0.5 - xc[1]) * ncoarse); - // xshift_[2] = (int)((0.5 - xc[2]) * ncoarse); + // xshift_[1] = (int)((0.5 - xc[1]) * ncoarse); + // xshift_[2] = (int)((0.5 - xc[2]) * ncoarse); } else { @@ -1458,15 +1475,15 @@ public: else { //... require alignment with coarser grid - music::ilog.Print("- Internal refinement bounding box: [%d,%d]x[%d,%d]x[%d,%d]", il, ir, jl, jr, kl, kr); + music::ilog.Print("- Internal refinement bounding box: [%d,%d]x[%d,%d]x[%d,%d]", il, ir, jl, jr, kl, kr); - il -= il % gridding_unit_; - jl -= jl % gridding_unit_; - kl -= kl % gridding_unit_; + il -= il % gridding_unit_; + jl -= jl % gridding_unit_; + kl -= kl % gridding_unit_; - ir = ((ir%gridding_unit_)!=0)? (ir/gridding_unit_ + 1)*gridding_unit_ : ir; - jr = ((jr%gridding_unit_)!=0)? (jr/gridding_unit_ + 1)*gridding_unit_ : jr; - kr = ((kr%gridding_unit_)!=0)? (kr/gridding_unit_ + 1)*gridding_unit_ : kr; + ir = ((ir % gridding_unit_) != 0) ? (ir / gridding_unit_ + 1) * gridding_unit_ : ir; + jr = ((jr % gridding_unit_) != 0) ? (jr / gridding_unit_ + 1) * gridding_unit_ : jr; + kr = ((kr % gridding_unit_) != 0) ? (kr / gridding_unit_ + 1) * gridding_unit_ : kr; } // if doing unigrid, set region to whole box @@ -1500,8 +1517,7 @@ public: { absoffsets_[levelmax_] = {(il + nresmax) % nresmax, (jl + nresmax) % nresmax, (kl + nresmax) % nresmax}; - len_[levelmax_] = { ir-il, jr-jl, kr-kl }; - + len_[levelmax_] = {ir - il, jr - jl, kr - kl}; if (equal_extent_) { @@ -1512,16 +1528,16 @@ public: throw std::runtime_error("Specified equal_extent=yes conflicting with ref_dims which are not equal."); } size_t ilevel = levelmax_; - index_t nmax = *std::max_element(len_[ilevel].begin(),len_[ilevel].end()); + index_t nmax = *std::max_element(len_[ilevel].begin(), len_[ilevel].end()); - index3_t dx = {0,0,0}; - for( int idim=0; idim<3; ++idim ) + index3_t dx = {0, 0, 0}; + for (int idim = 0; idim < 3; ++idim) { dx[idim] = (index_t)((double)(nmax - len_[ilevel][idim]) * 0.5); absoffsets_[ilevel][idim] -= dx[idim]; len_[ilevel][idim] = nmax; } - + il = absoffsets_[ilevel][0]; jl = absoffsets_[ilevel][1]; kl = absoffsets_[ilevel][2]; @@ -1575,12 +1591,12 @@ public: { //... require alignment with coarser grid il -= il % gridding_unit_; - jl -= jl % gridding_unit_; - kl -= kl % gridding_unit_; + jl -= jl % gridding_unit_; + kl -= kl % gridding_unit_; - ir = ((ir%gridding_unit_)!=0)? (ir/gridding_unit_ + 1)*gridding_unit_ : ir; - jr = ((jr%gridding_unit_)!=0)? (jr/gridding_unit_ + 1)*gridding_unit_ : jr; - kr = ((kr%gridding_unit_)!=0)? (kr/gridding_unit_ + 1)*gridding_unit_ : kr; + ir = ((ir % gridding_unit_) != 0) ? (ir / gridding_unit_ + 1) * gridding_unit_ : ir; + jr = ((jr % gridding_unit_) != 0) ? (jr / gridding_unit_ + 1) * gridding_unit_ : jr; + kr = ((kr % gridding_unit_) != 0) ? (kr / gridding_unit_ + 1) * gridding_unit_ : kr; } if (il >= ir || jl >= jr || kl >= kr || il < 0 || jl < 0 || kl < 0) @@ -1588,19 +1604,19 @@ public: music::elog.Print("Internal refinement bounding box error: [%d,%d]x[%d,%d]x[%d,%d], level=%d", il, ir, jl, jr, kl, kr, ilevel); throw std::runtime_error("refinement_hierarchy: Internal refinement bounding box error 2"); } - absoffsets_[ilevel] = { il, jl, kl }; - len_[ilevel] = { ir-il, jr-jl, kr-kl }; + absoffsets_[ilevel] = {il, jl, kl}; + len_[ilevel] = {ir - il, jr - jl, kr - kl}; if (blocking_factor_) { - for( int idim=0; idim<3; ++idim ) - len_[ilevel][idim] += len_[ilevel][idim] % blocking_factor_; + for (int idim = 0; idim < 3; ++idim) + len_[ilevel][idim] += len_[ilevel][idim] % blocking_factor_; } if (equal_extent_) { - index_t nmax = *std::max_element(len_[ilevel].begin(),len_[ilevel].end()); - for( int idim=0; idim<3; ++idim ) + index_t nmax = *std::max_element(len_[ilevel].begin(), len_[ilevel].end()); + for (int idim = 0; idim < 3; ++idim) { index_t dx = (int)((double)(nmax - len_[ilevel][idim]) * 0.5); absoffsets_[ilevel][idim] -= dx; @@ -1618,17 +1634,17 @@ public: //... determine relative offsets between grids for (unsigned ilevel = levelmax_; ilevel > levelmin_; --ilevel) - for(int idim=0; idim<3; ++idim ) + for (int idim = 0; idim < 3; ++idim) { - offsets_[ilevel][idim] = (absoffsets_[ilevel][idim]/2 - absoffsets_[ilevel-1][idim]); + offsets_[ilevel][idim] = (absoffsets_[ilevel][idim] / 2 - absoffsets_[ilevel - 1][idim]); } //... do a forward sweep to ensure that absolute offsets are also correct now for (unsigned ilevel = levelmin_ + 1; ilevel <= levelmax_; ++ilevel) - for(int idim=0; idim<3; ++idim ) - { - absoffsets_[ilevel][idim] = 2 * absoffsets_[ilevel-1][idim] + 2 * offsets_[ilevel][idim]; - } + for (int idim = 0; idim < 3; ++idim) + { + absoffsets_[ilevel][idim] = 2 * absoffsets_[ilevel - 1][idim] + 2 * offsets_[ilevel][idim]; + } for (unsigned ilevel = levelmin_ + 1; ilevel <= levelmax_; ++ilevel) { @@ -1751,7 +1767,7 @@ public: { unsigned n = 1 << i; - if ( absoffsets_[i] == index3_t({0,0,0}) && len_[i]== index3_t({n,n,n}) ) + if (absoffsets_[i] == index3_t({0, 0, 0}) && len_[i] == index3_t({n, n, n})) { levelmin_ = i; } @@ -1816,7 +1832,7 @@ public: if (xshift_[0] != 0 || xshift_[1] != 0 || xshift_[2] != 0) music::ilog << " - Domain will be shifted by (" << xshift_[0] << ", " << xshift_[1] << ", " << xshift_[2] << ")\n" - << std::endl; + << std::endl; music::ilog << " - Grid structure:\n"; @@ -1844,5 +1860,3 @@ public: typedef GridHierarchy grid_hierarchy; typedef MeshvarBnd meshvar_bnd; typedef Meshvar meshvar; - -#endif diff --git a/src/mg_solver.hh b/src/mg_solver.hh index 33a237f..89e4083 100644 --- a/src/mg_solver.hh +++ b/src/mg_solver.hh @@ -192,7 +192,7 @@ void solver::GaussSeidel(real_t h, MeshvarBnd *u, const Meshvar h2 = h * h; for (int color = 0; color < 2; ++color) -#pragma omp parallel for + #pragma omp parallel for collapse(3) for (int ix = 0; ix < nx; ++ix) for (int iy = 0; iy < ny; ++iy) for (int iz = 0; iz < nz; ++iz) @@ -279,7 +279,7 @@ void solver::twoGrid(unsigned ilevel) ozp = uf->offset(2); meshvar_bnd tLu(*uc, false); -#pragma omp parallel for + #pragma omp parallel for for (int ix = 0; ix < nx / 2; ++ix) { int iix = 2 * ix; @@ -297,7 +297,7 @@ void solver::twoGrid(unsigned ilevel) oj = ff->offset(1); ok = ff->offset(2); -#pragma omp parallel for +#pragma omp parallel for collapse(3) for (int ix = oi; ix < oi + (int)ff->size(0) / 2; ++ix) for (int iy = oj; iy < oj + (int)ff->size(1) / 2; ++iy) for (int iz = ok; iz < ok + (int)ff->size(2) / 2; ++iz) @@ -319,7 +319,7 @@ void solver::twoGrid(unsigned ilevel) meshvar_bnd cc(*uc, false); //... compute correction on coarse grid -#pragma omp parallel for +#pragma omp parallel for collapse(3) for (int ix = 0; ix < (int)cc.size(0); ++ix) for (int iy = 0; iy < (int)cc.size(1); ++iy) for (int iz = 0; iz < (int)cc.size(2); ++iz) @@ -372,8 +372,7 @@ double solver::compute_error(const MeshvarBnd &u, const Meshvar double h = 1.0 / (1ul << ilevel), h2 = h * h; -#pragma omp parallel for reduction(+ \ - : err, count) +#pragma omp parallel for reduction(+ : err, count) collapse(3) for (int ix = 0; ix < nx; ++ix) for (int iy = 0; iy < ny; ++iy) for (int iz = 0; iz < nz; ++iz) @@ -410,8 +409,7 @@ double solver::compute_error(const GridHierarchy &uh, const Gri double h = 1.0 / (1ul << ilevel), h2 = h * h; -#pragma omp parallel for reduction(+ \ - : err, count) + #pragma omp parallel for reduction(+ : err, count, mean_res) collapse(3) for (int ix = 0; ix < nx; ++ix) for (int iy = 0; iy < ny; ++iy) for (int iz = 0; iz < nz; ++iz) @@ -461,8 +459,7 @@ double solver::compute_RMS_resid(const GridHierarchy &uh, const double sum = 0.0, sumd2 = 0.0; size_t count = 0; -#pragma omp parallel for reduction(+ \ - : sum, sumd2, count) + #pragma omp parallel for reduction(+ : sum, sumd2, count) for (int ix = 0; ix < nx; ++ix) for (int iy = 0; iy < ny; ++iy) for (int iz = 0; iz < nz; ++iz) diff --git a/src/plugins/random_music.cc b/src/plugins/random_music.cc index 39a71bd..d935258 100644 --- a/src/plugins/random_music.cc +++ b/src/plugins/random_music.cc @@ -531,22 +531,16 @@ void RNG_music::fill_grid(int ilevel, DensityGrid &A) if (nx == (int)A.size(0) + 2 * margin[0] && ny == (int)A.size(1) + 2 * margin[1] && nz == (int)A.size(2) + 2 * margin[2]) { - int ox = margin[0], oy = margin[1], oz = margin[2]; std::vector slice(ny * nz, 0.0); for (int i = 0; i < nx; ++i) { ifs.read(reinterpret_cast(&slice[0]), ny * nz * sizeof(real_t)); - if (i < ox) - continue; - if (i >= 3 * ox) - break; - #pragma omp parallel for - for (int j = oy; j < 3 * oy; ++j) - for (int k = oz; k < 3 * oz; ++k) - A(i - ox, j - oy, k - oz) = slice[j * nz + k]; + for (int j = 0; j < ny; ++j) + for (int k = 0; k < nz; ++k) + A(i, j, k) = slice[j * nz + k]; } ifs.close(); @@ -577,7 +571,7 @@ void RNG_music::fill_grid(int ilevel, DensityGrid &A) } else { - music::ulog.Print("Copying white noise from memory cache..."); + music::ilog.Print("Copying white noise from memory cache..."); if (mem_cache_[ilevel - levelmin_] == NULL) music::elog.Print("Tried to access mem-cached random numbers for level %d. But these are not available!\n", ilevel); diff --git a/src/plugins/random_music_wnoise_generator.cc b/src/plugins/random_music_wnoise_generator.cc index 0e1a7b7..6d7f405 100644 --- a/src/plugins/random_music_wnoise_generator.cc +++ b/src/plugins/random_music_wnoise_generator.cc @@ -627,16 +627,18 @@ music_wnoise_generator::music_wnoise_generator(music_wnoise_generator &rc, val *= val_phas * sqrt8; - if (i != (int)nxc / 2 && j != (int)nyc / 2 && k != (int)nzc / 2) - { - RE(cfine[qf]) = val.real(); - IM(cfine[qf]) = val.imag(); - } - else - { - // RE(cfine[qf]) = val.real(); - // IM(cfine[qf]) = 0.0; - } + // if (x0_ == NULL || lx_ == NULL){ + if (i != (int)nxc / 2 && j != (int)nyc / 2 && k != (int)nzc / 2) + { + RE(cfine[qf]) = val.real(); + IM(cfine[qf]) = val.imag(); + } + else + { + RE(cfine[qf]) = val.real(); + IM(cfine[qf]) = 0.0; + } + // } } delete[] rcoarse; diff --git a/src/plugins/random_panphasia.cc b/src/plugins/random_panphasia.cc index 84dfd82..3a2c52a 100644 --- a/src/plugins/random_panphasia.cc +++ b/src/plugins/random_panphasia.cc @@ -223,10 +223,9 @@ public: incongruent_fields_ = true; ngrid_ = 2 * ratio * pdescriptor_->i_base; grid_rescale_fac_ = (double)ngrid_ / (1 << levelmin_); - music::ilog.Print("PANPHASIA: will use a higher resolution:\n" - " (%d -> %d) * 2**ref compatible with PANPHASIA\n" - " will Fourier interpolate after.", - grid_m_, grid_p_); + music::ilog << "PANPHASIA: will use a higher resolution:" << std::endl + << " (" << grid_m_ << " -> " << grid_p_ + << ") * 2**ref compatible with PANPHASIA (will Fourier interpolate after)" << std::endl; } } @@ -790,7 +789,9 @@ void RNG_panphasia::fill_grid(int level, DensityGrid &R) delete[] pr3; delete[] pr4; - music::ulog.Print("Copying random field data %d,%d,%d -> %d,%d,%d", nxremap[0], nxremap[1], nxremap[2], nx[0], nx[1], nx[2]); + music::ilog.Print("Copying random field data %d,%d,%d -> %d,%d,%d", nxremap[0], nxremap[1], nxremap[2], nx[0], nx[1], nx[2]); + music::ilog.Print("iexpand_levt = %d,%d,%d", iexpand_left[0], iexpand_left[1], iexpand_left[2]); + // n = 1< &R) music::ulog.Print("done with PANPHASIA for level %d:\n mean=%g, var=%g", level, sum, sum2); music::ulog.Print("Copying into R array: nx[0],nx[1],nx[2] %d %d %d \n", nx[0], nx[1], nx[2]); - music::ilog.Print("PANPHASIA level %d mean and variance are\n

= %g | var(p) = %g", level, sum, sum2); + music::ilog << "PANPHASIA level " << level << " mean and variance are" << std::endl + << "

= "<< sum << " | var(p) = " << sum2 << std::endl; } namespace