From 5c23348262b6feb2bf1bd93d231925fa7df3a9ee Mon Sep 17 00:00:00 2001
From: Cristian C Lalescu <Cristian.Lalescu@ds.mpg.de>
Date: Thu, 16 Jul 2015 15:57:05 +0200
Subject: [PATCH] fix divfree

---
 bfps/NavierStokes.py      | 63 ++++++++++++++++++++++------------
 bfps/resize.py            |  2 +-
 src/fluid_solver.cpp      | 71 ++++++++++++++-------------------------
 src/fluid_solver_base.cpp | 15 +++++----
 test.py                   |  4 ++-
 5 files changed, 78 insertions(+), 77 deletions(-)

diff --git a/bfps/NavierStokes.py b/bfps/NavierStokes.py
index 1f143114..a9217e65 100644
--- a/bfps/NavierStokes.py
+++ b/bfps/NavierStokes.py
@@ -91,6 +91,7 @@ class NavierStokes(bfps.code):
                         fwrite((void*)stats, sizeof(double), 3, stat_file);
                     }
                     fs->write_spectrum("velocity", fs->cvelocity);
+                    fs->write_spectrum("vorticity", fs->cvorticity);
                 }
                 //endcpp
                 """
@@ -127,14 +128,14 @@ class NavierStokes(bfps.code):
                 fs->iteration = iter0;
                 fs->read('v', 'c');
                 fs->low_pass_Fourier(fs->cvorticity, 3, fs->kM);
-                fs->force_divfree(fs->cvorticity);
+                //fs->force_divfree(fs->cvorticity);
                 fs->symmetrize(fs->cvorticity, 3);
                 if (myrank == 0)
                 {
                     sprintf(fname, "%s_stats.bin", simname);
-                    stat_file = fopen(fname, "wb");
+                    stat_file = fopen(fname, "ab");
                 }
-                t = 0.0;
+                t = dt*iter0;
                 do_stats(fs);
                 //endcpp
                 """
@@ -152,6 +153,7 @@ class NavierStokes(bfps.code):
                     fclose(stat_file);
                 }
                 fs->write('v', 'c');
+                fs->write('u', 'r');
                 delete fs;
                 //endcpp
                 """
@@ -335,10 +337,10 @@ def test(opt):
     c.parameters['nx'] = opt.n
     c.parameters['ny'] = opt.n
     c.parameters['nz'] = opt.n
-    c.parameters['nu'] = 1e-1
+    c.parameters['nu'] = 5.5*opt.n**(-4./3)
     c.parameters['dt'] = 2e-3
     c.parameters['niter_todo'] = opt.nsteps
-    c.parameters['famplitude'] = 1.
+    c.parameters['famplitude'] = 0.0
     c.parameters['nparticles'] = 32
     if opt.particles:
         c.add_particles()
@@ -348,30 +350,47 @@ def test(opt):
     c.write_par(simname = 'test')
     if opt.run:
         c.generate_vector_field(simname = 'test')
-        c.generate_tracer_state(simname = 'test', species = 0)
-        c.generate_tracer_state(simname = 'test', species = 1)
+        if opt.particles:
+            c.generate_tracer_state(simname = 'test', species = 0)
+            c.generate_tracer_state(simname = 'test', species = 1)
         c.run(ncpu = opt.ncpu,
               simname = 'test')
     stats = c.read_stats()
-    k, espec = c.read_spec()
+    k, enespec = c.read_spec()
+    k, ensspec = c.read_spec(field = 'vorticity')
+
+    # plot energy and enstrophy
+    fig = plt.figure(figsize = (12, 6))
+    a = fig.add_subplot(121)
+    etaK = (c.parameters['nu']**2 / (stats['enstrophy']*2))**.25
+    a.plot(stats['t'], k[-3]*etaK, label = '$k_M \eta_K$')
+    a.plot(stats['t'], c.parameters['dt']*stats['vel_max'] / (2*np.pi/c.parameters['nx']),
+            label = '$\\frac{\\Delta t \\| u \\|_\infty}{\\Delta x}$')
+    a.legend(loc = 'best')
+    a = fig.add_subplot(122)
+    a.plot(stats['t'], stats['energy'], label = 'energy', color = (0, 0, 1))
+    a.set_ylabel('energy')
+    a.set_xlabel('$t$')
+    b = a.twinx()
+    b.plot(stats['t'], stats['enstrophy'], label = 'enstrophy', color = (1, 0, 0))
+    b.set_ylabel('enstrophy')
+    fig.savefig('stats.pdf', format = 'pdf')
 
     # plot spectra
-    fig = plt.figure(figsize=(6,6))
-    a = fig.add_subplot(111)
-    for i in range(espec.shape[0]):
-        a.plot(k, espec[i]['val'])
+    fig = plt.figure(figsize=(12,6))
+    a = fig.add_subplot(121)
+    for i in range(0, enespec.shape[0]):
+        a.plot(k, enespec[i]['val'], color = (i*1./enespec.shape[0], 0, 1 - i*1./enespec.shape[0]))
+    a.set_xscale('log')
+    a.set_yscale('log')
+    a.set_title('velocity')
+    a = fig.add_subplot(122)
+    for i in range(0, ensspec.shape[0]):
+        a.plot(k, ensspec[i]['val'], color = (i*1./ensspec.shape[0], 0, 1 - i*1./ensspec.shape[0]))
+        a.plot(k, k**2*enespec[i]['val'], color = (0, 1 - i*1./ensspec.shape[0], i*1./ensspec.shape[0]))
     a.set_xscale('log')
     a.set_yscale('log')
+    a.set_title('vorticity')
     fig.savefig('spectrum.pdf', format = 'pdf')
-
-    #plot consistency checks
-    fig = plt.figure(figsize=(6,6))
-    a = fig.add_subplot(111)
-    etaK = (c.parameters['nu']**2 / (stats['enstrophy']*2))**.25
-    a.plot(k[-3]*etaK, label = '$k_M \eta_K$')
-    a.plot(c.parameters['dt']*stats['vel_max'] / (2*np.pi/c.parameters['nx']),
-            label = '$\\frac{\\Delta t \\| u \\|_\infty}{\\Delta x}$')
-    a.legend(loc = 'best')
-    fig.savefig('consistency_checks.pdf', format = 'pdf')
     return None
 
diff --git a/bfps/resize.py b/bfps/resize.py
index 3e9a267f..173c3d96 100644
--- a/bfps/resize.py
+++ b/bfps/resize.py
@@ -210,7 +210,7 @@ def double(opt):
     c.write_src()
     c.write_par(simname = 'test1')
     if opt.run:
-        c.generate_vector_field(simname = 'test1')
+        #c.generate_vector_field(simname = 'test1')
         c.run(ncpu = opt.ncpu,
               simname = 'test1')
 
diff --git a/src/fluid_solver.cpp b/src/fluid_solver.cpp
index bc1402c5..e1af0841 100644
--- a/src/fluid_solver.cpp
+++ b/src/fluid_solver.cpp
@@ -28,36 +28,13 @@
 template <class rnumber>
 void fluid_solver<rnumber>::impose_zero_modes()
 {
-    this->cu[0][0] = 0.0;
-    this->cu[1][0] = 0.0;
-    this->cu[2][0] = 0.0;
-    this->cu[0][1] = 0.0;
-    this->cu[1][1] = 0.0;
-    this->cu[2][1] = 0.0;
-
-    this->cv[0][0][0] = 0.0;
-    this->cv[0][1][0] = 0.0;
-    this->cv[0][2][0] = 0.0;
-
-    this->cv[1][0][0] = 0.0;
-    this->cv[1][1][0] = 0.0;
-    this->cv[1][2][0] = 0.0;
-
-    this->cv[2][0][0] = 0.0;
-    this->cv[2][1][0] = 0.0;
-    this->cv[2][2][0] = 0.0;
-
-    this->cv[0][0][1] = 0.0;
-    this->cv[0][1][1] = 0.0;
-    this->cv[0][2][1] = 0.0;
-
-    this->cv[1][0][1] = 0.0;
-    this->cv[1][1][1] = 0.0;
-    this->cv[1][2][1] = 0.0;
-
-    this->cv[2][0][1] = 0.0;
-    this->cv[2][1][1] = 0.0;
-    this->cv[2][2][1] = 0.0;
+    if (this->cd->myrank == this->cd->rank[0])
+    {
+        std::fill_n((rnumber*)(this->cu), 6, 0.0);
+        std::fill_n((rnumber*)(this->cv[0]), 6, 0.0);
+        std::fill_n((rnumber*)(this->cv[1]), 6, 0.0);
+        std::fill_n((rnumber*)(this->cv[2]), 6, 0.0);
+    }
 }
 /*****************************************************************************/
 /* macro for specializations to numeric types compatible with FFTW           */
@@ -196,30 +173,29 @@ fluid_solver<R>::~fluid_solver() \
 template<> \
 void fluid_solver<R>::compute_vorticity() \
 { \
-    this->low_pass_Fourier(this->cu, 3, this->kM); \
+    /*this->low_pass_Fourier(this->cu, 3, this->kM);*/ \
     CLOOP( \
-            this->cvorticity[3*cindex+0][0] = (this->ky[yindex]*this->cu[3*cindex+2][1] - this->kz[zindex]*this->cu[3*cindex+1][1]); \
-            this->cvorticity[3*cindex+1][0] = (this->kz[zindex]*this->cu[3*cindex+0][1] - this->kx[xindex]*this->cu[3*cindex+2][1]); \
-            this->cvorticity[3*cindex+2][0] = (this->kx[xindex]*this->cu[3*cindex+1][1] - this->ky[yindex]*this->cu[3*cindex+0][1]); \
-            this->cvorticity[3*cindex+0][1] = (this->ky[yindex]*this->cu[3*cindex+2][0] - this->kz[zindex]*this->cu[3*cindex+1][0]); \
-            this->cvorticity[3*cindex+1][1] = (this->kz[zindex]*this->cu[3*cindex+0][0] - this->kx[xindex]*this->cu[3*cindex+2][0]); \
-            this->cvorticity[3*cindex+2][1] = (this->kx[xindex]*this->cu[3*cindex+1][0] - this->ky[yindex]*this->cu[3*cindex+0][0]); \
+            this->cvorticity[3*cindex+0][0] = -(this->ky[yindex]*this->cu[3*cindex+2][1] - this->kz[zindex]*this->cu[3*cindex+1][1]); \
+            this->cvorticity[3*cindex+1][0] = -(this->kz[zindex]*this->cu[3*cindex+0][1] - this->kx[xindex]*this->cu[3*cindex+2][1]); \
+            this->cvorticity[3*cindex+2][0] = -(this->kx[xindex]*this->cu[3*cindex+1][1] - this->ky[yindex]*this->cu[3*cindex+0][1]); \
+            this->cvorticity[3*cindex+0][1] =  (this->ky[yindex]*this->cu[3*cindex+2][0] - this->kz[zindex]*this->cu[3*cindex+1][0]); \
+            this->cvorticity[3*cindex+1][1] =  (this->kz[zindex]*this->cu[3*cindex+0][0] - this->kx[xindex]*this->cu[3*cindex+2][0]); \
+            this->cvorticity[3*cindex+2][1] =  (this->kx[xindex]*this->cu[3*cindex+1][0] - this->ky[yindex]*this->cu[3*cindex+0][0]); \
             ); \
-    this->impose_zero_modes(); \
-    this->symmetrize(this->cvorticity, 3); \
+    /*this->symmetrize(this->cvorticity, 3);*/ \
 } \
  \
 template<> \
 void fluid_solver<R>::compute_velocity(C *vorticity) \
 { \
     double k2; \
-    this->low_pass_Fourier(vorticity, 3, this->kM); \
+    /*this->low_pass_Fourier(vorticity, 3, this->kM);*/ \
     std::fill_n((R*)this->cu, this->cd->local_size*2, 0.0); \
     CLOOP( \
             k2 = (this->kx[xindex]*this->kx[xindex] + \
                   this->ky[yindex]*this->ky[yindex] + \
                   this->kz[zindex]*this->kz[zindex]); \
-            if (k2 < this->kM2) \
+            /*if (k2 < this->kM2)*/ \
             { \
                 this->cu[3*cindex+0][0] = -(this->ky[yindex]*vorticity[3*cindex+2][1] - this->kz[zindex]*vorticity[3*cindex+1][1]) / k2; \
                 this->cu[3*cindex+1][0] = -(this->kz[zindex]*vorticity[3*cindex+0][1] - this->kx[xindex]*vorticity[3*cindex+2][1]) / k2; \
@@ -228,7 +204,7 @@ void fluid_solver<R>::compute_velocity(C *vorticity) \
                 this->cu[3*cindex+1][1] =  (this->kz[zindex]*vorticity[3*cindex+0][0] - this->kx[xindex]*vorticity[3*cindex+2][0]) / k2; \
                 this->cu[3*cindex+2][1] =  (this->kx[xindex]*vorticity[3*cindex+1][0] - this->ky[yindex]*vorticity[3*cindex+0][0]) / k2; \
             } \
-            else \
+            /*else \
             { \
                 this->cu[3*cindex+0][0] = 0.0; \
                 this->cu[3*cindex+1][0] = 0.0; \
@@ -236,10 +212,12 @@ void fluid_solver<R>::compute_velocity(C *vorticity) \
                 this->cu[3*cindex+0][1] = 0.0; \
                 this->cu[3*cindex+1][1] = 0.0; \
                 this->cu[3*cindex+2][1] = 0.0; \
-            } \
+            }*/ \
             ); \
-    this->impose_zero_modes(); \
-    this->symmetrize(this->cu, 3); \
+    if (this->cd->myrank == this->cd->rank[0]) \
+        std::fill_n((R*)(this->cu), 6, 0.0); \
+    /*this->impose_zero_modes();*/ \
+    /*this->symmetrize(this->cu, 3);*/ \
     /*this->low_pass_Fourier(this->cu, 3, this->kM);*/ \
 } \
  \
@@ -349,8 +327,9 @@ void fluid_solver<R>::omega_nonlin( \
             this->cu[3*cindex+1][1] = tmpy1; \
             this->cu[3*cindex+2][1] = tmpz1; \
             ); \
-    /*this->symmetrize(this->cu, 3);*/ \
+    this->symmetrize(this->cu, 3); \
     this->add_forcing(this->cu, 1.0); \
+    /*this->force_divfree(this->cu);*/ \
 } \
  \
 template<> \
diff --git a/src/fluid_solver_base.cpp b/src/fluid_solver_base.cpp
index 33860eeb..b3156a56 100644
--- a/src/fluid_solver_base.cpp
+++ b/src/fluid_solver_base.cpp
@@ -314,7 +314,6 @@ void fluid_solver_base<R>::force_divfree(C *a) \
             k2 = (this->kx[xindex]*this->kx[xindex] + \
                   this->ky[yindex]*this->ky[yindex] + \
                   this->kz[zindex]*this->kz[zindex]); \
-            if (k2 >= this->kM2) \
             { \
                 tval[0] = (this->kx[xindex]*((*(a + cindex*3  ))[0]) + \
                            this->ky[yindex]*((*(a + cindex*3+1))[0]) + \
@@ -322,14 +321,16 @@ void fluid_solver_base<R>::force_divfree(C *a) \
                 tval[1] = (this->kx[xindex]*((*(a + cindex*3  ))[1]) + \
                            this->ky[yindex]*((*(a + cindex*3+1))[1]) + \
                            this->kz[zindex]*((*(a + cindex*3+2))[1]) ) / k2; \
-                a[cindex*3  ][0] -= tval[0]*this->kx[xindex]; \
-                a[cindex*3+1][1] -= tval[1]*this->kx[xindex]; \
-                a[cindex*3+2][0] -= tval[0]*this->ky[yindex]; \
-                a[cindex*3  ][1] -= tval[1]*this->ky[yindex]; \
-                a[cindex*3+1][0] -= tval[0]*this->kz[zindex]; \
-                a[cindex*3+2][1] -= tval[1]*this->kz[zindex]; \
+                for (int imag_part=0; imag_part<2; imag_part++) \
+                { \
+                    a[cindex*3  ][imag_part] -= tval[imag_part]*this->kx[xindex]; \
+                    a[cindex*3+1][imag_part] -= tval[imag_part]*this->ky[yindex]; \
+                    a[cindex*3+2][imag_part] -= tval[imag_part]*this->kz[zindex]; \
+                } \
             } \
             );\
+    if (this->cd->myrank == this->cd->rank[0]) \
+        std::fill_n((R*)(a), 6, 0.0); \
 } \
  \
 template<> \
diff --git a/test.py b/test.py
index d3b414a3..997fe738 100755
--- a/test.py
+++ b/test.py
@@ -23,6 +23,7 @@
 from bfps.test import convergence_test
 from bfps.NavierStokes import test as NStest
 from bfps.resize import double as resize_test
+from bfps.test_curl import test as test_curl
 
 import numpy as np
 import subprocess
@@ -275,6 +276,7 @@ if __name__ == '__main__':
     parser.add_argument('-n', type = int, dest = 'n', default = 64)
     parser.add_argument('--wd', type = str, dest = 'work_dir', default = 'data')
     opt = parser.parse_args()
+    test_curl(opt)
     #NStest(opt)
-    resize_test(opt)
+    #resize_test(opt)
 
-- 
GitLab