diff --git a/src/fluid_solver.cpp b/src/fluid_solver.cpp index 5bb2e356cfbd5d0bf52b82d87674e181eff5b501..26ec1305d91b52506b81c9593dc05e03fec934c7 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 c82ec809027f665a89800519cd95b7a4d2930a3b..0fc059b560320f9fa52c8c692172c011b5fdb0d6 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))