Merge branch 'feature/Stokes-drag' into develop

TurTLE/test/test_turtle_NSVE_Stokes_particles.py

@@ -35,25 +35,23 @@ import matplotlib.pyplot as plt
 import TurTLE
 from TurTLE import DNS
 
-
+def theory_exponential(t,initial_value,drag_coeff):
+    return initial_value*np.exp(-t*drag_coeff)
 
 def main():
     """
-        Run test where the trajectory of initially moving particles in a
-	quiescent flow is compared to the analytical exponential  
+        Run test where the trajectory of initially moving particles with varying drag coefficient in a
+	quiescent flow is compared to the analytically expected exponential. The generated plot shows the accuracy of the corresponding particles for a certain fixed timestep.  
     """
-    drag_coefficients = ['30.','20.','10.','1.','0.1']
+    drag_coefficients = ['30','20','10','1']
     niterations = 100
     nparticles = 1
     njobs = 1
     dt=0.01
-    
-    
     fig = plt.figure()
     plt.ylabel('$v$')
     plt.xlabel('$t$')
     plt.yscale('log')
-
     j=0
     for drag_coeff in drag_coefficients:
         c = DNS()
@@ -68,7 +66,6 @@ def main():
              '--niter_todo', '{0}'.format(niterations),
              '--niter_out', '{0}'.format(1),
              '--niter_stat', '1',
-             '--checkpoints_per_file', '{0}'.format(3),
              '--nparticles', '{0}'.format(nparticles),
              '--initial_field_amplitude', '0',
              '--initial_particle_vel', '0.05',
@@ -77,29 +74,25 @@ def main():
              '--njobs', '{0}'.format(njobs),
              '--wd', './'] +
              sys.argv[1:])
-
-        
         f = h5py.File('quiescent_nsve_stokes_particles_drag'+str(j)+'_particles.h5', 'r')
         f2 = h5py.File('quiescent_nsve_stokes_particles_drag'+str(j)+'.h5', 'r')
-        
         particle_momentum = [np.sqrt(f['tracers0/momentum/'+'{}'.format(i)][0][0]**2
                                      +f['tracers0/momentum/'+'{}'.format(i)][0][1]**2
                                      +f['tracers0/momentum/'+'{}'.format(i)][0][2]**2)
                              for i in range(niterations)]
         time_values = np.arange(0, niterations*dt, dt)
-        def theory_exponential(t):
-            return f2['parameters/initial_particle_vel']*np.exp(-t*np.float(drag_coeff))
-        
-
-        numerics, = plt.plot(time_values, particle_momentum, alpha = 0.3, lw=2)
-        plt.plot(time_values, theory_exponential(time_values), color=numerics.get_color(), ls = '--')
-        plt.plot(0,0,color=numerics.get_color(),label=r'$\frac{1}{St}=$'+drag_coeff)
-        plt.legend(loc='lower left')
-        print(f['tracers0/momentum/0'][:])
-        print(f['tracers0/momentum/100'][:]-f['tracers0/momentum/0'][:])
-        print(f['tracers0/position/100'][:]-f['tracers0/position/0'][:])
-        plt.savefig('Stokes_test1.pdf')
+        numerics, = plt.plot(time_values, particle_momentum, alpha = 1, lw=2)
+        plt.plot(time_values, 
+            theory_exponential(time_values, f2['parameters/initial_particle_vel'],np.float(drag_coeff)),
+            color='black', ls = '--')
+        plt.plot(0,0,color=numerics.get_color(),label=r'drag coefficient '+drag_coeff)
         j+=1
+    plt.plot(0, 
+            0,
+            color='black', ls = '--',
+            label='analytic sol.')
+    plt.legend(loc='lower left')
+    plt.savefig('Stokes_quiescent_test.pdf')
     return None
 
 if __name__ == '__main__':

cpp/full_code/NSVE_Stokes_particles.cpp

@@ -46,8 +46,6 @@ int NSVE_Stokes_particles<rnumber>::initialize(void)
 
     particles_inner_computer_2nd_order_Stokes<double, long long int> current_particles_inner_computer;
     current_particles_inner_computer.set_drag_coefficient(this->drag_coefficient);
-    //DEBUG_MSG("drag coefficient is set to %f \n", current_particles_inner_computer.get_drag_coefficient());
-    //DEBUG_MSG_WAIT(MPI_COMM_WORLD, "before call to particles_system_builder\n");
     this->ps = particles_system_builder_with_p2p(
                 this->fs->cvelocity,                                                    // (field object)
                 this->fs->kk,                                                           // (kspace object, contains dkx, dky, dkz)
@@ -79,7 +77,6 @@ int NSVE_Stokes_particles<rnumber>::initialize(void)
                 "tracers0",
                 "position/0");
     this->particles_sample_writer_mpi->setParticleFileLayout(this->ps->getParticleFileLayout());
-    //DEBUG_MSG("drag coefficient is after initialization %f \n", current_particles_inner_computer.get_drag_coefficient());
     return EXIT_SUCCESS;
 }