model/detail/EquilibrateDensitiesSolver.hh

0001 //------------------------------- -*- C++ -*- -------------------------------//
0002 // Copyright Celeritas contributors: see top-level COPYRIGHT file for details
0003 // SPDX-License-Identifier: (Apache-2.0 OR MIT)
0004 //---------------------------------------------------------------------------//
0005 //! \file celeritas/mucf/model/detail/EquilibrateDensitiesSolver.hh
0006 //---------------------------------------------------------------------------//
0007 #pragma once
0008
0009 #include "corecel/cont/EnumArray.hh"
0010 #include "celeritas/Constants.hh"
0011 #include "celeritas/mucf/Types.hh"
0012
0013 namespace celeritas
0014 {
0015 namespace detail
0016 {
0017 //---------------------------------------------------------------------------//
0018 /*!
0019  * Calculate dt mixture densities after reaching thermodynamic
0020  * equilibrium based isotopic fraction, density, and material temperature.
0021  *
0022  * Based on the theory from https://www.osti.gov/biblio/6205719.
0023  *
0024  * The equilibrated densities are needed to correctly calculate the cycle time
0025  * of dd, dt, and tt fusion cycles.
0026  *
0027  * \sa MucfMaterialInserter
0028  */
0029 class EquilibrateDensitiesSolver
0030 {
0031   public:
0032     //! Enum for safely accessing hydrogen isoprotologue molecules
0033     enum class MucfIsoprotologueMolecule
0034     {
0035         protium_protium,
0036         protium_deuterium,
0037         protium_tritium,
0038         deuterium_deuterium,
0039         deuterium_tritium,
0040         tritium_tritium,
0041         size_
0042     };
0043
0044     //!@{
0045     //! \name Type aliases
0046     using LhdArray = EnumArray<MucfIsotope, real_type>;
0047     using EquilibriumArray = EnumArray<MucfIsoprotologueMolecule, real_type>;
0048     //!@}
0049
0050     // Construct with material information
0051     EquilibrateDensitiesSolver(LhdArray const& lhd_densities);
0052
0053     // Calculate equilibrated isoprotologue densities
0054     EquilibriumArray operator()(real_type const temperature);
0055
0056   private:
0057     //// DATA ////
0058
0059     LhdArray lhd_densities_;  //!< Number densities in units of LHD
0060     real_type total_density_;  //!< Total LHD density
0061     real_type inv_tot_density_;  //!< Inverse of total LHD density
0062     static constexpr Constant r_gas_ = constants::k_boltzmann
0063                                        * constants::na_avogadro;
0064
0065     //// HELPER FUNCTIONS ////
0066
0067     // {
0068     // Convergence limit parameters
0069     // Acceptance error between current and previous equilibrium iteration
0070     static real_type constexpr convergence_err() { return 1e-6; }
0071
0072     // Maximum number of iterations to reach convergence
0073     static size_type constexpr max_iterations() { return 1000; }
0074     // }
0075
0076     // Calculate equilibrium constant: \f$ H_2 + D_2 \rightleftharpoons 2HD \f$
0077     real_type calc_hd_equilibrium_constant(real_type temperature);
0078
0079     // Calculate equilibrium constant: \f$ H_2 + T_2 \rightleftharpoons 2HT \f$
0080     real_type calc_ht_equilibrium_constant(real_type temperature);
0081
0082     // Calculate equilibrium constant: \f$ D_2 + T_2 \rightleftharpoons 2DT \f$
0083     real_type calc_dt_equilibrium_constant(real_type temperature);
0084
0085     // Equilibrate a pair of isotopes and write the new density in the array
0086     void equilibrate_pair(MucfIsoprotologueMolecule molecule_aa,
0087                           MucfIsoprotologueMolecule molecule_bb,
0088                           MucfIsoprotologueMolecule molecule_ab,
0089                           real_type eq_constant_ab,
0090                           EquilibriumArray& input);
0091 };
0092
0093 //---------------------------------------------------------------------------//
0094 }  // namespace detail
0095 }  // namespace celeritas