diff --git a/src/fftw_tools.cpp b/src/fftw_tools.cpp
index 795982531bfde0e7a908b6547d943f35c5c7b9b5..02f83603e71931afebe0f92679475b13db597142 100644
--- a/src/fftw_tools.cpp
+++ b/src/fftw_tools.cpp
@@ -141,7 +141,7 @@ int clip_zero_padding(
rnumber *a,
int howmany)
{
- if (f->ndims != 3)
+ if (f->ndims < 3)
return EXIT_FAILURE;
rnumber *b = a;
ptrdiff_t copy_size = f->sizes[2] * howmany;
diff --git a/src/fluid_solver.cpp b/src/fluid_solver.cpp
index ad960fdd33b624b40149cf9ba68f957c4a8b95c3..89b3db0a0a330b9108190a50c0cae37acc887e5b 100644
--- a/src/fluid_solver.cpp
+++ b/src/fluid_solver.cpp
@@ -195,12 +195,12 @@ void fluid_solver<R>::compute_velocity(C *vorticity) \
k2 = (this->kx[xindex]*this->kx[xindex] + \
this->ky[yindex]*this->ky[yindex] + \
this->kz[zindex]*this->kz[zindex]); \
- this->cu[cindex*3+0][0] = (this->ky[yindex]*vorticity[cindex*3+2][1] - this->kz[zindex]*vorticity[cindex*3+1][1]) / k2; \
- this->cu[cindex*3+1][0] = (this->kz[zindex]*vorticity[cindex*3+0][1] - this->kx[xindex]*vorticity[cindex*3+2][1]) / k2; \
- this->cu[cindex*3+2][0] = (this->kx[xindex]*vorticity[cindex*3+1][1] - this->ky[yindex]*vorticity[cindex*3+0][1]) / k2; \
- this->cu[cindex*3+0][1] = (this->ky[yindex]*vorticity[cindex*3+2][0] - this->kz[zindex]*vorticity[cindex*3+1][0]) / k2; \
- this->cu[cindex*3+1][1] = (this->kz[zindex]*vorticity[cindex*3+0][0] - this->kx[xindex]*vorticity[cindex*3+2][0]) / k2; \
- this->cu[cindex*3+2][1] = (this->kx[xindex]*vorticity[cindex*3+1][0] - this->ky[yindex]*vorticity[cindex*3+0][0]) / k2; \
+ this->cu[cindex*3+0][0] = -(this->ky[yindex]*vorticity[cindex*3+2][1] - this->kz[zindex]*vorticity[cindex*3+1][1]) / k2; \
+ this->cu[cindex*3+1][0] = -(this->kz[zindex]*vorticity[cindex*3+0][1] - this->kx[xindex]*vorticity[cindex*3+2][1]) / k2; \
+ this->cu[cindex*3+2][0] = -(this->kx[xindex]*vorticity[cindex*3+1][1] - this->ky[yindex]*vorticity[cindex*3+0][1]) / k2; \
+ this->cu[cindex*3+0][1] = (this->ky[yindex]*vorticity[cindex*3+2][0] - this->kz[zindex]*vorticity[cindex*3+1][0]) / k2; \
+ this->cu[cindex*3+1][1] = (this->kz[zindex]*vorticity[cindex*3+0][0] - this->kx[xindex]*vorticity[cindex*3+2][0]) / k2; \
+ this->cu[cindex*3+2][1] = (this->kx[xindex]*vorticity[cindex*3+1][0] - this->ky[yindex]*vorticity[cindex*3+0][0]) / k2; \
); \
this->impose_zero_modes(); \
this->symmetrize(this->cu, 3); \
@@ -227,6 +227,7 @@ void fluid_solver<R>::omega_nonlin( \
this->ru[rindex*3+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); \
@@ -253,6 +254,8 @@ void fluid_solver<R>::omega_nonlin( \
this->cv[src][cindex*3+2][0] /= this->normalization_factor; \
this->cv[src][cindex*3+2][1] /= this->normalization_factor; \
); \
+ this->symmetrize(this->cu, 3); \
+ this->symmetrize(this->cv[src], 3); \
} \
\
template<> \
diff --git a/src/fluid_solver_base.cpp b/src/fluid_solver_base.cpp
index 57e85324e5daa127edde7ee288cde03d1d30b386..b38492b4ed4d94d327c01a9ad9de4e3bc87ff7da 100644
--- a/src/fluid_solver_base.cpp
+++ b/src/fluid_solver_base.cpp
@@ -157,7 +157,7 @@ R fluid_solver_base<R>::correl_vec(C *a, C *b) \
(void*)(&val_process), \
(void*)(&val), \
1, MPI_DOUBLE_PRECISION, MPI_SUM, MPI_COMM_WORLD); \
- return R(val / this->normalization_factor); \
+ return R(val); \
} \
\
template<> \
diff --git a/test.py b/test.py
index 3a3f538e1eeb25afecc0a7053b79c6e4ae2e4510..f5bfe888f477b02156fc737e3e3caddd6c83529a 100755
--- a/test.py
+++ b/test.py
@@ -75,6 +75,7 @@ class stat_test(code):
self.parameters['dt'] = 0.01
self.parameters['nu'] = 0.1
self.includes += '#include <cstring>\n'
+ self.includes += '#include "fftw_tools.hpp"\n'
self.variables += ('double t;\n' +
'FILE *stat_file;\n'
'double stats[2];\n')
@@ -127,9 +128,16 @@ class stat_test(code):
do_stats(fs);
fftwf_execute(*((fftwf_plan*)fs->c2r_velocity));
sprintf(fname, "%s_rvelocity_i%.5x", simname, fs->iteration);
+ clip_zero_padding<float>(fs->rd, fs->rvelocity, 3);
fs->rd->write(
fname,
(void*)fs->rvelocity);
+ fftwf_execute(*((fftwf_plan*)fs->c2r_vorticity));
+ sprintf(fname, "%s_rvorticity_i%.5x", simname, fs->iteration);
+ clip_zero_padding<float>(fs->rd, fs->rvorticity, 3);
+ fs->rd->write(
+ fname,
+ (void*)fs->rvorticity);
for (; fs->iteration < iter0 + niter_todo;)
{
fs->step(dt);
@@ -141,8 +149,15 @@ class stat_test(code):
fs->cd->write(
fname,
(void*)fs->cvorticity);
+ fftwf_execute(*((fftwf_plan*)fs->c2r_vorticity));
+ sprintf(fname, "%s_rvorticity_i%.5x", simname, fs->iteration);
+ clip_zero_padding<float>(fs->rd, fs->rvorticity, 3);
+ fs->rd->write(
+ fname,
+ (void*)fs->rvorticity);
fftwf_execute(*((fftwf_plan*)fs->c2r_velocity));
sprintf(fname, "%s_rvelocity_i%.5x", simname, fs->iteration);
+ clip_zero_padding<float>(fs->rd, fs->rvelocity, 3);
fs->rd->write(
fname,
(void*)fs->rvelocity);
@@ -154,45 +169,80 @@ class stat_test(code):
if __name__ == '__main__':
parser = argparse.ArgumentParser()
parser.add_argument('test_name', type = str)
+ parser.add_argument('--run', dest = 'run', action = 'store_true')
parser.add_argument('--ncpu', dest = 'ncpu', default = 2)
parser.add_argument('--nsteps', dest = 'nsteps', default = 8)
parser.add_argument('-n', type = int, dest = 'n', default = 32)
opt = parser.parse_args()
- np.random.seed(7547)
- Kdata00 = generate_data_3D(opt.n, p = 1.5).astype(np.complex64)
- Kdata01 = generate_data_3D(opt.n, p = 1.5).astype(np.complex64)
- Kdata02 = generate_data_3D(opt.n, p = 1.5).astype(np.complex64)
- Kdata0 = np.zeros(
- Kdata00.shape + (3,),
- Kdata00.dtype)
- Kdata0[..., 0] = Kdata00
- Kdata0[..., 1] = Kdata01
- Kdata0[..., 2] = Kdata02
- c = stat_test(name = opt.test_name)
- c.write_src()
- c.parameters['nx'] = opt.n
- c.parameters['ny'] = opt.n
- c.parameters['nz'] = opt.n
- c.parameters['nu'] = 0.1
- c.parameters['dt'] = 0.01
- c.parameters['niter_todo'] = opt.nsteps
- c.write_par(simname = 'test1')
- Kdata0.tofile("test1_kvorticity_i00000")
- c.run(ncpu = opt.ncpu, simname = 'test1')
- c.parameters['dt'] = 0.005
- c.parameters['niter_todo'] = opt.nsteps*2
- c.write_par(simname = 'test2')
- Kdata0.tofile("test2_kvorticity_i00000")
- c.run(ncpu = opt.ncpu, simname = 'test2')
+ if opt.run:
+ np.random.seed(7547)
+ Kdata00 = generate_data_3D(opt.n, p = 1.).astype(np.complex64)
+ Kdata01 = generate_data_3D(opt.n, p = 1.).astype(np.complex64)
+ Kdata02 = generate_data_3D(opt.n, p = 1.).astype(np.complex64)
+ Kdata0 = np.zeros(
+ Kdata00.shape + (3,),
+ Kdata00.dtype)
+ Kdata0[..., 0] = Kdata00
+ Kdata0[..., 1] = Kdata01
+ Kdata0[..., 2] = Kdata02
+ c = stat_test(name = opt.test_name)
+ c.write_src()
+ c.parameters['nx'] = opt.n
+ c.parameters['ny'] = opt.n
+ c.parameters['nz'] = opt.n
+ c.parameters['nu'] = 0.1
+ c.parameters['dt'] = 0.01
+ c.parameters['niter_todo'] = opt.nsteps
+ c.write_par(simname = 'test1')
+ Kdata0.tofile("test1_kvorticity_i00000")
+ c.run(ncpu = opt.ncpu, simname = 'test1')
+ c.parameters['dt'] = 0.005
+ c.parameters['niter_todo'] = opt.nsteps*2
+ c.write_par(simname = 'test2')
+ Kdata0.tofile("test2_kvorticity_i00000")
+ c.run(ncpu = opt.ncpu, simname = 'test2')
+ Rdata = np.fromfile(
+ 'test1_rvorticity_i00000',
+ dtype = np.float32).reshape(opt.n, opt.n, opt.n, 3)
+ tdata = Rdata.transpose(3, 0, 1, 2).copy()
+ tdata.tofile('input_for_vortex')
dtype = pickle.load(open(opt.test_name + '_dtype.pickle'))
stats1 = np.fromfile('test1_stats.bin', dtype = dtype)
stats2 = np.fromfile('test2_stats.bin', dtype = dtype)
- fig = plt.figure(figsize = (6,6))
- a = fig.add_subplot(111)
- print stats1['energy']
- print stats2['energy']
+ stats_vortex = np.loadtxt('../vortex/sim_000000.log')
+ 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 = 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)
fig.savefig('test.pdf', format = 'pdf')
+ def plot_vel_cut(fname, axis):
+ axis.set_axis_off()
+ Rdata0 = np.fromfile(
+ fname,
+ dtype = np.float32).reshape((opt.n, opt.n, opt.n, 3))
+ energy = np.sum(Rdata0[13, :, :, :]**2, axis = 2)*.5
+ axis.imshow(energy, interpolation='none')
+ axis.set_title('{0}'.format(np.average(Rdata0[..., 0]**2 +
+ Rdata0[..., 1]**2 +
+ Rdata0[..., 2]**2)*.5))
+ return None
+
+ fig = plt.figure(figsize=(12, 12))
+ a = fig.add_subplot(221)
+ plot_vel_cut('test1_rvelocity_i00000', a)
+ a = fig.add_subplot(222)
+ a.set_axis_off()
+ plot_vel_cut('test2_rvelocity_i00000', a)
+ a = fig.add_subplot(223)
+ plot_vel_cut('test1_rvelocity_i00008', a)
+ a = fig.add_subplot(224)
+ plot_vel_cut('test2_rvelocity_i00010', a)
+ fig.savefig('vel_cut.pdf', format = 'pdf')
+