Line data    Source code

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