diff --git a/TurTLE/DNS.py b/TurTLE/DNS.py index d7da0753d41aef9cd4edc08218155fe795bd3027..71c01187090e0bdc3478712d0b603111a9efd88e 100644 --- a/TurTLE/DNS.py +++ b/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, diff --git a/TurTLE/test/test_turtle_NSVE_Stokes_particles.py b/TurTLE/test/test_turtle_NSVE_Stokes_particles.py index d1372600fd63ffd0b96c59363e13ff85b908e288..6a51ee30ee07cbe4071c7c09a5229000474eb449 100644 --- a/TurTLE/test/test_turtle_NSVE_Stokes_particles.py +++ b/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__': diff --git a/cpp/full_code/NSVE_Stokes_particles.cpp b/cpp/full_code/NSVE_Stokes_particles.cpp index ef1c13ce5aaf78dc11285a790d0b29b34d892aca..19bc0f9151ada7becc10f4525352b0b14760d792 100644 --- a/cpp/full_code/NSVE_Stokes_particles.cpp +++ b/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; } diff --git a/cpp/full_code/NSVE_Stokes_particles.hpp b/cpp/full_code/NSVE_Stokes_particles.hpp index eec5247530a5b43647d847be262c99345495cd78..8e818041b966fe26a6163609e5f2fb50ee8f2489 100644 --- a/cpp/full_code/NSVE_Stokes_particles.hpp +++ b/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;