broken --- testing with Kolmogorov forcing

src/fluid_solver.cpp

@@ -211,6 +211,23 @@ void fluid_solver<R>::compute_velocity(C *vorticity) \
 } \
  \
 template<> \
+void fluid_solver<R>::add_forcing(\
+        C *field, R factor) \
+{ \
+    ptrdiff_t cindex; \
+    if (this->cd->myrank == this->cd->rank[this->fmode]) \
+    { \
+        cindex = ((this->fmode - this->cd->starts[0]) * this->cd->sizes[1])*this->cd->sizes[2]; \
+        field[cindex*3+2][1] -= this->famplitude*factor/2; \
+    } \
+    if (this->cd->myrank == this->cd->rank[this->cd->sizes[0] - this->fmode]) \
+    { \
+        cindex = ((this->cd->sizes[0] - this->fmode - this->cd->starts[0]) * this->cd->sizes[1])*this->cd->sizes[2]; \
+        field[cindex*3+2][1] += this->famplitude*factor/2; \
+    } \
+} \
+ \
+template<> \
 void fluid_solver<R>::omega_nonlin( \
         int src) \
 { \
@@ -251,8 +268,8 @@ void fluid_solver<R>::omega_nonlin( \
             this->cu[cindex*3+1][1] = tmpy1 / this->normalization_factor;\
             this->cu[cindex*3+2][1] = tmpz1 / this->normalization_factor;\
             ); \
+    this->add_forcing(this->cu, 1.0); \
     this->symmetrize(this->cu, 3); \
-    this->symmetrize(this->cv[src], 3); \
 } \
  \
 template<> \

src/fluid_solver.hpp

@@ -60,7 +60,9 @@ class fluid_solver:public fluid_solver_base<rnumber>
         void *vr2c[3], *vc2r[3];
 
         /* physical parameters */
-        rnumber nu;
+        double nu;
+        int fmode;
+        double famplitude;
 
         /* methods */
         fluid_solver(
@@ -78,6 +80,7 @@ class fluid_solver:public fluid_solver_base<rnumber>
         void omega_nonlin(int src);
         void step(double dt);
         void impose_zero_modes(void);
+        void add_forcing(rnumber (*field)[2], rnumber factor);
 
         int read(char field, char representation);
         int write(char field, char representation);

test.py

@@ -74,6 +74,8 @@ class stat_test(code):
         self.parameters['niter_todo'] = 8
         self.parameters['dt'] = 0.01
         self.parameters['nu'] = 0.1
+        self.parameters['famplitude'] = 1.0
+        self.parameters['fmode'] = 1
         self.includes += '#include <cstring>\n'
         self.includes += '#include "fftw_tools.hpp"\n'
         self.variables += ('double t;\n' +
@@ -110,6 +112,8 @@ class stat_test(code):
                 char fname[512];
                 fs = new fluid_solver<float>(simname, nx, ny, nz);
                 fs->nu = nu;
+                fs->fmode = fmode;
+                fs->famplitude = famplitude;
                 fs->iteration = iter0;
                 if (myrank == fs->cd->io_myrank)
                     {
@@ -143,8 +147,8 @@ 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 = 16)
+    parser.add_argument('--ncpu', type = int, dest = 'ncpu', default = 2)
+    parser.add_argument('--nsteps', type = int, dest = 'nsteps', default = 16)
     parser.add_argument('-n', type = int, dest = 'n', default = 32)
     opt = parser.parse_args()
 
@@ -164,8 +168,8 @@ if __name__ == '__main__':
         c.parameters['nx'] = opt.n
         c.parameters['ny'] = opt.n
         c.parameters['nz'] = opt.n
-        c.parameters['nu'] = 0.1
-        c.parameters['dt'] = 0.004
+        c.parameters['nu'] = 1.0
+        c.parameters['dt'] = 0.05
         c.parameters['niter_todo'] = opt.nsteps
         c.write_par(simname = 'test1')
         Kdata0.tofile("test1_cvorticity_i00000")
@@ -174,7 +178,7 @@ if __name__ == '__main__':
         c.parameters['niter_todo'] *= 2
         c.write_par(simname = 'test2')
         Kdata0.tofile("test2_cvorticity_i00000")
-        c.run(ncpu = opt.ncpu, simname = 'test2')
+        #c.run(ncpu = opt.ncpu, simname = 'test2')
         Rdata = np.fromfile(
                 'test1_rvorticity_i00000',
                 dtype = np.float32).reshape(opt.n, opt.n, opt.n, 3)
@@ -205,7 +209,7 @@ if __name__ == '__main__':
         axis.set_title('{0}'.format(np.average(Rdata0[..., 0]**2 +
                                                Rdata0[..., 1]**2 +
                                                Rdata0[..., 2]**2)*.5))
-        return None
+        return Rdata0
 
     fig = plt.figure(figsize=(12, 12))
     a = fig.add_subplot(221)
@@ -214,8 +218,14 @@ if __name__ == '__main__':
     a.set_axis_off()
     plot_vel_cut('test2_rvelocity_i00000', a)
     a = fig.add_subplot(223)
-    plot_vel_cut('test1_rvelocity_i{0:0>5x}'.format(stats1.shape[0]-1), a)
+    tmp = plot_vel_cut('test1_rvelocity_i{0:0>5x}'.format(stats1.shape[0]-1), a)
     a = fig.add_subplot(224)
     plot_vel_cut('test2_rvelocity_i{0:0>5x}'.format(stats2.shape[0]-1), a)
     fig.savefig('vel_cut.pdf', format = 'pdf')
 
+    fig = plt.figure()
+    a = fig.add_subplot(111)
+    a.plot(tmp[0, :, 0, 0])
+    a.plot(tmp[0, :, 0, 1])
+    a.plot(tmp[0, :, 0, 2])
+    fig.savefig('ux_vs_y.pdf', format = 'pdf')