#include #include #include #include "epot_bicgstabsolver.hpp" #include "meshvectorfield.hpp" #include "gtkplotter.hpp" #include "geomplotter.hpp" #include "geometry.hpp" #include "func_solid.hpp" #include "epot_efield.hpp" #include "random.hpp" #include "error.hpp" #include "ibsimu.hpp" #include "trajectorydiagnostics.hpp" #include "particledatabase.hpp" #include "particlediagplotter.hpp" using namespace std; void snapshot( ParticleDataBase3D &pdb, string fn, double t ) { ibsimu.message(1) << "Snapshot at t = " << t << "\n"; string fn2 = fn + "_" + to_string(t) + ".txt"; ofstream of( fn2.c_str() ); for( size_t i = 0; i < pdb.size(); i++ ) { Particle3D p = pdb.particle(i); if( p.get_status() != PARTICLE_OK ) continue; // When using particle coordinates, please note that as a // Leapfrog algorithm is used, the velocity coordinates are // 0.5*dt behind the location coordinates in time. of << p[0] << " " << p[1] << " " << p[2] << " " << p[3] << " " << p[4] << " " << p[5] << " " << p[6] << "\n"; } } void simu( int *argc, char ***argv ) { double h = 1e-3; Vec3D origo( -30e-3, -30e-3, -30e-3 ); double sizereq[3] = { 60e-3, 60e-3, 60e-3 }; Int3D meshsize( (int)floor(sizereq[0]/h)+1, (int)floor(sizereq[1]/h)+1, (int)floor(sizereq[2]/h)+1 ); Geometry geom( MODE_3D, meshsize, origo, h ); geom.set_boundary( 1, Bound(BOUND_DIRICHLET, 0.0) ); geom.set_boundary( 2, Bound(BOUND_DIRICHLET, 0.0) ); geom.set_boundary( 3, Bound(BOUND_DIRICHLET, 0.0) ); geom.set_boundary( 4, Bound(BOUND_DIRICHLET, 0.0) ); geom.set_boundary( 5, Bound(BOUND_DIRICHLET, 0.0) ); geom.set_boundary( 6, Bound(BOUND_DIRICHLET, 0.0) ); geom.build_mesh(); EpotBiCGSTABSolver solver( geom ); EpotField epot( geom ); MeshScalarField scharge( geom ); solver.solve( epot, scharge ); MeshVectorField bfield; EpotEfield efield( epot ); field_extrpl_e efldextrpl[6] = { FIELD_EXTRAPOLATE, FIELD_EXTRAPOLATE, FIELD_EXTRAPOLATE, FIELD_EXTRAPOLATE, FIELD_EXTRAPOLATE, FIELD_EXTRAPOLATE }; efield.set_extrapolation( efldextrpl ); ParticleDataBase3D pdb( geom ); bool pmirror[6] = { false, false, false, false, false, false }; pdb.set_mirror( pmirror ); pdb.set_save_trajectories( 1 ); ibsimu.message(1) << "Adding particles\n"; // Add a sphere of particles with zero starting velocity, each particle carries // a charge of 1000e. double q = 1.0; double m = 1.0; double r = 10e-3; double C = CHARGE_E; double N = 1000; size_t i = 0; QRandom rng(3); double x[3]; while( i < N ) { rng.get(x); x[0] = r*(-1+2*x[0]); x[1] = r*(-1+2*x[1]); x[2] = r*(-1+2*x[2]); if( x[0]*x[0] + x[1]*x[1] + x[2]*x[2] > r*r ) continue; // Strange, but critical: created particles have to be created // with time (first coordinate) zero. That is used to recognize new particles. pdb.add_particle( C, q, m, ParticleP3D(0,x[0],0,x[1],0,x[2],0) ); i++; } snapshot( pdb, "pout", 0.0 ); ofstream of( "field.txt" ); ibsimu.message(1) << "Particle-in-cell loop\n"; double t = 0.0; // This example has 1 us time step. Be careful with time step // size, Courant-Friedrichs-Lewy condition dt < h/v should be taken in account. double dt = 1e-6; i = 0; while( t < 0.3e-3 ) { pdb.step_particles( scharge, efield, bfield, dt ); solver.solve( epot, scharge ); efield.recalculate(); i++; t = i*dt; if( i%10 == 0 ) snapshot( pdb, "pout", t ); of << t << " " << epot(Vec3D(0,0,0)) << "\n"; } MeshScalarField tdens(geom); pdb.build_trajectory_density_field(tdens); GTKPlotter plotter( argc, argv ); plotter.set_geometry( &geom ); plotter.set_epot( &epot ); plotter.set_bfield( &bfield ); plotter.set_trajdens( &tdens ); plotter.set_scharge( &scharge ); plotter.set_particledatabase( &pdb ); plotter.new_geometry_plot_window(); plotter.run(); } int main( int argc, char **argv ) { try { //ibsimu.set_message_output( "ibsimu.txt" ); ibsimu.set_message_threshold( MSG_VERBOSE, 1 ); ibsimu.set_thread_count( 4 ); simu( &argc, &argv ); } catch( Error e ) { e.print_error_message( ibsimu.message( 0 ) ); exit( 1 ); } return( 0 ); }