diff --git a/src/collision/atomic.cc b/src/collision/atomic.cc new file mode 100644 index 00000000..fafe31e5 --- /dev/null +++ b/src/collision/atomic.cc @@ -0,0 +1,502 @@ +#include "atomic.h" + +// We are splitting the colllsion process up into discrete sections +// because some processes only use part, and others need to do work +// in-between. +#define CMOV(a,b) if(t0m * _spj->m)/(_spi->m + _spj->m); \ + const float mu_i = _spj->m/(_spi->m + _spj->m); \ + const float mu_j = _spi->m/(_spi->m + _spj->m); \ + \ + float E, dd, R, R1, ur, urx, ury, urz, tx, ty, tz, t0, t1, t2, stack[3]; \ + float cmx, cmy, cmz; \ + int d0, d1, d2; \ + \ + R = 1; \ + R1 = 0; \ + \ + /* Relative velocity. Needed for all collisions. */ \ + urx = _pi->ux - _pj->ux; \ + ury = _pi->uy - _pj->uy; \ + urz = _pi->uz - _pj->uz; \ + \ + /* Center of mass velocity. Needed for inelastic collisions. */ \ + cmx = mu_j*_pi->ux + mu_i*_pj->ux; \ + cmy = mu_j*_pi->uy + mu_i*_pj->uy; \ + cmz = mu_j*_pi->uz + mu_i*_pj->uz; \ + \ + /* There are lots of ways to formulate T vector formation */ \ + /* This has no branches (but uses L1 heavily) */ \ + \ + t0 = urx*urx; d0=0; d1=1; d2=2; t1=t0; ur = t0; \ + t0 = ury*ury; CMOV(d0,1); CMOV(d1,2); CMOV(d2,0); CMOV(t1,t0); ur += t0; \ + t0 = urz*urz; CMOV(d0,2); CMOV(d1,0); CMOV(d2,1); ur += t0; \ + E = (AP)->alpha_m2*mu*ur; \ + ur = sqrtf( ur ); \ + \ + stack[0] = urx; \ + stack[1] = ury; \ + stack[2] = urz; \ + t1 = stack[d1]; \ + t2 = stack[d2]; \ + t0 = 1 / sqrtf( t1*t1 + t2*t2 + FLT_MIN ); \ + stack[d0] = 0; \ + stack[d1] = t0*t2; \ + stack[d2] = -t0*t1; \ + tx = stack[0]; \ + ty = stack[1]; \ + tz = stack[2] + +#define SCATTER(DELTA,PHI) \ + dd = (DELTA); \ + if(!isfinite(dd) || dd > 1.30e19f) dd = 1.30e19f; \ + t0 = 2*dd/(1+dd*dd); \ + t1 = (PHI); \ + t2 = t0*sin(t1); \ + t1 = t0*ur*cos(t1); \ + t0 *= -dd; \ + \ + stack[0] = (t0*urx + t1*tx) + t2*( ury*tz - urz*ty ); \ + stack[1] = (t0*ury + t1*ty) + t2*( urz*tx - urx*tz ); \ + stack[2] = (t0*urz + t1*tz) + t2*( urx*ty - ury*tx ) + +#define REDUCE_ENERGY(GAMMA) \ + R = (GAMMA); \ + R1 = 1-R \ + +#define FINALIZE_COLLISION(TYPE) \ + d0 = (TYPE); \ + if( d0 & 1 ){ \ + _pi->ux = (_pi->ux + mu_i*stack[0])*R + R1*cmx; \ + _pi->uy = (_pi->uy + mu_i*stack[1])*R + R1*cmy; \ + _pi->uz = (_pi->uz + mu_i*stack[2])*R + R1*cmz; \ + } \ + if( d0 & 2 ){ \ + _pj->ux = (_pj->ux - mu_j*stack[0])*R + R1*cmx; \ + _pj->uy = (_pj->uy - mu_j*stack[1])*R + R1*cmy; \ + _pj->uz = (_pj->uz - mu_j*stack[2])*R + R1*cmz; \ + } \ +} while(0) + +#define ELASTIC_COLLISION(AP, SPI, SPJ, PI, PJ, RNG, TYPE, DELTA) \ + BEGIN_COLLISION(AP,SPI,SPJ,PI,PJ); \ + SCATTER( DELTA(RNG, E), 2*M_PI*frand_c0(RNG) ); \ + FINALIZE_COLLISION(TYPE) + +#define INELASTIC_COLLISION(AP, SPI, SPJ, PI, PJ, RNG, TYPE, DELTA, GAMMA) \ + BEGIN_COLLISION(AP,SPI,SPJ,PI,PJ); \ + SCATTER( DELTA(RNG, E), 2*M_PI*frand_c0(RNG) ); \ + REDUCE_ENERGY( GAMMA(RNG, E) ); \ + FINALIZE_COLLISION(TYPE) + +// Convenience factory for binary rate constants. +#define atomic_rate_constant(sigma) do { \ + float dux = pi->ux - pj->ux; \ + float duy = pi->uy - pj->uy; \ + float duz = pi->uz - pj->uz; \ + float du2 = dux*dux + duy*duy + duz*duz; \ + float mu = (spi->m*spj->m)/(spi->m+spj->m); \ + return ap->Nca02 * (*(sigma))(ap->alpha_m2*mu*du2) * sqrtf(du2); \ +} while(0); + +#define rate_factory(ABXX) \ +float \ +ABXX##_rate_constant( const atomic_properties_t * RESTRICT ap, \ + const species_t * RESTRICT spi, \ + const species_t * RESTRICT spj, \ + const particle_t * RESTRICT pi, \ + const particle_t * RESTRICT pj ) { \ + atomic_rate_constant(ap->sigma_##ABXX); \ +} + +// Convenience factory for binary collisions. We should _not_ use +// restricted pointers because technically they are aliased in the +// macros. +#define inelastic_collision_factory(ABXX, DELTA) \ +void \ +ABXX##_collision( const atomic_properties_t * ap, \ + const species_t * spi, \ + const species_t * spj, \ + /**/ particle_t * pi, \ + /**/ particle_t * pj, \ + /**/ rng_t * rng, \ + const int type ) { \ + INELASTIC_COLLISION(ap, spi, spj, pi, pj, rng, type, \ + ap->delta_##DELTA, ap->gamma_##ABXX ); \ +} + +#define elastic_collision_factory(ABXX) \ +void \ +ABXX##_collision( const atomic_properties_t * ap, \ + const species_t * spi, \ + const species_t * spj, \ + /**/ particle_t * pi, \ + /**/ particle_t * pj, \ + /**/ rng_t * rng, \ + const int type ) { \ + ELASTIC_COLLISION(ap, spi, spj, pi, pj, rng, type, ap->delta_##ABXX); \ +} + +// Now actually create the wrappers. +rate_factory(pt_el); +rate_factory(pt_cx); +rate_factory(pt_ex); +rate_factory(et_el); +rate_factory(et_ex); +rate_factory(tt_el); +rate_factory(tt_ex); + +elastic_collision_factory(pt_el); +elastic_collision_factory(pt_cx); +elastic_collision_factory(et_el); +elastic_collision_factory(tt_el); + +// Inelastic collisions use the scattering angle +// from their elastic counterpart +inelastic_collision_factory(pt_ex, pt_el); +inelastic_collision_factory(et_ex, et_el); +inelastic_collision_factory(tt_ex, tt_el); + +#undef rate_factory +#undef elastic_collision_factory +#undef inelastic_collision_factory + +// Ionization and recombination collisions are special +// and use the chemical collision operator +float +ep_rx_rate_constant( atomic_properties_t * RESTRICT ap, + species_t ** RESTRICT reactants, + particle_t ** RESTRICT particles ) { + particle_t * RESTRICT pi = particles[0]; + particle_t * RESTRICT pj = particles[1]; + species_t * RESTRICT spi = reactants[0]; + species_t * RESTRICT spj = reactants[1]; + atomic_rate_constant( ap->sigma_ep_rx ); +} + +void +ep_rx_collision( atomic_properties_t * RESTRICT ap, + species_t ** RESTRICT reactants, + particle_t ** RESTRICT reactant_particles, + species_t ** RESTRICT products, + particle_t ** RESTRICT product_particles, + rng_t * RESTRICT rng, + int update ) { + + // In free recombination, momentum conservation tells us exactly what the + // product velocity should be (assuming photon momentum is negligible). + particle_t * RESTRICT pi = reactant_particles[0]; + particle_t * RESTRICT pj = reactant_particles[1]; + particle_t * RESTRICT pk = product_particles[0]; + species_t * RESTRICT spi = reactants[0]; + species_t * RESTRICT spj = reactants[1]; + species_t * RESTRICT spk = products[0]; + + float mi = spi->m / spk->m; + float mj = spj->m / spk->m; + + pk->ux = mi * pi->ux + mj * pj->ux ; + pk->uy = mi * pi->uy + mj * pj->uy ; + pk->uz = mi * pi->uz + mj * pj->uz ; + +} + +float +et_in_rate_constant( atomic_properties_t * RESTRICT ap, + species_t ** RESTRICT reactants, + particle_t ** RESTRICT particles ) { + particle_t * RESTRICT pi = particles[0]; + particle_t * RESTRICT pj = particles[1]; + species_t * RESTRICT spi = reactants[0]; + species_t * RESTRICT spj = reactants[1]; + atomic_rate_constant( ap->sigma_et_in ); +} + +void +et_in_collision( atomic_properties_t * RESTRICT ap, + species_t ** RESTRICT reactants, + particle_t ** RESTRICT reactant_particles, + species_t ** RESTRICT products, + particle_t ** RESTRICT product_particles, + rng_t * RESTRICT rng, + int update ) { + + // FIXME: Designation of electron, target, parent are hardcoded + // and assuming the order specified in the public interface + // below. + + particle_t * pe1 = reactant_particles[0]; + particle_t * pt = reactant_particles[1]; + particle_t * pe2 = product_particles[0]; + particle_t * pp = product_particles[1]; + + const species_t * spe1 = reactants[0]; + const species_t * spt = reactants[1]; + const species_t * spe2 = products[0]; + const species_t * spp = products[1]; + + float me2, mt, dE, gamma; + int type; + + // Scatter the primary electron off of the target and reduce the energy, + // but only update the electron if we are told to. In the current + // implementation, we are assuming that we can write the triple differential + // cross-section (d^3 sigma(E1)/ dE2 dOmega1 dOmega2) as the seperable product + // of the singly differential cross section d sigma(E1)/ dE2 and the + // elastic scattering cross-sections. This is not a good approximation in + // general. It would be better to fully implement the triply differential + // cross-section, but data is sparse. Testing this for now. + + if( update & 1 ) type = 1; type += 2; + BEGIN_COLLISION(ap, spe1, spt, pe1, pt); + REDUCE_ENERGY(ap->gamma_et_in(rng, E, &dE)); + SCATTER( ap->delta_et_el(rng, dE), 2*M_PI*frand_c0(rng) ); // Scatter pe2 + + dE = sqrtf(dE*mu/(E*spe2->m)); // Secondary v/ur + pe2->ux = cmx - dE*stack[0]; // Initialize secondary ux moving with pt + pe2->uy = cmy - dE*stack[1]; // Initialize secondary uy moving with pt + pe2->uz = cmz - dE*stack[2]; // Initialize secondary uz moving with pt + + SCATTER( ap->delta_et_el(rng, E), 2*M_PI*frand_c0(rng) ); // Scatter pe1 + FINALIZE_COLLISION(type); + + // Finally, copy target momentum to parent and reduce it by the secondary + // momentum. Note that if implemented correctly at the user level, + // m_p + m_e = m_t in order to ensure mass conservation, but even if this + // is not true, we always conserve momentum. How this is specifically + // implemented ie debatble. + + me2 = spe2->m / spp->m; + mt = spt->m / spp->m; + pp->ux = mt * pt->ux - me2 * pe2->ux; + pp->uy = mt * pt->uy - me2 * pe2->uy; + pp->uz = mt * pt->uz - me2 * pe2->uz; + +} + +float +eep_rx_rate_constant( atomic_properties_t * RESTRICT ap, + species_t ** RESTRICT reactants, + particle_t ** RESTRICT particles ){ + + ERROR(("Not implemented.")); + return 0; + +} + +void +eep_rx_collision( atomic_properties_t * RESTRICT ap, + species_t ** RESTRICT reactants, + particle_t ** RESTRICT reactant_particles, + species_t ** RESTRICT products, + particle_t ** RESTRICT product_particles, + rng_t * RESTRICT rng, + int update ) { + + ERROR(("Not implemented.")); + +} + +#undef CMOV +#undef BEGIN_COLLISION +#undef SCATTER +#undef REDUCE_ENERGY +#undef FINALIZE_COLLISION +#undef ELASTIC_COLLISION +#undef INELASTIC_COLLISION +#undef atomic_rate_constant + + +// Helper routine to decide how to pick the sampling fraction for each process. +double +sampling_fraction( const atomic_properties_t * RESTRICT ap, + const species_t * RESTRICT spi, + const species_t * RESTRICT spj, + /**/ float (*const sigma)(float), + const int interval, + const double Emax, + const double Emin, + const int nsample){ + + double E, K, dE, Kmax = 0; + const float mu = (spi->m*spj->m)/(spi->m+spj->m); + + // Find the maxmium sigma v over the energy range using brute force. + E = Emin; + dE = (Emax-Emin)/nsample; + for(int i=0 ; i<=nsample ; ++i, E += dE){ + K = (*sigma)(ap->alpha_m2*E) * sqrt(E); + Kmax = K > Kmax ? K : Kmax; + } + Kmax *= ap->Nca02 * sqrt(2/mu) ; + + // Now figure out sampling fraction. The number of collisions between species + // i and j is equal to sample np = sample ni nj, which should be larger than + // max(wi,wj) ni nj Kmax interval dt / dV. So, + // sample > Kmax dtinterval_dV * max(wi,wj) + + if( spi->np && spj->np ){ + // Implicit assumption that all particles in each species have the same weight. + double wi = spi->p[0].w; + double wj = spj->p[0].w; + + // Add a factor of 2 for safety. + Kmax *= (wi > wj ? wi : wj) * spi->g->dt * interval / spi->g->dV; + return Kmax; + } + else { + WARNING(("Cannot determine an appropriate sampling fraction for collisions " + "between %s and %s because one (or more) species are empty. " + "Setting to 1\n", spi->name, spj->name)); + return 1; + } +} + +void +checkpt_atomic_properties( const atomic_properties_t * ap ){ + CHECKPT(ap, 1); + CHECKPT_SYM( ap->sigma_pt_el ); CHECKPT_SYM( ap->delta_pt_el ); + CHECKPT_SYM( ap->sigma_pt_cx ); CHECKPT_SYM( ap->delta_pt_cx ); + CHECKPT_SYM( ap->sigma_et_el ); CHECKPT_SYM( ap->delta_et_el ); + CHECKPT_SYM( ap->sigma_tt_el ); CHECKPT_SYM( ap->delta_tt_el ); + CHECKPT_SYM( ap->sigma_pt_ex ); CHECKPT_SYM( ap->gamma_pt_ex ); + CHECKPT_SYM( ap->sigma_et_ex ); CHECKPT_SYM( ap->gamma_et_ex ); + CHECKPT_SYM( ap->sigma_et_in ); CHECKPT_SYM( ap->gamma_et_in ); + CHECKPT_SYM( ap->sigma_tt_ex ); CHECKPT_SYM( ap->gamma_tt_ex ); + CHECKPT_SYM( ap->sigma_ep_rx ); +} + +atomic_properties_t * +restore_atomic_properties( void ){ + atomic_properties_t * ap; + RESTORE( ap ); + RESTORE_SYM( ap->sigma_pt_el ); RESTORE_SYM( ap->delta_pt_el ); + RESTORE_SYM( ap->sigma_pt_cx ); RESTORE_SYM( ap->delta_pt_cx ); + RESTORE_SYM( ap->sigma_et_el ); RESTORE_SYM( ap->delta_et_el ); + RESTORE_SYM( ap->sigma_tt_el ); RESTORE_SYM( ap->delta_tt_el ); + RESTORE_SYM( ap->sigma_pt_ex ); RESTORE_SYM( ap->gamma_pt_ex ); + RESTORE_SYM( ap->sigma_et_ex ); RESTORE_SYM( ap->gamma_et_ex ); + RESTORE_SYM( ap->sigma_et_in ); RESTORE_SYM( ap->gamma_et_in ); + RESTORE_SYM( ap->sigma_tt_ex ); RESTORE_SYM( ap->gamma_tt_ex ); + RESTORE_SYM( ap->sigma_ep_rx ); + return ap; +} + +/* Public interface **********************************************************/ + +#define add_binary_op(SPI,SPJ,ABXX) if(SPI&&SPJ&&ap->sigma_##ABXX) \ +append_collision_op( \ + binary_collision_model( \ + #ABXX, \ + (binary_rate_constant_func_t) ABXX##_rate_constant, \ + (binary_collision_func_t) ABXX##_collision, \ + ap, SPI, SPJ, rp, \ + sample > 0 ? sample : sampling_fraction(ap,SPI,SPJ,ap->sigma_##ABXX, \ + interval,Emax,0,100), \ + interval, \ + mass_action ), \ + cop_list \ + ) + +collision_op_t * +append_atomic_collision_operators( /**/ unsigned int processes, + /**/ species_t * RESTRICT electron, + /**/ species_t * RESTRICT target, + /**/ species_t * RESTRICT parent, + const atomic_properties_t * RESTRICT base, + /**/ field_array_t * RESTRICT fa, + /**/ rng_pool_t * RESTRICT rp, + /**/ collision_op_t ** cop_list, + const double Na02, + const double alpha, + const int interval, + const double sample, + const double Emax + ) { + + if( !electron || !electron->g || !target || !target->g || + !parent || !parent->g || !base || !rp ) ERROR(( "Bad args" )); + + if( sample <= 0 && Emax <= 0 ) ERROR(("Either sample or Emax must be given.")); + + atomic_properties_t *ap; + MALLOC( ap, 1 ); + COPY( ap, base, 1 ); + + // Set up properties. + ap->alpha_m2 = 1/(alpha*alpha); + ap->Nca02 = Na02 * electron->g->cvac; + + // Register the object. + REGISTER_OBJECT( ap, + checkpt_atomic_properties, + restore_atomic_properties, NULL ); + + // Create the collision operators. + if( processes & parent_target_elastic ) add_binary_op(parent, target, pt_el); + if( processes & parent_target_inelastic ) add_binary_op(parent, target, pt_ex); + if( processes & parent_target_charge_exchange ) add_binary_op(parent, target, pt_cx); + if( processes & electron_target_elastic ) add_binary_op(electron, target, et_el); + if( processes & electron_target_inelastic ) add_binary_op(electron, target, et_ex); + if( processes & target_target_elastic ) add_binary_op(target, target, tt_el); + if( processes & target_target_inelastic ) add_binary_op(target, target, tt_ex); + + // Non-number preserving operators. These use the chemical collision operator. + + // Ionization. + if( processes & electron_target_ionization && ap->sigma_et_in ){ + species_t *reactants[2] = { electron, target }; + species_t *products[2] = { electron, parent }; + int consumable[2] = {0, 1}; + double s = sample; + + if( sample <= 0 ) + s = sampling_fraction(ap, electron, target, ap->sigma_et_in, + interval, Emax, 0, 100); + + collision_op_t * et_in = chemical_collision_model("et_in", ap, + reactants, consumable, 2, + products, 2, + (chemical_rate_constant_func_t) et_in_rate_constant, + (chemical_collision_func_t) et_in_collision, + rp, fa, s, 1, interval); + append_collision_op(et_in, cop_list); + } + + // Free recombination. + if( processes & electron_parent_recombination && ap->sigma_ep_rx ){ + species_t *reactants[2] = { electron, parent }; + species_t *products[1] = { target }; + int consumable[2] = {1, 1}; + double s = sample; + + if( sample <= 0 ) + // Divide by 0 at E=0 ! + s = sampling_fraction(ap, electron, parent, ap->sigma_ep_rx, + interval, Emax, 1e-2*Emax, 99); + + collision_op_t * ep_rx = chemical_collision_model("ep_rx", ap, + reactants, consumable, 2, + products, 1, + (chemical_rate_constant_func_t) ep_rx_rate_constant, + (chemical_collision_func_t) ep_rx_collision, + rp, fa, s, 1, interval); + append_collision_op(ep_rx, cop_list); + } + + // Three-body recombination. + if( processes & three_body_recombination ){ + ERROR(("Three body recombination is not implemented.")); + } + + return *cop_list; +} + +#undef add_binary_op diff --git a/src/collision/atomic.h b/src/collision/atomic.h new file mode 100644 index 00000000..d5dae7ab --- /dev/null +++ b/src/collision/atomic.h @@ -0,0 +1,85 @@ +#ifndef _atomic_h_ +#define _atomic_h_ +#include "collision.h" + +typedef struct atomic_properties { + float alpha_m2; // alpha^-2 = 2*Ry/me c^2 + float Nca02; // n0 c a_0^2 has units of 1/time + + // Atomic data routines. They are denoted by _ab_xx where a and b are the + // species [ e: electrons, z: primary ions, z1: parent ions ] and xx are the processes: + // el : elastic scattering + // cx : charge exchange (implemented as backward-dominated scattering) + // ex : excitation (inelastic scattering) + // in : ionization + // rx : recombination + // Note that elastic scattering cross-sections only make sense when z=0, such + // that the primary ion is actually a neutral atom. Otherwise, elastic + // collisions will likely be dominated by Coulomb collisions and a different + // collisional model should be applied. + + // Cross sections. Each should take on input the collision energy in units of + // Ry and return the cross-section in units of a0^2. + float (*sigma_pt_el)(float E); + float (*sigma_pt_cx)(float E); + float (*sigma_pt_ex)(float E); + float (*sigma_et_el)(float E); + float (*sigma_et_ex)(float E); + float (*sigma_et_in)(float E); + float (*sigma_ep_rx)(float E); + float (*sigma_tt_el)(float E); + float (*sigma_tt_ex)(float E); + + // Delta = tan(theta/2). Each should take on input the collision energy. + float (*delta_pt_el)(rng_t* rng, float E); + float (*delta_pt_cx)(rng_t* rng, float E); + float (*delta_et_el)(rng_t* rng, float E); + float (*delta_tt_el)(rng_t* rng, float E); + + // Gamma is the coefficient of restitution for inelastic processes, defined + // by gamma = sqrt(KE'/KE) where KE' is the post-collision energy and KE + // is the collision energy. Each should take on input the collision energy. + float (*gamma_pt_ex)(rng_t* rng, float E); + float (*gamma_et_ex)(rng_t* rng, float E); + float (*gamma_tt_ex)(rng_t* rng, float E); + + // For ionization, we also set the secondary energy (in Ry) + float (*gamma_et_in)(rng_t* rng, float Eprimary, float *Esecondary); + +} atomic_properties_t; + +enum atomic_processes { + parent_target_elastic = 1 << 0, + parent_target_charge_exchange = 1 << 1, + parent_target_inelastic = 1 << 2, + electron_target_elastic = 1 << 3, + electron_target_inelastic = 1 << 4, + electron_target_ionization = 1 << 5, + electron_parent_recombination = 1 << 6, + target_target_elastic = 1 << 7, + target_target_inelastic = 1 << 8, + three_body_recombination = 1 << 9, + all_atomic_processes = 0xFFFF +}; + +BEGIN_C_DECLS + +collision_op_t * +append_atomic_collision_operators( /**/ unsigned int processes, + /**/ species_t * RESTRICT electron, + /**/ species_t * RESTRICT target, + /**/ species_t * RESTRICT parent, + const atomic_properties_t * RESTRICT base, + /**/ field_array_t * RESTRICT fa, + /**/ rng_pool_t * RESTRICT rp, + /**/ collision_op_t ** cop_list, + const double Na02, + const double alpha, + const int interval, + const double sample, + const double Emax + ); + +END_C_DECLS + +#endif /* _atomic_h_ */ diff --git a/src/collision/atomic_models/HeI.h b/src/collision/atomic_models/HeI.h new file mode 100644 index 00000000..55320fdf --- /dev/null +++ b/src/collision/atomic_models/HeI.h @@ -0,0 +1,341 @@ +#ifndef _HeI_h_ +#define _HeI_h_ +#include +#include "atom.h" + +/* Atomic data routines for collisions with Helium atoms. + + The coefficents used here are based on the methodology developed by + Khrabrov and Kaganovich and Wang et al. 2017. The basic idea is to + construct a plausible physical model for the differential scattering + cross-section and then use experimental measurements of the total, + momentum, and/or viscous cross-sections to determine free-parameters + in the model. Specifically, we use the following models: + + en_el : Scattering in a screened Coulomb potential with an energy + dependent screening-length. + et_ex : Cross-sections are a degenercy averaged cross section for + excitation from 11S -> nls with n <= 4. Raw cross-sections + are taken from the analytic fits of Ralchenko 2008. Gamma + is chosen in a consistent manner. Scattering angle uses + the elastic angles. + et_in : Cross-section is direct ionization from 11S taken from + Ralchenko 2008. + ep_rx : Recombination is tricky. At high energy, resonances due to + dielectronic recombination dominate the total cross-section, + but are difficult to include directly since they require very + high resolution in energy space to capture properly. Here, + raw data from Nahar 2010 is used and pre-processed to give + a correct resonace averaged structure. + in_el/cx : Forward scattering due to a polarization potential and + backwards scattering due to charge exchange. Forward and + backward scattering amplitudes are independent with energy + dependent amplitudes. + nn_el : Normally distributed forward scattering with variance + chosen to match total and viscous cross-sections. + + For details on how the scattering angles are chosen, see Khraborv and + Kaganovich and Wang et al. However, here we have changed some details + in the actual implementation. Specifically, the references cited + specifically compute cos(theta) of the scattering angle, however this + is not numerically stable for small theta and will not conserve + momentum. Here we compute delta = tan(theta/2) instead and use that to + then compute sin(theta) and cos(theta) in the actual collision algorithm. + This way minimizes the effect of roundoff, but the tradeoff is that we + have to worry about delta going to infinity for perfect backscattering. + Also, instead of using look-up tables, we are using polynomial + interpolants in log2 space. + + References: + Nahar, S.N., New Astronomy (15), 2010 + Ralchenko, Y. et al., Atomic Data and Nuclear Data Tables (94), 2008 + +*/ + + +float +HeI_sigma_pt_el(float E){ + E = log2(E + 7.349865e-04f); + return exp2(((((( 6.012956e-06 *E +2.815147e-05)*E + -3.985465e-04)*E -1.499984e-03)*E + +9.442074e-03)*E -2.464035e-01)*E +7.177672e+00); +} + +float +HeI_sigma_pt_cx(float E){ + E = log2(E + 7.349865e-04f); + return exp2(((((( 4.581837e-06 *E +1.009797e-05)*E + -3.056459e-04)*E -2.674257e-04)*E + +4.621865e-04)*E -1.343024e-01)*E +6.043904e+00); +} + +float +HeI_sigma_et_el(float E){ + E = log2(E + 7.349865e-04f); + return exp2(((((( 4.307671e-06 *E +1.273896e-04)*E + +8.749181e-04)*E -6.549400e-03)*E + -1.136261e-01)*E -4.852011e-01)*E +3.771030e+00); +} + +float +HeI_sigma_et_ex(float E){ + float x = E/1.4568; + if (E < 1) return 0; + + // This is a sum of all excitation cross-sections out of 11S. + // This may look very long, but is actually shorter than directly + // computing the first 4 directly (which is needed for 10% accuracy) + x = log2(x); + return exp2((((((((((( 6.555543e-05 *x -2.368052e-03)*x + +3.769421e-02)*x -3.476701e-01)*x + +2.054920e+00)*x -8.124464e+00)*x + +2.175721e+01)*x -3.904986e+01)*x + +4.534678e+01)*x -3.196479e+01)*x + +1.249841e+01)*x -3.129806e+00); + +} + +float +HeI_sigma_et_in(float E){ + if( E <= 1.807140 ) return 0; + + // Equation 9 from Ralchenko 2008. + float x = 1.807140/E; // I/E + float y = 1-x; // (1-I/E) + return 5.9071e+00*x*(-5.857e-01*log(x)+(((3.317e+00 *y -2.521e+00)*y + +7.680e-01)*y -4.457e-01)*y); +} + +float +HeI_sigma_ep_rx(float E){ + + // Radiative recombinaton. Simple but 5% accurate for E >~ 2e-3 + float sigma = 0.9e-5*pow(E+7.349865e-04f, -1.125); + if( E < 2 || E > 4 ) return sigma; + + // Resonance averaged dielectronic recombination. These resonances + // are very sharp, but I am averaging them over an energy interval + // of 0.05 Ry to smear them out. Note when doing averaging, we are + // really doing = / v so that the total rate is + // conserved (assuming that <.> denotes convolution with a window + // of unit area). This is a little overkill here, we could get + // away with only 3 peaks probably, but why not. + + float x; + x = E - 2.889545; + sigma += 2.01915e-04 * exp( -96.345*x*x ); + x = E - 2.994697; + sigma += 4.39237e-03 * exp( -195.95*x*x ); + x = E - 3.703070; + sigma += 1.55912e-05 * exp( -21.751*x*x ); + x = E - 2.584947; + sigma += 4.86425e-05 * exp( -70.603*x*x ); + x = E - 3.551628; + sigma += 1.89497e-04 * exp( -214.397*x*x ); + + return sigma; + +} + +float +HeI_sigma_tt_el(float E){ + E = log2(E + 7.349865e-04f); + return exp2(((( 1.354527e-05 *E +5.235340e-05)*E + -5.428023e-03)*E -1.122875e-01)*E +6.294592e+00); +} + +float +HeI_delta_pt_el(rng_t * rng, float E){ + const float sqrt_sqrt_half = 0.8408964152537146f; + const float one_8 = 0.125f; + const float five_128 = 0.0390625f; + const float fifteen_1024 = 0.0146484375f; + + float inv_sqrt_sqrt_a, a, x; + E = log2(E + 7.349865e-04f); + a = (((((( -5.692770e-06 *E -1.911994e-04)*E + +3.192144e-05)*E +2.383839e-02)*E + -1.522081e-01)*E -3.641591e+00)*E -2.128630e+01); + inv_sqrt_sqrt_a = exp2(-0.25f*a); + a = exp2(a); + x = inv_sqrt_sqrt_a - sqrt_sqrt_half*(1 - a*(one_8 - a*(five_128 - a*fifteen_1024))); + x = inv_sqrt_sqrt_a - x*frand(rng); + x *= x; + x *= x; + return sqrtf( (1-a*x) / ((2+a)*x - 1) ); +} + +float +HeI_delta_pt_cx(rng_t * rng, float E){ + const float sqrt_sqrt_half = 0.8408964152537146f; + const float one_8 = 0.125f; + const float five_128 = 0.0390625f; + const float fifteen_1024 = 0.0146484375f; + + float inv_sqrt_sqrt_b, b, x; + E = log2(E + 7.349865e-04f); + b = (((((( -1.070780e-05 *E -2.346134e-04)*E + +6.388344e-04)*E +2.652242e-02)*E + -1.743078e-01)*E -3.978732e+00)*E -1.685851e+01); + inv_sqrt_sqrt_b = exp2(-0.25f*b); + b = exp2(b); + x = inv_sqrt_sqrt_b - sqrt_sqrt_half*(1 - b*(one_8 - b*(five_128 - b*fifteen_1024))); + x = inv_sqrt_sqrt_b + x*(frand(rng)-1); + x *= x; + x *= x; + return sqrtf(2*x/(1-b*x) - 1); +} + +float +HeI_delta_et_el(rng_t * rng, float E){ + float x,z; + E *= 13.605693f; + z = 2.45f * sqrtf(E); + z = 1 + (z-19.9324f)/(E-2.302041f*z+19.9324f) - z/(E-4.171429f*z+90.1221f); + x = 1-2*frand(rng); + x = (z + x)/(1 + z*x); + return sqrtf((1-x)/(1+x)); +} + +float +HeI_delta_tt_el(rng_t * rng, float E){ + E = log2(E + 7.349865e-04f); + return frandn(rng) * exp2((((-8.664093e-05 *E -2.570293e-03)*E + -2.534166e-02)*E -2.399119e-01)*E + -2.575432e+00); +} + +float +HeI_gamma_et_ex(rng_t *rng, float E){ + // At each energy, we should choose the excitation energy E_i with + // probability sigma_i / sigma_tot, However, this requires evaluating + // many cross-sections. Instead, we just use (E), then compute + // gamma(E) = sqrt(1-/E) + float x = E/1.4568 - 1; + if(x <= 0) return x < 0 ; + x = log2(x); + x = exp2((2.327251e-02*x +8.498602e-01)*x +1.086187e+00); + return x/(1+x); +} + +float +HeI_gamma_et_in(rng_t *rng, float Eprimary, float *Esecondary){ + + // In the non-relativistic binary-encounter-dipole model, the singly + // differential cross-section is given by + // + // d sigma(W,T)/ dW = (S/B)/(t+u+1) [ ((Ni/N)-2)/(t+1) (1/(w+1) + 1/(t-w)) + // + (2-Ni/N) (1/(w+1)^2 + 1/(t-w)^2) + // + ln(t)/(N (w+1)) df/dw ] + // + // where : + // W = Energy of secondary electron + // B = Binding energy of electron + // T = Reduced kinetic energy of primary me v^2 / 2 + // N = Number of bound electrons in the subshell + // U = Average kinetic energy of bound electrons + // df/dw = differential oscillator strength + // Ni = int_0^inf df/dw dw = total oscillator strength + // and lower case letters denote normalization by B. + // + // For ground state He: + // u = 1.607 + // N = 2 + // df/dw = 8.24012e+00 y^3 -1.04769e+01 y^4 +3.96496e+00 y^5 +4.45976e-02 y^6 + // Ni/N - 2 = -1.18604 + // y = 1/(1+w) + // + // Now integrating the cross-section to obtain the CDF, we find + // + // CDF(w|t) ~ (Ni/N-2) / t+1 [ ln(w+1/t-w) + (t+1) (1/w+1 - 1/t-w) + // + ln(t) - (t+1) (1 - 1/t) ] + // - ln(t)/N sum_i a_i/i ( (1+w)^-i - 1 ) + // + // And, the upper bound on w is given by w <= (t-1)/2 so the proportionality + // coefficient is given by + // + // CDF_norm = (Ni/N-2) / t+1 [ ln(t) - (t+1) (1 - 1/t) ] + // - ln(t)/N sum_i a_i/i ( 2^i (t+1)^-i - 1 ) + // + // Now this CDF is not easy to invert since it depends on t and w in + // non-trivial ways, but we can fairly quickly do some Newton iterations. + + int max_iter = 100; + float CDF, dCDFdw, w, w1, tw, err; + const float a3 = 2.74671e+00; // a3 / 3 + const float a4 = -2.61923e+00; // a4 / 4 + const float a5 = 7.92992e-01; // a5 / 5 + const float a6 = 7.43293e-03; // a6 / 6 + const float t = Eprimary/1.807140; // t = E/I + const float t1 = t+1; + const float c2 = 0.5*log(t); // ln(t)/N + const float c1 = -1.18604 / t1; // (Ni/N - 2)/(t+1) + const float c0 = c1 * (2*c2 - t1*(1-1/t)) + c2 * (a3 + a4 + a5 + a6); + + w1 = 2/t1; + const float norm = c0 - c2*((((a6*w1 + a5)*w1 + a4)*w1 + a3)*w1*(w1*w1) ); + const float R = frand(rng)*norm; + const float wmax = 0.5*(t-1); + + // SHhould never happen, but test just in case. + if(t < 1){ + (*Esecondary) = 0; + return 0; + } + + // This is a good starting guess. + w = R*wmax/(t*norm); + for( ; max_iter ; --max_iter ) { + w1 = 1/(w+1); + tw = 1/(t-w); + CDF = c0 + c1*( -log(w1/tw) + t1*(w1 - tw) ) + - c2*((((a6*w1 + a5)*w1 + a4)*w1 + a3)*w1*(w1*w1) ); + dCDFdw = c1*( w1*(1-t1*w1) + tw*(1-t1*tw) ) + + c2*w1*((((6*a6*w1 + 5*a5)*w1 + 4*a4)*w1 + 3*a3)*w1*(w1*w1) ); + err = CDF-R; + w -= err/dCDFdw; + if(w < 0) w = 0; + if(w > wmax) w = wmax; + if( err*err <= 1e-6 ) break; + } + + // Set the secondary energy and compute gamma. + (*Esecondary) = w*1.807140; + return sqrtf(1 - (w+1)/t); +} + +/* Public interface **********************************************************/ + +const atomic_properties_t HeI = { + .alpha_m2 = 0, + .Nca02 = 0, + + // Cross-sections for the various processes. + .sigma_pt_el = &HeI_sigma_pt_el , // He + He+ -> He + He+ + .sigma_pt_cx = &HeI_sigma_pt_cx , // He + He+ -> He+ + He + .sigma_pt_ex = NULL , // He + He+ -> He* + He+ + .sigma_et_el = &HeI_sigma_et_el , // He + e -> He + e + .sigma_et_ex = &HeI_sigma_et_ex , // He + e -> He* + e + .sigma_et_in = &HeI_sigma_et_in , // He + e -> He+ + 2e + .sigma_ep_rx = &HeI_sigma_ep_rx , // He+ + e -> He + .sigma_tt_el = &HeI_sigma_tt_el , // He + He -> He + He + .sigma_tt_ex = NULL , // He + He -> He* + He + + // Scattering angles for the various elastic processes. + // Inelastic processes use the corresponding elastic delta. + .delta_pt_el = &HeI_delta_pt_el , + .delta_pt_cx = &HeI_delta_pt_cx , + .delta_et_el = &HeI_delta_et_el , + .delta_tt_el = &HeI_delta_tt_el , + + // Coefficient of restitution for the various processes. + .gamma_pt_ex = NULL , + .gamma_et_ex = &HeI_gamma_et_ex , + .gamma_tt_ex = NULL , + + // Ionization coefficient of restitution and secondary energy. + .gamma_et_in = &HeI_gamma_et_in , + +} ; + +#endif /* _HeI_h_ */ diff --git a/src/collision/binary.cc b/src/collision/binary.cc index 2fb761fc..5c5a8221 100644 --- a/src/collision/binary.cc +++ b/src/collision/binary.cc @@ -91,7 +91,8 @@ binary_collision_model( const char * RESTRICT name, species_t * RESTRICT spj, rng_pool_t * RESTRICT rp, double sample, - int interval ) + int interval, + int sample_strategy ) { binary_collision_model_t * cm; @@ -119,10 +120,16 @@ binary_collision_model( const char * RESTRICT name, ERROR( ( "collision model parameters must be checkpoint registered" ) ); } + if ( sample_strategy != per_particle || + sample_strategy != mass_action ) + { + ERROR( ( "Bad sampling strategy." ) ); + } + MALLOC( cm, 1 ); MALLOC( cm->name, len+1 ); - strcpy( cm->name, name ); + strcpy( cm->name, name ); cm->rate_constant = rate_constant; cm->collision = collision; @@ -132,6 +139,7 @@ binary_collision_model( const char * RESTRICT name, cm->rp = rp; cm->sample = sample; cm->interval = interval; + cm->strategy = sample_strategy; return new_collision_op_internal( cm, ( collision_op_func_t ) apply_binary_collision_model, diff --git a/src/collision/binary.h b/src/collision/binary.h index aa16dd37..6d8b8a92 100644 --- a/src/collision/binary.h +++ b/src/collision/binary.h @@ -15,6 +15,7 @@ typedef struct binary_collision_model double sample; int interval; int n_large_pr[ MAX_PIPELINE ]; + int strategy; } binary_collision_model_t; void diff --git a/src/collision/chemical.cc b/src/collision/chemical.cc new file mode 100644 index 00000000..892505a3 --- /dev/null +++ b/src/collision/chemical.cc @@ -0,0 +1,146 @@ +#define IN_collision + +#include "chemical.h" + +/* Private interface *********************************************************/ + +void +apply_chemical_collision_model( chemical_collision_model_t * cm ) +{ + + species_t * sp; + species_t ** reactants = cm->reactants; + grid_t * g = cm->reactants[0]->g; + + const int N = cm->n_reactants; + int n; + + if ( cm->interval < 1 || ( g->step % cm->interval ) ) + { + return; + } + + FOR_REACTANTS if( sp->last_sorted != g->step ) sort_p( sp ); + + // Conditionally execute this when more abstractions are available. + apply_chemical_collision_model_pipeline( cm ); + +} + +void +checkpt_chemical_collision_model( const collision_op_t * cop ) { + const chemical_collision_model_t * cm = + (const chemical_collision_model_t *)cop->params; + int i; + CHECKPT( cm, 1 ); + CHECKPT_STR( cm->name ); + CHECKPT_SYM( cm->rate_constant ); + CHECKPT_SYM( cm->collision ); + CHECKPT_PTR( cm->params ); + CHECKPT_PTR( cm->rp ); + CHECKPT_PTR( cm->fa ); + CHECKPT( cm->consumable, cm->n_reactants ); + for(i=0 ; i < cm->n_products ; ++i ) CHECKPT_PTR( cm->products[i] ); + for(i=0 ; i < cm->n_reactants ; ++i ) CHECKPT_PTR( cm->reactants[i] ); + + checkpt_collision_op_internal( cop ); +} + +collision_op_t * +restore_chemical_collision_model( void ) { + chemical_collision_model_t * cm; + int i; + RESTORE( cm ); + MALLOC( cm->reactants, cm->n_reactants ); + MALLOC( cm->products, cm->n_products ); + MALLOC( cm->n_produced, MAX_PIPELINE ); + MALLOC( cm->n_modified, cm->n_reactants * MAX_PIPELINE ); + MALLOC( cm->reactant_movers, cm->n_reactants * MAX_PIPELINE ); + MALLOC( cm->product_particles, cm->n_products * MAX_PIPELINE ); + RESTORE_STR( cm->name ); + RESTORE_SYM( cm->rate_constant ); + RESTORE_SYM( cm->collision ); + RESTORE_PTR( cm->params ); + RESTORE_PTR( cm->rp ); + RESTORE_PTR( cm->fa ); + RESTORE( cm->consumable ); + for(i=0 ; i < cm->n_products ; ++i ) RESTORE_PTR( cm->products[i] ); + for(i=0 ; i < cm->n_reactants ; ++i ) RESTORE_PTR( cm->reactants[i] ); + return restore_collision_op_internal( cm ); +} + +void +delete_chemical_collision_model( collision_op_t * cop ) { + chemical_collision_model_t * cm = (chemical_collision_model_t *)cop->params; + FREE( cm->name ); + FREE( cm->reactants ); FREE( cm->reactant_movers ); FREE( cm->n_modified ); + FREE( cm->products ); FREE( cm->product_particles ); FREE( cm->n_produced ); + FREE( cm->consumable ); + FREE( cm ); + delete_collision_op_internal( cop ); +} + +/* Public interface **********************************************************/ + +collision_op_t * +chemical_collision_model( const char * RESTRICT name, + /**/ void * RESTRICT params, + /**/ species_t ** RESTRICT reactants, + const int * RESTRICT consumable, + const int n_reactants, + /**/ species_t ** RESTRICT products, + const int n_products, + /**/ chemical_rate_constant_func_t rate_constant, + /**/ chemical_collision_func_t collision, + /**/ rng_pool_t * RESTRICT rp, + /**/ field_array_t * RESTRICT fa, + const double sample, + const int dynamic_sampling, + const int interval ) +{ + chemical_collision_model_t * cm; + size_t len = name ? strlen(name) : 0; + + // if( !rate_constant || !collision || !spi || !spj || spi->g!=spj->g || + // !rp || rp->n_rng= 1 for chemical collisions.")); + + const int N = n_reactants; + const int M = n_products; + + MALLOC( cm, 1 ); + MALLOC( cm->name, len+1 ); strcpy( cm->name, name ); + MALLOC( cm->reactants, N ); COPY( cm->reactants, reactants, N); + MALLOC( cm->consumable, N ); COPY( cm->consumable, consumable, N); + MALLOC( cm->products, M ); COPY( cm->products, products, M); + MALLOC( cm->n_produced, MAX_PIPELINE ); + MALLOC( cm->n_modified, N*MAX_PIPELINE ); + MALLOC( cm->reactant_movers, N*MAX_PIPELINE ); + MALLOC( cm->product_particles, M*MAX_PIPELINE ); + + cm->n_reactants = N; + cm->n_products = M; + cm->dynamic_sampling = dynamic_sampling; + cm->rate_constant = rate_constant; + cm->collision = collision; + cm->params = params; + cm->rp = rp; + cm->fa = fa; + cm->sample = sample; + cm->interval = interval; + + return new_collision_op_internal( cm, + (collision_op_func_t)apply_chemical_collision_model, + delete_chemical_collision_model, + (checkpt_func_t)checkpt_chemical_collision_model, + (restore_func_t)restore_chemical_collision_model, + NULL ); +} diff --git a/src/collision/chemical.h b/src/collision/chemical.h new file mode 100644 index 00000000..45cedd22 --- /dev/null +++ b/src/collision/chemical.h @@ -0,0 +1,44 @@ +#ifndef _chemical_h_ +#define _chemical_h_ + +#include "collision_private.h" + +typedef struct chemical_collision_model { + char * name; + chemical_rate_constant_func_t rate_constant; + chemical_collision_func_t collision; + void * params; + rng_pool_t * rp; + field_array_t * fa; + int interval; + + // Species information. + int n_reactants, n_products, *consumable; + species_t **reactants, **products; + + // Sampling parameters. + double sample; + int dynamic_sampling; + + // Information about the number of collisions and their probabilities. + int n_large_pr[ MAX_PIPELINE ]; + int n_tested[ MAX_PIPELINE ]; + int unaccumulated[ MAX_PIPELINE ]; + float pr_max[ MAX_PIPELINE ]; + + // Information on the number of particles produced and modified. + int *n_modified; // size = N * MAX_PIPELINE + particle_mover_t **reactant_movers; // size = N * MAX_PIPELINE + int *n_produced; // size = MAX_PIPELINE + particle_t **product_particles; // size = M * MAX_PIPELINE + +} chemical_collision_model_t; + +void +apply_chemical_collision_model_pipeline( chemical_collision_model_t * cm ); + +// Helper macros for looping. +#define FOR_PRODUCTS for( sp=products[0], m=0 ; m P_0 + P_1 + ... P_M + + The N reactants (R_i) may be consumed during the reaction, in which + case their statistical weight should decrease post-collision, or they + may act as catalysts. All products are always generated and are positioned + at the cetner of mass of the reactants. Charge conservation is maintained + by the collision algorithm, and so the specific implementation only needs + to take care of: + 1. Specifying the chemical reaction rate constant + 2. Adjusting momentum and statistical weight for each reactant + 3. Adjusting momentum for each product + + A chemical_rate_constant_func_t returns the lab-frame rate constant + for the collisions between monochromatic beams of reactant species sp[i], + 0m) and momentum p[i]->u{xyz} + (normalized to sp[i]->m sp[i]->g->cvac where cvac is the speed of light + in vacuum).bg + + The returned value has units of VOLUME^{N-1} / TIME. + + Parameters: + params : Pointer to the parameters for this collision_op_t + reactants : Pointers to the N reactant species + particles : Pointers to the N reacting particles +*/ + +typedef float +(*chemical_rate_constant_func_t)( void * RESTRICT params, + species_t ** RESTRICT reactants, + particle_t ** RESTRICT particles ); + +/* A chemical_collision_func_t implements the microscopic physics of a + collision between reactants. + + Since reacting particles may have different weights, not all particle + momenta may need to be updated during this collision. To determine if + a reacting particle should be updated, the control logic is : + + if( update_flag & (1 << n) ) update particles[n]->u{x,y,z} + + Product particles are generated with the appropriate statistical weight by + the calling function. A chemical_collision_func_t only needs to set the + momenta of the products. If any reactant particles are consumed during the + reaction, their statistical weight should be reduced by the statistical + weight of the product particles (which is equal to the minimum of the + reactant weights). Catalyzing reactants should not have their weight + adjusted. Any reactants with weight <= 0 are automatically garbage + collected. + + Parameters: + params : Pointer to the parameters for this collision_op_t + reactants : Pointers to the N reactant species + reactant_particles : Pointers to the N reacting particles + products : Pointers to the M product species + product_particles : Pointers to the M product particles + rng : Random number generator + update_flag : Specifies which reacting particles should be adjusted +*/ + +typedef void +(*chemical_collision_func_t)( void * RESTRICT params, + species_t ** RESTRICT reactants, + particle_t ** RESTRICT reactant_particles, + species_t ** RESTRICT products, + particle_t ** RESTRICT product_particles, + rng_t * RESTRICT rng, + int update_flag ); + +collision_op_t * +chemical_collision_model( const char * RESTRICT name, + /**/ void * RESTRICT params, + /**/ species_t ** RESTRICT reactants, + const int * RESTRICT consumable, + const int n_reactants, + /**/ species_t ** RESTRICT products, + const int n_products, + /**/ chemical_rate_constant_func_t rate_constant, + /**/ chemical_collision_func_t collision, + /**/ rng_pool_t * RESTRICT rp, + /**/ field_array_t * RESTRICT fa, + const double sample, + const int dynamic_sampling, + const int interval ); /* In hard_sphere.c */ @@ -301,7 +397,8 @@ hard_sphere( const char * RESTRICT name, /* Model name */ const float rj, /* Species-j p. radius (LENGTH) */ rng_pool_t * RESTRICT rp, /* Entropy pool */ const double sample, /* Sampling density */ - const int interval ); /* How often to apply this */ + const int interval, /* How often to apply this */ + const int strategy ); /* Sampling strategy */ /* In large_angle_coulomb.c */ @@ -331,7 +428,8 @@ large_angle_coulomb( const char * RESTRICT name, /* Model name */ const float bmax, /* Impact parameter cutoff */ rng_pool_t * RESTRICT rp, /* Entropy pool */ const double sample, /* Sampling density */ - const int interval ); /* How often to apply this */ + const int interval, /* How often to apply this */ + const int strategy ); /* Sampling strategy */ END_C_DECLS diff --git a/src/collision/hard_sphere.cc b/src/collision/hard_sphere.cc index 495d9220..6c210836 100644 --- a/src/collision/hard_sphere.cc +++ b/src/collision/hard_sphere.cc @@ -26,7 +26,7 @@ typedef struct hard_sphere and parallel, vjp, to vr. This yields: K = { [ nj pi R^2 ] / [ vtj^3 (2 pi)^(3/2) ] } integral ( vr^2 - 2|vr|vjp + vjp^2 + vjP^2 )^(1/2) - exp( -vjp^2 / ( 2 vtj^2 ) ) + exp( -vjp^2 / ( 2 vtj^2 ) ) exp( -vjP^2 / ( 2 vtj^2 ) ) 2 pi vjP dvjP dvjp where the vjP integration is from 0 to infinity and the vjp is @@ -34,7 +34,7 @@ typedef struct hard_sphere we have: K = { [ nj pi^2 R^2 ] / [ vtj^3 (2 pi)^(3/2) ] } integral ( vr^2 - 2|vr|vjp + vjp^2 + z )^(1/2) - exp( -vjp^2 / ( 2 vtj^2 ) ) + exp( -vjp^2 / ( 2 vtj^2 ) ) exp( -z / ( 2 vtj^2 ) ) dz dvjp Performing the z integral, with the help of Mathematica, we get: @@ -189,7 +189,7 @@ hard_sphere_rate_constant( const hard_sphere_t * RESTRICT hs, Yielding: dvr = -2 [ 1 - (b/R)^2 ] vr0 + |vr0| 2 (b/R) [ 1 - (b/R)^2 ]^(1/2) T - + To convert the scattering angle above into T, we need to determine to two vectors perpendicular to the incident velocity, vr0. Given the desire to not introduce any preferred directions in the @@ -204,14 +204,14 @@ hard_sphere_rate_constant( const hard_sphere_t * RESTRICT hs, formed as the cross product of VR and T1 (and normalizing). This given: - dvr = -2 [ 1 - (b/R)^2 ] vr0 + + dvr = -2 [ 1 - (b/R)^2 ] vr0 + |vr0| 2 (b/R) [ 1 - (b/R)^2 ]^(1/2) ( cs T1 + sn T2 ) Noting that (b/R) cs and (b/R) sn are just the coordinates of the point picked in the unit circle originally and that T2 = vr0 x T1 / |vr0|, we can simplify the above: - dvr = -2 [ 1 - (b/R)^2 ] vr0 + + dvr = -2 [ 1 - (b/R)^2 ] vr0 + 2 [ 1 - (b/R)^2 ]^(1/2) [ |vr0| bcs_R T1 + bsn_R ( vr0 x T1 ) ] Requiring conservation of momentum (which implies @@ -302,7 +302,7 @@ hard_sphere_fluid_collision( const hard_sphere_t * RESTRICT hs, ury -= w*frandn(rng); urz -= w*frandn(rng); } - + COMPUTE_MOMENTUM_TRANSFER(urx,ury,urz,ax,ay,az,rng); w = hs->twomu_mi; @@ -329,14 +329,14 @@ hard_sphere_collision( const hard_sphere_t * RESTRICT hs, COMPUTE_MOMENTUM_TRANSFER(urx,ury,urz,ax,ay,az,rng); if( type & 1 ) { - w = hs->twomu_mi; + w = hs->twomu_mi; pi->ux -= w*ax; pi->uy -= w*ay; pi->uz -= w*az; } if( type & 2 ) { - w = hs->twomu_mj; + w = hs->twomu_mj; pj->ux += w*ax; pj->uy += w*ay; pj->uz += w*az; @@ -410,7 +410,8 @@ hard_sphere( const char * RESTRICT name, /* Model name */ const float rj, /* Species-j p. radius (LENGTH) */ rng_pool_t * RESTRICT rp, /* Entropy pool */ const double sample, /* Sampling density */ - const int interval ) /* How often to apply this */ + const int interval, /* How often to apply this */ + const int strategy ) /* Sampling strategy */ { hard_sphere_t * hs; @@ -427,5 +428,5 @@ hard_sphere( const char * RESTRICT name, /* Model name */ return binary_collision_model( name, (binary_rate_constant_func_t)hard_sphere_rate_constant, (binary_collision_func_t) hard_sphere_collision, - hs, spi, spj, rp, sample, interval ); + hs, spi, spj, rp, sample, interval, strategy ); } diff --git a/src/collision/large_angle_coulomb.cc b/src/collision/large_angle_coulomb.cc index d2e2845e..230ac8b0 100644 --- a/src/collision/large_angle_coulomb.cc +++ b/src/collision/large_angle_coulomb.cc @@ -44,7 +44,7 @@ large_angle_coulomb_rate_constant( The only difference is that, in a large angle Coulomb collision, the relationship between the impact parameter and b is: - b = ( qi qj cot theta/2 ) / ( 4 pi eps0 mu vr0^2 ) + b = ( qi qj cot theta/2 ) / ( 4 pi eps0 mu vr0^2 ) or, given b: theta = 2 acot ( ( 4 pi eps0 mu vr0^2 b ) / ( qi qj ) ) = 2 acot B @@ -70,16 +70,16 @@ large_angle_coulomb_rate_constant( = 2 cos acot B sin acot B = 2 B / ( B^2+1 ) such that: - dvr = -2 (1 / (B^2+1)) vr0 + dvr = -2 (1 / (B^2+1)) vr0 + 2 (B / (B^2+1)) cos phi |vr0| T1 + 2 (B / (B^2+1)) sin phi vr0 x T1 Normalizing to -2 c: - dur = (1 / (B^2+1)) ur0 + dur = (1 / (B^2+1)) ur0 - (B / (B^2+1)) cos phi |ur0| T1 - (B / (B^2+1)) sin phi ur0 x T1 */ - + #define CMOV(a,b) if(t0twomu_mi; @@ -177,14 +177,14 @@ large_angle_coulomb_collision( COMPUTE_MOMENTUM_TRANSFER(urx,ury,urz,ax,ay,az,rng); if( type & 1 ) { - w = lac->twomu_mi; + w = lac->twomu_mi; pi->ux -= w*ax; pi->uy -= w*ay; pi->uz -= w*az; } if( type & 2 ) { - w = lac->twomu_mj; + w = lac->twomu_mj; pj->ux += w*ax; pj->uy += w*ay; pj->uz += w*az; @@ -230,7 +230,7 @@ large_angle_coulomb_fluid( ERROR(( "Bad args" )); MALLOC( lac, 1 ); - + lac->cc = (4.*M_PI*sp->g->eps0*sp->m*m0*sp->g->cvac*sp->g->cvac*bmax) / ((sp->m + m0)*sp->q*q0); lac->twomu_mi = 2.*m0 / (sp->m + m0); @@ -264,11 +264,12 @@ large_angle_coulomb( const char * RESTRICT name, /* Model name */ const float bmax, /* Impact parameter cutoff */ rng_pool_t * RESTRICT rp, /* Entropy pool */ const double sample, /* Sampling density */ - const int interval ) /* How often to apply this */ + const int interval, /* How often to apply this */ + const int strategy ) /* Sampling strategy */ { large_angle_coulomb_t * lac; - if( !spi || !spi->q || spi->m<=0 || + if( !spi || !spi->q || spi->m<=0 || !spj || !spj->q || spj->m<=0 || spi->g!=spj->g ) ERROR(( "Bad args" )); MALLOC( lac, 1 ); @@ -285,5 +286,5 @@ large_angle_coulomb( const char * RESTRICT name, /* Model name */ return binary_collision_model( name, (binary_rate_constant_func_t)large_angle_coulomb_rate_constant, (binary_collision_func_t) large_angle_coulomb_collision, - lac, spi, spj, rp, sample, interval ); + lac, spi, spj, rp, sample, interval, strategy ); } diff --git a/src/collision/pipeline/binary_pipeline.cc b/src/collision/pipeline/binary_pipeline.cc index e3595309..805f7619 100644 --- a/src/collision/pipeline/binary_pipeline.cc +++ b/src/collision/pipeline/binary_pipeline.cc @@ -37,9 +37,12 @@ binary_pipeline_scalar( binary_collision_model_t * RESTRICT cm, /**/ particle_t * RESTRICT spj_p = spj->p; const int * RESTRICT spj_partition = spj->partition; - const double sample = (spi_p==spj_p ? 0.5 : 1)*cm->sample; const float dtinterval_dV = ( g->dt * (float)cm->interval ) / g->dV; + double sample = cm->sample; + const int strategy = cm->strategy; + if( strategy == per_particle && spi == spj ) sample *= 0.5; + float pr_norm, pr_coll, wk, wl, w_max, w_min; int v, v1, k, k0, nk, rk, l, l0, nl, rl, np, nc, type, n_large_pr = 0; @@ -81,7 +84,13 @@ binary_pipeline_scalar( binary_collision_model_t * RESTRICT cm, nl = nk; rl = rk; np = nk*(nk+1) >> 1; - nc = (int)( 0.5 + sample*(double)nk ); + + switch(strategy) { + case per_particle: nc = (int)( 0.5 + sample*(double)nk ) ; break ; + case mass_action: nc = (int)( 0.5 + sample*(double)np ) ; break ; + default: nc = 0; + } + } else @@ -98,7 +107,13 @@ binary_pipeline_scalar( binary_collision_model_t * RESTRICT cm, if( !nl ) continue; /* Nothing to do */ rl = UINT_MAX / (unsigned)nl; np = nk*nl; - nc = (int)( 0.5 + sample*(double)(nk>nl ? nk : nl) ); + + switch(strategy) { + case per_particle: nc = (int)( 0.5 + sample*(double)(nk>nl ? nk : nl) ) ; break ; + case mass_action: nc = (int)( 0.5 + sample*(double)np ) ; break ; + default: nc = 0; + } + } /* Determine the collision rate to probability normalization: @@ -113,13 +128,13 @@ binary_pipeline_scalar( binary_collision_model_t * RESTRICT cm, /* Pick a pair of computational particles uniformly at random from all pairs of particles in the voxel. Note that the while test virtually always fails (this manner of - splitting up the nk, nl guarantees a uniform prob + splitting up the nk, nl guarantees a uniform prob of getting k on 0:nk-1 and l on 0:nl-1 and uses the preferred high order randgen bits). */ - + do { k = (int)(uirand(rng)/rk); } while( k==nk ); k += k0; do { l = (int)(uirand(rng)/rl); } while( l==nl ); l += l0; - + /* Compute the probability that a physical particle in the species whose candidate computational particle has the least weight (and comoving with this computational particle) will @@ -145,7 +160,7 @@ binary_pipeline_scalar( binary_collision_model_t * RESTRICT cm, The other particle should be updated with probability of w_min / w_max, such that, on average, detailed balance is preserved. */ - + w_min = (wk>wl) ? wl : wk; type = 1; if( wl==w_min ) type++; if( w_max==w_min || w_max*frand_c0(rng)i - ((particle_mover_t*)b)->i ; +} + +void +chemical_pipeline_scalar( chemical_collision_model_t * RESTRICT cm, + int pipeline_rank, + int n_pipeline ) +{ + if( pipeline_rank==n_pipeline ) return; /* No host straggler cleanup */ + + chemical_rate_constant_func_t rate_constant = cm->rate_constant; + chemical_collision_func_t collision = cm->collision; + + /**/ void * RESTRICT params = cm->params; + /**/ species_t ** RESTRICT reactants = cm->reactants; + /**/ species_t ** RESTRICT products = cm->products; + const int * RESTRICT consumable = cm->consumable; + /**/ rng_t * RESTRICT rng = cm->rp->rng[ pipeline_rank ]; + const grid_t * RESTRICT g = cm->reactants[0]->g; + + const int N = cm->n_reactants; + const int M = cm->n_products; + const double sample = cm->sample; + const float dtinterval_dV = ( g->dt * (float)cm->interval ) / g->dV; + + species_t * sp; + float pr_norm, pr_coll, w_prod, w_min, pr_max = 0; + int n, m, k, v, v1, nc, type, n_large_pr = 0, n_tested = 0, unaccumulated = 0; + unsigned long long np; + + // Reactant variables + // i : Particle index for reactant i + // i0 : Start of particles for reactant i in this voxel + // ni : Number of particles of reactant i in this voxel + // ri : Scale from random integer to particle index + // im : Index of next available reactant mover + // nm : Number of available reactant movers + // wi : Statistical weight of reactant particle i + // mass : Species mass of reactant i + // mpart : Temporary variable for particle mass (weight * species mass) + // mtot : Temporary variable for sum(mpart) + // com : center of mass coordinates + // n_consumed : Number of reactant i consumed + // reactant_particles : Pointers to the reactant particles. + int i[N], i0[N], ni[N], ri[N], im[N], nm[N], n_modified[N]; + float wi[N], mass[N], com[3], mpart, mtot; + particle_t *reactant_particles[N]; + + // Product variables + // j : Index of the next available particle for product j + // j0 : Index of first available particle for product j + // nj : How much free space is available for product j + // product_particles : Pointers to the product particles + int j[M], j0[M], nj[M], n_produced=0; + particle_t *product_particles[M]; + + // Determine the free space available for this pipeline to inject product + // particles. + FOR_PRODUCTS { + DISTRIBUTE(sp->max_np-sp->np, 16, pipeline_rank, n_pipeline, j0[m], nj[m]); + j0[m] += sp->np; + j[m] = j0[m]; + } + + // Determine which particle movers belong to this pipeline. We are using + // particle movers to temporarially hold deleted particles since they are + // already allocated and sized appropriately. + FOR_REACTANTS { + DISTRIBUTE(sp->max_nm, 8, pipeline_rank, n_pipeline, im[n], nm[n]); + n_modified[n] = 0; + mass[n] = sp->m; + } + + // Stripe the (mostly non-ghost) voxels over threads for load balance. + v = VOXEL( 0,0,0, g->nx,g->ny,g->nz ) + pipeline_rank; + v1 = VOXEL( g->nx,g->ny,g->nz, g->nx,g->ny,g->nz ) + 1; + for( ; vpartition[v ]; + ni[n] = sp->partition[v+1] - i0[n]; + if( !ni[n] ) goto next_voxel; + ri[n] = UINT_MAX / (unsigned)ni[n]; + + // Check to see if this species is repeated in the reactants, and + // if so, count the number of repetitions and update arrays. + for( k=1, m=n+1 ; mp + i[n]; + + // Compute prod(w) and find min(w) + wi[n] = reactant_particles[n]->w; + w_prod *= wi[n]; + if( n == 0 || wi[n] < w_min ) w_min = wi[n]; + } + + // There is a possibility that if a particle is picked for collisions + // that it's weight could be <= 0. In this case, we should really just + // pick another set of particles in order to maintain the correct + // probabilities, however this could lead to an infinite loop if _all_ + // particles are consumed. The correct way to do this would be to + // recompute probabilities after each collision, taking into account + // particle consumption, but this is a lot of overhead. For now, we + // will just skip this collision and take the loss in accuracy. + if( w_min <= 0 ) continue; + + + /* Compute the probability that a physical particle in the + species whose candidate computational particle has the least + weight (and comoving with this computational particle) will + collide off a beam of physical particles of the density and + momentum of the computational particle in the other species. + If this probability is bigger than one, make a note for + diagnostic use. */ + + pr_coll = w_prod / w_min * pr_norm * + rate_constant( params, reactants, reactant_particles ); + + if( pr_coll > 1 ) n_large_pr++; + if( pr_coll > pr_max ) pr_max = pr_coll; + + /* Yes, >= so that 0 rate constants guarantee no collision and + yes, _c0, so that 1 probabilities guarantee a collision */ + if( frand_c0(rng)>=pr_coll ) continue; /* Didn't collide */ + + // Reactants had a collision. Determine which computational particles + // should be updated by the collision process. We should always update the + // particle of least weight. The other particles should be updated with + // probability of w_min / w[i] such that, on average, detailed balance is + // preserved. + + type = 0; + com[0] = com[1] = com[2] = 0; + FOR_REACTANTS{ + if( wi[n] == w_min || wi[n]*frand_c0(rng)dx; + com[1] += mpart * reactant_particles[n]->dy; + com[2] += mpart * reactant_particles[n]->dz; + mtot += mpart; + } + com[0] /= mtot; + com[1] /= mtot; + com[2] /= mtot; + + // Now create the product particles. They should all be created with + // a statistical weight of w_min and located at the center of mass. + // NOTE: If cells are concave, then COM may be outside of the cell! + + FOR_PRODUCTS { + product_particles[m] = sp->p + j[m]; + product_particles[m]->i = v; + product_particles[m]->dx = com[0]; + product_particles[m]->dy = com[1]; + product_particles[m]->dz = com[2]; + product_particles[m]->w = w_min; + ++j[m]; + } + + // Do the collision. The collision function should set u{x,y,z} and d{x,y,z} + // for the product particles, and may also modify u{x,y,z} and w for the + // reactant particles. + collision( params, reactants, reactant_particles, + products, product_particles, rng, type ); + + // Check for reactant consumption which is flagged by a change in particle + // statistical weight post-collison. These modified particles (along with + // product particles) need to undergo rhob accumulation, however doing this + // from within the pipeline is significant effort since pipelines would want + // to write to the same location. To do this, we would either need to use + // blocking constructs or have each pipeline allocate a field array and then + // reduce across pipelines. It is much easier to offload this work to the host, + // although this raises an issue if more particles are modified than exist + // particle movers. This should not be a common event, and can be avoided by + // either allocating more movers or decreasing the collisional timestep. + FOR_REACTANTS if( consumable[n] ) { + reactant_particles[n]->w -= w_min; + if( nm[n] <= n_modified[n] ) { + ++unaccumulated; + } + else { + sp->pm[im[n]].i = i[n]; + sp->pm[im[n]].dispx = w_min; + ++im[n]; + ++n_modified[n]; + } + } + + // Check that we have enough product storage left to continue to process + // collisions. The use of goto is not ideal, but cleaner than 3x break. + FOR_PRODUCTS if( j[m] - j0[m] >= nj[m] ) goto final; + + } + + next_voxel:; + + } + + final: + + // Save pipeline-produced information to the model structure. + FOR_REACTANTS { + cm->n_modified[N*pipeline_rank + n] = n_modified[n]; + cm->reactant_movers[N*pipeline_rank + n] = sp->pm + im[n] - n_modified[n]; + } + + cm->n_produced[pipeline_rank] = j[0]-j0[0]; + FOR_PRODUCTS cm->product_particles[M*pipeline_rank + m] = sp->p + j0[m]; + + cm->unaccumulated[pipeline_rank] = unaccumulated; + cm->n_large_pr[pipeline_rank] = n_large_pr; + cm->pr_max[pipeline_rank] = pr_max; + cm->n_tested[pipeline_rank] = n_tested; +} + +void +apply_chemical_collision_model_pipeline( chemical_collision_model_t * cm ) +{ + species_t *sp; + double pr_max = 0; + int n, m, p, np, nm, n_large_pr = 0, n_tested = 0, unaccumulated = 0; + particle_t *p0; + particle_mover_t *m0; + float dw, w0; + /**/ species_t ** RESTRICT reactants = cm->reactants; + /**/ species_t ** RESTRICT products = cm->products; + const grid_t * RESTRICT g = reactants[0]->g; + /**/ field_t * RESTRICT f = cm->fa->f; + const int N = cm->n_reactants; + const int M = cm->n_products; + + if( cm->interval<1 || (g->step % cm->interval) ) return; + FOR_REACTANTS if( sp->last_sorted != g->step ) sort_p( sp ); + EXEC_PIPELINES( chemical, cm, 0 ); + WAIT_PIPELINES(); + + // Combine elements from all pipelines and do some reductions. + for( p=0; pn_large_pr[p]; + n_tested += cm->n_tested[p]; + unaccumulated += cm->unaccumulated[p]; + if( cm->pr_max[p] > pr_max ) pr_max = cm->pr_max[p]; + + // Eliminate holes in product particle arrays and accumulate bound charge. + FOR_PRODUCTS { + np = cm->n_produced[p]; + p0 = cm->product_particles[p*M+m]; + if( p0 != sp->p + sp->np ) MOVE(sp->p + sp->np, p0, np); + for( ; np ; --np, ++(sp->np) ) accumulate_rhob(f, sp->p+sp->np, g, -sp->q); + } + + // Eliminate holes in reactant particle movers. + FOR_REACTANTS { + nm = cm->n_modified[p*N+n]; + m0 = cm->reactant_movers[p*N+n]; + if( m0 != sp->pm + sp->nm ) MOVE(sp->pm + sp->nm, m0, nm); + sp->nm += nm; + } + + } + + // Since we striped voxels for better load balancing, the reactant particle + // movers are out of order. We need to sort them before we can backfill. This + // could be optimized significantly, but we expect len(pm) << len(p) so + // complex sorting functions are more effort than they are worth here. + FOR_REACTANTS { + nm = 0; // Number of particles deleted from this species + qsort(sp->pm, sp->nm, sizeof(particle_mover_t), compare_particle_mover); + for( ; sp->nm ; --(sp->nm) ){ + p0 = sp->p + sp->pm[sp->nm-1].i; + + // Fractional accumulation of rhob. Charge = dw q + dw = sp->pm[sp->nm-1].dispx; + w0 = p0->w; + p0->w = dw; + accumulate_rhob(f, p0, g, sp->q); + p0->w = w0; + + // Check for complete consumption and backfill. + if(w0 <= 0) { + *p0 = sp->p[ --(sp->np) ]; + ++nm; + } + + } + if( nm ) sort_p(sp); + } + + if( n_large_pr ) { + WARNING(( "%i/%i computational particle groupings encountered a large " + "collision probability in collision model \"%s\". The " + "collision rate for such groupings will be lower than it should " + "be physically. Consider lowering the collision operator " + "interval, increasing the sampling or reducing the timestep.", + n_large_pr, n_tested, cm->name )); + if( cm->dynamic_sampling ) cm->sample *= 2*pr_max; + } + + if( unaccumulated ) + WARNING(( "%i reactant particles were consumed without proper rhob " + "accumulation in collision model \"%s\". Local charge " + "conservation may be broken.", unaccumulated, cm->name )); + +}