From 9d1f196bbef214d980a27b464779137543f106c2 Mon Sep 17 00:00:00 2001 From: itayhorin Date: Tue, 25 Oct 2022 13:42:24 +0300 Subject: [PATCH] implemented u_com and tried making a ghost particle for each collision and failed (core dumped) --- include/openmc/particle.h | 2 +- prints | 36 +++++++++++++++++ src/particle.cpp | 16 ++++++++ src/tallies/tally_scoring.cpp | 73 +++++++++++++++++++++++------------ 4 files changed, 102 insertions(+), 25 deletions(-) create mode 100644 prints diff --git a/include/openmc/particle.h b/include/openmc/particle.h index da2f6a61fc8..71c2cef581e 100644 --- a/include/openmc/particle.h +++ b/include/openmc/particle.h @@ -54,7 +54,7 @@ class Particle : public ParticleData { //! simply as a secondary particle. //! \param src Source site data void from_source(const SourceSite* src); - + void initilze_ghost_particle(Particle& p,Direction u, double E); // Coarse-grained particle events void event_calculate_xs(); void event_advance(); diff --git a/prints b/prints new file mode 100644 index 00000000000..b72c9356daf --- /dev/null +++ b/prints @@ -0,0 +1,36 @@ + fmt::print("post energy = {0} pre energy = {1}\n",p.E(),p.E_last()); + fmt::print("event mt = {}\n",p.event_mt()); + fmt::print("post weight = {0} pre weight = {1}\n",p.wgt(),p.wgt_last()); + fmt::print("pos = {0} , {1} , {2}\n",p.r().x,p.r().y,p.r().z); + fmt::print("u = {0} , {1} , {2}\n",p.u().x,p.u().y,p.u().z); + fmt::print("u_last = {0} , {1} , {2}\n",p.u_last().x,p.u_last().y,p.u_last().z); + fmt::print("distance = {}\n",distance); + fmt::print("xs = {}\n",p.macro_xs().total); + fmt::print("number of nuclear = {}\n",n); + fmt::print("name of 1# nuclicde = {}\n",nuc->name_); + fmt::print("name of 2# nuclicde = {}\n",nuc2->name_); + fmt::print("reaction MT = {}\n",rx->mt_); + fmt::print("d type = {}\n",typeid(d.get()).name()); + fmt::print("energy distribution at 0 {}\n",d_->angle().get_energy(0)); + fmt::print("energy distribution at 1 {}\n",d_->angle().get_energy(1)); + fmt::print("post energy = {}\n",p.E); + fmt::print(p.E); + + // Neutron velocity in LAB + Direction v_n = vel * p.u_last(); + + // // Sample velocity of target nucleus + // Direction v_t {}; + // if (!p.neutron_xs(i_nuclide).use_ptable) { + // v_t = sample_target_velocity(*nuc, p.E_last(), p.u_last(), v_n, + // p.neutron_xs(i_nuclide).elastic, kT, p.current_seed()); + // } + + // Velocity of center-of-mass + Direction v_cm = (v_n ) / (awr + 1.0); //(v_n + awr * v_t) / (awr + 1.0); + + // Transform to CM frame + v_n -= v_cm; + + // Find speed of neutron in CM + vel = v_n.norm(); diff --git a/src/particle.cpp b/src/particle.cpp index 4552bf6fdff..3f174b695a8 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -60,6 +60,22 @@ double Particle::speed() const return C_LIGHT * std::sqrt(1 - inv_gamma * inv_gamma); } +void Particle::initilze_ghost_particle(Particle& p,Direction u_new, double E_new) +{ + type() = p.type(); + wgt() = p.wgt_last(); + wgt_last() = p.wgt_last(); + r() = p.r(); + r_last() = p.r(); + u() = u_new; + u_last() = u_new; + //u_last() = p.u_last(); + E() = E_new; //settings::run_CE ? p.E_last() : p.g_last(); + E_last() = E_new; + time() = p.time_last(); + time_last() = p.time_last(); +} + void Particle::create_secondary( double wgt, Direction u, double E, ParticleType type) { diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index 316c0149699..7b8f295935b 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -2367,44 +2367,69 @@ void score_collision_tally(Particle& p) double distance; distance = p.r().x*p.r().x+p.r().y*p.r().y+p.r().z*p.r().z; - //fmt::print("collison happened\n"); - if(p.E_last()>2e5) - { - //fmt::print("post energy = {0} pre energy = {1}\n",p.E(),p.E_last()); - //fmt::print("event mt = {}\n",p.event_mt()); - //fmt::print("post weight = {0} pre weight = {1}\n",p.wgt(),p.wgt_last()); - //fmt::print("pos = {0} , {1} , {2}\n",p.r().x,p.r().y,p.r().z); - //fmt::print("u = {0} , {1} , {2}\n",p.u().x,p.u().y,p.u().z); - //fmt::print("u_last = {0} , {1} , {2}\n",p.u_last().x,p.u_last().y,p.u_last().z); - //fmt::print("distance = {}\n",distance); - //fmt::print("xs = {}\n",p.macro_xs().total); - const auto& mat {model::materials[p.material()]}; int n = mat->nuclide_.size(); const auto& nuc {data::nuclides[0]}; - const auto& nuc2 {data::nuclides[1]}; - //fmt::print("number of nuclear = {}\n",n); - //fmt::print("name of 1# nuclicde = {}\n",nuc->name_); - //fmt::print("name of 2# nuclicde = {}\n",nuc2->name_); - int i_nuclide=0; + double awr = nuc->awr_; const auto& rx {nuc->reactions_[0]}; - //fmt::print("reaction MT = {}\n",rx->mt_); + auto& d = rx->products_[0].distribution_[0]; auto d_ = dynamic_cast(d.get()); - //fmt::print("d type = {}\n",typeid(d.get()).name()); - //fmt::print("energy distribution at 0 {}\n",d_->angle().get_energy(0)); - //fmt::print("energy distribution at 1 {}\n",d_->angle().get_energy(1)); + + Direction u_lab {0-p.r().x,0-p.r().y,0-p.r().z}; + u_lab = u_lab/u_lab.norm(); + double cos_lab = u_lab.dot(p.u_last()); + double theta_lab = std::acos(cos_lab); + double sin_lab = std::sin(theta_lab); + double mu_COM = (std::sqrt(awr*awr - sin_lab*sin_lab)*cos_lab - sin_lab*sin_lab)/awr; + double vel = std::sqrt(p.E_last()); + double theta_pdf = d_->angle().get_pdf_value(p.E_last(),mu_COM,p.current_seed()); + double E_ghost = p.E_last()*(1+awr*awr+2*awr*mu_COM)/(1+awr)/(1+awr); + Particle ghost_particle=Particle(); + ghost_particle.initilze_ghost_particle(p,u_lab,E_ghost); + ghost_particle.event_calculate_xs(); + + //fmt::print("collison happened\n"); + if(p.E_last()>2e5) + { + fmt::print("pos = {0} , {1} , {2}\n",p.r().x,p.r().y,p.r().z); + fmt::print("u_last = {0} , {1} , {2}\n",p.u_last().x,p.u_last().y,p.u_last().z); + fmt::print("u_lab = {0} , {1} , {2}\n",u_lab.x,u_lab.y,u_lab.z); + + + if(cos_lab > 1) + { + cos_lab = 1; + if(cos_lab>1.0001) + fmt::print("-----------------------error cos_lab = {}-----------------------\n",cos_lab); } - //fmt::print("post energy = {}\n",p.E); - //fmt::print(p.E); + fmt::print("cos_lab = {}\n",cos_lab); + fmt::print("awr = {}\n",awr); + + + fmt::print("mu_COM = {}\n",mu_COM); + + + fmt::print("theta pdf = {}\n",theta_pdf); + + + fmt::print("ghost_particle E = {}\n",ghost_particle.E()); + fmt::print("ghost_particle u = {0} , {1} , {2}\n",ghost_particle.u().x,ghost_particle.u().y,ghost_particle.u().z); + + fmt::print("ghost_particle XS = {}\n",ghost_particle.macro_xs().total); + + } + + if (p.type() == ParticleType::neutron || p.type() == ParticleType::photon) { flux = 0;//p.wgt_last() / p.macro_xs().total; if(p.event_mt() == 2) { - flux = 0.5*p.wgt()*exp(-(distance-1)*p.macro_xs().total)/(2*3.14*distance*distance); + flux = 0.5*ghost_particle.wgt()*exp(-(distance-1)*ghost_particle.macro_xs().total)/(2*3.14*distance*distance); } } + for (auto i_tally : model::active_collision_tallies) { const Tally& tally {*model::tallies[i_tally]};