diff --git a/.DS_Store b/.DS_Store deleted file mode 100644 index 254ba9333..000000000 Binary files a/.DS_Store and /dev/null differ diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 5ba5074c6..3791701ae 100755 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -6,7 +6,7 @@ compile: stage: compile only: - develop - - smilei_Poisson_AM + - DDfusion script: # Force workdir cleaning in case of retried @@ -28,7 +28,7 @@ makerun1D: stage: makerun only: - develop - - smilei_Poisson_AM + - DDfusion script: # Move in test dir @@ -40,7 +40,7 @@ makerun2D: stage: makerun only: - develop - - smilei_Poisson_AM + - DDfusion script: # Move in test dir @@ -52,7 +52,7 @@ makerun3D: stage: makerun only: - develop - - smilei_Poisson_AM + - DDfusion script: # Move in test dir @@ -64,7 +64,7 @@ makerunAM: stage: makerun only: - develop - - smilei_Poisson_AM + - DDfusion script: # Move in test dir @@ -76,7 +76,7 @@ makerunV: stage: makerun only: - develop - - smilei_Poisson_AM + - DDfusion script: # Move in test dir @@ -88,7 +88,7 @@ makerunCollisions: stage: makerun only: - develop - - smilei_Poisson_AM + - DDfusion script: # Move in test dir diff --git a/benchmarks/tst_collisions4_DD_fusion_rate.py b/benchmarks/tst_collisions4_DD_fusion_rate.py new file mode 100644 index 000000000..0d438eafb --- /dev/null +++ b/benchmarks/tst_collisions4_DD_fusion_rate.py @@ -0,0 +1,133 @@ +# --------------------------------------------- +# SIMULATION PARAMETERS FOR THE PIC-CODE SMILEI +# --------------------------------------------- + +import math +L0 = 2.*math.pi # conversion from normalization length to wavelength + + +Main( + geometry = "1Dcartesian", + + number_of_patches = [ 8 ], + + interpolation_order = 2, + + timestep = 0.05 * L0, + simulation_time = 10 * L0, + + + time_fields_frozen = 100000000000., + + cell_length = [5.*L0], + grid_length = [1600.*L0], + + EM_boundary_conditions = [ ["periodic"] ], + + reference_angular_frequency_SI = L0 * 3e8 /1.e-6, + + random_seed = smilei_mpi_rank +) + + +i = 0 +for Da_nppc, Db_nppc, v in [ + [100, 100, 0.1], + [100, 10, 0.1], + [10 , 100, 0.1], + [100, 100, 0.00082], + [100, 100, 0.00201], + [100, 100, 0.00493], + [100, 100, 0.01209], + [100, 100, 0.02960], + [100, 100, 0.07236], + [100, 100, 0.17545], + [100, 100, 0.40641], + [100, 100, 0.76966], + [100, 100, 0.97376], + ]: + + Da = "Da_"+str(i) + Db = "Db_"+str(i) + He = "He_"+str(i) + + Species( + name = Da, + position_initialization = "regular", + momentum_initialization = "maxwell-juettner", + particles_per_cell= Da_nppc, + atomic_number = 1, + mass = 3870.5, + charge = 0., + number_density = 100., + mean_velocity = [v, 0., 0.], + temperature = [0.0000001]*3, + time_frozen = 100000000.0, + boundary_conditions = [ + ["periodic", "periodic"], + ], + ) + + Species( + name = Db, + position_initialization = "regular", + momentum_initialization = "maxwell-juettner", + particles_per_cell= Db_nppc, + atomic_number = 1, + mass = 3870.5, + charge = 0., + number_density = 100., + mean_velocity = [-v, 0., 0.], + temperature = [0.00000001]*3, + time_frozen = 100000000.0, + boundary_conditions = [ + ["periodic", "periodic"], + ], + ) + + Species( + name = He, + position_initialization = "regular", + momentum_initialization = "maxwell-juettner", + particles_per_cell= 0, + atomic_number = 2, + mass = 5497.9, + charge = 0., + number_density = 0., + time_frozen = 100000000.0, + boundary_conditions = [ + ["periodic", "periodic"], + ], + ) + + Collisions( + species1 = [Da], + species2 = [Db], + coulomb_log = 0.001, + nuclear_reaction = [He], + debug_every = 10 + ) + + DiagParticleBinning( + deposited_quantity = "weight", + every = 20, + species = [He], + axes = [ + ["x", 0, Main.grid_length[0], 1] + ] + ) + + DiagParticleBinning( + deposited_quantity = "weight", + every = 20, + species = [Da, Db], + axes = [ + ["x", 0, Main.grid_length[0], 1] + ] + ) + + i+=1 + +DiagScalar( + every = 50, +) \ No newline at end of file diff --git a/doc/Sphinx/namelist.rst b/doc/Sphinx/namelist.rst index 94408fa43..6bb2a7642 100755 --- a/doc/Sphinx/namelist.rst +++ b/doc/Sphinx/namelist.rst @@ -1961,6 +1961,7 @@ Collisions coulomb_log = 5., debug_every = 1000, ionizing = False, + # nuclear_reaction = [], ) @@ -2013,7 +2014,7 @@ Collisions .. py:data:: ionizing - :default: False + :default: ``False`` :ref:`Collisional ionization ` is set when this parameter is not ``False``. It can either be set to the name of a pre-existing electron species (where the ionized @@ -2023,7 +2024,36 @@ Collisions One of the species groups must be all electrons (:py:data:`mass` = 1), and the other one all ions of the same :py:data:`atomic_number`. +.. rst-class:: experimental + +.. py:data:: nuclear_reaction + + :type: a list of strings + :default: ``None`` (no nuclear reaction) + + A list of the species names for the products of nuclear reactions + that may occur during collisions. You may omit product species if they are not necessary + for the simulation. + + All members of :py:data:`species1` must be the same type of atoms, which is automatically + recognized by their :py:data:`mass` and :py:data:`atomic_number`. The same applies for + all members of :py:data:`species2`. + + In the current version, only the reaction D(d,n)He³ is available. + +.. rst-class:: experimental + +.. py:data:: nuclear_reaction_multiplier + + :type: a float + :default: 0. (automatically adjusted) + The rate multiplier for nuclear reactions. It is a positive number that artificially + increases the occurence of reactions so that a good statistics is obtained. The number + of actual reaction products is adjusted by changing their weights in order to provide + a physically correct number of reactions. Leave this number to ``0.`` for an automatic + rate multiplier: the final number of produced macro-particles will be of the same order + as that of reactants. diff --git a/src/Checkpoint/Checkpoint.cpp b/src/Checkpoint/Checkpoint.cpp index 8f6b2ce16..b04d81db1 100755 --- a/src/Checkpoint/Checkpoint.cpp +++ b/src/Checkpoint/Checkpoint.cpp @@ -26,6 +26,7 @@ #include "DiagnosticScreen.h" #include "DiagnosticTrack.h" #include "LaserEnvelope.h" +#include "Collisions.h" using namespace std; @@ -255,7 +256,7 @@ void Checkpoint::dumpAll( VectorPatch &vecPatches, unsigned int itime, SmileiMP string patchName=Tools::merge( "patch-", patch_name.str() ); hid_t patch_gid = H5::group( fid, patchName.c_str() ); - dumpPatch( vecPatches( ipatch )->EMfields, vecPatches( ipatch )->vecSpecies, params, patch_gid ); + dumpPatch( vecPatches( ipatch )->EMfields, vecPatches( ipatch )->vecSpecies, vecPatches( ipatch )->vecCollisions, params, patch_gid ); // Random number generator state H5::attr( patch_gid, "xorshift32_state", vecPatches( ipatch )->xorshift32_state ); @@ -283,7 +284,7 @@ void Checkpoint::dumpAll( VectorPatch &vecPatches, unsigned int itime, SmileiMP } -void Checkpoint::dumpPatch( ElectroMagn *EMfields, std::vector vecSpecies, Params ¶ms, hid_t patch_gid ) +void Checkpoint::dumpPatch( ElectroMagn *EMfields, std::vector vecSpecies, std::vector &vecCollisions, Params ¶ms, hid_t patch_gid ) { if ( params.geometry != "AMcylindrical" ) { dumpFieldsPerProc( patch_gid, EMfields->Ex_ ); @@ -442,6 +443,13 @@ void Checkpoint::dumpPatch( ElectroMagn *EMfields, std::vector vecSpe H5Gclose( gid ); } // End for ispec + + // Manage some collisions parameters + std::vector rate_multiplier( vecCollisions.size() ); + for( unsigned int icoll = 0; icollNuclearReaction->rate_multiplier_; + } + H5::vect( patch_gid, "collisions_rate_multiplier", rate_multiplier ); }; @@ -539,7 +547,7 @@ void Checkpoint::restartAll( VectorPatch &vecPatches, SmileiMPI *smpi, SimWindo string patchName=Tools::merge( "patch-", patch_name.str() ); hid_t patch_gid = H5Gopen( fid, patchName.c_str(), H5P_DEFAULT ); - restartPatch( vecPatches( ipatch )->EMfields, vecPatches( ipatch )->vecSpecies, params, patch_gid ); + restartPatch( vecPatches( ipatch )->EMfields, vecPatches( ipatch )->vecSpecies, vecPatches( ipatch )->vecCollisions, params, patch_gid ); // Random number generator state H5::getAttr( patch_gid, "xorshift32_state", vecPatches( ipatch )->xorshift32_state ); @@ -566,7 +574,7 @@ void Checkpoint::restartAll( VectorPatch &vecPatches, SmileiMPI *smpi, SimWindo } -void Checkpoint::restartPatch( ElectroMagn *EMfields, std::vector &vecSpecies, Params ¶ms, hid_t patch_gid ) +void Checkpoint::restartPatch( ElectroMagn *EMfields, std::vector &vecSpecies, std::vector &vecCollisions, Params ¶ms, hid_t patch_gid ) { if ( params.geometry != "AMcylindrical" ) { restartFieldsPerProc( patch_gid, EMfields->Ex_ ); @@ -754,6 +762,14 @@ void Checkpoint::restartPatch( ElectroMagn *EMfields, std::vector &ve H5Gclose( gid ); } + // Manage some collisions parameters + if( H5::getVectSize( patch_gid, "collisions_rate_multiplier" ) > 0 ) { + std::vector rate_multiplier; + H5::getVect( patch_gid, "collisions_rate_multiplier", rate_multiplier, true ); + for( unsigned int icoll = 0; icollNuclearReaction->rate_multiplier_ = rate_multiplier[icoll]; + } + } } void Checkpoint::dumpFieldsPerProc( hid_t fid, Field *field ) diff --git a/src/Checkpoint/Checkpoint.h b/src/Checkpoint/Checkpoint.h index afa222bb5..4bb95e32f 100755 --- a/src/Checkpoint/Checkpoint.h +++ b/src/Checkpoint/Checkpoint.h @@ -23,6 +23,7 @@ class Field; class cField; class Species; class VectorPatch; +class Collisions; #include @@ -42,7 +43,7 @@ class Checkpoint //! restart everything to file per processor void readPatchDistribution( SmileiMPI *smpi, SimWindow *simWin ); void restartAll( VectorPatch &vecPatches, SmileiMPI *smpi, SimWindow *simWin, Params ¶ms, OpenPMDparams &openPMD ); - void restartPatch( ElectroMagn *EMfields, std::vector &vecSpecies, Params ¶ms, hid_t patch_gid ); + void restartPatch( ElectroMagn *EMfields, std::vector &vecSpecies, std::vector &vecCollisions, Params ¶ms, hid_t patch_gid ); //! restart field per proc void restartFieldsPerProc( hid_t fid, Field *field ); @@ -58,7 +59,7 @@ class Checkpoint //! dump everything to file per processor void dumpAll( VectorPatch &vecPatches, unsigned int itime, SmileiMPI *smpi, SimWindow *simWin, Params ¶ms ); - void dumpPatch( ElectroMagn *EMfields, std::vector vecSpecies, Params ¶ms, hid_t patch_gid ); + void dumpPatch( ElectroMagn *EMfields, std::vector vecSpecies, std::vector &vecCollisions, Params ¶ms, hid_t patch_gid ); //! incremental number of times we've done a dump unsigned int dump_number; diff --git a/src/Collisions/CollisionalFusionDD.cpp b/src/Collisions/CollisionalFusionDD.cpp new file mode 100644 index 000000000..e19fdbd19 --- /dev/null +++ b/src/Collisions/CollisionalFusionDD.cpp @@ -0,0 +1,110 @@ +#include "CollisionalFusionDD.h" + +#include "Collisions.h" +#include "Species.h" +#include "Patch.h" + +#include + +using namespace std; + +// Coefficients used for energy interpolation +// The list of energies is in logarithmic scale, +// with Emin=1 keV, Emax=631 MeV and npoints=50. +const int CollisionalFusionDD::npoints = 50; +const double CollisionalFusionDD::npointsm1 = ( double )( npoints-1 ); +const double CollisionalFusionDD::a1 = log(511./(2.*2.013553)); // = ln(me*c^2 / Emin / n_nucleons) +const double CollisionalFusionDD::a2 = 3.669039; // = (npoints-1) / ln( Emax/Emin ) +const double CollisionalFusionDD::a3 = log(511000./(2.*2.013553));; // = ln(me*c^2 / eV / n_nucleons) +// Log of cross-section in units of 4 pi re^2 +const double CollisionalFusionDD::DB_log_crossSection[50] = { + -27.307, -23.595, -20.383, -17.607, -15.216, -13.167, -11.418, -9.930, -8.666, -7.593, + -6.682, -5.911, -5.260, -4.710, -4.248, -3.858, -3.530, -3.252, -3.015, -2.814, -2.645, + -2.507, -2.398, -2.321, -2.273, -2.252, -2.255, -2.276, -2.310, -2.352, -2.402, -2.464, + -2.547, -2.664, -2.831, -3.064, -3.377, -3.773, -4.249, -4.794, -5.390, -6.021, -6.677, + -7.357, -8.072, -8.835, -9.648, -10.498, -11.387, -12.459 +}; + +// Constructor +CollisionalFusionDD::CollisionalFusionDD( + Params *params, + vector product_particles, + vector product_species, + double rate_multiplier +) +: CollisionalNuclearReaction(params, &product_particles, &product_species, rate_multiplier) +{ +} + +// Cloning constructor +CollisionalFusionDD::CollisionalFusionDD( CollisionalNuclearReaction *NR ) +: CollisionalNuclearReaction( NR ) +{ +} + + +bool CollisionalFusionDD::occurs( double U, double coeff, double m1, double m2, double g1, double g2, double &ekin, double &log_ekin, double &W ) +{ + // Interpolate the total cross-section at some value of ekin = m1(g1-1) + m2(g2-1) + ekin = m1 * (g1-1.) + m2 * (g2-1.); + log_ekin = log( ekin ); + double x = a2*( a1 + log_ekin ); + double cs; + // if energy below Emin, approximate to 0. + if( x < 0. ) { + cs = 0.; + } + // if energy within table range, interpolate + else if( x < npointsm1 ) { + int i = int( x ); + double a = x - ( double )i; + cs = exp( + ( DB_log_crossSection[i+1]-DB_log_crossSection[i] )*a + DB_log_crossSection[i] + ); + } + // if energy above table range, extrapolate + else { + double a = x - npointsm1; + cs = exp( + ( DB_log_crossSection[npoints-1]-DB_log_crossSection[npoints-2] )*a + DB_log_crossSection[npoints-1] + ); + } + + // Calculate probability for fusion + double prob = coeff * cs * rate_multiplier_; + tot_probability_ += prob; + if( U > exp( -prob ) ) { + W /= rate_multiplier_; + return true; + } else { + return false; + } +} + +void CollisionalFusionDD::makeProducts( + double U, double ekin, double log_ekin, double q, + Particles *&p3, Particles *&p4, + double &p3_COM, double &p4_COM, + double &q3, double &q4, + double &cosX +) { + // Sample the products angle from empirical fits + double A = 0.083 * (a3 + log_ekin); + A *= A*A; // ^3 + A *= A; // ^6 + A = 1. / min(10., 1. + A); + double B = 0.06 * (a3 + log_ekin); + B = 2. + B*B; + cosX = 1. - A*U - (1.-A)*pow(U, B); + + // Calculate the resulting momenta + const double Q = 6.397; // Qvalue + const double m_n = 1838.7; + const double m_He = 5497.9; + p3_COM = sqrt( (ekin+Q) * (ekin+Q+2.*m_n) * (ekin+Q+2.*m_He) * (ekin+Q+2.*m_n+2.*m_He) ) + / ( ( ekin+Q+m_n+m_He ) * (2.*m_He) ); + + // Other properties + q3 = q; + p3 = product_particles_[0]; // helium3 +} diff --git a/src/Collisions/CollisionalFusionDD.h b/src/Collisions/CollisionalFusionDD.h new file mode 100644 index 000000000..52e77d6fb --- /dev/null +++ b/src/Collisions/CollisionalFusionDD.h @@ -0,0 +1,37 @@ + +#ifndef COLLISIONALFUSIONDD_H +#define COLLISIONALFUSIONDD_H + +#include + +#include "CollisionalNuclearReaction.h" + +class CollisionalFusionDD : public CollisionalNuclearReaction +{ + +public: + //! Constructor + CollisionalFusionDD( Params*, std::vector, std::vector, double ); + //! Cloning Constructor + CollisionalFusionDD( CollisionalNuclearReaction * ); + //! Destructor + ~CollisionalFusionDD() {}; + + //! Method to apply the nuclear reaction + bool occurs( double U, double coeff, double m1, double m2, double g1, double g2, double &ekin, double &log_ekin, double &W ) override; + //! Method to prepare the products of the reaction + void makeProducts( double U, double etot, double log_ekin, double q, Particles *&p3, Particles *&p4, double &p3_COM, double &p4_COM, double &q3, double &q4, double &cosX ) override; + + std::string name() override { return "D-D fusion"; }; + + //! static parameters to read the database + static const double a1, a2, a3, npointsm1; + static const int npoints; + static const double DB_log_crossSection[50]; + +private: + + +}; + +#endif diff --git a/src/Collisions/CollisionalIonization.cpp b/src/Collisions/CollisionalIonization.cpp index 81897035e..4abf3899f 100755 --- a/src/Collisions/CollisionalIonization.cpp +++ b/src/Collisions/CollisionalIonization.cpp @@ -19,23 +19,18 @@ const double CollisionalIonization::a1 = 510998.9 ; // = me*c^2/Emin const double CollisionalIonization::a2 = 6.142165 ; // = (npoints-1) / ln( Emax/Emin ) // Constructor -CollisionalIonization::CollisionalIonization( int Z, int nDim_, double reference_angular_frequency_SI, int ionization_electrons, Particles* particles ) +CollisionalIonization::CollisionalIonization( int Z, Params *params, int ionization_electrons, Particles* particles ) { - nDim = nDim_; atomic_number = Z; rate .resize( Z ); irate.resize( Z ); prob .resize( Z ); ionization_electrons_ = ionization_electrons; - if( particles ) { - new_electrons.tracked = particles->tracked; - new_electrons.isQuantumParameter = particles->isQuantumParameter; - new_electrons.isMonteCarlo = particles->isMonteCarlo; + if( params ) { + new_electrons.initialize( 0, *particles ); } - new_electrons.initialize( 0, nDim ); - if( Z>0 ) { - dataBaseIndex = createDatabase( reference_angular_frequency_SI ); + dataBaseIndex = createDatabase( params->reference_angular_frequency_SI ); assignDatabase( dataBaseIndex ); } } @@ -44,16 +39,12 @@ CollisionalIonization::CollisionalIonization( int Z, int nDim_, double reference CollisionalIonization::CollisionalIonization( CollisionalIonization *CI ) { - nDim = CI->nDim; atomic_number = CI->atomic_number; rate .resize( atomic_number ); irate.resize( atomic_number ); prob .resize( atomic_number ); - ionization_electrons_ = CI->ionization_electrons_; - new_electrons.tracked = CI->new_electrons.tracked; - new_electrons.isQuantumParameter = CI->new_electrons.isQuantumParameter; - new_electrons.isMonteCarlo = CI->new_electrons.isMonteCarlo; - new_electrons.initialize( 0, nDim ); + ionization_electrons_ = CI->ionization_electrons_; + new_electrons.initialize( 0, CI->new_electrons ); dataBaseIndex = CI->dataBaseIndex; assignDatabase( dataBaseIndex ); @@ -285,9 +276,9 @@ void CollisionalIonization::calculate( double gamma_s, double gammae, double gam e = ( ( *lostEnergy )[Zstar][i+1]-( *lostEnergy )[Zstar][i] )*a + ( *lostEnergy )[Zstar][i]; } else { // if energy above table range, extrapolate a = x - npointsm1; - cs = ( ( *crossSection )[Zstar][npointsm1]-( *crossSection )[Zstar][npointsm1-1] )*a + ( *crossSection )[Zstar][npointsm1]; - w = ( *transferredEnergy )[Zstar][npointsm1]; - e = ( *lostEnergy )[Zstar][npointsm1]; + cs = ( ( *crossSection )[Zstar][npoints-1]-( *crossSection )[Zstar][npoints-2] )*a + ( *crossSection )[Zstar][npoints-1]; + w = ( *transferredEnergy )[Zstar][npoints-1]; + e = ( *lostEnergy )[Zstar][npoints-1]; } if( e > gamma_s-1. ) { break; diff --git a/src/Collisions/CollisionalIonization.h b/src/Collisions/CollisionalIonization.h index 41c9287ce..74613cbe6 100755 --- a/src/Collisions/CollisionalIonization.h +++ b/src/Collisions/CollisionalIonization.h @@ -15,7 +15,7 @@ class CollisionalIonization public: //! Constructor - CollisionalIonization( int, int, double, int, Particles* ); + CollisionalIonization( int, Params*, int, Particles* ); //! Cloning Constructor CollisionalIonization( CollisionalIonization * ); //! Destructor @@ -62,9 +62,6 @@ class CollisionalIonization unsigned int dataBaseIndex; private: - - //! Simulation dimension - int nDim; //! Atomic number int atomic_number; @@ -107,7 +104,7 @@ class CollisionalIonization class CollisionalNoIonization : public CollisionalIonization { public: - CollisionalNoIonization() : CollisionalIonization( 0, 0, 0., -1, NULL ) {}; + CollisionalNoIonization() : CollisionalIonization( 0, NULL, -1, NULL ) {}; ~CollisionalNoIonization() {}; unsigned int createDatabase( double ) override @@ -124,6 +121,4 @@ class CollisionalNoIonization : public CollisionalIonization }; -//! Class to hold the databases - #endif diff --git a/src/Collisions/CollisionalNuclearReaction.cpp b/src/Collisions/CollisionalNuclearReaction.cpp new file mode 100644 index 000000000..f8694c2de --- /dev/null +++ b/src/Collisions/CollisionalNuclearReaction.cpp @@ -0,0 +1,106 @@ +#include "CollisionalNuclearReaction.h" + +#include "Collisions.h" +#include "Species.h" +#include "Particles.h" +#include "Patch.h" + +#include + + +using namespace std; + +// Constructor +CollisionalNuclearReaction::CollisionalNuclearReaction( + Params *params, + vector *product_particles, + vector *product_species, + double rate_multiplier + ) +{ + if( rate_multiplier == 0. ) { + auto_multiplier_ = true; + rate_multiplier_ = 1.; + } else { + auto_multiplier_ = false; + rate_multiplier_ = rate_multiplier; + } + product_particles_.resize(0); + product_species_.resize(0); + if( params ) { + for( unsigned int i=0; isize(); i++ ) { + if( product_particles->at(i) != NULL ) { + product_particles_.push_back( new Particles() ); + product_particles_.back()->initialize( 0, *product_particles->at(i) ); + product_species_.push_back(product_species->at(i)); + } + } + } +} + +// Cloning Constructor +CollisionalNuclearReaction::CollisionalNuclearReaction( CollisionalNuclearReaction *CNR ) +{ + product_species_ = CNR->product_species_; + rate_multiplier_ = CNR->rate_multiplier_; + auto_multiplier_ = CNR->auto_multiplier_; + product_particles_.resize( CNR->product_particles_.size(), NULL ); + for( unsigned int i=0; iproduct_particles_.size(); i++ ) { + product_particles_[i] = new Particles(); + product_particles_[i]->initialize( 0, *CNR->product_particles_[i] ); + } +} + +// Cloning Constructor +CollisionalNuclearReaction::~CollisionalNuclearReaction() +{ + for( unsigned int i=0; i &localDiags, + bool intra_collisions, vector sg1, vector sg2, + double npairs, int itime +) { + // Move new particles in place + for( unsigned int i=0; ivecSpecies[product_species_[i]]->importParticles( params, patch, *product_particles_[i], localDiags ); + } + + // Remove reactants that have been (very rare) + for( unsigned int is = 0; is < sg1.size(); is++ ) { + patch->vecSpecies[sg1[is]]->eraseWeightlessParticles(); + } + if( ! intra_collisions ) { + for( unsigned int is = 0; is < sg2.size(); is++ ) { + patch->vecSpecies[sg2[is]]->eraseWeightlessParticles(); + } + } + + if( auto_multiplier_ ) { + // Update the rate multiplier + double goal = tot_probability_ * (double) params.n_time / npairs; + if( goal > 0. ) { + if( goal > 2. ) { + rate_multiplier_ *= 0.01; + } else if( goal > 1.5 ) { + rate_multiplier_ *= 0.1; + } else if( goal > 1. ) { + rate_multiplier_ *= 0.3; + } else if( rate_multiplier_ < 1.e14 ) { + if( goal < 0.01 ) { + rate_multiplier_ *= 8.; + } else if( goal < 0.1 ) { + rate_multiplier_ *= 2.; + } + } + if( rate_multiplier_ < 1. ) { + rate_multiplier_ = 1.; + } + } + } +} diff --git a/src/Collisions/CollisionalNuclearReaction.h b/src/Collisions/CollisionalNuclearReaction.h new file mode 100644 index 000000000..5a76b4a5e --- /dev/null +++ b/src/Collisions/CollisionalNuclearReaction.h @@ -0,0 +1,73 @@ +#ifndef COLLISIONALNUCLEARREACTION_H +#define COLLISIONALNUCLEARREACTION_H + +#include + +#include "Tools.h" +#include "Species.h" +#include "Params.h" + +class Patch; + +class CollisionalNuclearReaction +{ + +public: + //! Constructor + CollisionalNuclearReaction( Params *, std::vector*, std::vector*, double ); + //! Cloning Constructor + CollisionalNuclearReaction( CollisionalNuclearReaction * ); + //! Destructor + virtual ~CollisionalNuclearReaction(); + + //! Prepare the nuclear reaction + virtual void prepare() { + tot_probability_ = 0.; + }; + //! Test the occurence of the nuclear reaction + virtual bool occurs( double U, double coeff, double m1, double m2, double g1, double g2, double &etot, double &log_ekin, double &W ) = 0; + //! Prepare the products of the reaction + virtual void makeProducts( double U, double ekin, double log_ekin, double q, Particles *&p3, Particles *&p4, double &p3_COM, double &p4_COM, double &q3, double &q4, double &cosX ) = 0; + //! Finish the nuclear reaction and put new electrons in place + virtual void finish( Params &, Patch *, std::vector &, bool, std::vector, std::vector, double npairs, int itime ); + + //! Get the name of the reaction + virtual std::string name() = 0; + + //! New electrons temporary species + std::vector product_particles_; + + //! Coefficient to adapt the rate of create of new particles + double rate_multiplier_; + + //! True if rate multiplier isn't automatically adjusted + double auto_multiplier_; + + //! sum of probabilities in this patch + double tot_probability_; + +private: + + //! Species where products are sent + std::vector product_species_; + +}; + +//! Class to make empty reaction objects +class CollisionalNoNuclearReaction : public CollisionalNuclearReaction +{ +public: + CollisionalNoNuclearReaction() : CollisionalNuclearReaction( NULL, NULL, NULL, 0. ) {}; + ~CollisionalNoNuclearReaction() {}; + + void prepare() override {}; + bool occurs( double U, double coeff, double m1, double m2, double g1, double g2, double &etot, double &log_ekin, double &W ) override { + return false; + }; + void makeProducts( double U, double ekin, double log_ekin, double q, Particles *&p3, Particles *&p4, double &p3_COM, double &p4_COM, double &q3, double &q4, double &cosX ) override {}; + void finish( Params ¶ms, Patch *patch, std::vector &diags, bool, std::vector, std::vector, double npairs, int itime ) override {}; + std::string name() { return ""; } +}; + + +#endif diff --git a/src/Collisions/Collisions.cpp b/src/Collisions/Collisions.cpp index 26b3434c7..178aa23b8 100755 --- a/src/Collisions/Collisions.cpp +++ b/src/Collisions/Collisions.cpp @@ -24,34 +24,27 @@ Collisions::Collisions( double coulomb_log, bool intra_collisions, int debug_every, - int Z, - int ionization_electrons, - Particles * ionization_particles, - int nDim, + CollisionalIonization *ionization, + CollisionalNuclearReaction *nuclear_reaction, string filename ) : + Ionization( ionization ), + NuclearReaction( nuclear_reaction ), n_collisions_( n_collisions ), species_group1_( species_group1 ), species_group2_( species_group2 ), coulomb_log_( coulomb_log ), intra_collisions_( intra_collisions ), debug_every_( debug_every ), - atomic_number( Z ), filename_( filename ) { - // Create the ionization object - if( ionization_electrons >= 0 ) { - Ionization = new CollisionalIonization( Z, nDim, params.reference_angular_frequency_SI, ionization_electrons, ionization_particles ); - } else { - Ionization = new CollisionalNoIonization(); - } coeff1_ = 4.046650232e-21*params.reference_angular_frequency_SI; // h*omega/(2*me*c^2) coeff2_ = 2.817940327e-15*params.reference_angular_frequency_SI/299792458.; // re omega / c } // Cloning Constructor -Collisions::Collisions( Collisions *coll, int nDim ) +Collisions::Collisions( Collisions *coll ) { n_collisions_ = coll->n_collisions_ ; @@ -60,15 +53,20 @@ Collisions::Collisions( Collisions *coll, int nDim ) coulomb_log_ = coll->coulomb_log_ ; intra_collisions_ = coll->intra_collisions_; debug_every_ = coll->debug_every_ ; - atomic_number = coll->atomic_number ; filename_ = coll->filename_ ; coeff1_ = coll->coeff1_ ; coeff2_ = coll->coeff2_ ; - if( atomic_number>0 ) { - Ionization = new CollisionalIonization( coll->Ionization ); - } else { + if( dynamic_cast( coll->Ionization ) ) { Ionization = new CollisionalNoIonization(); + } else { + Ionization = new CollisionalIonization( coll->Ionization ); + } + + if( dynamic_cast( coll->NuclearReaction ) ) { + NuclearReaction = new CollisionalNoNuclearReaction(); + } else if ( dynamic_cast( coll->NuclearReaction ) ) { + NuclearReaction = new CollisionalFusionDD( coll->NuclearReaction ); } } @@ -76,6 +74,7 @@ Collisions::Collisions( Collisions *coll, int nDim ) Collisions::~Collisions() { delete Ionization; + delete NuclearReaction; } // Declare other static variables here @@ -179,7 +178,7 @@ void Collisions::collide( Params ¶ms, Patch *patch, int itime, vector 0 && itime % debug_every_ == 0 ); // debug only every N timesteps - if( debug ) { - ncol = 0.; - smean_ = 0.; - logLmean_ = 0.; - //temperature = 0.; - } + ncol = 0.; + smean_ = 0.; + logLmean_ = 0.; + //temperature = 0.; + + NuclearReaction->prepare(); // Loop bins of particles (typically, cells, but may also be clusters) unsigned int nbin = patch->vecSpecies[0]->first_index.size(); @@ -318,7 +317,7 @@ void Collisions::collide( Params ¶ms, Patch *patch, int itime, vectorprepare3( params.timestep, inv_cell_volume ); // Now start the real loop on pairs of particles @@ -344,19 +343,17 @@ void Collisions::collide( Params ¶ms, Patch *patch, int itime, vectorparticles; p2 = s2->particles; - m12 = s1->mass_ / s2->mass_; // mass ratio - logL = coulomb_log_; double U1 = patch->xorshift32() * patch->xorshift32_invmax; double U2 = patch->xorshift32() * patch->xorshift32_invmax; double phi = patch->xorshift32() * patch->xorshift32_invmax * twoPi; - s = one_collision( p1, i1, s1->mass_, p2, i2, m12, coeff1_, coeff2_, coeff3, coeff4, n123, n223, debye2, logL, U1, U2, phi ); + s = one_collision( p1, i1, s1->mass_, p2, i2, s2->mass_, coeff1_, coeff2_, coeff3, coeff4, n123, n223, debye2, logL, U1, U2, phi ); - // Handle ionization + // Handle ionization & nuclear reaction Ionization->apply( patch, p1, i1, p2, i2 ); + ncol ++; if( debug ) { - ncol += 1; smean_ += s; logLmean_ += logL; //temperature += m1 * (sqrt(1.+pow(p1->momentum(0,i1),2)+pow(p1->momentum(1,i1),2)+pow(p1->momentum(2,i1),2))-1.); @@ -367,6 +364,7 @@ void Collisions::collide( Params ¶ms, Patch *patch, int itime, vectorfinish( params, patch, localDiags ); + NuclearReaction->finish( params, patch, localDiags, intra_collisions_, species_group1_, species_group2_, ncol, itime ); if( debug && ncol>0. ) { smean_ /= ncol; @@ -388,7 +386,8 @@ void Collisions::debug( Params ¶ms, int itime, unsigned int icoll, VectorPat vector smean( npatch, 0. ); vector logLmean( npatch, 0. ); //vector temperature=(npatch, 0.); - vector debye_length_squared( npatch, 0. ); + vector debye_length( npatch, 0. ); + vector nuclear_reaction_multiplier( npatch, 0. ); // Collect info for all patches for( unsigned int ipatch=0; ipatchvecCollisions[icoll]->smean_ ; logLmean [ipatch] = vecPatches( ipatch )->vecCollisions[icoll]->logLmean_ ; //temperature[ipatch] = vecPatches(ipatch)->vecCollisions[icoll]->temperature; - debye_length_squared[ipatch] = 0.; unsigned int nbin = vecPatches( ipatch )->debye_length_squared.size(); for( unsigned int ibin=0; ibindebye_length_squared[ibin]; + debye_length[ipatch] += vecPatches( ipatch )->debye_length_squared[ibin]; } + debye_length[ipatch] = sqrt( debye_length[ipatch] / nbin ); + nuclear_reaction_multiplier[ipatch] = vecPatches( ipatch )->vecCollisions[icoll]->NuclearReaction->rate_multiplier_; } // Open the HDF5 file @@ -438,7 +438,10 @@ void Collisions::debug( Params ¶ms, int itime, unsigned int icoll, VectorPat //H5Dwrite( dset_id, H5T_NATIVE_DOUBLE, memspace, filespace, transfer, &temperature[0] ); //H5Dclose(dset_id); dset_id = H5Dcreate( group, "debyelength", H5T_NATIVE_DOUBLE, filespace, H5P_DEFAULT, plist_id, H5P_DEFAULT ); - H5Dwrite( dset_id, H5T_NATIVE_DOUBLE, memspace, filespace, transfer, &debye_length_squared[0] ); + H5Dwrite( dset_id, H5T_NATIVE_DOUBLE, memspace, filespace, transfer, &debye_length[0] ); + H5Dclose( dset_id ); + dset_id = H5Dcreate( group, "nuclear_reaction_multiplier", H5T_NATIVE_DOUBLE, filespace, H5P_DEFAULT, plist_id, H5P_DEFAULT ); + H5Dwrite( dset_id, H5T_NATIVE_DOUBLE, memspace, filespace, transfer, &nuclear_reaction_multiplier[0] ); H5Dclose( dset_id ); // Close all H5Pclose( plist_id ); diff --git a/src/Collisions/Collisions.h b/src/Collisions/Collisions.h index d77e8c8ef..706f64e03 100755 --- a/src/Collisions/Collisions.h +++ b/src/Collisions/Collisions.h @@ -7,6 +7,8 @@ #include "Tools.h" #include "H5.h" #include "CollisionalIonization.h" +#include "CollisionalNuclearReaction.h" +#include "CollisionalFusionDD.h" class Patch; class Params; @@ -18,12 +20,20 @@ class Collisions public: //! Constructor for Collisions between two species - Collisions( Params ¶ms, unsigned int n_collisions, std::vector, - std::vector, double coulomb_log, bool intra_collisions, - int debug_every, int Z, int ionization_electrons, Particles * ionization_particles, - int nDim, std::string ); + Collisions( + Params ¶ms, + unsigned int n_collisions, + std::vector, + std::vector, + double coulomb_log, + bool intra_collisions, + int debug_every, + CollisionalIonization *ionization, + CollisionalNuclearReaction *nuclear_reaction, + std::string + ); //! Cloning Constructor - Collisions( Collisions *, int ); + Collisions( Collisions * ); //! destructor virtual ~Collisions(); @@ -42,6 +52,8 @@ class Collisions //! CollisionalIonization object, created if ionization required CollisionalIonization *Ionization; + CollisionalNuclearReaction *NuclearReaction; + protected: //! Identification number of the Collisions object @@ -59,9 +71,6 @@ class Collisions //! Number of timesteps between each dump of collisions debugging int debug_every_; - //! Species atomic number, in case of ionization - int atomic_number; - //! Hdf5 file name std::string filename_; @@ -79,7 +88,7 @@ class Collisions double m1, Particles *p2, unsigned int i2, - double m12, + double m2, double coeff1, double coeff2, double coeff3, @@ -93,34 +102,31 @@ class Collisions double phi ) { - double qqm, qqm2, gamma1, gamma2, gamma12, gamma12_inv, - COM_vx, COM_vy, COM_vz, COM_vsquare, COM_gamma, - term1, term2, term3, term4, term5, term6, - vcv1, vcv2, px_COM, py_COM, pz_COM, p2_COM, p_COM, gamma1_COM, gamma2_COM, - bmin, s, vrel, smax, - cosX, sinX, sinXcosPhi, sinXsinPhi, p_perp, inv_p_perp, + double term6, cosX, sinX, sinXcosPhi, sinXsinPhi, p_perp, inv_p_perp, newpx_COM, newpy_COM, newpz_COM, vcp; - - // Calculate stuff - qqm = p1->charge( i1 ) * p2->charge( i2 ) / m1; - qqm2 = qqm * qqm; + + // If one weight is zero, then skip. Can happen after nuclear reaction + double minW = std::min( p1->weight(i1), p2->weight(i2) ); + if( minW <= 0. ) return 0.; // Get momenta and calculate gammas - gamma1 = sqrt( 1. + p1->momentum( 0, i1 )*p1->momentum( 0, i1 ) + p1->momentum( 1, i1 )*p1->momentum( 1, i1 ) + p1->momentum( 2, i1 )*p1->momentum( 2, i1 ) ); - gamma2 = sqrt( 1. + p2->momentum( 0, i2 )*p2->momentum( 0, i2 ) + p2->momentum( 1, i2 )*p2->momentum( 1, i2 ) + p2->momentum( 2, i2 )*p2->momentum( 2, i2 ) ); - gamma12 = m12 * gamma1 + gamma2; - gamma12_inv = 1./gamma12; + double m12 = m1 / m2; + double gamma1 = sqrt( 1. + p1->momentum( 0, i1 )*p1->momentum( 0, i1 ) + p1->momentum( 1, i1 )*p1->momentum( 1, i1 ) + p1->momentum( 2, i1 )*p1->momentum( 2, i1 ) ); + double gamma2 = sqrt( 1. + p2->momentum( 0, i2 )*p2->momentum( 0, i2 ) + p2->momentum( 1, i2 )*p2->momentum( 1, i2 ) + p2->momentum( 2, i2 )*p2->momentum( 2, i2 ) ); + double gamma12 = m12 * gamma1 + gamma2; + double gamma12_inv = 1./gamma12; // Calculate the center-of-mass (COM) frame // Quantities starting with "COM" are those of the COM itself, expressed in the lab frame. // They are NOT quantities relative to the COM. - COM_vx = ( m12 * ( p1->momentum( 0, i1 ) ) + p2->momentum( 0, i2 ) ) * gamma12_inv; - COM_vy = ( m12 * ( p1->momentum( 1, i1 ) ) + p2->momentum( 1, i2 ) ) * gamma12_inv; - COM_vz = ( m12 * ( p1->momentum( 2, i1 ) ) + p2->momentum( 2, i2 ) ) * gamma12_inv; - COM_vsquare = COM_vx*COM_vx + COM_vy*COM_vy + COM_vz*COM_vz; + double COM_vx = ( m12 * ( p1->momentum( 0, i1 ) ) + p2->momentum( 0, i2 ) ) * gamma12_inv; + double COM_vy = ( m12 * ( p1->momentum( 1, i1 ) ) + p2->momentum( 1, i2 ) ) * gamma12_inv; + double COM_vz = ( m12 * ( p1->momentum( 2, i1 ) ) + p2->momentum( 2, i2 ) ) * gamma12_inv; + double COM_vsquare = COM_vx*COM_vx + COM_vy*COM_vy + COM_vz*COM_vz; // Change the momentum to the COM frame (we work only on particle 1) // Quantities ending with "COM" are quantities of the particle expressed in the COM frame. + double COM_gamma, term1, term2, vcv1, vcv2, px_COM, py_COM, pz_COM, gamma1_COM, gamma2_COM; if( COM_vsquare != 0. ) { COM_gamma = 1./sqrt( 1.-COM_vsquare ); term1 = ( COM_gamma - 1. ) / COM_vsquare; @@ -142,17 +148,128 @@ class Collisions gamma1_COM = gamma1; gamma2_COM = gamma2; } - p2_COM = px_COM*px_COM + py_COM*py_COM + pz_COM*pz_COM; - p_COM = sqrt( p2_COM ); + double p2_COM = px_COM*px_COM + py_COM*py_COM + pz_COM*pz_COM; + double p_COM = sqrt( p2_COM ); // Calculate some intermediate quantities - term3 = COM_gamma * gamma12_inv; - term4 = gamma1_COM * gamma2_COM; - term5 = term4/p2_COM + m12; + double term3 = COM_gamma * gamma12_inv; + double term4 = gamma1_COM * gamma2_COM; + double term5 = term4/p2_COM + m12; + double vrel = p_COM/term3/term4; // relative velocity + + // We first try to do a nuclear reaction + // If succesful, then no need to do a collision + double E, logE; + if( NuclearReaction->occurs( U1, vrel*coeff3, m1, m2, gamma1_COM, gamma2_COM, E, logE, minW ) ) { + // Reduce the weight of both reactants + // If becomes zero, then the particle will be discarded later + p1->weight(i1) -= minW; + p2->weight(i2) -= minW; + + // Get the magnitude and the angle of the outgoing products in the COM frame + Particles *p3 = NULL, *p4 = NULL; + double p3_COM, p4_COM, q3, q4; + double tot_charge = p1->charge( i1 ) + p2->charge( i2 ); + U2 = 2*U2 - 1.; + if( U2 > 0. ) { + NuclearReaction->makeProducts( abs(U2), E, logE, tot_charge, p3, p4, p3_COM, p4_COM, q3, q4, cosX ); + } else { + NuclearReaction->makeProducts( abs(U2), E, logE, tot_charge, p4, p3, p4_COM, p3_COM, q4, q3, cosX ); + } + + // Calculate combination of angles + sinX = sqrt( 1. - cosX*cosX ); + sinXcosPhi = sinX*cos( phi ); + sinXsinPhi = sinX*sin( phi ); + + // Calculate the deflection in the COM frame + p_perp = sqrt( px_COM*px_COM + py_COM*py_COM ); + if( p_perp > 1.e-10*p_COM ) { // make sure p_perp is not too small + inv_p_perp = 1./p_perp; + newpx_COM = ( px_COM * pz_COM * sinXcosPhi - py_COM * p_COM * sinXsinPhi ) * inv_p_perp + px_COM * cosX; + newpy_COM = ( py_COM * pz_COM * sinXcosPhi + px_COM * p_COM * sinXsinPhi ) * inv_p_perp + py_COM * cosX; + newpz_COM = -p_perp * sinXcosPhi + pz_COM * cosX; + } else { // if p_perp is too small, we use the limit px->0, py=0 + newpx_COM = p_COM * sinXcosPhi; + newpy_COM = p_COM * sinXsinPhi; + newpz_COM = p_COM * cosX; + } + + // Go back to the lab frame and store the results in the particle array + vcp = COM_vx * newpx_COM + COM_vy * newpy_COM + COM_vz * newpz_COM; + double newW1, newW2; + if( p1->charge(i1) != 0. || p2->charge(i2) != 0. ) { + double weight_factor = minW / tot_charge; + newW1 = p1->charge( i1 ) * weight_factor; + newW2 = p2->charge( i2 ) * weight_factor; + } else { + newW1 = minW; + newW2 = 0.; + } + if( p3 ) { + double momentum_ratio = p3_COM / p_COM; + term6 = momentum_ratio*term1*vcp + sqrt( p3_COM*p3_COM + 1. ) * COM_gamma; + double newpx = momentum_ratio * newpx_COM + COM_vx * term6; + double newpy = momentum_ratio * newpy_COM + COM_vy * term6; + double newpz = momentum_ratio * newpz_COM + COM_vz * term6; + // Make new particle at position of particle 1 + if( newW1 > 0. ) { + p1->cp_particle_safe( i1, *p3 ); + p3->Weight.back() = newW1; + p3->Charge.back() = q3; + p3->Momentum[0].back() = newpx; + p3->Momentum[1].back() = newpy; + p3->Momentum[2].back() = newpz; + } + // Make new particle at position of particle 2 + if( newW2 > 0. ) { + p2->cp_particle_safe( i2, *p3 ); + p3->Weight.back() = newW2; + p3->Charge.back() = q3; + p3->Momentum[0].back() = newpx; + p3->Momentum[1].back() = newpy; + p3->Momentum[2].back() = newpz; + } + } + if( p4 ) { + double momentum_ratio = p4_COM / p_COM; + term6 = -momentum_ratio*term1*vcp + sqrt( p4_COM*p4_COM + 1. ) * COM_gamma; + double newpx = -momentum_ratio * newpx_COM + COM_vx * term6; + double newpy = -momentum_ratio * newpy_COM + COM_vy * term6; + double newpz = -momentum_ratio * newpz_COM + COM_vz * term6; + // Make new particle at position of particle 1 + if( newW1 > 0. ) { + p1->cp_particle_safe( i1, *p4 ); + p4->Weight.back() = newW1; + p4->Charge.back() = q4; + p4->Momentum[0].back() = newpx; + p4->Momentum[1].back() = newpy; + p4->Momentum[2].back() = newpz; + } + // Make new particle at position of particle 2 + if( newW2 > 0. ) { + p2->cp_particle_safe( i2, *p4 ); + p4->Weight.back() = newW2; + p4->Charge.back() = q4; + p4->Momentum[0].back() = newpx; + p4->Momentum[1].back() = newpy; + p4->Momentum[2].back() = newpz; + } + } + + if( p1->weight(i1) == 0. || p2->weight(i2) == 0. ) { + return 0.; // no collision + } + + } // end nuclear reaction + + // Calculate stuff + double qqm = p1->charge( i1 ) * p2->charge( i2 ) / m1; + double qqm2 = qqm * qqm; // Calculate coulomb log if necessary if( logL <= 0. ) { // if auto-calculation requested - bmin = std::max( coeff1/m1/p_COM, std::abs( coeff2*qqm*term3*term5 ) ); // min impact parameter + double bmin = std::max( coeff1/m1/p_COM, std::abs( coeff2*qqm*term3*term5 ) ); // min impact parameter logL = 0.5*log( 1.+debye2/( bmin*bmin ) ); if( logL < 2. ) { logL = 2.; @@ -160,11 +277,10 @@ class Collisions } // Calculate the collision parameter s12 (similar to number of real collisions) - s = coeff3 * logL * qqm2 * term3 * p_COM * term5*term5 / ( gamma1*gamma2 ); + double s = coeff3 * logL * qqm2 * term3 * p_COM * term5*term5 / ( gamma1*gamma2 ); // Low-temperature correction - vrel = p_COM/term3/term4; // relative velocity - smax = coeff4 * ( m12+1. ) * vrel / std::max( m12*n123, n223 ); + double smax = coeff4 * ( m12+1. ) * vrel / std::max( m12*n123, n223 ); if( s>smax ) { s = smax; } @@ -224,7 +340,8 @@ class Collisions } return s; - } + + }; }; diff --git a/src/Collisions/CollisionsFactory.h b/src/Collisions/CollisionsFactory.h index fbb091c5a..d155a1d5d 100755 --- a/src/Collisions/CollisionsFactory.h +++ b/src/Collisions/CollisionsFactory.h @@ -1,4 +1,3 @@ - #ifndef COLLISIONSFACTORY_H #define COLLISIONSFACTORY_H @@ -104,6 +103,7 @@ class CollisionsFactory ERROR( "In collisions #" << n_collisions << ": `ionizing` must be True, False, or the name of an electron species" ); } + CollisionalIonization *Ionization; if( ionization ) { if( intra ) { @@ -148,8 +148,92 @@ class CollisionsFactory } } ionization_particles = vecSpecies[ionization_electrons]->particles; + + // Create the ionization object + Ionization = new CollisionalIonization( Z, ¶ms, ionization_electrons, ionization_particles ); + + } else { + + Ionization = new CollisionalNoIonization(); + } + // D-D fusion + std::vector nuclear_reaction(0); + PyObject * py_nuclear_reaction = PyTools::extract_py( "nuclear_reaction", "Collisions", n_collisions ); + CollisionalNuclearReaction *NuclearReaction; + // If fusion, verify parameters + if( py_nuclear_reaction == Py_None ) { + + NuclearReaction = new CollisionalNoNuclearReaction(); + + } else { + + // Extract the content of the list nuclear_reaction + std::vector py_list; + if( ! PyTools::convert( py_nuclear_reaction, py_list ) ) { + ERROR( "In collisions #" << n_collisions << ": nuclear_reaction should be a list" ); + } + std::vector nuclear_reaction( 0 ); + if( py_list.size() ) { + PyTools::convert( py_list, nuclear_reaction ); + for( unsigned int i=0; imass_ != s0->mass_ ) { + ERROR( "In collisions #" << n_collisions << ": nuclear_reaction requires all `species" + << g+1 << "` to have equal masses" ); + } + if( s->atomic_number_ != s0->atomic_number_ ) { + ERROR( "In collisions #" << n_collisions << ": nuclear_reaction requires all `species" + << g+1 << "` to have equal atomic_number" ); + } + } + } + + // Rate multiplier + double rate_multiplier = 0.; + PyTools::extract( "nuclear_reaction_multiplier", rate_multiplier, "Collisions", n_collisions ); + + // Find products + std::vector products = params.FindSpecies( vecSpecies, nuclear_reaction ); + + // Type of reaction + unsigned int A0 = round( vecSpecies[sgroup[0][0]]->mass_ / 1822.89 ); + unsigned int A1 = round( vecSpecies[sgroup[1][0]]->mass_ / 1822.89 ); + std::vector product_particles; + std::vector product_species; + // D-D fusion + if( Z0 == 1 && Z1 == 1 && A0 == 2 && A1 == 2 ) { + + std::vector Z(2); Z[0] = 2; Z[1] = 0; + std::vector A(2); A[0] = 3; A[1] = 1; + std::vector name(2); name[0] = "helium3"; name[1] = "neutron"; + findProducts( vecSpecies, products, Z, A, name, product_particles, product_species, n_collisions ); + + NuclearReaction = new CollisionalFusionDD( ¶ms, product_particles, product_species, rate_multiplier ); + + // Unknown types + } else { + ERROR( "In collisions #" << n_collisions << ": nuclear_reaction not available between (Z="<0 ) { MESSAGE( 2, "Collisional ionization with atomic number "<name_ << "`" ); } + if( py_nuclear_reaction != Py_None ) { + MESSAGE( 2, "Collisional nuclear reaction "<< NuclearReaction->name() ); + } // If debugging log requested if( debug_every>0 ) { @@ -218,10 +305,8 @@ class CollisionsFactory sgroup[1], clog, intra, debug_every, - Z, - ionization_electrons, - ionization_particles, - params.nDim_particle, + Ionization, + NuclearReaction, filename ); } else { @@ -232,10 +317,8 @@ class CollisionsFactory sgroup[1], clog, intra, debug_every, - Z, - ionization_electrons, - ionization_particles, - params.nDim_particle, + Ionization, + NuclearReaction, filename ); } @@ -273,21 +356,64 @@ class CollisionsFactory //! Clone a vector of Collisions objects - static std::vector clone( std::vector vecCollisions, Params ¶ms ) + static std::vector clone( std::vector vecCollisions ) { std::vector newVecCollisions( 0 ); for( unsigned int i=0; i( vecCollisions[i] ) ) { - newVecCollisions.push_back( new CollisionsSingle( vecCollisions[i], params.nDim_particle ) ); + newVecCollisions.push_back( new CollisionsSingle( vecCollisions[i] ) ); } else { - newVecCollisions.push_back( new Collisions( vecCollisions[i], params.nDim_particle ) ); + newVecCollisions.push_back( new Collisions( vecCollisions[i] ) ); } } return newVecCollisions; } + //! Utility for nuclear reactions: find products in the provided list + static void findProducts( + std::vector vecSpecies, + std::vector products, + std::vector Z, + std::vector A, + std::vector name, + std::vector &product_particles, + std::vector &product_species, + unsigned int n_coll + ) + { + unsigned int n = Z.size(); + product_particles.resize( n, NULL ); + product_species.resize( n, 0 ); + + std::ostringstream list(""); + list << name[0]; + for( unsigned int i=1; iatomic_number_; + unsigned int Aj = round( s->mass_ / 1822.89 ); + + bool product_found = false; + for( unsigned int i=0; iparticles; + product_found = true; + } + } + if( ! product_found ) { + ERROR( "In collisions #" << n_coll << ", nuclear_reaction : species `"<name_<<"` is not one of "<vecSpecies[species_group1_[0]]; s2 = patch->vecSpecies[species_group2_[0]]; bool debug = ( debug_every_ > 0 && itime % debug_every_ == 0 ); // debug only every N timesteps + ncol = 0; if( debug ) { - ncol = 0.; smean_ = 0.; logLmean_ = 0.; //temperature = 0.; } + NuclearReaction->prepare(); + // Loop bins of particles (typically, cells, but may also be clusters) unsigned int nbin = patch->vecSpecies[0]->first_index.size(); for( unsigned int ibin = 0 ; ibin < nbin ; ibin++ ) { @@ -128,7 +130,6 @@ void CollisionsSingle::collide( Params ¶ms, Patch *patch, int itime, vector< coeff3 = params.timestep * n1*n2/n12; coeff4 = pow( 3.*coeff2_, -1./3. ) * coeff3; coeff3 *= coeff2_; - m12 = s1->mass_ / s2->mass_; // mass ratio // Prepare the ionization Ionization->prepare3( params.timestep, inv_cell_volume ); @@ -144,13 +145,13 @@ void CollisionsSingle::collide( Params ¶ms, Patch *patch, int itime, vector< double U2 = patch->xorshift32() * patch->xorshift32_invmax; double phi = patch->xorshift32() * patch->xorshift32_invmax * twoPi; - s = one_collision( p1, i1, s1->mass_, p2, i2, m12, coeff1_, coeff2_, coeff3, coeff4, n123, n223, debye2, logL, U1, U2, phi ); + s = one_collision( p1, i1, s1->mass_, p2, i2, s2->mass_, coeff1_, coeff2_, coeff3, coeff4, n123, n223, debye2, logL, U1, U2, phi ); // Handle ionization Ionization->apply( patch, p1, i1, p2, i2 ); + ncol ++; if( debug ) { - ncol += 1; smean_ += s; logLmean_ += logL; //temperature += m1 * (sqrt(1.+pow(p1->momentum(0,i1),2)+pow(p1->momentum(1,i1),2)+pow(p1->momentum(2,i1),2))-1.); @@ -161,6 +162,7 @@ void CollisionsSingle::collide( Params ¶ms, Patch *patch, int itime, vector< } // end loop on bins Ionization->finish( params, patch, localDiags ); + NuclearReaction->finish( params, patch, localDiags, intra_collisions_, species_group1_, species_group2_, ncol, itime ); if( debug && ncol>0. ) { smean_ /= ncol; diff --git a/src/Collisions/CollisionsSingle.h b/src/Collisions/CollisionsSingle.h index 39d659b01..0f5df3852 100755 --- a/src/Collisions/CollisionsSingle.h +++ b/src/Collisions/CollisionsSingle.h @@ -18,33 +18,32 @@ class CollisionsSingle : public Collisions public: //! Constructor for Collisions between two species - CollisionsSingle( Params ¶ms, - unsigned int n_collisions, - std::vector sg1, - std::vector sg2, - double coulomb_log, - bool intra_collisions, - int debug_every, - int Z, - int ionization_electrons, - Particles * ionization_particles, - int nDim, - std::string fname - ) : Collisions( params, - n_collisions, - sg1, - sg2, - coulomb_log, - intra_collisions, - debug_every, - Z, - ionization_electrons, - ionization_particles, - nDim, - fname - ) {} ; + CollisionsSingle( + Params ¶ms, + unsigned int n_collisions, + std::vector sg1, + std::vector sg2, + double coulomb_log, + bool intra_collisions, + int debug_every, + CollisionalIonization *ionization, + CollisionalNuclearReaction *nuclear_reaction, + std::string fname + ) : Collisions( + params, + n_collisions, + sg1, + sg2, + coulomb_log, + intra_collisions, + debug_every, + ionization, + nuclear_reaction, + fname + ) {} ; + //! Cloning Constructor - CollisionsSingle( Collisions *coll, int ndim ) : Collisions( coll, ndim ) {}; + CollisionsSingle( Collisions *coll ) : Collisions( coll ) {}; //! destructor ~CollisionsSingle() {}; diff --git a/src/Patch/Patch.cpp b/src/Patch/Patch.cpp index f83000cd1..5a4fe2c65 100755 --- a/src/Patch/Patch.cpp +++ b/src/Patch/Patch.cpp @@ -192,7 +192,7 @@ void Patch::finishCloning( Patch *patch, Params ¶ms, SmileiMPI *smpi, unsign EMfields = ElectroMagnFactory::clone( patch->EMfields, params, vecSpecies, this, n_moved ); // clone the collisions - vecCollisions = CollisionsFactory::clone( patch->vecCollisions, params ); + vecCollisions = CollisionsFactory::clone( patch->vecCollisions ); // Clone the particle injector particle_injector_vector_ = ParticleInjectorFactory::cloneVector( patch->particle_injector_vector_, params, patch); diff --git a/src/Patch/VectorPatch.cpp b/src/Patch/VectorPatch.cpp index 3082dcf35..68b3cfe40 100755 --- a/src/Patch/VectorPatch.cpp +++ b/src/Patch/VectorPatch.cpp @@ -3424,20 +3424,22 @@ void VectorPatch::applyCollisions( Params ¶ms, int itime, Timers &timers ) { timers.collisions.restart(); - if( Collisions::debye_length_required ) + if( Collisions::debye_length_required ) { #pragma omp for schedule(runtime) for( unsigned int ipatch=0 ; ipatchvecCollisions.size(); - + #pragma omp for schedule(runtime) - for( unsigned int ipatch=0 ; ipatchvecCollisions[icoll]->collide( params, patches_[ipatch], itime, localDiags ); } - + } + #pragma omp single for( unsigned int icoll=0 ; icollpush_back( ( *double_prop[iprop] )[ipart] ); } @@ -331,7 +343,7 @@ void Particles::cp_particle_safe( unsigned int ipart, Particles &dest_parts ) dest_parts.short_prop[iprop]->push_back( ( *short_prop[iprop] )[ipart] ); } - for( unsigned int iprop=0 ; iproppush_back( ( *uint64_prop[iprop] )[ipart] ); } } diff --git a/src/Species/Particles.h b/src/Species/Particles.h index 27638384c..cf3ae0e1a 100755 --- a/src/Species/Particles.h +++ b/src/Species/Particles.h @@ -38,6 +38,9 @@ class Particles //! Set capacity of Particles vectors void reserve( unsigned int n_part_max, unsigned int nDim ); + + //! Initialize like another particle, but only reserve space + void initialize_reserve( unsigned int n_part_max, Particles &part ); //! Resize Particles vectors void resize( unsigned int nParticles, unsigned int nDim ); diff --git a/src/Species/Species.cpp b/src/Species/Species.cpp index 5fc2a4f76..04e8ee996 100755 --- a/src/Species/Species.cpp +++ b/src/Species/Species.cpp @@ -1456,4 +1456,29 @@ void Species::check( Patch *patch, std::string title ) << " - pz: " << sum_pz << " - ck: " << sum_ck << '\n'; -}; +} + +void Species::eraseWeightlessParticles() +{ + unsigned int nbins = first_index.size(); + unsigned int i = 0, available_i = 0; + + // Loop all particles, bin per bin + // Overwrite over earlier particles to erase them + for( unsigned int ibin = 0; ibin < nbins; ibin++ ) { + first_index[ibin] = available_i; + while( i < (unsigned int) last_index[ibin] ) { + if( particles->weight(i) > 0. ) { + if( i > available_i ) { + particles->overwrite_part( i, available_i ); + } + available_i ++; + } + i++; + } + last_index[ibin] = available_i; + } + + // Remove trailing particles + particles->erase_particle_trail( available_i ); +} diff --git a/src/Species/Species.h b/src/Species/Species.h index 759e852b8..54d3f748b 100755 --- a/src/Species/Species.h +++ b/src/Species/Species.h @@ -568,7 +568,9 @@ class Species return s_gamma; } - + + //! Erase all particles with zero weight + void eraseWeightlessParticles(); protected: diff --git a/src/Species/SpeciesFactory.h b/src/Species/SpeciesFactory.h index 6d3f2bf64..1dc194f93 100755 --- a/src/Species/SpeciesFactory.h +++ b/src/Species/SpeciesFactory.h @@ -82,7 +82,7 @@ class SpeciesFactory } } else { // Cancelation of the letter case for `radiation_model` - std::transform( radiation_model.begin(), radiation_model.end(), radiation_model.begin(), tolower ); + std::transform( radiation_model.begin(), radiation_model.end(), radiation_model.begin(), ::tolower ); // Name simplification if( radiation_model=="monte-carlo" ) { @@ -308,7 +308,7 @@ class SpeciesFactory // Cancelation of the letter case for `merging_method_` std::transform( this_species->merging_method_.begin(), this_species->merging_method_.end(), - this_species->merging_method_.begin(), tolower ); + this_species->merging_method_.begin(), ::tolower ); if( ( this_species->merging_method_ != "vranic_spherical" ) && ( this_species->merging_method_ != "vranic_cartesian" ) && @@ -1105,97 +1105,77 @@ class SpeciesFactory } } - // Loop species to find the electron species for ionizable species + // Loop species to find related species for( unsigned int ispec1 = 0; ispec1Ionize ) { - continue; - } - - // Loop all other species - for( unsigned int ispec2 = 0; ispec2ionization_electrons == returned_species[ispec2]->name_ ) { - if( ispec1==ispec2 ) { - ERROR( "For species '"<name_<<"' ionization_electrons must be a distinct species" ); - } - if( returned_species[ispec2]->mass_!=1 ) { - ERROR( "For species '"<name_<<"' ionization_electrons must be a species with mass==1" ); - } - returned_species[ispec1]->electron_species_index = ispec2; - returned_species[ispec1]->electron_species = returned_species[ispec2]; - returned_species[ispec1]->Ionize->new_electrons.tracked = returned_species[ispec1]->electron_species->particles->tracked; - returned_species[ispec1]->Ionize->new_electrons.isQuantumParameter = returned_species[ispec1]->electron_species->particles->isQuantumParameter; - returned_species[ispec1]->Ionize->new_electrons.isMonteCarlo = returned_species[ispec1]->electron_species->particles->isMonteCarlo; - returned_species[ispec1]->Ionize->new_electrons.initialize( 0, params.nDim_particle ); - if( ( !returned_species[ispec1]->getNbrOfParticles() ) && ( !returned_species[ispec2]->getNbrOfParticles() ) ) { - if( returned_species[ispec1]->atomic_number_!=0 ) { - int max_eon_number = returned_species[ispec1]->getNbrOfParticles() * returned_species[ispec1]->atomic_number_; - returned_species[ispec2]->particles->reserve( max_eon_number, returned_species[ispec2]->particles->dimension() ); - } else if( returned_species[ispec1]->maximum_charge_state_!=0 ) { - int max_eon_number = returned_species[ispec1]->getNbrOfParticles() * returned_species[ispec1]->maximum_charge_state_; - returned_species[ispec2]->particles->reserve( max_eon_number, returned_species[ispec2]->particles->dimension() ); - } - } - break; - } - } - if( returned_species[ispec1]->electron_species_index==-1 ) { - ERROR( "For species '"<name_<<"' ionization_electrons named " << returned_species[ispec1]->ionization_electrons << " could not be found" ); - } - } - - // Loop species to find the photon species for radiation species - for( unsigned int ispec1 = 0; ispec1Radiate ) { - continue; - } - - // No emission of discrete photon, only scalar diagnostics are updated - if( returned_species[ispec1]->radiation_photon_species.empty() ) { - returned_species[ispec1]->photon_species_index = -1; - returned_species[ispec1]->photon_species = NULL; - } - // Else, there will be emission of macro-photons. - else { - unsigned int ispec2 = 0; - for( ispec2 = 0; ispec2radiation_photon_species == returned_species[ispec2]->name_ ) { + + // Ionizable species + if( returned_species[ispec1]->Ionize ) { + // Loop all other species + for( unsigned int ispec2 = 0; ispec2ionization_electrons == returned_species[ispec2]->name_ ) { if( ispec1==ispec2 ) { - ERROR( "For species '"<name_<<"' radiation_photon_species must be a distinct photon species" ); + ERROR( "For species '"<name_<<"' ionization_electrons must be a distinct species" ); } - if( returned_species[ispec2]->mass_!=0 ) { - ERROR( "For species '"<name_<<"' radiation_photon_species must be a photon species with mass==0" ); + if( returned_species[ispec2]->mass_!=1 ) { + ERROR( "For species '"<name_<<"' ionization_electrons must be a species with mass==1" ); } - returned_species[ispec1]->photon_species_index = ispec2; - returned_species[ispec1]->photon_species = returned_species[ispec2]; - returned_species[ispec1]->Radiate->new_photons_.tracked = returned_species[ispec1]->photon_species->particles->tracked; - returned_species[ispec1]->Radiate->new_photons_.isQuantumParameter = returned_species[ispec1]->photon_species->particles->isQuantumParameter; - returned_species[ispec1]->Radiate->new_photons_.isMonteCarlo = returned_species[ispec1]->photon_species->particles->isMonteCarlo; - returned_species[ispec1]->Radiate->new_photons_.initialize( 0, - params.nDim_particle ); - //returned_species[ispec1]->Radiate->new_photons_.initialize(returned_species[ispec1]->getNbrOfParticles(), - // params.nDim_particle ); - returned_species[ispec2]->particles->reserve( returned_species[ispec1]->getNbrOfParticles(), - returned_species[ispec2]->particles->dimension() ); + returned_species[ispec1]->electron_species_index = ispec2; + returned_species[ispec1]->electron_species = returned_species[ispec2]; + + int max_eon_number = + returned_species[ispec1]->getNbrOfParticles() + * ( returned_species[ispec1]->atomic_number_ || returned_species[ispec1]->maximum_charge_state_ ); + returned_species[ispec1]->Ionize->new_electrons.initialize_reserve( + max_eon_number, *returned_species[ispec1]->electron_species->particles + ); break; } } - if( ispec2 == returned_species.size() ) { - ERROR( "Species '" << returned_species[ispec1]->radiation_photon_species << "' does not exist." ) + if( returned_species[ispec1]->electron_species_index==-1 ) { + ERROR( "For species '"<name_<<"' ionization_electrons named " << returned_species[ispec1]->ionization_electrons << " could not be found" ); } } - } - - // Loop species to find the electron and positron - // species for multiphoton Breit-wheeler - for( unsigned int ispec1 = 0; ispec1Multiphoton_Breit_Wheeler_process ) { - continue; - } else { + + // Radiating species + if( returned_species[ispec1]->Radiate ) { + // No emission of discrete photon, only scalar diagnostics are updated + if( returned_species[ispec1]->radiation_photon_species.empty() ) { + returned_species[ispec1]->photon_species_index = -1; + returned_species[ispec1]->photon_species = NULL; + } + // Else, there will be emission of macro-photons. + else { + unsigned int ispec2 = 0; + for( ispec2 = 0; ispec2radiation_photon_species == returned_species[ispec2]->name_ ) { + if( ispec1==ispec2 ) { + ERROR( "For species '"<name_<<"' radiation_photon_species must be a distinct photon species" ); + } + if( returned_species[ispec2]->mass_!=0 ) { + ERROR( "For species '"<name_<<"' radiation_photon_species must be a photon species with mass==0" ); + } + returned_species[ispec1]->photon_species_index = ispec2; + returned_species[ispec1]->photon_species = returned_species[ispec2]; + returned_species[ispec1]->Radiate->new_photons_.initialize_reserve( + returned_species[ispec1]->getNbrOfParticles(), + *returned_species[ispec1]->photon_species->particles + ); + break; + } + } + if( ispec2 == returned_species.size() ) { + ERROR( "Species '" << returned_species[ispec1]->radiation_photon_species << "' does not exist." ) + } + } + } + + // Breit-Wheeler species + if( returned_species[ispec1]->Multiphoton_Breit_Wheeler_process ) { unsigned int ispec2; for( int k=0; k<2; k++ ) { ispec2 = 0; while( ispec2multiphoton_Breit_Wheeler_[k] == returned_species[ispec2]->name_ ) { if( ispec1==ispec2 ) { ERROR( "For species '" << returned_species[ispec1]->name_ @@ -1207,13 +1187,10 @@ class SpeciesFactory } returned_species[ispec1]->mBW_pair_species_index[k] = ispec2; returned_species[ispec1]->mBW_pair_species[k] = returned_species[ispec2]; - returned_species[ispec1]->Multiphoton_Breit_Wheeler_process->new_pair[k].tracked = returned_species[ispec1]->mBW_pair_species[k]->particles->tracked; - returned_species[ispec1]->Multiphoton_Breit_Wheeler_process->new_pair[k].isQuantumParameter = returned_species[ispec1]->mBW_pair_species[k]->particles->isQuantumParameter; - returned_species[ispec1]->Multiphoton_Breit_Wheeler_process->new_pair[k].isMonteCarlo = returned_species[ispec1]->mBW_pair_species[k]->particles->isMonteCarlo; - returned_species[ispec1]->Multiphoton_Breit_Wheeler_process->new_pair[k].initialize( 0, - params.nDim_particle ); - returned_species[ispec2]->particles->reserve( returned_species[ispec1]->getNbrOfParticles(), - returned_species[ispec2]->particles->dimension() ); + returned_species[ispec1]->Multiphoton_Breit_Wheeler_process->new_pair[k].initialize_reserve( + returned_species[ispec1]->getNbrOfParticles(), + *returned_species[ispec1]->mBW_pair_species[k]->particles + ); ispec2 = returned_species.size() + 1; } ispec2++ ; diff --git a/validation/analyses/validate_tst_collisions4_DD_fusion_rate.py b/validation/analyses/validate_tst_collisions4_DD_fusion_rate.py new file mode 100644 index 000000000..a037150da --- /dev/null +++ b/validation/analyses/validate_tst_collisions4_DD_fusion_rate.py @@ -0,0 +1,52 @@ +import happi +from numpy import array, log10 + +S = happi.Open(["./restart*"], verbose=False) + +v = [] +N = [] +Nref = [] +Npart = [] +for i in range(0,len(S.namelist.DiagParticleBinning),2): + v += [S.namelist.Species["Da_%d"%(i//2)].mean_velocity[0]] + N += [S.ParticleBinning(i ,sum={"x":"all"}).getData()[-1]] + Nref += [S.ParticleBinning(i+1,sum={"x":"all"}).getData()[-1]] + Npart += [S.Scalar("Ntot_He_%d"%(i//2)).getData()[-1]] + +Validate("log10 of created helium3 density", log10(array(N)+1e-40), 0.2) +Validate("number of created helium3", Npart, 300.) + + +# # Some ploting to compare w theory +# E_theory = array([ +# 1.00000000e+03, 1.49328996e+03, 2.22991489e+03, 3.32990952e+03, +# 4.97252044e+03, 7.42541482e+03, 1.10882974e+04, 1.65580431e+04, +# 2.47259595e+04, 3.69230270e+04, 5.51367853e+04, 8.23352077e+04, +# 1.22950339e+05, 1.83600506e+05, 2.74168792e+05, 4.09413503e+05, +# 6.11373072e+05, 9.12957268e+05, 1.36330992e+06, 2.03581701e+06, +# 3.04006509e+06, 4.53969867e+06, 6.77908643e+06, 1.01231417e+07, +# 1.51167858e+07, 2.25737444e+07, 3.37091458e+07, 5.03375289e+07, +# 7.51685263e+07, 1.12248405e+08, 1.67619416e+08, 2.50304391e+08, +# 3.73777033e+08, 5.58157490e+08 +# ]) * 2.0141 # in eV +# crosssection_theory = array([ +# 1.37959339e-09, 2.71808933e-07, 1.93984386e-05, 5.96656956e-04, +# 9.09434767e-03, 7.80510801e-02, 4.23198290e-01, 1.59539840e+00, +# 4.50717898e+00, 1.01191544e+01, 1.89534351e+01, 3.08794670e+01, +# 4.52831325e+01, 6.12170894e+01, 7.72386840e+01, 9.12459059e+01, +# 1.00937123e+02, 1.04962615e+02, 1.03790887e+02, 9.92555901e+01, +# 9.30691435e+01, 8.55168643e+01, 7.52550684e+01, 6.07107292e+01, +# 4.26722519e+01, 2.52075010e+01, 1.24955923e+01, 5.39117181e+00, +# 2.13561911e+00, 8.02190038e-01, 2.81582182e-01, 8.91232341e-02, +# 2.57170089e-02, 6.56861761e-03 +# ]) * 1e-31 # in m^2 +# mass_deuterium = 3671.46*511000. +# Wr = S.namelist.Main.reference_angular_frequency_SI # in Hz +# density = S.namelist.Species[0].number_density * 0.0003142 * Wr**2 # in m^-3 +# time = S.namelist.Main.simulation_time / Wr # in s +# vrel = 3e8 * sqrt(1. - (E_theory/mass_deuterium + 1.)**-2) # in m/s +# N_theory = vrel * density * crosssection_theory * time +# E = (1./sqrt(1.-array(v)**2)-1.) * mass_deuterium +# clf() +# loglog(E_theory,N_theory,'-') +# loglog(E, array(N) / array(Nref), '+') diff --git a/validation/references/tst_collisions4_DD_fusion_rate.py.txt b/validation/references/tst_collisions4_DD_fusion_rate.py.txt new file mode 100644 index 000000000..751424f60 Binary files /dev/null and b/validation/references/tst_collisions4_DD_fusion_rate.py.txt differ