From 06b92ff37e166f09522a6904217fcfc4f332620b Mon Sep 17 00:00:00 2001
From: Chichi Lalescu <clalesc1@jhu.edu>
Date: Fri, 19 Jun 2015 21:50:49 +0200
Subject: [PATCH] add impose_zero_modes. fixes NANs
---
src/fluid_solver.cpp | 32 ++++++++++++++++++++++++++------
src/fluid_solver.hpp | 1 +
src/fluid_solver_base.cpp | 17 ++++++++++++++---
3 files changed, 41 insertions(+), 9 deletions(-)
diff --git a/src/fluid_solver.cpp b/src/fluid_solver.cpp
index fbfe6ce5..a8991ec1 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 12f3a824..a26db8ad 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 645aa675..7ab4a730 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); \
}
--
GitLab