diff --git a/src/fluid_solver.cpp b/src/fluid_solver.cpp index fbfe6ce5e245b0a211b1e2aabe41d1b75509085b..a8991ec173afb3e086fc778ea6a1873bf3e68c4b 100644 --- a/src/fluid_solver.cpp +++ b/src/fluid_solver.cpp @@ -25,6 +25,25 @@ #include "fftw_tools.hpp" +template <class rnumber> +void fluid_solver<rnumber>::impose_zero_modes() +{ + this->cu[0][0] = 0.0; + this->cu[1][0] = 0.0; + this->cu[2][0] = 0.0; + + this->cv[0][0][1] = 0.0; + this->cv[0][1][1] = 0.0; + this->cv[0][2][1] = 0.0; + + this->cv[1][0][1] = 0.0; + this->cv[1][1][1] = 0.0; + this->cv[1][2][1] = 0.0; + + this->cv[2][0][1] = 0.0; + this->cv[2][1][1] = 0.0; + this->cv[2][2][1] = 0.0; +} /*****************************************************************************/ /* macro for specializations to numeric types compatible with FFTW */ @@ -169,6 +188,7 @@ void fluid_solver<R>::omega_nonlin( \ this->cu[cindex*3+1][1] = (this->kz[zindex]*this->cv[src][cindex*3+0][0] - this->kx[xindex]*this->cv[src][cindex*3+2][0]) / k2; \ this->cu[cindex*3+2][1] = (this->kx[xindex]*this->cv[src][cindex*3+1][0] - this->ky[yindex]*this->cv[src][cindex*3+0][0]) / k2; \ ); \ + this->impose_zero_modes(); \ /* get fields from Fourier space to real space */ \ FFTW(execute)(*((FFTW(plan)*)this->c2r_velocity )); \ FFTW(execute)(*((FFTW(plan)*)this->vc2r[src])); \ @@ -187,12 +207,12 @@ void fluid_solver<R>::omega_nonlin( \ /* $\imath k \times DFT(u \times \omega)$ */ \ R tmpx1, tmpy1, tmpz1; \ CLOOP( \ - tmpx0 = (this->ky[yindex]*this->cu[cindex*3+2][1] - this->kz[zindex]*this->cu[cindex*3+1][1]); \ - tmpy0 = (this->kz[zindex]*this->cu[cindex*3+0][1] - this->kx[xindex]*this->cu[cindex*3+2][1]); \ - tmpz0 = (this->kx[xindex]*this->cu[cindex*3+1][1] - this->ky[yindex]*this->cu[cindex*3+0][1]); \ - tmpx1 = (this->ky[yindex]*this->cu[cindex*3+2][0] - this->kz[zindex]*this->cu[cindex*3+1][0]); \ - tmpy1 = (this->kz[zindex]*this->cu[cindex*3+0][0] - this->kx[xindex]*this->cu[cindex*3+2][0]); \ - tmpz1 = (this->kx[xindex]*this->cu[cindex*3+1][0] - this->ky[yindex]*this->cu[cindex*3+0][0]); \ + tmpx0 = -(this->ky[yindex]*this->cu[cindex*3+2][1] - this->kz[zindex]*this->cu[cindex*3+1][1]); \ + tmpy0 = -(this->kz[zindex]*this->cu[cindex*3+0][1] - this->kx[xindex]*this->cu[cindex*3+2][1]); \ + tmpz0 = -(this->kx[xindex]*this->cu[cindex*3+1][1] - this->ky[yindex]*this->cu[cindex*3+0][1]); \ + tmpx1 = (this->ky[yindex]*this->cu[cindex*3+2][0] - this->kz[zindex]*this->cu[cindex*3+1][0]); \ + tmpy1 = (this->kz[zindex]*this->cu[cindex*3+0][0] - this->kx[xindex]*this->cu[cindex*3+2][0]); \ + tmpz1 = (this->kx[xindex]*this->cu[cindex*3+1][0] - this->ky[yindex]*this->cu[cindex*3+0][0]); \ this->cu[cindex*3+0][0] = tmpx0;\ this->cu[cindex*3+1][0] = tmpy0;\ this->cu[cindex*3+2][0] = tmpz0;\ diff --git a/src/fluid_solver.hpp b/src/fluid_solver.hpp index 12f3a824dfa4b18cccb945670d6b7447bfd84e1e..a26db8ada8c984b871d79c5c929bc732e1661d42 100644 --- a/src/fluid_solver.hpp +++ b/src/fluid_solver.hpp @@ -74,6 +74,7 @@ class fluid_solver:public fluid_solver_base<rnumber> void omega_nonlin(int src); void step(double dt); + void impose_zero_modes(); }; #endif//FLUID_SOLVER diff --git a/src/fluid_solver_base.cpp b/src/fluid_solver_base.cpp index 645aa67591e25a1936d19d434804f193b66b05d9..7ab4a730993e7e1fe294e92aab3788d36e17ac7a 100644 --- a/src/fluid_solver_base.cpp +++ b/src/fluid_solver_base.cpp @@ -106,6 +106,13 @@ fluid_solver_base<R>::fluid_solver_base( \ else \ this->knullz[i] = true; \ } \ + this->kM = this->kMx; \ + if (this->kM < this->kMy) this->kM = this->kMy; \ + if (this->kM < this->kMz) this->kM = this->kMz; \ + this->kM2 = this->kM * this->kM; \ + this->dk = this->dkx; \ + if (this->dk > this->dky) this->dk = this->dky; \ + if (this->dk > this->dkz) this->dk = this->dkz; \ } \ \ template<> \ @@ -133,14 +140,18 @@ R fluid_solver_base<R>::correl_vec(C *a, C *b) \ k2 = (this->kx[xindex]*this->kx[xindex] + \ this->ky[yindex]*this->ky[yindex] + \ this->kz[zindex]*this->kz[zindex]); \ - factor = (xindex == 0) ? 1 : 2; \ - val_process += factor * ((*(a + cindex))[0] * (*(b + cindex))[0] + \ - (*(a + cindex))[1] * (*(b + cindex))[1]); \ + if (k2 < this->kM2) \ + { \ + factor = (xindex == 0) ? 1 : 2; \ + val_process += factor * ((*(a + cindex))[0] * (*(b + cindex))[0] + \ + (*(a + cindex))[1] * (*(b + cindex))[1]); \ + } \ );\ MPI_Allreduce( \ (void*)(&val_process), \ (void*)(&val), \ 1, MPI_DOUBLE_PRECISION, MPI_SUM, MPI_COMM_WORLD); \ + /*return R(val / (this->rd->full_size / 3)) / (this->rd->full_size / 3);*/ \ return R(val); \ }