debugging

bfps/cpp/fluid_solver.cpp

@@ -692,9 +692,9 @@ void fluid_solver<R>::compute_Lagrangian_acceleration(R *acceleration) \
                 this->cv[1][tindex+2][1] = k2*pressure[cindex][1]; \
             } \
             ); \
-    char fname[256]; \
+    /*char fname[256]; \
     this->fill_up_filename("pressure_trS2_enstrophy", fname); \
-    this->cd->write(fname, this->cv[1]); \
+    this->cd->write(fname, this->cv[1]);*/ \
     /*double difference = 0.0; \
     double itval, rtval; \
     CLOOP_K2_NXMODES( \

tests/test_base.py

@@ -128,6 +128,31 @@ def launch(
               njobs = opt.njobs)
     return c
 
+def acceleration_test(c):
+    import numpy as np
+    import matplotlib.pyplot as plt
+    d = c.get_data_file()
+    pos = d['particles/tracers4/state'].value
+    vel = d['particles/tracers4/velocity'].value
+    acc = d['particles/tracers4/acceleration'].value
+
+    num_acc1 = (- vel[ :-2] + vel[2:])/(2*d['parameters/dt'].value*d['parameters/niter_part'].value)
+    num_vel = (- pos[ :-2] + pos[2:])/(2*d['parameters/dt'].value*d['parameters/niter_part'].value)
+    num_acc2 = (pos[ :-2] - 2*pos[1:-1] + pos[2:])/((d['parameters/dt'].value*d['parameters/niter_part'].value)**2)
+    num_acc3 = (-vel[4:] + 8*vel[3:-1] - 8*vel[1:-3] + vel[:-4])/(12*d['parameters/dt'].value*d['parameters/niter_part'].value)
+
+    pid = np.unravel_index(np.argmax(np.abs(num_acc3 - acc[2:-2])), dims = num_acc3.shape)[1]
+    fig = plt.figure()
+    a = fig.add_subplot(111)
+    col = ['red', 'green', 'blue']
+    for cc in range(3):
+        a.plot(num_acc1[1:, pid, cc], color = col[cc])
+        a.plot(num_acc2[1:, pid, cc], color = col[cc], dashes = (2, 2))
+        a.plot(num_acc3[1:, pid, cc], color = col[cc], dashes = (3, 4))
+        a.plot(acc[2:, pid, cc], color = col[cc], dashes = (1, 1))
+    fig.savefig(os.path.join(c.work_dir, 'acc_test_{0}.pdf'.format(c.simname)))
+    return None
+
 if __name__ == '__main__':
     print('this file doesn\'t do anything')
 

tests/test_convergence.py

@@ -93,6 +93,7 @@ def convergence_test(
     c2.compute_statistics()
     c2.set_plt_style({'dashes': (3, 4)})
     for c in [c0, c1, c2]:
+        acceleration_test(c)
         c.style.update({'label' : '${0}\\times {1} \\times {2}$'.format(c.parameters['nx'],
                                                                         c.parameters['ny'],
                                                                         c.parameters['nz'])})

tests/test_plain.py

@@ -38,6 +38,7 @@ def plain(opt):
     opt.work_dir = wd + '/N{0:0>3x}_1'.format(opt.n)
     c0 = launch(opt, dt = 0.2/opt.n)
     c0.compute_statistics()
+    acceleration_test(c0)
     if not opt.multiplejob:
         return None
     assert(opt.niter_todo % 3 == 0)

tox_plain.ini

@@ -10,7 +10,10 @@ changedir =
 commands =
     cp -r {toxinidir}/tests {envtmpdir}
     #python tests/test_plain.py -n 256 --run --initialize --ncpu 8 --niter_todo 8 --precision single --wd "data/single"
-    python tests/test_plain.py -n 32 --run --initialize --ncpu 2 --nparticles 16 --niter_todo 64 --precision single --wd "data/single"
+    python tests/test_plain.py \
+        -n 32 --run --initialize --ncpu 2 \
+        --nparticles 256 --niter_todo 64 \
+        --precision single --wd "data/single"
     #python tests/test_plain.py -n 32 --run --initialize --multiplejob --ncpu 2 --nparticles 16 --niter_todo 64 --precision single --wd "data/single"
     #python tests/test_plain.py -n 32 --run --initialize --multiplejob --ncpu 2 --nparticles 16 --niter_todo 64 --precision double --wd "data/double"