From 19e5faa49a5cb67a0b30519e9905eb9a6af58091 Mon Sep 17 00:00:00 2001
From: Chichi Lalescu <clalesc1@jhu.edu>
Date: Thu, 2 Jul 2015 10:39:51 +0200
Subject: [PATCH] still broken.
still debugging. at the moment, fields are set to zero a bunch of times,
but it doesn't seem to help
---
src/fluid_solver.cpp | 15 ++++++++++-----
src/fluid_solver_base.cpp | 18 ++++++++++++------
test.py | 1 +
3 files changed, 23 insertions(+), 11 deletions(-)
diff --git a/src/fluid_solver.cpp b/src/fluid_solver.cpp
index b721a6c9..a0dbff11 100644
--- a/src/fluid_solver.cpp
+++ b/src/fluid_solver.cpp
@@ -262,11 +262,15 @@ void fluid_solver<R>::omega_nonlin( \
{ \
assert(src >= 0 && src < 3); \
/* compute velocity */ \
- /*std::fill_n((R*)this->cu, 2*this->cd->local_size, 0);*/ \
+ std::fill_n((R*)this->cu, 2*this->cd->local_size, 0); \
+ this->low_pass_Fourier(this->cv[src], 3, this->kM); \
this->compute_velocity(this->cv[src]); \
+ this->low_pass_Fourier(this->cu, 3, this->kM); \
/* get fields from Fourier space to real space */ \
- FFTW(execute)(*((FFTW(plan)*)this->c2r_velocity )); \
- FFTW(execute)(*((FFTW(plan)*)this->vc2r[src])); \
+ /*std::fill_n(this->ru, 2*this->cd->local_size, 0); */ \
+ /*std::fill_n(this->rv[src], 2*this->cd->local_size, 0);*/ \
+ FFTW(execute)(*((FFTW(plan)*)this->c2r_velocity )); \
+ FFTW(execute)(*((FFTW(plan)*)this->vc2r[src])); \
/* compute cross product $u \times \omega$, and normalize */ \
R tmpx0, tmpy0, tmpz0; \
RLOOP ( \
@@ -279,9 +283,10 @@ void fluid_solver<R>::omega_nonlin( \
); \
/* go back to Fourier space */ \
/* TODO: is 0 padding needed here? */ \
+ std::fill_n((R*)this->cu, 2*this->cd->local_size, 0); \
FFTW(execute)(*((FFTW(plan)*)this->r2c_velocity )); \
this->low_pass_Fourier(this->cu, 3, this->kM); \
- this->symmetrize(this->cu, 3); \
+ /*this->symmetrize(this->cu, 3);*/ \
/* $\imath k \times Fourier(u \times \omega)$ */ \
R tmpx1, tmpy1, tmpz1; \
CLOOP( \
@@ -298,7 +303,7 @@ void fluid_solver<R>::omega_nonlin( \
this->cu[3*cindex+1][1] = tmpy1 / this->normalization_factor;\
this->cu[3*cindex+2][1] = tmpz1 / this->normalization_factor;\
); \
- this->symmetrize(this->cu, 3); \
+ /*this->symmetrize(this->cu, 3);*/ \
this->add_forcing(this->cu, 1.0); \
} \
\
diff --git a/src/fluid_solver_base.cpp b/src/fluid_solver_base.cpp
index 2c6c329d..e552451d 100644
--- a/src/fluid_solver_base.cpp
+++ b/src/fluid_solver_base.cpp
@@ -245,6 +245,7 @@ void fluid_solver_base<R>::symmetrize(C *data, const int howmany) \
C *buffer; \
buffer = FFTW(alloc_complex)(howmany*this->cd->sizes[1]); \
ptrdiff_t yy; \
+ ptrdiff_t tindex; \
int ranksrc, rankdst; \
for (yy = 1; yy < this->cd->sizes[0]/2; yy++) { \
ranksrc = this->cd->rank[yy]; \
@@ -271,13 +272,15 @@ void fluid_solver_base<R>::symmetrize(C *data, const int howmany) \
if (this->cd->myrank == rankdst) \
{ \
for (ii = 1; ii < this->cd->sizes[1]; ii++) \
- for (cc = 0; cc < howmany; cc++) { \
+ for (cc = 0; cc < howmany; cc++) \
+ { \
(*((data + howmany*((this->cd->sizes[0] - yy - this->cd->starts[0])*this->cd->sizes[1] + ii)*this->cd->sizes[2]) + cc))[0] = \
(*(buffer + howmany*(this->cd->sizes[1]-ii)+cc))[0]; \
(*((data + howmany*((this->cd->sizes[0] - yy - this->cd->starts[0])*this->cd->sizes[1] + ii)*this->cd->sizes[2]) + cc))[1] = \
-(*(buffer + howmany*(this->cd->sizes[1]-ii)+cc))[1]; \
} \
- for (cc = 0; cc < howmany; cc++) { \
+ for (cc = 0; cc < howmany; cc++) \
+ { \
(*((data + cc + howmany*(this->cd->sizes[0] - yy - this->cd->starts[0])*this->cd->sizes[1]*this->cd->sizes[2])))[0] = (*(buffer + cc))[0]; \
(*((data + cc + howmany*(this->cd->sizes[0] - yy - this->cd->starts[0])*this->cd->sizes[1]*this->cd->sizes[2])))[1] = -(*(buffer + cc))[1]; \
} \
@@ -288,12 +291,15 @@ void fluid_solver_base<R>::symmetrize(C *data, const int howmany) \
/* put asymmetric data to 0 */\
/*if (this->cd->myrank == this->cd->rank[this->cd->sizes[0]/2]) \
{ \
- for (cc = 0; cc < howmany; cc++) \
+ tindex = howmany*(this->cd->sizes[0]/2 - this->cd->starts[0])*this->cd->sizes[1]*this->cd->sizes[2]; \
+ for (ii = 0; ii < this->cd->sizes[1]; ii++) \
{ \
- data[cc + howmany*(this->cd->sizes[0]/2 - this->cd->starts[0])*this->cd->sizes[1]*this->cd->sizes[2]][0] = 0.0; \
- data[cc + howmany*(this->cd->sizes[0]/2 - this->cd->starts[0])*this->cd->sizes[1]*this->cd->sizes[2]][1] = 0.0; \
+ std::fill_n((R*)(data + tindex), howmany*2*this->cd->sizes[2], 0.0); \
+ tindex += howmany*this->cd->sizes[2]; \
} \
- }*/ \
+ } \
+ tindex = howmany*(); \
+ std::fill_n((R*)(data + tindex), howmany*2, 0.0);*/ \
} \
\
template<> \
diff --git a/test.py b/test.py
index ba3a3ae7..67888784 100755
--- a/test.py
+++ b/test.py
@@ -301,5 +301,6 @@ if __name__ == '__main__':
a.plot(ycoord, amp*np.sin(ycoord), dashes = (2, 2))
a.plot(ycoord, vfin[0, :, 0, 2])
a.plot(ycoord, -np.cos(ycoord), dashes = (2, 2))
+ a.grid()
fig.savefig('ux_vs_y.pdf', format = 'pdf')
--
GitLab