Line data    Source code

       1             : //* This file is part of SALAMANDER: Software for Advanced Large-scale Analysis of MAgnetic confinement for Numerical Design, Engineering & Research,
       2             : //* A multiphysics application for modeling plasma facing components
       3             : //* https://github.com/idaholab/salamander
       4             : //* https://mooseframework.inl.gov/salamander
       5             : //*
       6             : //* SALAMANDER is powered by the MOOSE Framework
       7             : //* https://www.mooseframework.inl.gov
       8             : //*
       9             : //* Licensed under LGPL 2.1, please see LICENSE for details
      10             : //* https://www.gnu.org/licenses/lgpl-2.1.html
      11             : //*
      12             : //* Copyright 2025, Battelle Energy Alliance, LLC
      13             : //* ALL RIGHTS RESERVED
      14             : //*
      15             : 
      16             : #include "UniformGridParticleInitializer.h"
      17             : #include "MooseRandom.h"
      18             : #include <limits>
      19             : 
      20             : registerMooseObject("SalamanderApp", UniformGridParticleInitializer);
      21             : 
      22             : InputParameters
      23          86 : UniformGridParticleInitializer::validParams()
      24             : {
      25          86 :   auto params = ParticleInitializerBase::validParams();
      26          86 :   params.addClassDescription("Particle initializer that places particles along a line with "
      27             :                              "approximate uniform spacing between particles");
      28         172 :   params.addRangeCheckedParam<unsigned int>(
      29             :       "total_particles",
      30             :       "total_particles != 0",
      31             :       "The number of computational particles that should be placed along the line");
      32         172 :   params.addRangeCheckedParam<Real>("number_density",
      33             :                                     "number_density > 0.0",
      34             :                                     "The number density to represent with computational particles");
      35             : 
      36          86 :   return params;
      37           0 : }
      38             : 
      39          44 : UniformGridParticleInitializer::UniformGridParticleInitializer(const InputParameters & parameters)
      40             :   : ParticleInitializerBase(parameters),
      41          44 :     _number_density(getParam<Real>("number_density")),
      42         132 :     _total_particles(getParam<unsigned int>("total_particles"))
      43             : {
      44          44 :   if (_mesh_dimension != 1)
      45           2 :     mooseError("The simulation must be in 1D in order to use the UniformGridParticleInitializer");
      46          42 : }
      47             : 
      48             : std::vector<InitialParticleData>
      49          28 : UniformGridParticleInitializer::getParticleData() const
      50             : {
      51             :   Real local_xmin = std::numeric_limits<float>::max();
      52          28 :   Real global_xmax = std::numeric_limits<float>::lowest();
      53             :   Real local_volume = 0;
      54         860 :   for (const auto elem : *_fe_problem.mesh().getActiveLocalElementRange())
      55             :   {
      56         832 :     local_volume += elem->volume();
      57        2496 :     for (const auto & node : elem->node_ref_range())
      58             :     {
      59        1664 :       if (node(0) < local_xmin)
      60             :         local_xmin = node(0);
      61             : 
      62        1664 :       if (node(0) > global_xmax)
      63         860 :         global_xmax = node(0);
      64             :     }
      65             :   }
      66          28 :   Real global_volume = local_volume;
      67          28 :   Real global_xmin = local_xmin;
      68             : 
      69          28 :   comm().sum(global_volume);
      70          28 :   comm().min(global_xmin);
      71          28 :   comm().max(global_xmax);
      72             : 
      73          28 :   double fraction = local_volume / global_volume;
      74          28 :   double min_frac = fraction;
      75          28 :   double max_frac = fraction;
      76             : 
      77          28 :   comm().min(min_frac);
      78          28 :   comm().max(max_frac);
      79             : 
      80             :   // Doing some rounding here to help reduce the cases where the total number of requested particles
      81             :   // does not match the total number to be created.
      82             :   // Without this rounding, even in cases where the total number of particles requested is divided
      83             :   // evenly by the number of processors, the number of particles created does not match the
      84             :   // requested number.
      85          28 :   uint local_particle_count = std::round(double(_total_particles) * local_volume / global_volume);
      86          28 :   uint global_particle_count = local_particle_count;
      87             : 
      88          28 :   comm().sum(global_particle_count);
      89             : 
      90          28 :   if (global_particle_count != _total_particles)
      91             :   {
      92           0 :     std::ostringstream oss;
      93           0 :     oss << _total_particles << " particles across " << comm().size() << " processes were requested."
      94             :         << std::endl;
      95           0 :     oss << "But " << global_particle_count << " will be created because of the mesh partition.";
      96           0 :     mooseWarning(oss.str());
      97           0 :   }
      98             : 
      99          28 :   std::vector<InitialParticleData> data = std::vector<InitialParticleData>(local_particle_count);
     100             : 
     101          28 :   Real dx = (global_xmax - global_xmin) / (_total_particles);
     102             : 
     103             :   MooseRandom generator;
     104          28 :   generator.seed(_seed);
     105             : 
     106          28 :   uint particle_count = 0;
     107          28 :   Point curr_point = Point((particle_count + 0.5) * dx + local_xmin);
     108         856 :   for (const auto elem : *_fe_problem.mesh().getActiveLocalElementRange())
     109             :   {
     110             :     // the particles that are currently in the element
     111             :     auto particle_idxs = std::vector<uint>();
     112        1664 :     while (elem->contains_point(curr_point) && particle_count < local_particle_count)
     113             :     {
     114         832 :       particle_idxs.push_back(particle_count);
     115         832 :       data[particle_count].elem = elem;
     116         832 :       data[particle_count].species = _species;
     117         832 :       data[particle_count].mass = _mass;
     118         832 :       data[particle_count].charge = _charge;
     119         832 :       data[particle_count].position = curr_point;
     120         832 :       data[particle_count].velocity = Point();
     121        3328 :       for (const auto j : make_range(uint(3)))
     122        2496 :         data[particle_count].velocity(j) = _velocity_distributions[j]->quantile(generator.rand());
     123         832 :       particle_count++;
     124         832 :       curr_point = Point((particle_count + 0.5) * dx + local_xmin);
     125             :     }
     126             : 
     127        1664 :     for (const auto idx : particle_idxs)
     128             :     {
     129         832 :       data[idx].weight = _number_density * elem->volume() / (particle_idxs.size());
     130             :     }
     131             : 
     132             :     particle_idxs.clear();
     133             : 
     134         832 :     if (particle_count == _total_particles)
     135             :       break;
     136             :   }
     137             : 
     138          28 :   return data;
     139           0 : }