1 /* ---------------------------------------------------------------------- 2 SPARTA - Stochastic PArallel Rarefied-gas Time-accurate Analyzer 3 http://sparta.sandia.gov 4 Steve Plimpton, sjplimp@sandia.gov, Michael Gallis, magalli@sandia.gov 5 Sandia National Laboratories 6 7 Copyright (2014) Sandia Corporation. Under the terms of Contract 8 DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains 9 certain rights in this software. This software is distributed under 10 the GNU General Public License. 11 12 See the README file in the top-level SPARTA directory. 13 ------------------------------------------------------------------------- */ 14 15 #ifdef REACT_CLASS 16 17 ReactStyle(tce/kk,ReactTCEKokkos) 18 19 #else 20 21 #ifndef SPARTA_REACT_TCE_KOKKOS_H 22 #define SPARTA_REACT_TCE_KOKKOS_H 23 24 #include "react_bird_kokkos.h" 25 #include "kokkos_type.h" 26 27 namespace SPARTA_NS { 28 29 class ReactTCEKokkos : public ReactBirdKokkos { 30 public: 31 ReactTCEKokkos(class SPARTA *, int, char **); 32 ReactTCEKokkos(class SPARTA* sparta) : ReactBirdKokkos(sparta) {}; 33 void init(); 34 int attempt(Particle::OnePart *, Particle::OnePart *, 35 double, double, double, double &, int &) { return 0; } 36 37 /* ---------------------------------------------------------------------- */ 38 39 enum{DISSOCIATION,EXCHANGE,IONIZATION,RECOMBINATION}; // other files 40 41 KOKKOS_INLINE_FUNCTION 42 int attempt_kk(Particle::OnePart *ip, Particle::OnePart *jp, 43 double pre_etrans, double pre_erot, double pre_evib, 44 double &post_etotal, int &kspecies, 45 int &recomb_species, double &recomb_density, const t_species_1d_const &d_species) const 46 { 47 OneReactionKokkos *r; 48 49 const int isp = ip->ispecies; 50 const int jsp = jp->ispecies; 51 52 const double pre_ave_rotdof = (d_species[isp].rotdof + d_species[jsp].rotdof)/2.0; 53 54 const int n = d_reactions(isp,jsp).n; 55 if (n == 0) return 0; 56 auto& d_list = d_reactions(isp,jsp).d_list; 57 58 // probablity to compare to reaction probability 59 60 double react_prob = 0.0; 61 rand_type rand_gen = rand_pool.get_state(); 62 const double random_prob = rand_gen.drand(); 63 rand_pool.free_state(rand_gen); 64 65 // loop over possible reactions for these 2 species 66 67 for (int i = 0; i < n; i++) { 68 r = &d_rlist[d_list[i]]; 69 70 // ignore energetically impossible reactions 71 72 const double pre_etotal = pre_etrans + pre_erot + pre_evib; 73 74 double ecc = pre_etrans; 75 if (pre_ave_rotdof > 0.1) ecc += pre_erot*r->d_coeff[0]/pre_ave_rotdof; 76 77 const double e_excess = ecc - r->d_coeff[1]; 78 if (e_excess <= 0.0) continue; 79 80 // compute probability of reaction 81 82 switch (r->type) { 83 case DISSOCIATION: 84 case IONIZATION: 85 case EXCHANGE: 86 { 87 react_prob += r->d_coeff[2] * 88 pow(ecc-r->d_coeff[1],r->d_coeff[3]) * 89 pow(1.0-r->d_coeff[1]/ecc,r->d_coeff[5]); 90 break; 91 } 92 93 case RECOMBINATION: 94 { 95 // skip if no 3rd particle chosen by Collide::collisions() 96 // this includes effect of boost factor to skip recomb reactions 97 // check if this recomb reaction is the same one 98 // that the 3rd particle species maps to, else skip it 99 // this effectively skips all recombinations reactions 100 // if selected a 3rd particle species that matches none of them 101 // scale probability by boost factor to restore correct stats 102 103 if (recomb_species < 0) continue; 104 auto& d_sp2recomb = d_reactions(isp,jsp).d_sp2recomb; 105 if (d_sp2recomb[recomb_species] != d_list[i]) continue; 106 107 react_prob += recomb_boost * recomb_density * r->d_coeff[2] * 108 pow(ecc,r->d_coeff[3]) * 109 pow(1.0-r->d_coeff[1]/ecc,r->d_coeff[5]); 110 break; 111 } 112 113 default: 114 //error->one(FLERR,"Unknown outcome in reaction"); 115 //d_error_flag() = 1; 116 Kokkos::abort("ReactTCEKokkos: Unknown outcome in reaction\n"); 117 break; 118 } 119 120 // test against random number to see if this reaction occurs 121 // if it does, reset species of I,J and optional K to product species 122 // J particle can be destroyed in recombination reaction, set species = -1 123 // K particle can be created in a dissociation or ionization reaction, 124 // set its kspecies, parent will create it 125 // important NOTE: 126 // does not matter what order I,J reactants are in compared 127 // to order the reactants are listed in the reaction file 128 // for two reasons: 129 // a) list of N possible reactions above includes all reactions 130 // that I,J species are in, regardless of order 131 // b) properties of pre-reaction state, stored in precoln, 132 // as computed by setup_collision(), 133 // and used by perform_collision() after reaction has taken place, 134 // only store combined properties of I,J, 135 // nothing that is I-specific or J-specific 136 137 if (react_prob > random_prob) { 138 Kokkos::atomic_fetch_add(&d_tally_reactions[d_list[i]],1); 139 ip->ispecies = r->d_products[0]; 140 141 // Previous statment did not destroy the 2nd species (B) if 142 // recombination was specified as A+B->AB+M (which has nproductus=2) 143 // but only for the A+B->AB specication form (which has nproductus=1) 144 145 switch (r->type) { 146 case DISSOCIATION: 147 case IONIZATION: 148 case EXCHANGE: 149 { 150 jp->ispecies = r->d_products[1]; 151 break; 152 } 153 case RECOMBINATION: 154 { 155 // always destroy 2nd reactant species 156 157 jp->ispecies = -1; 158 break; 159 } 160 } 161 162 if (r->nproduct > 2) kspecies = r->d_products[2]; 163 else kspecies = -1; 164 165 post_etotal = pre_etotal + r->d_coeff[4]; 166 167 return 1; 168 } 169 } 170 171 return 0; 172 } 173 174 /* ---------------------------------------------------------------------- */ 175 176 protected: 177 DAT::tdual_int_scalar k_error_flag; 178 DAT::t_int_scalar d_error_flag; 179 HAT::t_int_scalar h_error_flag; 180 }; 181 182 } 183 184 #endif 185 #endif 186 187 /* ERROR/WARNING messages: 188 189 */ 190