From 901d6bfbe6b282a57477723745b0db729a3d8af4 Mon Sep 17 00:00:00 2001
From: Cristian C Lalescu <Cristian.Lalescu@ds.mpg.de>
Date: Thu, 2 Jul 2015 15:35:14 +0200
Subject: [PATCH] fix nonlinear term computation
I was normalizing many many times because of the previous error, so now
everything seems to make sense, and I get collapse of energy and
enstrophy with those given by vortex
---
src/fluid_solver.cpp | 19 +++++++++----------
test.py | 10 ++++++----
2 files changed, 15 insertions(+), 14 deletions(-)
diff --git a/src/fluid_solver.cpp b/src/fluid_solver.cpp
index 5bb2e356..26ec1305 100644
--- a/src/fluid_solver.cpp
+++ b/src/fluid_solver.cpp
@@ -271,15 +271,14 @@ void fluid_solver<R>::omega_nonlin( \
/* compute cross product $u \times \omega$, and normalize */ \
R tmpx0, tmpy0, tmpz0; \
RLOOP ( \
- tmpx0 = (this->ru[(3*rindex)+1]*this->rv[src][(3*rindex)+2] - this->ru[(3*rindex)+2]*this->rv[src][(3*rindex)+1]) / this->normalization_factor; \
- tmpy0 = (this->ru[(3*rindex)+2]*this->rv[src][(3*rindex)+0] - this->ru[(3*rindex)+0]*this->rv[src][(3*rindex)+2]) / this->normalization_factor; \
- tmpz0 = (this->ru[(3*rindex)+0]*this->rv[src][(3*rindex)+1] - this->ru[(3*rindex)+1]*this->rv[src][(3*rindex)+0]) / this->normalization_factor; \
+ tmpx0 = (this->ru[(3*rindex)+1]*this->rv[src][(3*rindex)+2] - this->ru[(3*rindex)+2]*this->rv[src][(3*rindex)+1]); \
+ tmpy0 = (this->ru[(3*rindex)+2]*this->rv[src][(3*rindex)+0] - this->ru[(3*rindex)+0]*this->rv[src][(3*rindex)+2]); \
+ tmpz0 = (this->ru[(3*rindex)+0]*this->rv[src][(3*rindex)+1] - this->ru[(3*rindex)+1]*this->rv[src][(3*rindex)+0]); \
this->ru[(3*rindex)+0] = tmpx0 / this->normalization_factor; \
this->ru[(3*rindex)+1] = tmpy0 / this->normalization_factor; \
this->ru[(3*rindex)+2] = tmpz0 / this->normalization_factor; \
); \
/* go back to Fourier space */ \
- /* TODO: is 0 padding needed here? */ \
FFTW(execute)(*((FFTW(plan)*)this->r2c_velocity )); \
this->low_pass_Fourier(this->cu, 3, this->kM); \
this->symmetrize(this->cu, 3); \
@@ -292,12 +291,12 @@ void fluid_solver<R>::omega_nonlin( \
tmpx1 = (this->ky[yindex]*this->cu[3*cindex+2][0] - this->kz[zindex]*this->cu[3*cindex+1][0]); \
tmpy1 = (this->kz[zindex]*this->cu[3*cindex+0][0] - this->kx[xindex]*this->cu[3*cindex+2][0]); \
tmpz1 = (this->kx[xindex]*this->cu[3*cindex+1][0] - this->ky[yindex]*this->cu[3*cindex+0][0]); \
- this->cu[3*cindex+0][0] = tmpx0 / this->normalization_factor;\
- this->cu[3*cindex+1][0] = tmpy0 / this->normalization_factor;\
- this->cu[3*cindex+2][0] = tmpz0 / this->normalization_factor;\
- this->cu[3*cindex+0][1] = tmpx1 / this->normalization_factor;\
- this->cu[3*cindex+1][1] = tmpy1 / this->normalization_factor;\
- this->cu[3*cindex+2][1] = tmpz1 / this->normalization_factor;\
+ this->cu[3*cindex+0][0] = tmpx0; \
+ this->cu[3*cindex+1][0] = tmpy0; \
+ this->cu[3*cindex+2][0] = tmpz0; \
+ this->cu[3*cindex+0][1] = tmpx1; \
+ this->cu[3*cindex+1][1] = tmpy1; \
+ this->cu[3*cindex+2][1] = tmpz1; \
); \
/*this->symmetrize(this->cu, 3);*/ \
this->add_forcing(this->cu, 1.0); \
diff --git a/test.py b/test.py
index c82ec809..0fc059b5 100755
--- a/test.py
+++ b/test.py
@@ -223,15 +223,17 @@ def convergence_test(opt):
stats1 = np.fromfile('test1_stats.bin', dtype = dtype)
stats2 = np.fromfile('test2_stats.bin', dtype = dtype)
stats_vortex = np.loadtxt('../vortex/sim_000000.log')
+ print stats1
+ print stats2
fig = plt.figure(figsize = (12,6))
a = fig.add_subplot(121)
a.plot(stats1['t'], stats1['energy'])
- a.plot(stats2['t'], stats2['energy'])
- a.plot(stats_vortex[:, 2], stats_vortex[:, 3])
+ a.plot(stats2['t'], stats2['energy'], dashes = (3, 3))
+ a.plot(stats_vortex[:, 2], stats_vortex[:, 3], dashes = (2, 4))
a = fig.add_subplot(122)
a.plot(stats1['t'], stats1['enstrophy'])
- a.plot(stats2['t'], stats2['enstrophy'])
- a.plot(stats_vortex[:, 2], stats_vortex[:, 9]/2)
+ a.plot(stats2['t'], stats2['enstrophy'], dashes = (3, 3))
+ a.plot(stats_vortex[:, 2], stats_vortex[:, 9]/2, dashes = (2, 4))
fig.savefig('test.pdf', format = 'pdf')
#fig = plt.figure(figsize=(12, 12))
--
GitLab