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 "ThreadedRadialGreensConvolutionLoop.h"
      10             : #include "Function.h"
      11             : 
      12             : // Remove after idaholab/moose#26102
      13             : #include "MagpieNanoflann.h"
      14             : 
      15          24 : ThreadedRadialGreensConvolutionLoop::ThreadedRadialGreensConvolutionLoop(
      16          24 :     RadialGreensConvolution & green)
      17          24 :   : _green(green), _function_name(_green.parameters().get<FunctionName>("function"))
      18             : {
      19          24 : }
      20             : 
      21             : // Splitting Constructor
      22           6 : ThreadedRadialGreensConvolutionLoop::ThreadedRadialGreensConvolutionLoop(
      23           6 :     const ThreadedRadialGreensConvolutionLoop & x, Threads::split /*split*/)
      24           6 :   : _green(x._green), _function_name(x._function_name)
      25             : {
      26           6 : }
      27             : 
      28             : void
      29          30 : ThreadedRadialGreensConvolutionLoop::operator()(const QPDataRange & qpdata_range)
      30             : {
      31             :   // fetch data from parent
      32          30 :   const auto r_cut = _green._r_cut;
      33             :   const auto & qp_data = _green._qp_data;
      34             :   const auto & correction_integral = _green._correction_integral;
      35             :   const auto & kd_tree = _green._kd_tree;
      36          30 :   const auto dim = _green._dim;
      37          30 :   const auto t = _green._t;
      38             : 
      39             :   // fetch thread copy of the function for safe evaluation
      40             :   ParallelUniqueId puid;
      41          30 :   const Function & function = _green._fe_problem.getFunction(_function_name, puid.id);
      42             : 
      43             :   // radial bin size
      44          30 :   const Real dr = r_cut / correction_integral.size();
      45             : 
      46             :   // integral of the convolution (used if normalization is requested)
      47          30 :   _convolution_integral = 0.0;
      48             : 
      49             :   // tree search data structures
      50             :   std::vector<nanoflann::ResultItem<std::size_t, Real>> ret_matches;
      51             :   nanoflann::SearchParameters search_params;
      52             : 
      53             :   // result map entry
      54          30 :   const auto end_it = _green._convolution.end();
      55             :   auto it = end_it;
      56             : 
      57             :   // iterate over qp range
      58      142950 :   for (auto && local_qp : qpdata_range)
      59             :   {
      60             :     // Look up result map iterator only if we enter a new element. this saves a bunch
      61             :     // of map lookups because same element entries are consecutive in the qp_data vector.
      62      142920 :     if (it == end_it || it->first != local_qp._elem_id)
      63       34260 :       it = _green._convolution.find(local_qp._elem_id);
      64             : 
      65             :     // initialize result entry
      66             :     mooseAssert(it != end_it, "Current element id not found in result set.");
      67      142920 :     auto & sum = it->second[local_qp._qp];
      68      142920 :     sum = 0.0;
      69             : 
      70             :     // if the variable is periodic we need to perform extra searches translated onto
      71             :     // the periodic neighbors
      72      142920 :     std::list<Point> cell_vector = {Point()};
      73      440640 :     for (unsigned int j = 0; j < dim; ++j)
      74      297720 :       if (_green._periodic[j])
      75             :       {
      76             :         std::list<Point> new_cell_vector;
      77             : 
      78      641280 :         for (const auto & cell : cell_vector)
      79             :         {
      80      343560 :           if (local_qp._q_point(j) + _green._periodic_vector[j](j) - r_cut <
      81             :               _green._periodic_max[j])
      82       50814 :             new_cell_vector.push_back(cell + _green._periodic_vector[j]);
      83             : 
      84      343560 :           if (local_qp._q_point(j) - _green._periodic_vector[j](j) + r_cut >
      85             :               _green._periodic_min[j])
      86       50814 :             new_cell_vector.push_back(cell - _green._periodic_vector[j]);
      87             :         }
      88             : 
      89      297720 :         cell_vector.splice(cell_vector.end(), new_cell_vector);
      90             :       }
      91             : 
      92             :     // perform radius search and aggregate data considering potential periodicity
      93             :     Point center;
      94      387468 :     for (const auto & cell : cell_vector)
      95             :     {
      96             :       ret_matches.clear();
      97      244548 :       center = local_qp._q_point + cell;
      98             :       std::size_t n_result =
      99      244548 :           kd_tree->radiusSearch(&(center(0)), r_cut * r_cut, ret_matches, search_params);
     100    54405582 :       for (std::size_t j = 0; j < n_result; ++j)
     101             :       {
     102    54161034 :         const auto & other_qp = qp_data[ret_matches[j].first];
     103    54161034 :         const Real r = std::sqrt(ret_matches[j].second);
     104             : 
     105             :         // R is the equivalent sphere radius for the quadrature point. The spherical integral
     106             :         // integral_0^R 1/(4pi*r^2) *4pi*r^2 = R
     107             :         // So R is the integral over the geometric attenuation at the center quadrature point
     108    54161034 :         switch (dim)
     109             :         {
     110        1824 :           case 1:
     111             :           {
     112             :             // correction integral is the integral over the geometric attenuation
     113             :             // times the Green's function over a 2D disc inscribed into the r_cut
     114             :             // sphere and perpendicular to the mesh.
     115        1824 :             Real add = correction_integral[std::floor(r / dr)] * other_qp._volume;
     116             : 
     117        1824 :             if (r == 0)
     118             :             {
     119         120 :               const Real R = 0.5 * other_qp._volume;
     120         120 :               add += function.value(t, Point(0.0, 0.0, 0.0)) *
     121             :                      (
     122             :                          // add the center sphere attenuation integral
     123         120 :                          R +
     124             :                          // add the section missing or overlapping between correction_integral and
     125             :                          // center sphere
     126         120 :                          _green.attenuationIntegral(R, _green._zero_dh, 0.0, 1) * other_qp._volume);
     127             :             }
     128        1824 :             sum += add * other_qp._value;
     129        1824 :             break;
     130             :           }
     131             : 
     132    44229735 :           case 2:
     133             :           {
     134             :             // correction integral is the integral over the geometric attenuation
     135             :             // times the Green's function over a 1D line segment inscribed into the r_cut
     136             :             // sphere and perpendicular to the mesh.
     137    44229735 :             Real add = correction_integral[std::floor(r / dr)] * other_qp._volume;
     138    44229735 :             if (r == 0)
     139             :             {
     140      130800 :               const Real R = std::sqrt(other_qp._volume / (2.0 * libMesh::pi));
     141      130800 :               add += function.value(t, Point(0.0, 0.0, 0.0)) *
     142             :                      (
     143             :                          // add the center sphere attenuation integral
     144      130800 :                          R +
     145             :                          // add the section missing or overlapping between correction_integral and
     146             :                          // center sphere
     147      130800 :                          _green.attenuationIntegral(R, _green._zero_dh, 0.0, 2) * other_qp._volume);
     148             :             }
     149    44229735 :             sum += add * other_qp._volume * other_qp._value;
     150    44229735 :             break;
     151             :           }
     152             : 
     153     9929475 :           case 3:
     154             :           {
     155     9929475 :             const Real G = function.value(t, Point(r, 0.0, 0.0));
     156     9929475 :             if (r == 0)
     157             :             {
     158             :               // R is the integral over the geometric attenuation in a sphere around the origin
     159       12000 :               const Real R = std::cbrt(3.0 / 4.0 * other_qp._volume / libMesh::pi);
     160       12000 :               sum += G * R * other_qp._value;
     161             :             }
     162             :             else
     163     9917475 :               sum += G * 0.25 / (libMesh::pi * r * r) * other_qp._volume * other_qp._value;
     164             :           }
     165             :         }
     166             :       }
     167             :     }
     168             : 
     169             :     // integrate the convolution result
     170      142920 :     _convolution_integral += sum;
     171             :   }
     172          30 : }