diff --git a/doc/Sphinx/namelist.rst b/doc/Sphinx/namelist.rst index 8ad7d54e8..39b9aafb1 100644 --- a/doc/Sphinx/namelist.rst +++ b/doc/Sphinx/namelist.rst @@ -180,10 +180,10 @@ The block ``Main`` is **mandatory** and has the following syntax:: .. py:data:: patch_arrangement :default: ``"hilbertian"`` - + Determines the ordering of patches and the way they are separated into the various MPI processes. Options are: - + * ``"hilbertian"``: following the Hilbert curve (see :ref:`this explanation`). * ``"linearized_XY"`` in 2D or ``"linearized_XYZ"`` in 3D: following the row-major (C-style) ordering. @@ -1731,13 +1731,13 @@ tables. :default: 128 - Discretization of the *chimin* and *xip* tables in the *chipa* direction. + Discretization of the *chimin* and *xip* tables in the *particle_chi* direction. .. py:data:: xip_chiph_dim :default: 128 - Discretization of the *xip* tables in the *chiph* direction. + Discretization of the *xip* tables in the *photon_chi* direction. .. py:data:: output_format @@ -1749,7 +1749,7 @@ tables. :default: 1e-3 - Threshold on the particle quantum parameter *chipa*. When a particle has a + Threshold on the particle quantum parameter *particle_chi*. When a particle has a quantum parameter below this threshold, radiation reaction is not taken into account. @@ -1757,7 +1757,7 @@ tables. :default: 1e-2 - Threshold on the particle quantum parameter *chipa* between the continuous + Threshold on the particle quantum parameter *particle_chi* between the continuous and the discontinuous radiation model. .. py:data:: table_path @@ -1871,13 +1871,13 @@ There are two tables used for the multiphoton Breit-Wheeler refers to as the :default: 128 - Discretization of the *chimin* and *xip* tables in the *chiph* direction. + Discretization of the *chimin* and *xip* tables in the *photon_chi* direction. .. py:data:: xip_chipa_dim :default: 128 - Discretization of the *xip* tables in the *chipa* direction. + Discretization of the *xip* tables in the *particle_chi* direction. -------------------------------------------------------------------------------- @@ -2173,16 +2173,16 @@ To add one probe diagnostic, include the block ``DiagProbe``:: .. py:data:: fields :default: ``[]``, which means ``["Ex", "Ey", "Ez", "Bx", "By", "Bz", "Jx", "Jy", "Jz", "Rho"]`` - + A list of fields among ``"Ex"``, ``"Ey"``, ``"Ez"``, ``"Bx"``, ``"By"``, ``"Bz"``, ``"Jx"``, ``"Jy"``, ``"Jz"`` and ``"Rho"``. Only listed fields will be saved although they are all calculated. - + The contributions of each species to the currents and the density are also available, although they are not included by default. They may be added to the list as ``"Jx_abc"``, ``"Jy_abc"``, ``"Jz_abc"`` or ``"Rho_abc"``, where ``abc`` is the species name. - + In the case of an envelope model for the laser (see :doc:`laser_envelope`), the following fields are also available: ``"Env_A_abs"``, ``"Env_Chi"``, ``"Env_E_abs"``. diff --git a/src/MultiphotonBreitWheeler/MultiphotonBreitWheeler.cpp b/src/MultiphotonBreitWheeler/MultiphotonBreitWheeler.cpp index dd2b8fd24..4fb76df47 100644 --- a/src/MultiphotonBreitWheeler/MultiphotonBreitWheeler.cpp +++ b/src/MultiphotonBreitWheeler/MultiphotonBreitWheeler.cpp @@ -19,17 +19,17 @@ MultiphotonBreitWheeler::MultiphotonBreitWheeler(Params& params, Species * species) { // Dimension position - nDim_ = params.nDim_particle; + n_dimensions_ = params.nDim_particle; // Time step dt = params.timestep; // Normalized Schwinger Electric Field - norm_E_Schwinger = params.electron_mass*params.c_vacuum*params.c_vacuum + norm_E_Schwinger_ = params.electron_mass*params.c_vacuum*params.c_vacuum / (params.red_planck_cst*params.reference_angular_frequency_SI); - // Inverse of norm_E_Schwinger - inv_norm_E_Schwinger = 1./norm_E_Schwinger; + // Inverse of norm_E_Schwinger_ + inv_norm_E_Schwinger_ = 1./norm_E_Schwinger_; // Number of positrons and electrons generated per event this->mBW_pair_creation_sampling[0] = species->mBW_pair_creation_sampling[0]; @@ -161,7 +161,7 @@ void MultiphotonBreitWheeler::operator() (Particles &particles, // for now particles could be created outside of the local domain // without been subject do boundary conditions (including domain exchange) //double* position[3]; - //for ( int i = 0 ; i chiph_threshold, - // else chiph is too low to induce a decay - if (((*gamma)[ipart] > 2.) && (chiph[ipart] > chiph_threashold)) + // We also check that photon_chi > chiph_threshold, + // else photon_chi is too low to induce a decay + if (((*gamma)[ipart] > 2.) && (photon_chi[ipart] > chiph_threashold)) { // Init local variables event_time = 0; // New even // If tau[ipart] <= 0, this is a new process - if (tau[ipart] <= epsilon_tau) + if (tau[ipart] <= epsilon_tau_) { // New final optical depth to reach for emision - while (tau[ipart] <= epsilon_tau) + while (tau[ipart] <= epsilon_tau_) tau[ipart] = -log(1.-Rand::uniform()); } // Photon decay: emission under progress - // If epsilon_tau > 0 - else if (tau[ipart] > epsilon_tau) + // If epsilon_tau_ > 0 + else if (tau[ipart] > epsilon_tau_) { // from the cross section - temp = MultiphotonBreitWheelerTables.compute_dNBWdt(chiph[ipart],(*gamma)[ipart]); + temp = MultiphotonBreitWheelerTables.compute_dNBWdt(photon_chi[ipart],(*gamma)[ipart]); // Time to decay // If this time is above the remaining iteration time, @@ -241,17 +241,17 @@ void MultiphotonBreitWheeler::operator() (Particles &particles, // If the final optical depth is reached // The photon decays into pairs - if (tau[ipart] <= epsilon_tau) + if (tau[ipart] <= epsilon_tau_) { // Update of the position // Move the photons //#ifdef __DEBUG -// for ( int i = 0 ; i= T_dim-1) { ichiph = T_dim-2; - dNBWdt = 0.38*pow(chiph,-1./3.); + dNBWdt = 0.38*pow(photon_chi,-1./3.); } else { @@ -189,39 +189,39 @@ double MultiphotonBreitWheelerTables::compute_dNBWdt(double chiph, double gamma) dNBWdt = (T_table[ichiph+1]*fabs(logchiph-logchiphm) + T_table[ichiph]*fabs(logchiphp - logchiph))*T_chiph_inv_delta; } - return factor_dNBWdt*dNBWdt*chiph/gamma; + return factor_dNBWdt*dNBWdt*photon_chi/gamma; } // ----------------------------------------------------------------------------- -//! Computation of the value T(chiph) using the approximated +//! Computation of the value T(photon_chi) using the approximated //! formula of Erber -//! \param chiph photon quantum parameter +//! \param photon_chi photon quantum parameter //! \param nbit number of iteration for the Bessel evaluation //! \param eps epsilon for the Bessel evaluation // ----------------------------------------------------------------------------- -double MultiphotonBreitWheelerTables::compute_Erber_T(double chiph,int nbit, +double MultiphotonBreitWheelerTables::compute_Erber_T(double photon_chi,int nbit, double eps) { // Values for Bessel results //double I,dI; double K; - //userFunctions::modified_bessel_IK(1./3.,4./(3.*chiph),I,dI,K,dK,nbit,eps); - K = userFunctions::modified_bessel_K(1./3.,4./(3.*chiph),nbit,eps, false); + //userFunctions::modified_bessel_IK(1./3.,4./(3.*photon_chi),I,dI,K,dK,nbit,eps); + K = userFunctions::modified_bessel_K(1./3.,4./(3.*photon_chi),nbit,eps, false); - return 0.16*K*K/chiph; + return 0.16*K*K/photon_chi; } // ----------------------------------------------------------------------------- -//! Computation of the value of the integration of dT(chiph)/dhicpa +//! Computation of the value of the integration of dT(photon_chi)/dhicpa //! using the formula of Ritus. -//! \param chiph photon quantum parameter -//! \param chipa particle quantum parameter for integration (=0.5*chiph for full integration) +//! \param photon_chi photon quantum parameter +//! \param particle_chi particle quantum parameter for integration (=0.5*photon_chi for full integration) //! \param nbit number of iteration for the Gauss-Legendre integration //! \param eps epsilon for the Bessel evaluation // ----------------------------------------------------------------------------- -double MultiphotonBreitWheelerTables::compute_integration_Ritus_dTdchi(double chiph, - double chipa, +double MultiphotonBreitWheelerTables::compute_integration_Ritus_dTdchi(double photon_chi, + double particle_chi, int nbit, double eps) { @@ -235,14 +235,14 @@ double MultiphotonBreitWheelerTables::compute_integration_Ritus_dTdchi(double ch double T; // gauss Legendre coefficients - userFunctions::gauss_legendre_coef(log10(chipa)-50.,log10(chipa), + userFunctions::gauss_legendre_coef(log10(particle_chi)-50.,log10(particle_chi), gauleg_x, gauleg_w, nbit, eps); T = 0; for(i=0 ; i< nbit ; i++) { u = pow(10.,gauleg_x[i]); - T += gauleg_w[i]*compute_Ritus_dTdchi(chiph,u,nbit,eps)*u*log(10.); + T += gauleg_w[i]*compute_Ritus_dTdchi(photon_chi,u,nbit,eps)*u*log(10.); } return T; @@ -252,9 +252,9 @@ double MultiphotonBreitWheelerTables::compute_integration_Ritus_dTdchi(double ch //! Computation of the electron and positron quantum parameters for //! the multiphoton Breit-Wheeler pair creation // -//! \param chiph photon quantum parameter +//! \param photon_chi photon quantum parameter // ----------------------------------------------------------------------------- -double * MultiphotonBreitWheelerTables::compute_pair_chi(double chiph) +double * MultiphotonBreitWheelerTables::compute_pair_chi(double photon_chi) { // Parameters double * chi = new double[2]; @@ -268,30 +268,30 @@ double * MultiphotonBreitWheelerTables::compute_pair_chi(double chiph) int ixip; // ----------------------------------------------------------- - // Computation of the index associated to the given chiph + // Computation of the index associated to the given photon_chi // ----------------------------------------------------------- - logchiph = log10(chiph); + logchiph = log10(photon_chi); // Lower boundary of the table - if (chiph < xip_chiph_min) + if (photon_chi < xip_chiph_min) { ichiph = 0; } // Upper boundary of the table - else if (chiph >= xip_chiph_max) + else if (photon_chi >= xip_chiph_max) { ichiph = xip_chiph_dim-1; } // Inside the table else { - // Use floor so that chiph corresponding to ichiph is <= given chiph + // Use floor so that photon_chi corresponding to ichiph is <= given photon_chi ichiph = int(floor((logchiph-xip_log10_chiph_min)*(xip_chiph_inv_delta))); } // --------------------------------------- - // Search of the index ichiph for chiph + // Search of the index ichiph for photon_chi // --------------------------------------- // First, we compute a random xip in [0,1[ @@ -326,8 +326,8 @@ double * MultiphotonBreitWheelerTables::compute_pair_chi(double chiph) &xip_table[ichiph*xip_chipa_dim],xipp,xip_chipa_dim); } - // Delta for the chipa dimension - delta_chipa = (log10(0.5*chiph)-xip_chipamin_table[ichiph]) + // Delta for the particle_chi dimension + delta_chipa = (log10(0.5*photon_chi)-xip_chipamin_table[ichiph]) * xip_inv_chipa_dim_minus_one; ixip = ichiph*xip_chipa_dim + ichipa; @@ -345,7 +345,7 @@ double * MultiphotonBreitWheelerTables::compute_pair_chi(double chiph) chi[1] = pow(10,log10_chipam*(1.0-d) + log10_chipap*(d)); // Electron quantum parameter - chi[0] = chiph - chi[1]; + chi[0] = photon_chi - chi[1]; } // If xip <= 0.5, the positron will bring more energy than the electron else @@ -354,21 +354,21 @@ double * MultiphotonBreitWheelerTables::compute_pair_chi(double chiph) chi[0] = pow(10,log10_chipam*(1.0-d) + log10_chipap*(d)); // Positron quantum parameter - chi[1] = chiph - chi[0]; + chi[1] = photon_chi - chi[0]; } return chi; } // ----------------------------------------------------------------------------- -//! Computation of the value dT/dchipa(chiph) using the formula of Ritus -//! \param chiph photon quantum parameter -//! \param chipa particle quantum parameter +//! Computation of the value dT/dchipa(photon_chi) using the formula of Ritus +//! \param photon_chi photon quantum parameter +//! \param particle_chi particle quantum parameter //! \param nbit number of iteration for the Gauss-Legendre integration //! \param eps epsilon for the Bessel evaluation // ----------------------------------------------------------------------------- -double MultiphotonBreitWheelerTables::compute_Ritus_dTdchi(double chiph, - double chipa,int nbit,double eps) +double MultiphotonBreitWheelerTables::compute_Ritus_dTdchi(double photon_chi, + double particle_chi,int nbit,double eps) { // Parameters: @@ -381,12 +381,12 @@ double MultiphotonBreitWheelerTables::compute_Ritus_dTdchi(double chiph, //double I,dI; //double K,dK; - y = (chiph/(3.*chipa*(chiph-chipa))); + y = (photon_chi/(3.*particle_chi*(photon_chi-particle_chi))); //userFunctions::modified_bessel_IK(2./3.,2.*y,I,dI,K,dK,500,1e-16); - //p1 = (2. - 3.*chiph*y)*K; + //p1 = (2. - 3.*photon_chi*y)*K; - p1 = (2. - 3.*chiph*y)*userFunctions::modified_bessel_K(2./3.,2*y,500,1e-16,false); + p1 = (2. - 3.*photon_chi*y)*userFunctions::modified_bessel_K(2./3.,2*y,500,1e-16,false); // gauss Legendre coefficients userFunctions::gauss_legendre_coef(log10(2*y),log10(2*y)+50., @@ -403,7 +403,7 @@ double MultiphotonBreitWheelerTables::compute_Ritus_dTdchi(double chiph, //p2 += gauleg_w[i]*K*u*log(10.); } - return (p2 - p1)/(M_PI*sqrt(3.)*pow(chiph,2)); + return (p2 - p1)/(M_PI*sqrt(3.)*pow(photon_chi,2)); } // ----------------------------------------------------------------------------- @@ -438,7 +438,7 @@ void MultiphotonBreitWheelerTables::compute_T_table(SmileiMPI *smpi) { // Temporary photon chi value - double chiph; + double photon_chi; // For percentages double pct = 0.; double dpct = 10.; @@ -494,11 +494,11 @@ void MultiphotonBreitWheelerTables::compute_T_table(SmileiMPI *smpi) // Loop over the table values for(int i = 0 ; i < length_table[rank] ; i++) { - chiph = pow(10.,(imin_table[rank] + i)*T_chiph_delta + photon_chi = pow(10.,(imin_table[rank] + i)*T_chiph_delta + T_log10_chiph_min); - buffer[i] = 2.*MultiphotonBreitWheelerTables::compute_integration_Ritus_dTdchi(chiph,0.5*chiph,200,1e-15); - //buffer[i] = MultiphotonBreitWheelerTables::compute_Erber_T(chiph,500000,1e-10); + buffer[i] = 2.*MultiphotonBreitWheelerTables::compute_integration_Ritus_dTdchi(photon_chi,0.5*photon_chi,200,1e-15); + //buffer[i] = MultiphotonBreitWheelerTables::compute_Erber_T(photon_chi,500000,1e-10); if (100.*i >= length_table[rank]*pct) { @@ -533,7 +533,7 @@ void MultiphotonBreitWheelerTables::compute_T_table(SmileiMPI *smpi) //! Computation of the minimum particle quantum parameter chipamin //! for the photon xip array and computation of the photon xip array. // -//! \details Under the minimum chipa value, the particle kinetic energy is +//! \details Under the minimum particle_chi value, the particle kinetic energy is //! considered negligible. All energy goes to the other. // //! \param smpi Object of class SmileiMPI containing MPI properties @@ -562,9 +562,9 @@ void MultiphotonBreitWheelerTables::compute_xip_table(SmileiMPI *smpi) if (!table_exists) { // Parameters: - double chipa; // Temporary particle chi value - double chiph; // Temporary photon chi value - double chipa_delta; // Temporary delta for chiph + double particle_chi; // Temporary particle chi value + double photon_chi; // Temporary photon chi value + double chipa_delta; // Temporary delta for photon_chi double logchipamin; // Temporary log10 of photon chi value double xip; // Temporary xip double numerator; @@ -637,21 +637,21 @@ void MultiphotonBreitWheelerTables::compute_xip_table(SmileiMPI *smpi) logchipamin = (imin_table[rank] + ichiph)*xip_chiph_delta + xip_log10_chiph_min; - chiph = pow(10.,logchipamin); + photon_chi = pow(10.,logchipamin); - logchipamin = log10(0.5*chiph); + logchipamin = log10(0.5*photon_chi); // Denominator of xip - denominator = MultiphotonBreitWheelerTables::compute_integration_Ritus_dTdchi(chiph, - 0.5*chiph,200,1e-15); + denominator = MultiphotonBreitWheelerTables::compute_integration_Ritus_dTdchi(photon_chi, + 0.5*photon_chi,200,1e-15); k = 0; while(k < xip_power) { logchipamin -= pow(0.1,k); - chipa = pow(10.,logchipamin); - numerator = MultiphotonBreitWheelerTables::compute_integration_Ritus_dTdchi(chiph, - chipa,200,1e-15); + particle_chi = pow(10.,logchipamin); + numerator = MultiphotonBreitWheelerTables::compute_integration_Ritus_dTdchi(photon_chi, + particle_chi,200,1e-15); if (numerator == 0 || denominator == 0) { @@ -691,37 +691,37 @@ void MultiphotonBreitWheelerTables::compute_xip_table(SmileiMPI *smpi) // Allocation of the local buffer buffer = new double [length_table[rank]*xip_chipa_dim]; - // Loop for xip in the chiph dimension + // Loop for xip in the photon_chi dimension dpct = std::max(10.,100./(length_table[rank]*xip_chipa_dim)); pct = dpct; for(int ichiph = 0 ; ichiph < length_table[rank] ; ichiph++) { - // log10(chiph) for the current ichiph - chiph = (imin_table[rank] + ichiph)*xip_chiph_delta + // log10(photon_chi) for the current ichiph + photon_chi = (imin_table[rank] + ichiph)*xip_chiph_delta + xip_log10_chiph_min; - // chiph for the current ichiph - chiph = pow(10.,chiph); + // photon_chi for the current ichiph + photon_chi = pow(10.,photon_chi); - // Delta on chipa - chipa_delta = (log10(0.5*chiph) - xip_chipamin_table[imin_table[rank] + ichiph]) + // Delta on particle_chi + chipa_delta = (log10(0.5*photon_chi) - xip_chipamin_table[imin_table[rank] + ichiph]) / (xip_chipa_dim - 1); // Denominator of xip - denominator = MultiphotonBreitWheelerTables::compute_integration_Ritus_dTdchi(chiph, - 0.5*chiph,200,1e-15); + denominator = MultiphotonBreitWheelerTables::compute_integration_Ritus_dTdchi(photon_chi, + 0.5*photon_chi,200,1e-15); - // Loop in the chipa dimension + // Loop in the particle_chi dimension for (int ichipa = 0 ; ichipa < xip_chipa_dim ; ichipa ++) { - // Local chipa value - chipa = pow(10.,ichipa*chipa_delta + + // Local particle_chi value + particle_chi = pow(10.,ichipa*chipa_delta + xip_chipamin_table[imin_table[rank] + ichiph]); // Numerator of xip - numerator = MultiphotonBreitWheelerTables::compute_integration_Ritus_dTdchi(chiph, - chipa,200,1e-15); + numerator = MultiphotonBreitWheelerTables::compute_integration_Ritus_dTdchi(photon_chi, + particle_chi,200,1e-15); // Update local buffer value // buffer[ichiph*xip_chipa_dim + ichipa] = std::min(1.,numerator / 2*denominator); @@ -826,7 +826,7 @@ void MultiphotonBreitWheelerTables::output_T_table() file << "Multiphoton Breit-Wheeler T table\n"; - file << "Dimension chiph - chiph min - chiph max \n"; + file << "Dimension photon_chi - photon_chi min - photon_chi max \n"; file << T_dim; file << " " @@ -945,7 +945,7 @@ void MultiphotonBreitWheelerTables::output_xip_table() file << "Table xip_chipamin and xip for multiphoton Breit-Wheeler \n"; - file << "Dimension chiph - Dimension chipa - chiph min - chiph max \n"; + file << "Dimension photon_chi - Dimension particle_chi - photon_chi min - photon_chi max \n"; file << xip_chiph_dim << " " << xip_chipa_dim << " " @@ -1534,7 +1534,7 @@ void MultiphotonBreitWheelerTables::bcast_xip_table(SmileiMPI *smpi) // Inverse of delta xip_chiph_inv_delta = 1./xip_chiph_delta; - // Inverse chipa discetization (regularly used) + // Inverse particle_chi discetization (regularly used) xip_inv_chipa_dim_minus_one = 1./(xip_chipa_dim - 1.); } diff --git a/src/MultiphotonBreitWheeler/MultiphotonBreitWheelerTables.h b/src/MultiphotonBreitWheeler/MultiphotonBreitWheelerTables.h index 8a360d1d2..1886bb32e 100644 --- a/src/MultiphotonBreitWheeler/MultiphotonBreitWheelerTables.h +++ b/src/MultiphotonBreitWheeler/MultiphotonBreitWheelerTables.h @@ -48,39 +48,39 @@ class MultiphotonBreitWheelerTables // --------------------------------------------------------------------- //! Computation of the production rate of pairs per photon - //! \param chiph photon quantum parameter + //! \param photon_chi photon quantum parameter //! \param gamma photon normalized energy - double compute_dNBWdt(double chiph, double gamma); + double compute_dNBWdt(double photon_chi, double gamma); - //! Computation of the value T(chiph) using the approximated + //! Computation of the value T(photon_chi) using the approximated //! formula of Erber - //! \param chiph photon quantum parameter + //! \param photon_chi photon quantum parameter //! \param nbit number of iteration for the Bessel evaluation //! \param eps epsilon for the Bessel evaluation - double compute_Erber_T(double chiph,int nbit, + double compute_Erber_T(double photon_chi,int nbit, double eps); - //! Computation of the value T(chiph) using the formula of Ritus - //! \param chiph photon quantum parameter - //! \param chipa particle quantum parameter for integration (=0.5*chiph for full integration) + //! Computation of the value T(photon_chi) using the formula of Ritus + //! \param photon_chi photon quantum parameter + //! \param particle_chi particle quantum parameter for integration (=0.5*photon_chi for full integration) //! \param nbit number of iteration for the Gauss-Legendre integration //! \param eps epsilon for the Bessel evaluation - double compute_integration_Ritus_dTdchi(double chiph, - double chipa, + double compute_integration_Ritus_dTdchi(double photon_chi, + double particle_chi, int nbit, double eps); - //! Computation of the value T(chiph) using the formula of Ritus - //! \param chiph photon quantum parameter + //! Computation of the value T(photon_chi) using the formula of Ritus + //! \param photon_chi photon quantum parameter //! \param nbit number of iteration for the Gauss-Legendre integration //! \param eps epsilon for the Bessel evaluation - double compute_Ritus_dTdchi(double chiph, - double chipa,int nbit,double eps); + double compute_Ritus_dTdchi(double photon_chi, + double particle_chi,int nbit,double eps); //! Computation of the electron and positron quantum parameters for //! the multiphoton Breit-Wheeler pair creation - //! \param chiph photon quantum parameter - double * compute_pair_chi(double chiph); + //! \param photon_chi photon quantum parameter + double * compute_pair_chi(double photon_chi); // --------------------------------------------------------------------- // TABLE COMPUTATION @@ -94,7 +94,7 @@ class MultiphotonBreitWheelerTables //! Computation of the minimum particle quantum parameter chipamin //! for the photon xip array and computation of the photon xip array. - //! \details Under the minimum chipa value, the particle kinetic energy is + //! \details Under the minimum particle_chi value, the particle kinetic energy is //! considered negligible. All energy goes to the other. //! \param smpi Object of class SmileiMPI containing MPI properties void compute_xip_table(SmileiMPI *smpi); @@ -187,10 +187,10 @@ class MultiphotonBreitWheelerTables bool T_computed; // --------------------------------------------- - // Table chipa min for xip table + // Table particle_chi min for xip table // --------------------------------------------- - //! Table containing the chipa min values + //! Table containing the particle_chi min values //! Under this value, electron kinetic energy of the pair is //! considered negligible std::vector xip_chipamin_table; @@ -205,27 +205,27 @@ class MultiphotonBreitWheelerTables //! This enables to compute the energy repartition between the electron and the positron std::vector xip_table; - //! Minimum boundary for chiph in the table xip and xip_chipamin + //! Minimum boundary for photon_chi in the table xip and xip_chipamin double xip_chiph_min; - //! Logarithm of the minimum boundary for chiph in the table xip + //! Logarithm of the minimum boundary for photon_chi in the table xip //! and xip_chipamin double xip_log10_chiph_min; - //! Maximum boundary for chiph in the table xip and xip_chipamin + //! Maximum boundary for photon_chi in the table xip and xip_chipamin double xip_chiph_max; - //! Delta for the chiph discretization in the table xip and xip_chipamin + //! Delta for the photon_chi discretization in the table xip and xip_chipamin double xip_chiph_delta; - //! Inverse of the delta for the chiph discretization + //! Inverse of the delta for the photon_chi discretization //! in the table xip and xip_chipamin double xip_chiph_inv_delta; - //! Dimension of the discretized parameter chiph + //! Dimension of the discretized parameter photon_chi int xip_chiph_dim; - //! Dimension of the discretized parameter chipa + //! Dimension of the discretized parameter particle_chi int xip_chipa_dim; //! 1/(xip_chipa_dim - 1) diff --git a/src/Python/pyinit.py b/src/Python/pyinit.py index 41ff78600..9e4e37ae1 100644 --- a/src/Python/pyinit.py +++ b/src/Python/pyinit.py @@ -505,10 +505,10 @@ class RadiationReaction(SmileiComponent): xip_chiph_dim = 128 # Output format, can be "ascii", "binary", "hdf5" output_format = "hdf5" - # Threshold on chipa between the continuous and + # Threshold on particle_chi between the continuous and # the discontinuous approaches chipa_disc_min_threshold = 1e-2 - # Threshold on chipa: if chipa < 1E-3 no radiation reaction + # Threshold on particle_chi: if particle_chi < 1E-3 no radiation reaction chipa_radiation_threshold = 1e-3 # Path the tables/databases table_path = "./" diff --git a/src/Radiation/Radiation.cpp b/src/Radiation/Radiation.cpp index 5d02d18bd..247a95e14 100644 --- a/src/Radiation/Radiation.cpp +++ b/src/Radiation/Radiation.cpp @@ -16,24 +16,24 @@ // ----------------------------------------------------------------------------- Radiation::Radiation(Params& params, Species * species) { - // Dimension position - nDim_ = params.nDim_particle; + // Number of dimensions for the positions and momentums + n_dimensions_ = params.nDim_particle; // Time step - dt = params.timestep; + dt_ = params.timestep; // Inverse of the species mass one_over_mass_ = 1./species->mass; // Normalized Schwinger Electric Field - norm_E_Schwinger = params.electron_mass*params.c_vacuum*params.c_vacuum + norm_E_Schwinger_ = params.electron_mass*params.c_vacuum*params.c_vacuum / (params.red_planck_cst* params.reference_angular_frequency_SI); - // Inverse of norm_E_Schwinger - inv_norm_E_Schwinger = 1./norm_E_Schwinger; + // Inverse of norm_E_Schwinger_ + inv_norm_E_Schwinger_ = 1./norm_E_Schwinger_; // The thread radiated energy is initially null - radiated_energy = 0; + radiated_energy_ = 0; } // ----------------------------------------------------------------------------- @@ -52,7 +52,7 @@ Radiation::~Radiation() //! \param iend Index of the last particle //! \param ithread Thread index // ----------------------------------------------------------------------------- -void Radiation::compute_thread_chipa(Particles &particles, +void Radiation::computeParticlesChi(Particles &particles, SmileiMPI* smpi, int istart, int iend, @@ -105,7 +105,7 @@ void Radiation::compute_thread_chipa(Particles &particles, + momentum[2][ipart]*momentum[2][ipart]); // Computation of the Lorentz invariant quantum parameter - chi[ipart] = Radiation::compute_chipa(charge_over_mass2, + chi[ipart] = Radiation::computeParticleChi(charge_over_mass2, momentum[0][ipart],momentum[1][ipart],momentum[2][ipart], gamma, (*(Ex+ipart-ipart_ref)),(*(Ey+ipart-ipart_ref)),(*(Ez+ipart-ipart_ref)), diff --git a/src/Radiation/Radiation.h b/src/Radiation/Radiation.h index d45a484c6..10aa977c9 100644 --- a/src/Radiation/Radiation.h +++ b/src/Radiation/Radiation.h @@ -62,14 +62,14 @@ class Radiation //! \param By y component of the particle magnetic field //! \param Bz z component of the particle magnetic field //#pragma omp declare simd - double inline compute_chipa(double & charge_over_mass2, + double inline computeParticleChi(double & charge_over_mass2, double & px, double & py, double & pz, double & gamma, double & Ex, double & Ey, double & Ez, double & Bx, double & By, double & Bz) { - return fabs(charge_over_mass2)*inv_norm_E_Schwinger + return fabs(charge_over_mass2)*inv_norm_E_Schwinger_ * sqrt( fabs( pow(Ex*px + Ey*py + Ez*pz,2) - pow(gamma*Ex - By*pz + Bz*py,2) - pow(gamma*Ey - Bz*px + Bx*pz,2) @@ -79,14 +79,14 @@ class Radiation //! Return the total normalized radiated energy double inline getRadiatedEnergy() { - return radiated_energy; + return radiated_energy_; }; //! Set the total normalized radiated energy of the path //! \param value value of the radiated energy to be assigned void setRadiatedEnergy(double value) { - radiated_energy = value; + radiated_energy_ = value; }; //! Computation of the quantum parameter for the given @@ -96,7 +96,7 @@ class Radiation //! \param istart Index of the first particle //! \param iend Index of the last particle //! \param ithread Thread index - void compute_thread_chipa(Particles &particles, + void computeParticlesChi(Particles &particles, SmileiMPI* smpi, int istart, int iend, @@ -111,25 +111,25 @@ class Radiation // General parameters //! Dimension of position - int nDim_; + int n_dimensions_; - //! Inverse species mass + //! Inversed species mass double one_over_mass_; //! Time step - double dt; + double dt_; //! Radiated energy of the total thread - double radiated_energy; + double radiated_energy_; // _________________________________________ // Factors //! Normalized Schwinger Electric field - double norm_E_Schwinger; + double norm_E_Schwinger_; - //! Inverse Normalized Schwinger Electric field - double inv_norm_E_Schwinger; + //! Inversed Normalized Schwinger Electric field + double inv_norm_E_Schwinger_; private: diff --git a/src/Radiation/RadiationCorrLandauLifshitz.cpp b/src/Radiation/RadiationCorrLandauLifshitz.cpp index 78082a333..47d20a260 100644 --- a/src/Radiation/RadiationCorrLandauLifshitz.cpp +++ b/src/Radiation/RadiationCorrLandauLifshitz.cpp @@ -76,7 +76,7 @@ void RadiationCorrLandauLifshitz::operator() ( const double one_over_mass_2 = pow(one_over_mass_,2.); // Temporary quantum parameter - double chipa; + double particle_chi; // Temporary Lorentz factor double gamma; @@ -102,7 +102,7 @@ void RadiationCorrLandauLifshitz::operator() ( std::vector rad_norm_energy (iend-istart,0); // Reinitialize the cumulative radiated energy for the current thread - this->radiated_energy = 0.; + radiated_energy_ = 0.; // _______________________________________________________________ // Computation @@ -117,7 +117,7 @@ void RadiationCorrLandauLifshitz::operator() ( + momentum[2][ipart]*momentum[2][ipart]); // Computation of the Lorentz invariant quantum parameter - chipa = Radiation::compute_chipa(charge_over_mass2, + particle_chi = Radiation::computeParticleChi(charge_over_mass2, momentum[0][ipart],momentum[1][ipart],momentum[2][ipart], gamma, (*(Ex+ipart-ipart_ref)),(*(Ey+ipart-ipart_ref)),(*(Ez+ipart-ipart_ref)), @@ -125,12 +125,12 @@ void RadiationCorrLandauLifshitz::operator() ( // Effect on the momentum // (Should be vectorized with masked instructions) - if (chipa >= RadiationTables.get_chipa_radiation_threshold()) + if (particle_chi >= RadiationTables.get_chipa_radiation_threshold()) { // Radiated energy during the time step temp = - RadiationTables.get_corrected_cont_rad_energy_Ridgers(chipa,dt); + RadiationTables.get_corrected_cont_rad_energy_Ridgers(particle_chi,dt_); // Temporary factor temp *= gamma/(gamma*gamma - 1); @@ -159,5 +159,5 @@ void RadiationCorrLandauLifshitz::operator() ( { radiated_energy_loc += weight[ipart]*rad_norm_energy[ipart] ; } - radiated_energy += radiated_energy_loc; + radiated_energy_ += radiated_energy_loc; } diff --git a/src/Radiation/RadiationLandauLifshitz.cpp b/src/Radiation/RadiationLandauLifshitz.cpp index c0ba2c3c5..b0efaf3cc 100644 --- a/src/Radiation/RadiationLandauLifshitz.cpp +++ b/src/Radiation/RadiationLandauLifshitz.cpp @@ -73,7 +73,7 @@ void RadiationLandauLifshitz::operator() ( const double one_over_mass_2 = pow(one_over_mass_,2.); // Temporary quantum parameter - double chipa; + double particle_chi; // Temporary Lorentz factor double gamma; @@ -99,7 +99,7 @@ void RadiationLandauLifshitz::operator() ( std::vector rad_norm_energy (iend-istart,0); // Reinitialize the cumulative radiated energy for the current thread - this->radiated_energy = 0.; + radiated_energy_ = 0.; // _______________________________________________________________ // Computation @@ -114,19 +114,19 @@ void RadiationLandauLifshitz::operator() ( + momentum[2][ipart]*momentum[2][ipart]); // Computation of the Lorentz invariant quantum parameter - chipa = Radiation::compute_chipa(charge_over_mass2, + particle_chi = Radiation::computeParticleChi(charge_over_mass2, momentum[0][ipart],momentum[1][ipart],momentum[2][ipart], gamma, (*(Ex+ipart-ipart_ref)),(*(Ey+ipart-ipart_ref)),(*(Ez+ipart-ipart_ref)), (*(Bx+ipart-ipart_ref)),(*(By+ipart-ipart_ref)),(*(Bz+ipart-ipart_ref)) ); // Effect on the momentum - if (chipa >= RadiationTables.get_chipa_radiation_threshold()) + if (particle_chi >= RadiationTables.get_chipa_radiation_threshold()) { // Radiated energy during the time step temp = - RadiationTables.get_classical_cont_rad_energy(chipa,dt); + RadiationTables.get_classical_cont_rad_energy(particle_chi,dt_); // Effect on the momentum // Temporary factor @@ -155,5 +155,5 @@ void RadiationLandauLifshitz::operator() ( { radiated_energy_loc += weight[ipart]*rad_norm_energy[ipart - istart] ; } - radiated_energy += radiated_energy_loc; + radiated_energy_ += radiated_energy_loc; } diff --git a/src/Radiation/RadiationMonteCarlo.cpp b/src/Radiation/RadiationMonteCarlo.cpp index 8d8c31b32..cae52dd59 100644 --- a/src/Radiation/RadiationMonteCarlo.cpp +++ b/src/Radiation/RadiationMonteCarlo.cpp @@ -23,9 +23,9 @@ RadiationMonteCarlo::RadiationMonteCarlo(Params& params, Species * species) : Radiation(params, species) { - this->radiation_photon_sampling = species->radiation_photon_sampling; - this->radiation_photon_gamma_threshold = species->radiation_photon_gamma_threshold; - this->inv_radiation_photon_sampling = 1. / this->radiation_photon_sampling; + this->radiation_photon_sampling_ = species->radiation_photon_sampling_; + this->radiation_photon_gamma_threshold_ = species->radiation_photon_gamma_threshold_; + this->inv_radiation_photon_sampling_ = 1. / this->radiation_photon_sampling_; } // --------------------------------------------------------------------------------------------------------------------- @@ -78,7 +78,7 @@ void RadiationMonteCarlo::operator() ( const double one_over_mass_2 = pow(one_over_mass_,2.); // Temporary quantum parameter - double chipa; + double particle_chi; // Temporary Lorentz factor double gamma; @@ -105,7 +105,7 @@ void RadiationMonteCarlo::operator() ( // Position shortcut double* position[3]; - for ( int i = 0 ; iradiated_energy = 0.; + radiated_energy_ = 0.; // _______________________________________________________________ // Computation @@ -135,8 +135,8 @@ void RadiationMonteCarlo::operator() ( mc_it_nb = 0; // Monte-Carlo Manager inside the time step - while ((local_it_time < dt) - &&(mc_it_nb < mc_it_nb_max)) + while ((local_it_time < dt_) + &&(mc_it_nb < max_monte_carlo_iterations_)) { // Gamma @@ -145,51 +145,51 @@ void RadiationMonteCarlo::operator() ( + momentum[2][ipart]*momentum[2][ipart]); // Computation of the Lorentz invariant quantum parameter - chipa = Radiation::compute_chipa(charge_over_mass2, + particle_chi = Radiation::computeParticleChi(charge_over_mass2, momentum[0][ipart],momentum[1][ipart],momentum[2][ipart], gamma, (*(Ex+ipart-ipart_ref)),(*(Ey+ipart-ipart_ref)),(*(Ez+ipart-ipart_ref)), (*(Bx+ipart-ipart_ref)),(*(By+ipart-ipart_ref)),(*(Bz+ipart-ipart_ref)) ); // Update the quantum parameter in species - // chi[ipart] = chipa; + // chi[ipart] = particle_chi; // Discontinuous emission: New emission // If tau[ipart] <= 0, this is a new emission - // We also check that chipa > chipa_threshold, - // else chipa is too low to induce a discontinuous emission - if ((chipa > RadiationTables.get_chipa_disc_min_threshold()) - && (tau[ipart] <= epsilon_tau) ) + // We also check that particle_chi > chipa_threshold, + // else particle_chi is too low to induce a discontinuous emission + if ((particle_chi > RadiationTables.get_chipa_disc_min_threshold()) + && (tau[ipart] <= epsilon_tau_) ) { // New final optical depth to reach for emision - while (tau[ipart] <= epsilon_tau) + while (tau[ipart] <= epsilon_tau_) tau[ipart] = -log(1.-Rand::uniform()); } // Discontinuous emission: emission under progress - // If epsilon_tau > 0 - if (tau[ipart] > epsilon_tau) + // If epsilon_tau_ > 0 + if (tau[ipart] > epsilon_tau_) { // from the cross section - temp = RadiationTables.compute_dNphdt(chipa,gamma); + temp = RadiationTables.compute_dNphdt(particle_chi,gamma); // Time to discontinuous emission // If this time is > the remaining iteration time, // we have a synchronization - emission_time = std::min(tau[ipart]/temp, dt - local_it_time); + emission_time = std::min(tau[ipart]/temp, dt_ - local_it_time); // Update of the optical depth tau[ipart] -= temp*emission_time; // If the final optical depth is reached - if (tau[ipart] <= epsilon_tau) + if (tau[ipart] <= epsilon_tau_) { // Emission of a photon RadiationMonteCarlo::photon_emission(ipart, - chipa,gamma, + particle_chi,gamma, position, momentum, weight, @@ -210,41 +210,44 @@ void RadiationMonteCarlo::operator() ( } // Continuous emission - // chipa needs to be below the discontinuous threshold - // chipa needs to be above the continuous threshold + // particle_chi needs to be below the discontinuous threshold + // particle_chi needs to be above the continuous threshold // No discontunous emission is in progress: - // tau[ipart] <= epsilon_tau - else if ((chipa <= RadiationTables.get_chipa_disc_min_threshold()) - && (tau[ipart] <= epsilon_tau) - && (chipa > RadiationTables.get_chipa_radiation_threshold())) + // tau[ipart] <= epsilon_tau_ + else if ((particle_chi <= RadiationTables.get_chipa_disc_min_threshold()) + && (tau[ipart] <= epsilon_tau_) + && (particle_chi > RadiationTables.get_chipa_radiation_threshold()) + && (gamma > 1.)) { // Remaining time of the iteration - emission_time = dt - local_it_time; + emission_time = dt_ - local_it_time; // Radiated energy during emission_time cont_rad_energy = - RadiationTables.get_corrected_cont_rad_energy_Ridgers(chipa, + RadiationTables.get_corrected_cont_rad_energy_Ridgers(particle_chi, emission_time); // Effect on the momentum temp = cont_rad_energy*gamma/(gamma*gamma-1.); for ( int i = 0 ; i<3 ; i++ ) + { momentum[i][ipart] -= temp*momentum[i][ipart]; + } // Incrementation of the radiated energy cumulative parameter - radiated_energy += weight[ipart]*(gamma - sqrt(1.0 + radiated_energy_ += weight[ipart]*(gamma - sqrt(1.0 + momentum[0][ipart]*momentum[0][ipart] + momentum[1][ipart]*momentum[1][ipart] + momentum[2][ipart]*momentum[2][ipart])); // End for this particle - local_it_time = dt; + local_it_time = dt_; } - // No emission since chipa is too low - else if (chipa < RadiationTables.get_chipa_radiation_threshold()) + // No emission since particle_chi is too low + else // if (particle_chi < RadiationTables.get_chipa_radiation_threshold()) { - local_it_time = dt; + local_it_time = dt_; } } @@ -257,7 +260,7 @@ void RadiationMonteCarlo::operator() ( //! Perform the photon emission (creation of a super-photon //! and slow down of the emitting particle) //! \param ipart particle index -//! \param chipa particle quantum parameter +//! \param particle_chi particle quantum parameter //! \param gammapa particle gamma factor //! \param position particle position //! \param momentum particle momentum @@ -265,7 +268,7 @@ void RadiationMonteCarlo::operator() ( // for nonlinear inverse Compton scattering // --------------------------------------------------------------------------------------------------------------------- void RadiationMonteCarlo::photon_emission(int ipart, - double &chipa, + double & particle_chi, double & gammapa, double * position[3], double * momentum[3], @@ -275,16 +278,16 @@ void RadiationMonteCarlo::photon_emission(int ipart, { // ____________________________________________________ // Parameters - double chiph; // Photon quantum parameter + double photon_chi; // Photon quantum parameter double gammaph; // Photon gamma factor double inv_old_norm_p; //double new_norm_p; // Get the photon quantum parameter from the table xip - chiph = RadiationTables.compute_chiph_emission(chipa); + photon_chi = RadiationTables.compute_chiph_emission(particle_chi); // compute the photon gamma factor - gammaph = chiph/chipa*(gammapa-1.0); + gammaph = photon_chi/particle_chi*(gammapa-1.0); // ____________________________________________________ // Creation of the new photon @@ -309,7 +312,7 @@ void RadiationMonteCarlo::photon_emission(int ipart, // Creation of macro-photons if requested // Check that the photon_species is defined and the threshold on the energy if (photon_species - && (gammaph >= radiation_photon_gamma_threshold)) + && (gammaph >= radiation_photon_gamma_threshold_)) { /* --------------------------------------------------------------------- // First method: emission of a single photon @@ -319,7 +322,7 @@ void RadiationMonteCarlo::photon_emission(int ipart, int idNew = new_photons.size() - 1; - for (int i=0; iradiated_energy = 0.; + radiated_energy_ = 0.; const double chipa_radiation_threshold = RadiationTables.get_chipa_radiation_threshold(); const double factor_cla_rad_power = RadiationTables.get_factor_cla_rad_power(); @@ -137,7 +137,7 @@ void RadiationNiel::operator() ( + momentum[2][ipart]*momentum[2][ipart]); // Computation of the Lorentz invariant quantum parameter - chipa[ipart] = Radiation::compute_chipa(charge_over_mass2, + particle_chi[ipart] = Radiation::computeParticleChi(charge_over_mass2, momentum[0][ipart],momentum[1][ipart],momentum[2][ipart], gamma[ipart], (*(Ex+ipart-ipart_ref)),(*(Ey+ipart-ipart_ref)),(*(Ez+ipart-ipart_ref)), @@ -150,12 +150,12 @@ void RadiationNiel::operator() ( /*for (ipart=0 ; ipart < nbparticles; ipart++ ) { - // Below chipa = chipa_radiation_threshold, radiation losses are negligible - if (chipa[ipart] > chipa_radiation_threshold) + // Below particle_chi = chipa_radiation_threshold, radiation losses are negligible + if (particle_chi[ipart] > chipa_radiation_threshold) { // Pick a random number in the normal distribution of standard - // deviation sqrt(dt) (variance dt) + // deviation sqrt(dt_) (variance dt_) random_numbers[ipart] = Rand::normal(sqrtdt); } }*/ @@ -165,12 +165,12 @@ void RadiationNiel::operator() ( for (ipart=0 ; ipart < nbparticles; ipart++ ) { - // Below chipa = chipa_radiation_threshold, radiation losses are negligible - if (chipa[ipart] > chipa_radiation_threshold) + // Below particle_chi = chipa_radiation_threshold, radiation losses are negligible + if (particle_chi[ipart] > chipa_radiation_threshold) { // Pick a random number in the normal distribution of standard - // deviation sqrt(dt) (variance dt) + // deviation sqrt(dt_) (variance dt_) //random_numbers[ipart] = 2.*Rand::uniform() -1.; random_numbers[ipart] = 2.*drand48() -1.; } @@ -181,8 +181,8 @@ void RadiationNiel::operator() ( #pragma omp simd private(p,temp) for (ipart=0 ; ipart < nbparticles; ipart++ ) { - // Below chipa = chipa_radiation_threshold, radiation losses are negligible - if (chipa[ipart] > chipa_radiation_threshold) + // Below particle_chi = chipa_radiation_threshold, radiation losses are negligible + if (particle_chi[ipart] > chipa_radiation_threshold) { temp = -log((1.0-random_numbers[ipart])*(1.0+random_numbers[ipart])); @@ -226,13 +226,13 @@ void RadiationNiel::operator() ( for (ipart=0 ; ipart < nbparticles; ipart++ ) { - // Below chipa = chipa_radiation_threshold, radiation losses are negligible - if (chipa[ipart] > chipa_radiation_threshold) + // Below particle_chi = chipa_radiation_threshold, radiation losses are negligible + if (particle_chi[ipart] > chipa_radiation_threshold) { - //h = RadiationTables.get_h_Niel_from_fit_order10(chipa[ipart]); - //h = RadiationTables.get_h_Niel_from_fit_order5(chipa[ipart]); - temp = RadiationTables.get_h_Niel_from_table(chipa[ipart]); + //h = RadiationTables.get_h_Niel_from_fit_order10(particle_chi[ipart]); + //h = RadiationTables.get_h_Niel_from_fit_order5(particle_chi[ipart]); + temp = RadiationTables.get_h_Niel_from_table(particle_chi[ipart]); diffusion[ipart] = sqrt(factor_cla_rad_power*gamma[ipart]*temp)*random_numbers[ipart]; } @@ -245,11 +245,11 @@ void RadiationNiel::operator() ( for (ipart=0 ; ipart < nbparticles; ipart++ ) { - // Below chipa = chipa_radiation_threshold, radiation losses are negligible - if (chipa[ipart] > chipa_radiation_threshold) + // Below particle_chi = chipa_radiation_threshold, radiation losses are negligible + if (particle_chi[ipart] > chipa_radiation_threshold) { - temp = RadiationTables.get_h_Niel_from_fit_order5(chipa[ipart]); + temp = RadiationTables.get_h_Niel_from_fit_order5(particle_chi[ipart]); diffusion[ipart] = sqrt(factor_cla_rad_power*gamma[ipart]*temp)*random_numbers[ipart]; } @@ -262,11 +262,11 @@ void RadiationNiel::operator() ( for (ipart=0 ; ipart < nbparticles; ipart++ ) { - // Below chipa = chipa_radiation_threshold, radiation losses are negligible - if (chipa[ipart] > chipa_radiation_threshold) + // Below particle_chi = chipa_radiation_threshold, radiation losses are negligible + if (particle_chi[ipart] > chipa_radiation_threshold) { - temp = RadiationTables.get_h_Niel_from_fit_order10(chipa[ipart]); + temp = RadiationTables.get_h_Niel_from_fit_order10(particle_chi[ipart]); diffusion[ipart] = sqrt(factor_cla_rad_power*gamma[ipart]*temp)*random_numbers[ipart]; } @@ -280,11 +280,11 @@ void RadiationNiel::operator() ( for (ipart=0 ; ipart < nbparticles; ipart++ ) { - // Below chipa = chipa_radiation_threshold, radiation losses are negligible - if (chipa[ipart] > chipa_radiation_threshold) + // Below particle_chi = chipa_radiation_threshold, radiation losses are negligible + if (particle_chi[ipart] > chipa_radiation_threshold) { - temp = RadiationTables.get_h_Niel_from_fit_Ridgers(chipa[ipart]); + temp = RadiationTables.get_h_Niel_from_fit_Ridgers(particle_chi[ipart]); diffusion[ipart] = sqrt(factor_cla_rad_power*gamma[ipart]*temp)*random_numbers[ipart]; } @@ -296,13 +296,13 @@ void RadiationNiel::operator() ( #pragma omp simd private(temp,rad_energy) for (ipart=0 ; ipart chipa_radiation_threshold) + // Below particle_chi = chipa_radiation_threshold, radiation losses are negligible + if (particle_chi[ipart] > chipa_radiation_threshold) { // Radiated energy during the time step rad_energy = - RadiationTables.get_corrected_cont_rad_energy_Ridgers(chipa[ipart],dt); + RadiationTables.get_corrected_cont_rad_energy_Ridgers(particle_chi[ipart],dt_); // Effect on the momentum // Temporary factor @@ -330,7 +330,7 @@ void RadiationNiel::operator() ( + momentum[1][ipart]*momentum[1][ipart] + momentum[2][ipart]*momentum[2][ipart])); } - radiated_energy += radiated_energy_loc; + radiated_energy_ += radiated_energy_loc; //double t5 = MPI_Wtime(); diff --git a/src/Radiation/RadiationTables.cpp b/src/Radiation/RadiationTables.cpp index bb2ce4967..47d7d3e27 100644 --- a/src/Radiation/RadiationTables.cpp +++ b/src/Radiation/RadiationTables.cpp @@ -134,7 +134,7 @@ void RadiationTables::initParams(Params& params) // Path to the databases PyTools::extract("table_path", table_path, "RadiationReaction"); - // Radiation threshold on the quantum parameter chipa + // Radiation threshold on the quantum parameter particle_chi PyTools::extract("chipa_radiation_threshold", chipa_radiation_threshold,"RadiationReaction"); @@ -255,7 +255,7 @@ void RadiationTables::compute_h_table(SmileiMPI *smpi) { // Temporary particle chi value - double chipa; + double particle_chi; // For percentages double pct = 0.; double dpct = 10.; @@ -324,12 +324,12 @@ void RadiationTables::compute_h_table(SmileiMPI *smpi) // Loop over the table values for(int i = 0 ; i < length_table[rank] ; i++) { - chipa = pow(10.,(imin_table[rank] + i)*h_chipa_delta + particle_chi = pow(10.,(imin_table[rank] + i)*h_chipa_delta + h_log10_chipa_min); // 100 iterations is in theory sufficient to get the convergence // at 1e-15 - buffer[i] = RadiationTables::compute_h_Niel(chipa,400,1e-15); + buffer[i] = RadiationTables::compute_h_Niel(particle_chi,400,1e-15); if (100.*i >= length_table[rank]*pct) { @@ -378,7 +378,7 @@ void RadiationTables::compute_integfochi_table(SmileiMPI *smpi) // Get the MPI rank rank = smpi->getRank(); - MESSAGE(" --- Integration F/chipa table:"); + MESSAGE(" --- Integration F/particle_chi table:"); // Initial timer t0 = MPI_Wtime(); @@ -391,7 +391,7 @@ void RadiationTables::compute_integfochi_table(SmileiMPI *smpi) { // Temporary particle chi value - double chipa; + double particle_chi; // For percentages double pct = 0.; double dpct = 10.; @@ -447,11 +447,11 @@ void RadiationTables::compute_integfochi_table(SmileiMPI *smpi) // Loop over the table values for(int i = 0 ; i < length_table[rank] ; i++) { - chipa = pow(10.,(imin_table[rank] + i)*integfochi_chipa_delta + particle_chi = pow(10.,(imin_table[rank] + i)*integfochi_chipa_delta + integfochi_log10_chipa_min); - buffer[i] = RadiationTables::compute_integfochi(chipa, - 0.98e-40*chipa,0.98*chipa,400,1e-15); + buffer[i] = RadiationTables::compute_integfochi(particle_chi, + 0.98e-40*particle_chi,0.98*particle_chi,400,1e-15); //std::cout << rank << " " << buffer[i] << std::endl; @@ -486,7 +486,7 @@ void RadiationTables::compute_integfochi_table(SmileiMPI *smpi) //! Computation of the minimum photon quantum parameter chiphmin //! for the xip array and computation of the xip array. // -//! \details Under the minimum chiph value, photon energy is +//! \details Under the minimum photon_chi value, photon energy is //! considered negligible. // //! \param smpi Object of class SmileiMPI containing MPI properties @@ -515,9 +515,9 @@ void RadiationTables::compute_xip_table(SmileiMPI *smpi) if (!table_exists) { // Parameters: - double chipa; // Temporary particle chi value - double chiph; // Temporary photon chi value - double chiph_delta; // Temporary delta for chiph + double particle_chi; // Temporary particle chi value + double photon_chi; // Temporary photon chi value + double chiph_delta; // Temporary delta for photon_chi double logchiphmin; // Temporary log10 of photon chi value double xip; // Temporary xip double numerator; @@ -587,19 +587,19 @@ void RadiationTables::compute_xip_table(SmileiMPI *smpi) xip = 1; logchiphmin = (imin_table[rank] + ichipa)*xip_chipa_delta + xip_log10_chipa_min; - chipa = pow(10.,logchiphmin); + particle_chi = pow(10.,logchiphmin); // Denominator of xip - denominator = RadiationTables::compute_integfochi(chipa, - 0.99e-40*chipa,0.99*chipa,200,1e-13); + denominator = RadiationTables::compute_integfochi(particle_chi, + 0.99e-40*particle_chi,0.99*particle_chi,200,1e-13); k = 0; while(k < xip_power) { logchiphmin -= pow(0.1,k); - chiph = pow(10.,logchiphmin); - numerator = RadiationTables::compute_integfochi(chipa, - 0.99e-40*chiph,0.99*chiph,200,1e-13); + photon_chi = pow(10.,logchiphmin); + numerator = RadiationTables::compute_integfochi(particle_chi, + 0.99e-40*photon_chi,0.99*photon_chi,200,1e-13); if (numerator == 0) { @@ -638,34 +638,34 @@ void RadiationTables::compute_xip_table(SmileiMPI *smpi) // Allocation of the local buffer buffer = new double [length_table[rank]*xip_chiph_dim]; - // Loop for xip in the chipa dimension + // Loop for xip in the particle_chi dimension pct = 0; dpct = std::max(10.,100./(length_table[rank]*xip_chiph_dim)); for(int ichipa = 0 ; ichipa < length_table[rank] ; ichipa++) { - chipa = (imin_table[rank] + ichipa)*xip_chipa_delta + particle_chi = (imin_table[rank] + ichipa)*xip_chipa_delta + xip_log10_chipa_min; - chiph_delta = (chipa - xip_chiphmin_table[imin_table[rank] + ichipa]) + chiph_delta = (particle_chi - xip_chiphmin_table[imin_table[rank] + ichipa]) / (xip_chiph_dim - 1); - chipa = pow(10.,chipa); + particle_chi = pow(10.,particle_chi); // Denominator of xip - denominator = RadiationTables::compute_integfochi(chipa, - 1e-40*chipa,chipa,300,1e-15); + denominator = RadiationTables::compute_integfochi(particle_chi, + 1e-40*particle_chi,particle_chi,300,1e-15); - // Loop in the chiph dimension + // Loop in the photon_chi dimension for (int ichiph = 0 ; ichiph < xip_chiph_dim ; ichiph ++) { - // Local chiph value - chiph = pow(10.,ichiph*chiph_delta + + // Local photon_chi value + photon_chi = pow(10.,ichiph*chiph_delta + xip_chiphmin_table[imin_table[rank] + ichipa]); // Numerator of xip - numerator = RadiationTables::compute_integfochi(chipa, - 0.99e-40*chiph,0.99*chiph,300,1e-15); + numerator = RadiationTables::compute_integfochi(particle_chi, + 0.99e-40*photon_chi,0.99*photon_chi,300,1e-15); // Update local buffer value buffer[ichipa*xip_chiph_dim + ichiph] = std::min(1.,numerator / denominator); @@ -774,7 +774,7 @@ void RadiationTables::output_h_table() file << "Stochastic synchrotron-like radiation model of Niel \n"; - file << "Dimension chipa - chipa min - chipa max \n"; + file << "Dimension particle_chi - particle_chi min - particle_chi max \n"; file << h_dim; file << " " @@ -893,7 +893,7 @@ void RadiationTables::output_integfochi_table() file << "Table Integration F(CHI)/CHI for Nonlinear Compton Scattering \n"; - file << "Dimension chipa - chipa min - chipa max \n"; + file << "Dimension particle_chi - particle_chi min - particle_chi max \n"; file << integfochi_dim ; file << " " @@ -1017,7 +1017,7 @@ void RadiationTables::output_xip_table() file << "Table xip_chiphmin and xip for Nonlinear Compton Scattering \n"; - file << "Dimension chipa - Dimension chiph - chipa min - chipa max \n"; + file << "Dimension particle_chi - Dimension photon_chi - particle_chi min - particle_chi max \n"; file << xip_chipa_dim << " " << xip_chiph_dim << " " @@ -1181,16 +1181,16 @@ void RadiationTables::output_tables(SmileiMPI *smpi) // ----------------------------------------------------------------------------- // ----------------------------------------------------------------------------- -//! Computation of the photon quantum parameter chiph for emission +//! Computation of the photon quantum parameter photon_chi for emission //! ramdomly and using the tables xip and chiphmin // -//! \param chipa particle quantum parameter +//! \param particle_chi particle quantum parameter // ----------------------------------------------------------------------------- -double RadiationTables::compute_chiph_emission(double chipa) +double RadiationTables::compute_chiph_emission(double particle_chi) { - // Log10 of chipa + // Log10 of particle_chi double logchipa; - double chiph; + double photon_chi; double chiph_xip_delta; // Random xip double xip; @@ -1202,12 +1202,12 @@ double RadiationTables::compute_chiph_emission(double chipa) double d; int ixip; - logchipa = log10(chipa); + logchipa = log10(particle_chi); // --------------------------------------- - // index of chipa in xip_table + // index of particle_chi in xip_table // --------------------------------------- - // Use floor so that chipa corresponding to ichipa is <= given chipa + // Use floor so that particle_chi corresponding to ichipa is <= given particle_chi ichipa = int(floor((logchipa-xip_log10_chipa_min)*(xip_chipa_inv_delta))); // Checking that ichipa is in the range of the tables @@ -1222,21 +1222,21 @@ double RadiationTables::compute_chiph_emission(double chipa) } // --------------------------------------- - // Search of the index ichiph for chiph + // Search of the index ichiph for photon_chi // --------------------------------------- // First, we compute a random xip in [0,1[ xip = Rand::uniform(); // If the randomly computed xip if below the first one of the row, - // we take the first one which corresponds to the minimal photon chiph + // we take the first one which corresponds to the minimal photon photon_chi if (xip <= xip_table[ichipa*xip_chiph_dim]) { ichiph = 0; xip = xip_table[ichipa*xip_chiph_dim]; } // Above the last xip of the row, the last one corresponds - // to the maximal photon chiph + // to the maximal photon photon_chi else if (xip > xip_table[(ichipa+1)*xip_chiph_dim-2]) { ichiph = xip_chiph_dim-2; @@ -1250,21 +1250,21 @@ double RadiationTables::compute_chiph_emission(double chipa) &xip_table[ichipa*xip_chiph_dim],xip,xip_chiph_dim); } - // Corresponding chipa for ichipa + // Corresponding particle_chi for ichipa logchipa = ichipa*xip_chipa_delta+xip_log10_chipa_min; - // Delta for the corresponding chipa + // Delta for the corresponding particle_chi chiph_xip_delta = (logchipa - xip_chiphmin_table[ichipa]) *xip_inv_chiph_dim_minus_one; // -------------------------------------------------------------------- - // Compute chiph + // Compute photon_chi // This method is slow but more accurate than taking the nearest point // -------------------------------------------------------------------- ixip = ichipa*xip_chiph_dim + ichiph; - // Computation of the final chiph by interpolation + // Computation of the final photon_chi by interpolation if (xip_table[ixip+1] - xip_table[ixip] > 1e-15) { log10_chiphm = ichiph*chiph_xip_delta @@ -1274,39 +1274,39 @@ double RadiationTables::compute_chiph_emission(double chipa) d = (xip - xip_table[ixip]) / (xip_table[ixip+1] - xip_table[ixip]); // Chiph after linear interpolation in the logarithmic scale - chiph = pow(10.,log10_chiphm*(1.0-d) + log10_chiphp*(d)); + photon_chi = pow(10.,log10_chiphm*(1.0-d) + log10_chiphp*(d)); } else // For integration reasons, we can have xip_table[ixip+1] = xip_table[ixip] // In this case, no interpolation { - chiph = pow(10.,ichiph*chiph_xip_delta + photon_chi = pow(10.,ichiph*chiph_xip_delta + xip_chiphmin_table[ichipa]); } // ------------------------------------------------------------ - // Compute chiph + // Compute photon_chi // Fastest method using the nearest point but less accurate // ------------------------------------------------------------ - /*chiph = pow(10.,ichiph*chiph_xip_delta + /*photon_chi = pow(10.,ichiph*chiph_xip_delta + xip_chiphmin_table[ichipa]);*/ - return chiph; + return photon_chi; } // --------------------------------------------------------------------------------------------------------------------- //! Computation of the Cross Section dNph/dt which is also //! the number of photons generated per time unit. // -//! \param chipa particle quantum parameter +//! \param particle_chi particle quantum parameter //! \param gfpa particle gamma factor // --------------------------------------------------------------------------------------------------------------------- -double RadiationTables::compute_dNphdt(double chipa,double gfpa) +double RadiationTables::compute_dNphdt(double particle_chi,double gfpa) { - // Log of the particle quantum parameter chipa + // Log of the particle quantum parameter particle_chi double logchipa; double logchipam; double logchipap; @@ -1315,7 +1315,7 @@ double RadiationTables::compute_dNphdt(double chipa,double gfpa) // final value double dNphdt; - logchipa = log10(chipa); + logchipa = log10(particle_chi); // Lower index for interpolation in the table integfochi ichipa = int(floor((logchipa-integfochi_log10_chipa_min) @@ -1343,23 +1343,23 @@ double RadiationTables::compute_dNphdt(double chipa,double gfpa) integfochi_table[ichipa]*fabs(logchipap - logchipa))*integfochi_chipa_inv_delta; } - return factor_dNphdt*dNphdt*chipa/gfpa; + return factor_dNphdt*dNphdt*particle_chi/gfpa; } // --------------------------------------------------------------------------------------------------------------------- //! \brief //! Compute integration of F/chi between -//! using Gauss-Legendre for a given chipa value +//! using Gauss-Legendre for a given particle_chi value // // -//! \param chipa particle (electron for instance) quantum parameter +//! \param particle_chi particle (electron for instance) quantum parameter //! \param chipmin Minimal integration value (photon quantum parameter) //! \param chipmax Maximal integration value (photon quantum parameter) //! \param nbit number of points for integration //! \param eps integration accuracy // --------------------------------------------------------------------------------------------------------------------- -double RadiationTables::compute_integfochi(double chipa, +double RadiationTables::compute_integfochi(double particle_chi, double chipmin, double chipmax, int nbit, @@ -1372,7 +1372,7 @@ double RadiationTables::compute_integfochi(double chipa, double * gauleg_x = new double[nbit]; double * gauleg_w = new double[nbit]; // Photon quantum parameter - double chiph; + double photon_chi; // Integration result double integ; // Synchrotron emissivity @@ -1386,11 +1386,11 @@ double RadiationTables::compute_integfochi(double chipa, // Integration loop integ = 0; - #pragma omp parallel for reduction(+:integ) private(chiph,sync_emi) shared(chipa,gauleg_w,gauleg_x) + #pragma omp parallel for reduction(+:integ) private(photon_chi,sync_emi) shared(particle_chi,gauleg_w,gauleg_x) for(i=0 ; i< nbit ; i++) { - chiph = pow(10.,gauleg_x[i]); - sync_emi = RadiationTables::compute_sync_emissivity_ritus(chipa,chiph,200,1e-15); + photon_chi = pow(10.,gauleg_x[i]); + sync_emi = RadiationTables::compute_sync_emissivity_ritus(particle_chi,photon_chi,200,1e-15); integ += gauleg_w[i]*sync_emi*log(10); } @@ -1401,19 +1401,19 @@ double RadiationTables::compute_integfochi(double chipa, // --------------------------------------------------------------------------------------------------------------------- //! Computation of the synchrotron emissivity following the formulae of Ritus // -//! \param chipa particle quantum parameter -//! \param chiph photon quantum parameter +//! \param particle_chi particle quantum parameter +//! \param photon_chi photon quantum parameter //! \param nbit number of iterations for the Gauss-Legendre integration //! \param eps epsilon for the modified bessel function // --------------------------------------------------------------------------------------------------------------------- -double RadiationTables::compute_sync_emissivity_ritus(double chipa, - double chiph, int nbit, double eps) +double RadiationTables::compute_sync_emissivity_ritus(double particle_chi, + double photon_chi, int nbit, double eps) { //std::cout << "RadiationTables::compute_sync_emissivity_ritus" << std::endl; // The photon quantum parameter should be below the electron one - if (chipa > chiph) + if (particle_chi > photon_chi) { // Arrays for Gauss-Legendre integration double * gauleg_x = new double[nbit]; @@ -1426,13 +1426,13 @@ double RadiationTables::compute_sync_emissivity_ritus(double chipa, // Iterator int i; - double y = chiph/(3.*chipa*(chipa-chiph)); + double y = photon_chi/(3.*particle_chi*(particle_chi-photon_chi)); // Computation of Part. 1 // Call the modified Bessel function to get K userFunctions::modified_bessel_IK(2./3.,2*y,I,dI,K,dK,50000,eps,false); - part1 = (2. + 3.*chiph*y)*(K); + part1 = (2. + 3.*photon_chi*y)*(K); // Computation of Part. 2 // Using Gauss Legendre integration @@ -1449,34 +1449,34 @@ double RadiationTables::compute_sync_emissivity_ritus(double chipa, } // Factor for final result - y = 2*chiph/(3*pow(chipa,2.)); + y = 2*photon_chi/(3*pow(particle_chi,2.)); return (part1 - part2)*y; } - else if (chipa == chiph) + else if (particle_chi == photon_chi) { return 0; } else { - ERROR("In compute_sync_emissivity_ritus: chipa " << chipa - << " < chiph " - << chiph); + ERROR("In compute_sync_emissivity_ritus: particle_chi " << particle_chi + << " < photon_chi " + << photon_chi); return -1.; } } // ----------------------------------------------------------------------------- -//! Return the value of the function h(chipa) of Niel et al. +//! Return the value of the function h(particle_chi) of Niel et al. //! Use an integration of Gauss-Legendre // -//! \param chipa particle quantum parameter +//! \param particle_chi particle quantum parameter //! \param nbit number of iterations for the Gauss-Legendre integration //! \param eps epsilon for the modified bessel function // ----------------------------------------------------------------------------- -double RadiationTables::compute_h_Niel(double chipa, +double RadiationTables::compute_h_Niel(double particle_chi, int nbit, double eps) { // Arrays for Gauss-Legendre integration @@ -1503,8 +1503,8 @@ double RadiationTables::compute_h_Niel(double chipa, userFunctions::modified_bessel_IK(2./3.,nu,I,dI,K23,dK,50000,eps,false); h += gauleg_w[i]*log(10.)*nu - * (( 2.*pow(chipa*nu,3) / pow(2. + 3.*nu*chipa,3)*K53 ) - + ( 54.*pow(chipa,5)*pow(nu,4) / pow(2. + 3.*nu*chipa,5)*K23 )); + * (( 2.*pow(particle_chi*nu,3) / pow(2. + 3.*nu*particle_chi,3)*K53 ) + + ( 54.*pow(particle_chi,5)*pow(nu,4) / pow(2. + 3.*nu*particle_chi,5)*K23 )); } @@ -1513,17 +1513,17 @@ double RadiationTables::compute_h_Niel(double chipa, } // ----------------------------------------------------------------------------- -//! Return the value of the function h(chipa) of Niel et al. +//! Return the value of the function h(particle_chi) of Niel et al. //! from the computed table h_table -//! \param chipa particle quantum parameter +//! \param particle_chi particle quantum parameter // ----------------------------------------------------------------------------- -double RadiationTables::get_h_Niel_from_table(double chipa) +double RadiationTables::get_h_Niel_from_table(double particle_chi) { int ichipa; double d; // Position in the h_table - d = (log10(chipa)-h_log10_chipa_min)*h_chipa_inv_delta; + d = (log10(particle_chi)-h_log10_chipa_min)*h_chipa_inv_delta; ichipa = int(floor(d)); // distance for interpolation @@ -1537,16 +1537,16 @@ double RadiationTables::get_h_Niel_from_table(double chipa) //! Return the stochastic diffusive component of the pusher //! of Niel et al. //! \param gamma particle Lorentz factor -//! \param chipa particle quantum parameter +//! \param particle_chi particle quantum parameter // ----------------------------------------------------------------------------- double RadiationTables::get_Niel_stochastic_term(double gamma, - double chipa, + double particle_chi, double sqrtdt) { - // Get the value of h for the corresponding chipa + // Get the value of h for the corresponding particle_chi double h,r; - h = RadiationTables::get_h_Niel_from_table(chipa); + h = RadiationTables::get_h_Niel_from_table(particle_chi); // Pick a random number in the normal distribution of standard // deviation sqrt(dt) (variance dt) @@ -2205,7 +2205,7 @@ void RadiationTables::bcast_xip_table(SmileiMPI *smpi) // Inverse of delta xip_chipa_inv_delta = 1./xip_chipa_delta; - // Inverse chiph discetization (regularly used) + // Inverse photon_chi discetization (regularly used) xip_inv_chiph_dim_minus_one = 1./(xip_chiph_dim - 1.); } diff --git a/src/Radiation/RadiationTables.h b/src/Radiation/RadiationTables.h index 75851d2cf..cab58b684 100644 --- a/src/Radiation/RadiationTables.h +++ b/src/Radiation/RadiationTables.h @@ -46,17 +46,17 @@ class RadiationTables // --------------------------------------------------------------------- //! Synchrotron emissivity from Ritus - //! \param chipa particle quantum parameter - //! \param chiph photon quantum parameter + //! \param particle_chi particle quantum parameter + //! \param photon_chi photon quantum parameter //! \param nbit number of iterations for the Gauss-Legendre integration //! \param eps epsilon for the modified bessel function static double compute_sync_emissivity_ritus(double chie, - double chiph, + double photon_chi, int nbit, double eps); //! Computation of the cross-section dNph/dt - double compute_dNphdt(double chipa,double gfpa); + double compute_dNphdt(double particle_chi,double gfpa); //! Compute integration of F/chi between //! using Gauss-Legendre for a given chie value @@ -66,51 +66,51 @@ class RadiationTables int nbit, double eps); - //! Computation of the photon quantum parameter chiph for emission + //! Computation of the photon quantum parameter photon_chi for emission //! ramdomly and using the tables xip and chiphmin - //! \param chipa particle quantum parameter - double compute_chiph_emission(double chipa); + //! \param particle_chi particle quantum parameter + double compute_chiph_emission(double particle_chi); - //! Return the value of the function h(chipa) of Niel et al. + //! Return the value of the function h(particle_chi) of Niel et al. //! Use an integration of Gauss-Legendre // - //! \param chipa particle quantum parameter + //! \param particle_chi particle quantum parameter //! \param nbit number of iterations for the Gauss-Legendre integration //! \param eps epsilon for the modified bessel function - double compute_h_Niel(double chipa,int nbit, double eps); + double compute_h_Niel(double particle_chi,int nbit, double eps); - //! Return the value of the function h(chipa) of Niel et al. + //! Return the value of the function h(particle_chi) of Niel et al. //! from the computed table h_table - //! \param chipa particle quantum parameter - double get_h_Niel_from_table(double chipa); + //! \param particle_chi particle quantum parameter + double get_h_Niel_from_table(double particle_chi); //! Return the stochastic diffusive component of the pusher //! of Niel et al. //! \param gamma particle Lorentz factor - //! \param chipa particle quantum parameter + //! \param particle_chi particle quantum parameter //! \param dt time step double get_Niel_stochastic_term(double gamma, - double chipa, + double particle_chi, double dt); //! Computation of the corrected continuous quantum radiated energy - //! during dt from the quantum parameter chipa using the Ridgers + //! during dt from the quantum parameter particle_chi using the Ridgers //! formulae. - //! \param chipa particle quantum parameter + //! \param particle_chi particle quantum parameter //! \param dt time step //#pragma omp declare simd - double inline get_corrected_cont_rad_energy_Ridgers(double chipa, + double inline get_corrected_cont_rad_energy_Ridgers(double particle_chi, double dt) { - return compute_g_Ridgers(chipa)*dt*chipa*chipa*factor_cla_rad_power; + return compute_g_Ridgers(particle_chi)*dt*particle_chi*particle_chi*factor_cla_rad_power; }; //! Get of the classical continuous radiated energy during dt - //! \param chipa particle quantum parameter + //! \param particle_chi particle quantum parameter //! \param dt time step - double inline get_classical_cont_rad_energy(double chipa, double dt) + double inline get_classical_cont_rad_energy(double particle_chi, double dt) { - return dt*chipa*chipa*factor_cla_rad_power; + return dt*particle_chi*particle_chi*factor_cla_rad_power; }; //! Return the chipa_disc_min_threshold value @@ -129,12 +129,12 @@ class RadiationTables //! Computation of the function g of Erber using the Ridgers //! approximation formulae - //! \param chipa particle quantum parameter + //! \param particle_chi particle quantum parameter //#pragma omp declare simd - double inline compute_g_Ridgers(double chipa) + double inline compute_g_Ridgers(double particle_chi) { - return pow(1. + 4.8*(1.+chipa)*log(1. + 1.7*chipa) - + 2.44*chipa*chipa,-2./3.); + return pow(1. + 4.8*(1.+particle_chi)*log(1. + 1.7*particle_chi) + + 2.44*particle_chi*particle_chi,-2./3.); }; std::string inline get_h_computation_method() @@ -143,54 +143,54 @@ class RadiationTables } // ----------------------------------------------------------------------------- - //! Return the value of the function h(chipa) of Niel et al. + //! Return the value of the function h(particle_chi) of Niel et al. //! from a polynomial numerical fit at order 10 - //! Valid between chipa in 1E-3 and 1E1 - //! \param chipa particle quantum parameter + //! Valid between particle_chi in 1E-3 and 1E1 + //! \param particle_chi particle quantum parameter // ----------------------------------------------------------------------------- - double inline get_h_Niel_from_fit_order10(double chipa) + double inline get_h_Niel_from_fit_order10(double particle_chi) { // Max relative error ~2E-4 - return exp(-3.231764974833856e-08 * pow(log(chipa),10) - -7.574417415366786e-07 * pow(log(chipa),9) - -5.437005218419013e-06 * pow(log(chipa),8) - -4.359062260446135e-06 * pow(log(chipa),7) - + 5.417842511821415e-05 * pow(log(chipa),6) - -1.263905701127627e-04 * pow(log(chipa),5) - + 9.899812622393002e-04 * pow(log(chipa),4) - + 1.076648497464146e-02 * pow(log(chipa),3) - -1.624860613422593e-01 * pow(log(chipa),2) - + 1.496340836237785e+00 * log(chipa) + return exp(-3.231764974833856e-08 * pow(log(particle_chi),10) + -7.574417415366786e-07 * pow(log(particle_chi),9) + -5.437005218419013e-06 * pow(log(particle_chi),8) + -4.359062260446135e-06 * pow(log(particle_chi),7) + + 5.417842511821415e-05 * pow(log(particle_chi),6) + -1.263905701127627e-04 * pow(log(particle_chi),5) + + 9.899812622393002e-04 * pow(log(particle_chi),4) + + 1.076648497464146e-02 * pow(log(particle_chi),3) + -1.624860613422593e-01 * pow(log(particle_chi),2) + + 1.496340836237785e+00 * log(particle_chi) -2.756744141581370e+00); } // ----------------------------------------------------------------------------- - //! Return the value of the function h(chipa) of Niel et al. + //! Return the value of the function h(particle_chi) of Niel et al. //! from a polynomial numerical fit at order 5 - //! Valid between chipa in 1E-3 and 1E1 - //! \param chipa particle quantum parameter + //! Valid between particle_chi in 1E-3 and 1E1 + //! \param particle_chi particle quantum parameter // ----------------------------------------------------------------------------- - double inline get_h_Niel_from_fit_order5(double chipa) + double inline get_h_Niel_from_fit_order5(double particle_chi) { // Max relative error ~0.02 - return exp(1.399937206900322e-04 * pow(log(chipa),5) - + 3.123718241260330e-03 * pow(log(chipa),4) - + 1.096559086628964e-02 * pow(log(chipa),3) - -1.733977278199592e-01 * pow(log(chipa),2) - + 1.492675770100125e+00 * log(chipa) + return exp(1.399937206900322e-04 * pow(log(particle_chi),5) + + 3.123718241260330e-03 * pow(log(particle_chi),4) + + 1.096559086628964e-02 * pow(log(particle_chi),3) + -1.733977278199592e-01 * pow(log(particle_chi),2) + + 1.492675770100125e+00 * log(particle_chi) -2.748991631516466e+00 ); } // ----------------------------------------------------------------------------- - //! Return the value of the function h(chipa) of Niel et al. + //! Return the value of the function h(particle_chi) of Niel et al. //! using the numerical fit of Ridgers in //! Ridgers et al., ArXiv 1708.04511 (2017) - //! \param chipa particle quantum parameter + //! \param particle_chi particle quantum parameter // ----------------------------------------------------------------------------- - double inline get_h_Niel_from_fit_Ridgers(double chipa) + double inline get_h_Niel_from_fit_Ridgers(double particle_chi) { - return pow(chipa,3)*1.9846415503393384*pow(1. + - (1. + 4.528*chipa)*log(1.+12.29*chipa) + 4.632*pow(chipa,2),-7./6.); + return pow(particle_chi,3)*1.9846415503393384*pow(1. + + (1. + 4.528*particle_chi)*log(1.+12.29*particle_chi) + 4.632*pow(particle_chi,2),-7./6.); } // ----------------------------------------------------------------------------- @@ -287,7 +287,7 @@ class RadiationTables std::string table_path; //! Minimum threshold above which the Monte-Carlo algorithm is working - //! This avoids using the Monte-Carlo algorithm when chipa is too low + //! This avoids using the Monte-Carlo algorithm when particle_chi is too low double chipa_disc_min_threshold; //! Under this value, no radiation loss @@ -359,10 +359,10 @@ class RadiationTables bool integfochi_computed; // --------------------------------------------- - // Table chiph min for xip table + // Table photon_chi min for xip table // --------------------------------------------- - //! Table containing the chiph min values + //! Table containing the photon_chi min values //! Under this value, photon energy is //! considered negligible std::vector xip_chiphmin_table; @@ -375,27 +375,27 @@ class RadiationTables //! that gives gives the probability for a photon emission in the range \f$[0, \chi_{\gamma}]\f$ std::vector xip_table; - //! Minimum boundary for chipa in the table xip and xip_chiphmin + //! Minimum boundary for particle_chi in the table xip and xip_chiphmin double xip_chipa_min; - //! Logarithm of the minimum boundary for chipa in the table xip + //! Logarithm of the minimum boundary for particle_chi in the table xip //! and xip_chiphmin double xip_log10_chipa_min; - //! Maximum boundary for chipa in the table xip and xip_chiphmin + //! Maximum boundary for particle_chi in the table xip and xip_chiphmin double xip_chipa_max; - //! Delta for the chipa discretization in the table xip and xip_chiphmin + //! Delta for the particle_chi discretization in the table xip and xip_chiphmin double xip_chipa_delta; - //! Inverse of the delta for the chipa discretization + //! Inverse of the delta for the particle_chi discretization //! in the table xip and xip_chiphmin double xip_chipa_inv_delta; - //! Dimension of the discretized parameter chipa + //! Dimension of the discretized parameter particle_chi int xip_chipa_dim; - //! Dimension of the discretized parameter chiph + //! Dimension of the discretized parameter photon_chi int xip_chiph_dim; //! 1/(xip_chiph_dim - 1) diff --git a/src/Species/Species.cpp b/src/Species/Species.cpp index 6263349d5..d2e268cde 100644 --- a/src/Species/Species.cpp +++ b/src/Species/Species.cpp @@ -603,7 +603,7 @@ void Species::dynamics(double time_dual, unsigned int ispec, //Point to local thread dedicated buffers //Still needed for ionization vector *Epart = &(smpi->dynamics_Epart[ithread]); - + for (unsigned int ibin = 0 ; ibin < first_index.size() ; ibin++) { #ifdef __DETAILED_TIMERS @@ -649,7 +649,7 @@ void Species::dynamics(double time_dual, unsigned int ispec, nrj_radiation += Radiate->getRadiatedEnergy(); // Update the quantum parameter chi - Radiate->compute_thread_chipa(*particles, + Radiate->computeParticlesChi(*particles, smpi, first_index[ibin], last_index[ibin], @@ -884,14 +884,14 @@ void Species::projection_for_diags(double time_dual, unsigned int ispec, int n_species = patch->vecSpecies.size(); for ( unsigned int imode = 0; imoderho_AM_s[ifield] ? &(*emAM->rho_AM_s[ifield])(0) : &(*emAM->rho_AM_[imode])(0) ; buf[1] = emAM->Jl_s [ifield] ? &(*emAM->Jl_s [ifield])(0) : &(*emAM->Jl_[imode])(0) ; buf[2] = emAM->Jr_s [ifield] ? &(*emAM->Jr_s [ifield])(0) : &(*emAM->Jr_[imode])(0) ; buf[3] = emAM->Jt_s [ifield] ? &(*emAM->Jt_s [ifield])(0) : &(*emAM->Jt_[imode])(0) ; - + for (int iPart=first_index[ibin] ; iPartdensityFrozenComplex(buf[quantity], (*particles), iPart, quantity, imode); @@ -1517,7 +1517,7 @@ int Species::createParticles(vector n_space_to_create, Params& par double x = position[0][ippy]-min_loc ; unsigned int ibin = int(x * one_ov_dbin) ; int ip = indices[ibin] ; //Indice of the position of the particle in the particles array. - + unsigned int ijk[3] = {0,0,0}; for(unsigned int idim=0; idimposition(idim,ip) = position[idim][ippy]; @@ -1535,7 +1535,7 @@ int Species::createParticles(vector n_space_to_create, Params& par for(unsigned int idim=0; idim < 3; idim++) particles->momentum(idim,ip) = momentum[idim][ippy] ; } - + particles->weight(ip) = weight_arr[ippy] ; initCharge(1, ip, charge(ijk[0], ijk[1], ijk[2])); indices[ibin]++; @@ -1894,9 +1894,9 @@ void Species::ponderomotive_update_position_and_currents(double time_dual, unsig // calculate the particle updated position // ------------------------------- if (time_dual>time_frozen) { // moving particle - + smpi->dynamics_resize(ithread, nDim_field, last_index.back(), params.geometry=="AMcylindrical"); - + for (unsigned int ibin = 0 ; ibin < first_index.size() ; ibin++) { // Interpolate the ponderomotive potential and its gradient at the particle position, present and previous timestep diff --git a/src/Species/Species.h b/src/Species/Species.h index e579cfd7a..8a610e3f3 100644 --- a/src/Species/Species.h +++ b/src/Species/Species.h @@ -169,11 +169,11 @@ class Species //! Name of the species where radiated photons go std::string radiation_photon_species; //! Number of photons emitted per particle and per event - int radiation_photon_sampling; + int radiation_photon_sampling_; //! Threshold on the photon Lorentz factor under which the macro-photon //! is not generated but directly added to the energy scalar diags //! This enable to limit emission of useless low-energy photons - double radiation_photon_gamma_threshold; + double radiation_photon_gamma_threshold_; //! Pointer to the species where electron-positron pairs //! from the multiphoton Breit-Wheeler go diff --git a/src/Species/SpeciesAdaptiveV2.cpp b/src/Species/SpeciesAdaptiveV2.cpp index 2e02ff78f..6868435a7 100644 --- a/src/Species/SpeciesAdaptiveV2.cpp +++ b/src/Species/SpeciesAdaptiveV2.cpp @@ -151,7 +151,7 @@ void SpeciesAdaptiveV2::scalar_dynamics(double time_dual, unsigned int ispec, nrj_radiation += Radiate->getRadiatedEnergy(); // Update the quantum parameter chi - Radiate->compute_thread_chipa(*particles, + Radiate->computeParticlesChi(*particles, smpi, first_index[scell], last_index[scell], diff --git a/src/Species/SpeciesFactory.h b/src/Species/SpeciesFactory.h index 1f6fe2bc2..e1ff34eea 100644 --- a/src/Species/SpeciesFactory.h +++ b/src/Species/SpeciesFactory.h @@ -254,9 +254,9 @@ class SpeciesFactory { // Number of photons emitted per Monte-Carlo event if (PyTools::extract("radiation_photon_sampling", - thisSpecies->radiation_photon_sampling, "Species",ispec)) + thisSpecies->radiation_photon_sampling_, "Species",ispec)) { - if (thisSpecies->radiation_photon_sampling < 1) + if (thisSpecies->radiation_photon_sampling_ < 1) { ERROR("For species '" << species_name << "' radiation_photon_sampling should be > 1"); @@ -264,19 +264,19 @@ class SpeciesFactory { } else { - thisSpecies->radiation_photon_sampling = 1; + thisSpecies->radiation_photon_sampling_ = 1; } MESSAGE(2,"> Number of macro-photons emitted per MC event: " - << thisSpecies->radiation_photon_sampling); + << thisSpecies->radiation_photon_sampling_); // Photon energy threshold if (!PyTools::extract("radiation_photon_gamma_threshold", - thisSpecies->radiation_photon_gamma_threshold, "Species",ispec)) + thisSpecies->radiation_photon_gamma_threshold_, "Species",ispec)) { - thisSpecies->radiation_photon_gamma_threshold = 2.; + thisSpecies->radiation_photon_gamma_threshold_ = 2.; } MESSAGE(2,"> Photon energy threshold for macro-photon emission: " - << thisSpecies->radiation_photon_gamma_threshold); + << thisSpecies->radiation_photon_gamma_threshold_); } else { @@ -714,8 +714,8 @@ class SpeciesFactory { newSpecies->pusher = species->pusher; newSpecies->radiation_model = species->radiation_model; newSpecies->radiation_photon_species = species->radiation_photon_species; - newSpecies->radiation_photon_sampling = species->radiation_photon_sampling; - newSpecies->radiation_photon_gamma_threshold = species->radiation_photon_gamma_threshold; + newSpecies->radiation_photon_sampling_ = species->radiation_photon_sampling_; + newSpecies->radiation_photon_gamma_threshold_ = species->radiation_photon_gamma_threshold_; newSpecies->photon_species = species->photon_species; newSpecies->speciesNumber = species->speciesNumber; newSpecies->position_initialization_on_species = species->position_initialization_on_species; diff --git a/src/Species/SpeciesV.cpp b/src/Species/SpeciesV.cpp index 308fee1e4..c405f00c9 100644 --- a/src/Species/SpeciesV.cpp +++ b/src/Species/SpeciesV.cpp @@ -213,7 +213,7 @@ int ithread; nrj_radiation += Radiate->getRadiatedEnergy(); // Update the quantum parameter chi - Radiate->compute_thread_chipa(*particles, + Radiate->computeParticlesChi(*particles, smpi, first_index[ibin], last_index[ibin],