Line data    Source code

       1             : /**********************************************************************/
       2             : /*                     DO NOT MODIFY THIS HEADER                      */
       3             : /* MAGPIE - Mesoscale Atomistic Glue Program for Integrated Execution */
       4             : /*                                                                    */
       5             : /*            Copyright 2017 Battelle Energy Alliance, LLC            */
       6             : /*                        ALL RIGHTS RESERVED                         */
       7             : /**********************************************************************/
       8             : 
       9             : #include "DiscretePKAPDF.h"
      10             : #include "MooseError.h"
      11             : #include "MooseRandom.h"
      12             : #include "MagpieUtils.h"
      13             : 
      14             : #include "libmesh/vector_value.h"
      15             : 
      16           0 : DiscretePKAPDF::DiscretePKAPDF()
      17             :   : DiscretePKAPDFBase(),
      18           0 :     _probability_density_function(MultiIndex<Real>({1})),
      19           0 :     _marginal_cdf_mu(MultiIndex<Real>({1})),
      20           0 :     _marginal_cdf_phi(MultiIndex<Real>({1})),
      21           0 :     _marginal_cdf_energy(MultiIndex<Real>({1})),
      22           0 :     _marginal_cdf_zaid(MultiIndex<Real>({1}))
      23             : {
      24           0 : }
      25             : 
      26           1 : DiscretePKAPDF::DiscretePKAPDF(const std::vector<unsigned int> & ZAID,
      27             :                                const std::vector<Real> & energies,
      28           1 :                                const MultiIndex<Real> & probabilities)
      29             :   : DiscretePKAPDFBase(ZAID, energies),
      30           1 :     _na(probabilities.size()[2]),
      31           1 :     _dphi(2.0 * libMesh::pi / _na),
      32           1 :     _np(probabilities.size()[3]),
      33           1 :     _dmu(2.0 / _np),
      34           1 :     _probability_density_function(probabilities),
      35           1 :     _marginal_cdf_mu(probabilities),
      36           1 :     _marginal_cdf_phi(probabilities),
      37           1 :     _marginal_cdf_energy(probabilities),
      38           2 :     _marginal_cdf_zaid(probabilities)
      39             : {
      40             :   // check the size of _probabilities
      41           1 :   if (probabilities.dim() != 4)
      42           0 :     mooseError("probabilities MultiIndex object has wrong dimensions.");
      43             : 
      44           1 :   MultiIndex<Real>::size_type shape = probabilities.size();
      45           1 :   if (shape[0] != _nZA || shape[1] != _ng || shape[2] != _na || shape[3] != _np)
      46           0 :     mooseError("Size of probabilities is inconsistent with random variable input.");
      47             : 
      48           1 :   precomputeCDF(probabilities);
      49           1 :   computeMagnitude(probabilities);
      50           1 : }
      51             : 
      52             : void
      53           1 : DiscretePKAPDF::precomputeCDF(MultiIndex<Real> probabilities)
      54             : {
      55             :   /**
      56             :    * Compute the weighted pdf: pdf value * width of the bin
      57             :    */
      58           1 :   MultiIndex<Real>::size_type shape(4);
      59           1 :   MultiIndex<Real>::size_type index(4);
      60          48 :   for (auto it : probabilities)
      61             :   {
      62          48 :     index = it.first;
      63          48 :     Real wt = _dmu * _dphi * (_energies[index[1]] - _energies[index[1] + 1]);
      64          48 :     it.second *= wt;
      65             :   }
      66             : 
      67             :   /**
      68             :    *  precompute _marginal_cdf_mu
      69             :    */
      70             :   // Step 1: Marginalize the pdf
      71           1 :   shape.assign(1, _nZA);
      72           1 :   _marginal_cdf_zaid = MultiIndex<Real>(shape);
      73           3 :   for (unsigned int jZA = 0; jZA < _nZA; ++jZA)
      74             :   {
      75             :     Real integral_value = 0.0;
      76           6 :     for (unsigned int jE = 0; jE < _ng; ++jE)
      77          20 :       for (unsigned int jP = 0; jP < _na; ++jP)
      78          64 :         for (unsigned int jM = 0; jM < _np; ++jM)
      79          48 :           integral_value += probabilities({jZA, jE, jP, jM});
      80             : 
      81           2 :     _marginal_cdf_zaid({jZA}) = integral_value;
      82             :   }
      83             : 
      84             :   // Step 2: Compute cdf
      85           2 :   for (auto zaid : _marginal_cdf_zaid)
      86             :   {
      87           2 :     index = zaid.first;
      88           2 :     index[0] -= 1;
      89           2 :     if (zaid.first[0] > 0)
      90           1 :       zaid.second += _marginal_cdf_zaid(index);
      91             :   }
      92             : 
      93             :   // Step 3: Renormalize to ensure that cdf[-1] == 1
      94           2 :   for (auto zaid : _marginal_cdf_zaid)
      95             :   {
      96           2 :     index = zaid.first;
      97           2 :     index[0] = _nZA - 1;
      98           2 :     zaid.second /= _marginal_cdf_zaid(index);
      99             :   }
     100             : 
     101             :   /**
     102             :    *  precompute _marginal_cdf_energy
     103             :    */
     104             :   // Step 1: Marginalize the pdf
     105           1 :   shape.resize(2);
     106           1 :   index.resize(4);
     107           1 :   shape[0] = _nZA;
     108           1 :   shape[1] = _ng;
     109           1 :   _marginal_cdf_energy = MultiIndex<Real>(shape);
     110           3 :   for (unsigned int jZA = 0; jZA < _nZA; ++jZA)
     111           6 :     for (unsigned int jE = 0; jE < _ng; ++jE)
     112             :     {
     113             :       Real integral_value = 0.0;
     114          20 :       for (unsigned int jP = 0; jP < _na; ++jP)
     115          64 :         for (unsigned int jM = 0; jM < _np; ++jM)
     116          48 :           integral_value += probabilities({jZA, jE, jP, jM});
     117           4 :       _marginal_cdf_energy({jZA, jE}) = integral_value;
     118             :     }
     119             : 
     120             :   // Step 2: Compute cdf
     121           4 :   for (auto energy : _marginal_cdf_energy)
     122             :   {
     123           4 :     index = energy.first;
     124           4 :     index[1] -= 1;
     125           4 :     if (energy.first[1] > 0)
     126           2 :       energy.second += _marginal_cdf_energy(index);
     127             :   }
     128             : 
     129             :   // Step 3: Renormalize to ensure that cdf[-1] == 1
     130           4 :   for (auto energy : _marginal_cdf_energy)
     131             :   {
     132           4 :     index = energy.first;
     133           4 :     index[1] = _ng - 1;
     134           4 :     energy.second /= _marginal_cdf_energy(index);
     135             :   }
     136             : 
     137             :   /**
     138             :    *  precompute _marginal_cdf_phi
     139             :    */
     140             :   // Step 1: Marginalize the pdf
     141           1 :   shape.resize(3);
     142           1 :   index.resize(4);
     143           1 :   shape[0] = _nZA;
     144           1 :   shape[1] = _ng;
     145           1 :   shape[2] = _na;
     146           1 :   _marginal_cdf_phi = MultiIndex<Real>(shape);
     147           3 :   for (unsigned int jZA = 0; jZA < _nZA; ++jZA)
     148           6 :     for (unsigned int jE = 0; jE < _ng; ++jE)
     149          20 :       for (unsigned int jP = 0; jP < _na; ++jP)
     150             :       {
     151             :         Real integral_value = 0.0;
     152          64 :         for (unsigned int jM = 0; jM < _np; ++jM)
     153          48 :           integral_value += probabilities({jZA, jE, jP, jM});
     154          16 :         _marginal_cdf_phi({jZA, jE, jP}) = integral_value;
     155             :       }
     156             : 
     157             :   // Step 2: Compute cdf
     158          16 :   for (auto phi : _marginal_cdf_phi)
     159             :   {
     160          16 :     index = phi.first;
     161          16 :     index[2] -= 1;
     162          16 :     if (phi.first[2] > 0)
     163          12 :       phi.second += _marginal_cdf_phi(index);
     164             :   }
     165             : 
     166             :   // Step 3: Renormalize to ensure that cdf[-1] == 1
     167          16 :   for (auto phi : _marginal_cdf_phi)
     168             :   {
     169          16 :     index = phi.first;
     170          16 :     index[2] = _na - 1;
     171          16 :     phi.second /= _marginal_cdf_phi(index);
     172             :   }
     173             : 
     174             :   /**
     175             :    *  precompute _marginal_cdf_mu
     176             :    */
     177             :   // Step 1: No need to marginalize pdf for mu
     178           1 :   _marginal_cdf_mu = probabilities;
     179             : 
     180             :   // Step 2: Compute cdf
     181          48 :   for (auto mu : _marginal_cdf_mu)
     182             :   {
     183          48 :     index = mu.first;
     184          48 :     index[3] -= 1;
     185          48 :     if (mu.first[3] > 0)
     186          32 :       mu.second += _marginal_cdf_mu(index);
     187             :   }
     188             : 
     189             :   // Step 3: Renormalize to ensure that cdf[-1] == 1
     190          48 :   for (auto mu : _marginal_cdf_mu)
     191             :   {
     192          48 :     index = mu.first;
     193          48 :     index[3] = _np - 1;
     194          48 :     mu.second /= _marginal_cdf_mu(index);
     195             :   }
     196           1 : }
     197             : 
     198             : void
     199      200000 : DiscretePKAPDF::drawSample(std::vector<MyTRIM_NS::IonBase> & initial_state) const
     200             : {
     201             :   // resize initial_state
     202      200000 :   initial_state.resize(1);
     203             : 
     204             :   /**
     205             :    * Same the discrete pdfs by first sampling from the "most marginal"
     206             :    * _marginal_cdf_zaid. Then get the conditional _marginal_cdf_energy(zaid)
     207             :    * and sample the energy bin from it, then continue on sampling phi...
     208             :    */
     209             :   MultiIndex<Real>::size_type sampled_indices;
     210      200000 :   sampled_indices.push_back(sampleHelper(_marginal_cdf_zaid));
     211      200000 :   sampled_indices.push_back(sampleHelper(_marginal_cdf_energy, sampled_indices));
     212      200000 :   sampled_indices.push_back(sampleHelper(_marginal_cdf_phi, sampled_indices));
     213      400000 :   sampled_indices.push_back(sampleHelper(_marginal_cdf_mu, sampled_indices));
     214             : 
     215             :   // first we need to compute Z and m from ZAID
     216      200000 :   initial_state[0]._Z = MagpieUtils::getZFromZAID(_zaids[sampled_indices[0]]);
     217      200000 :   initial_state[0]._m = MagpieUtils::getAFromZAID(_zaids[sampled_indices[0]]);
     218             : 
     219             :   // the real random variables also need to be resampled uniformly
     220             :   // within bin index[j]
     221      200000 :   initial_state[0]._E =
     222      200000 :       (_energies[sampled_indices[1]] - _energies[sampled_indices[1] + 1]) * MooseRandom::rand() +
     223      200000 :       _energies[sampled_indices[1] + 1];
     224             : 
     225             :   // NOTE: we need to sample the direction in this class because the direction is anisotropic
     226      200000 :   Real sampled_phi = _dphi * MooseRandom::rand() + _dphi * sampled_indices[2];
     227      200000 :   Real sampled_mu = _dmu * MooseRandom::rand() + _dmu * sampled_indices[3] - 1.0;
     228             :   // NOTE: this is consistent with Rattlesnake's definition of omega
     229      200000 :   initial_state[0]._dir(0) = sampled_mu;
     230      200000 :   initial_state[0]._dir(1) = std::sqrt(1.0 - sampled_mu * sampled_mu) * std::cos(sampled_phi);
     231      200000 :   initial_state[0]._dir(2) = std::sqrt(1.0 - sampled_mu * sampled_mu) * std::sin(sampled_phi);
     232      200000 : }
     233             : 
     234             : void
     235           1 : DiscretePKAPDF::computeMagnitude(MultiIndex<Real> probabilities)
     236             : {
     237           1 :   _magnitude = 0.0;
     238          48 :   for (auto it : probabilities)
     239             :   {
     240          48 :     MultiIndex<Real>::size_type index = it.first;
     241          48 :     Real delE = _energies[index[1]] - _energies[index[1] + 1];
     242          48 :     _magnitude += delE * _dphi * _dmu * it.second;
     243             :   }
     244           1 : }
     245             : 
     246             : std::ostream &
     247           0 : operator<<(std::ostream & out, const DiscretePKAPDF & pdf)
     248             : {
     249           0 :   out << "Magnitude " << pdf._magnitude << "\n";
     250           0 :   out << "ZAIDs ";
     251           0 :   for (auto & z : pdf._zaids)
     252           0 :     out << z << " ";
     253           0 :   out << "\nEnergy group boundaries ";
     254           0 :   for (auto & e : pdf._energies)
     255           0 :     out << e << " ";
     256           0 :   out << "\nPolar cosine boundaries spacing " << pdf._dmu;
     257           0 :   out << "\nAzimuthal angle spacing " << pdf._dphi;
     258           0 :   out << "\nPDF values:\n";
     259           0 :   for (unsigned int jZA = 0; jZA < pdf._nZA; ++jZA)
     260           0 :     for (unsigned int jE = 0; jE < pdf._ng; ++jE)
     261           0 :       for (unsigned int jP = 0; jP < pdf._na; ++jP)
     262           0 :         for (unsigned int jM = 0; jM < pdf._np; ++jM)
     263             :         {
     264           0 :           Real mu = (Real(jM) + 0.5) * pdf._dmu - 1.0;
     265           0 :           Real phi = (Real(jP) + 0.5) * pdf._dphi;
     266             :           RealVectorValue omega(
     267           0 :               mu, std::cos(phi) * std::sqrt(1 - mu * mu), std::sin(phi) * std::sqrt(1 - mu * mu));
     268           0 :           out << "Zaid: " << pdf._zaids[jZA] << " Energies: [" << pdf._energies[jE + 1] << ", "
     269           0 :               << pdf._energies[jE] << "] "
     270           0 :               << " phi [" << jP * pdf._dphi << ", " << (jP + 1) * pdf._dphi << "] "
     271           0 :               << " mu: [" << jM * pdf._dmu - 1.0 << ", " << (jM + 1) * pdf._dmu - 1.0 << "] "
     272           0 :               << " omega " << omega
     273           0 :               << " pdf: " << pdf._probability_density_function({jZA, jE, jP, jM}) << "\n";
     274             :         }
     275           0 :   return out;
     276             : }