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