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 "BoundingBoxParticleInitializer.h"
      17             : #include "MooseRandom.h"
      18             : #include "ElementSampler.h"
      19             : #include "Distribution.h"
      20             : 
      21             : registerMooseObject("SalamanderApp", BoundingBoxParticleInitializer);
      22             : 
      23             : InputParameters
      24         275 : BoundingBoxParticleInitializer::validParams()
      25             : {
      26         275 :   auto params = PerElementParticleInitializer::validParams();
      27         275 :   params.addClassDescription(
      28             :       "This initializer performs the same task as [PerElementParticleInitializer.md], however, "
      29             :       "only particles that exist within the bounding box are created.");
      30         550 :   params.addRequiredParam<Point>("bottom_left", "The bottom left corner of the bounding box");
      31         550 :   params.addRequiredParam<Point>("top_right", "The top right corner of the bounding box");
      32             : 
      33         275 :   return params;
      34           0 : }
      35             : 
      36         139 : BoundingBoxParticleInitializer::BoundingBoxParticleInitializer(const InputParameters & parameters)
      37             :   : PerElementParticleInitializer(parameters),
      38         139 :     _bottom_left(getParam<Point>("bottom_left")),
      39         278 :     _top_right(getParam<Point>("top_right")),
      40        1112 :     _planes({{-1, 0, 0, _bottom_left(0)},
      41         139 :              {1, 0, 0, -_top_right(0)},
      42         139 :              {0, -1, 0, _bottom_left(1)},
      43         139 :              {0, 1, 0, -_top_right(1)},
      44         139 :              {0, 0, -1, _bottom_left(2)},
      45         139 :              {0, 0, 1, -_top_right(2)}})
      46             : {
      47         464 :   for (const auto i : make_range(_mesh_dimension))
      48         326 :     if (_top_right(i) <= _bottom_left(i))
      49           1 :       paramError("top_right",
      50           2 :                  "Component " + std::to_string(uint(i)) + " of 'top_right' is <= component " +
      51           1 :                      std::to_string(uint(i)) + " of 'bottom_left'");
      52             : 
      53         138 :   if (_mesh_dimension != Moose::dim)
      54         196 :     mooseWarning(std::to_string(uint(Moose::dim)) +
      55             :                  " components are required for libMesh::Point input.\n"
      56          66 :                  "However, your simulation is only " +
      57         130 :                  std::to_string(uint(_mesh_dimension)) +
      58             :                  " dimensional.\n"
      59          66 :                  "The extra component" +
      60         153 :                  std::string(_mesh_dimension == uint(2) ? "s" : "") +
      61             :                  " of the libMesh::Point input will be ignored.\n");
      62         414 : }
      63             : 
      64             : std::vector<InitialParticleData>
      65          91 : BoundingBoxParticleInitializer::getParticleData() const
      66             : {
      67             : 
      68             :   // first collect the elements that we will place particles in
      69             :   std::vector<const Elem *> valid_elems;
      70             : 
      71       28072 :   for (const auto elem : *_fe_problem.mesh().getActiveLocalElementRange())
      72             :   {
      73             :     unsigned int intersection_count = 0;
      74             :     unsigned int nodes_within = 0;
      75             :     bool elem_added = false;
      76             :     // if one of the nodes is in our box or on the edge
      77             :     // of the box then we need to add this to the elements we want
      78             :     // to place particles in
      79      140091 :     for (const auto & node : elem->node_ref_range())
      80             :     {
      81             :       unsigned int dim_valid = 0;
      82             :       // first we can check to see if one node is within the bounding box
      83             :       // this is sufficient for 1D
      84      452706 :       for (const auto i : make_range(_mesh_dimension))
      85             :       {
      86      337571 :         if (node(i) >= _bottom_left(i) && node(i) <= _top_right(i))
      87      119577 :           dim_valid++;
      88             :       }
      89             : 
      90      115135 :       if (dim_valid == _mesh_dimension)
      91             :       {
      92        3025 :         valid_elems.push_back(elem);
      93             :         elem_added = true;
      94             :         break;
      95             :       }
      96             : 
      97             :       // if we are in a higher dimension then we need to check for
      98             :       // elements which intersect the bounding box but do not have
      99             :       // any nodes within the bounding box
     100      112110 :       if (_mesh_dimension > 1)
     101             :       {
     102             :         // checking to see which side of each plane each node is on
     103      769466 :         for (const auto i : make_range(2 * _mesh_dimension))
     104             :         {
     105             :           Real plane_value = 0;
     106     2600688 :           for (const auto j : make_range(_mesh_dimension))
     107     1943228 :             plane_value += node(j) * _planes[i][j];
     108             : 
     109      657460 :           plane_value += _planes[i][3];
     110             : 
     111      657460 :           if (plane_value >= 0)
     112      217890 :             nodes_within++;
     113             :         }
     114      112006 :         if (nodes_within > 1)
     115      103546 :           intersection_count++;
     116             :       }
     117             :     }
     118        3025 :     if (elem_added)
     119        3025 :       continue;
     120             :     // if more than two planes are intersecting the
     121             :     // element then we need to put some particles within it as well
     122       24956 :     if (intersection_count > 1)
     123       24908 :       valid_elems.push_back(elem);
     124             :   }
     125             : 
     126             :   // if there are no elements for this processor: do nothing
     127          91 :   if (valid_elems.size() == 0)
     128           6 :     return {};
     129             : 
     130             :   std::vector<InitialParticleData> data;
     131             :   MooseRandom generator;
     132             :   // this objects allows us to uniformly sample space in elements
     133          85 :   SALAMANDER::ElementSampler sampler = SALAMANDER::ElementSampler(_fe_problem, _seed, generator);
     134       28017 :   for (const auto elem : valid_elems)
     135             :   {
     136             :     // first sample points in the element like we would if this were a uniform initialization
     137             :     // across the whole domain
     138       27933 :     const auto & physical_points = sampler.sampleElement(elem, _particles_per_element);
     139       27932 :     Real weight = _number_density * elem->volume() / (_particles_per_element);
     140      586572 :     for (const auto point : physical_points)
     141             :     {
     142             :       // we need to check to make sure that every point we have
     143             :       // created is actually in the bounding box as well
     144             :       unsigned int dim_valid = 0;
     145     2185280 :       for (const auto i : make_range(_mesh_dimension))
     146             :       {
     147     1626640 :         if (point(i) >= _bottom_left(i) && point(i) <= _top_right(i))
     148      583548 :           dim_valid++;
     149             :       }
     150             :       // if the point is in the box then we can add a particle for this point
     151      558640 :       if (dim_valid == _mesh_dimension)
     152             :       {
     153       22558 :         auto & particle = data.emplace_back();
     154       22558 :         particle.elem = elem;
     155       22558 :         particle.weight = weight;
     156       22558 :         particle.species = _species;
     157       22558 :         particle.mass = _mass;
     158       22558 :         particle.charge = _charge;
     159       22558 :         particle.position = point;
     160       90232 :         for (const auto i : make_range(uint(3)))
     161       67674 :           particle.velocity(i) = _velocity_distributions[0]->quantile(generator.rand());
     162             :       }
     163             :     }
     164             :   }
     165             :   return data;
     166         168 : }