test for different drag coefficients in quiescent flow

TurTLE/DNS.py

@@ -704,8 +704,9 @@ class DNS(_code):
             pars['spectrum_k_cutoff'] = float(16)
             pars['spectrum_coefficient'] = float(1)
         if dns_type == 'NSVE_Stokes_particles':
-            pars['initial_field_amplitude'] = float(0.05)
+            pars['initial_field_amplitude'] = float(0.0)
             pars['initial_particle_vel'] = float(0.05)
+            pars['drag_coefficient'] = float(0.1)
         return pars
     def prepare_launch(
             self,

TurTLE/test/test_turtle_NSVE_Stokes_particles.py

@@ -29,42 +29,77 @@ import os
 import numpy as np
 import h5py
 import sys
+import matplotlib.pyplot as plt
+
 
 import TurTLE
 from TurTLE import DNS
 
 
+
 def main():
     """
         Run test where the trajectory of initially moving particles in a
 	quiescent flow is compared to the analytical exponential  
     """
+    drag_coefficients = ['30.','20.','10.','1.','0.1']
     niterations = 100
-    nparticles = 10
+    nparticles = 1
     njobs = 1
-    c = DNS()
-    c.launch(
+    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()
+        c.launch(
             ['NSVE_Stokes_particles',
              '-n', '32',
-             '--simname', 'quiescent_nsve_stokes_particles',
+             '--simname', 'quiescent_nsve_stokes_particles_drag'+str(j),
              '--np', '4',
              '--ntpp', '1',
+             '--dt', '{0}'.format(dt),
              '--fftw_plan_rigor', 'FFTW_PATIENT',
              '--niter_todo', '{0}'.format(niterations),
-             '--niter_out', '{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',
+             '--drag_coefficient', drag_coeff, 
              '--famplitude', '0',
              '--njobs', '{0}'.format(njobs),
              '--wd', './'] +
              sys.argv[1:])
-    f = h5py.File('quiescent_nsve_stokes_particles_particles.h5', 'r')
-    print(f['tracers0/momentum/0'][:])
-    print(f['tracers0/momentum/100'][:]-f['tracers0/momentum/0'][:])
-    print(f['tracers0/position/100'][:]-f['tracers0/position/0'][:])
+
+        
+        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')
+        j+=1
     return None
 
 if __name__ == '__main__':

cpp/full_code/NSVE_Stokes_particles.cpp

@@ -45,8 +45,8 @@ int NSVE_Stokes_particles<rnumber>::initialize(void)
             this->fs->cvelocity->fftw_plan_rigor);
 
     particles_inner_computer_2nd_order_Stokes<double, long long int> current_particles_inner_computer;
-    current_particles_inner_computer.set_drag_coefficient(0.1);
-    DEBUG_MSG("drag coefficient is set to %f \n", current_particles_inner_computer.get_drag_coefficient());
+    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)
@@ -79,7 +79,7 @@ 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());
+    //DEBUG_MSG("drag coefficient is after initialization %f \n", current_particles_inner_computer.get_drag_coefficient());
     return EXIT_SUCCESS;
 }
 
@@ -207,6 +207,7 @@ int NSVE_Stokes_particles<rnumber>::read_parameters(void)
     this->tracers0_neighbours = hdf5_tools::read_value<int>(parameter_file, "parameters/tracers0_neighbours");
     this->tracers0_smoothness = hdf5_tools::read_value<int>(parameter_file, "parameters/tracers0_smoothness");
     this->tracers0_cutoff = hdf5_tools::read_value<double>(parameter_file, "parameters/tracers0_cutoff");
+    this->drag_coefficient = hdf5_tools::read_value<double>(parameter_file, "parameters/drag_coefficient");
     H5Fclose(parameter_file);
     return EXIT_SUCCESS;
 }

cpp/full_code/NSVE_Stokes_particles.hpp

@@ -58,7 +58,7 @@ class NSVE_Stokes_particles: public NSVE<rnumber>
         int tracers0_neighbours;
         int tracers0_smoothness;
         double tracers0_cutoff;
-
+        double drag_coefficient;
 
         /* other stuff */
         std::unique_ptr<abstract_particles_system<long long int, double>> ps;