diff --git a/TurTLE/test/test_turtle_NSVE_Stokes_particles.py b/TurTLE/test/test_turtle_NSVE_Stokes_particles.py index 6a51ee30ee07cbe4071c7c09a5229000474eb449..400a0cdb12ba0790343489a8570f9383a16daa80 100644 --- a/TurTLE/test/test_turtle_NSVE_Stokes_particles.py +++ b/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__': diff --git a/cpp/full_code/NSVE_Stokes_particles.cpp b/cpp/full_code/NSVE_Stokes_particles.cpp index 19bc0f9151ada7becc10f4525352b0b14760d792..a4b3c41e8908a853333d6041b073d8742021b87a 100644 --- a/cpp/full_code/NSVE_Stokes_particles.cpp +++ b/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; }