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 "BorisStepper.h"
      17             : 
      18             : registerMooseObject("SalamanderApp", BorisStepper);
      19             : 
      20             : InputParameters
      21         170 : BorisStepper::validParams()
      22             : {
      23         170 :   auto params = ParticleStepperBase::validParams();
      24         170 :   params.addClassDescription(
      25             :       "Electromagnetic particle stepper method implementing the Boris Algorithm.");
      26         340 :   params.addRequiredParam<std::vector<VariableName>>(
      27             :       "efield_components",
      28             :       "A list of 3 variables which represent the 3 components of the electric field");
      29         340 :   params.addRequiredParam<std::vector<VariableName>>(
      30             :       "bfield_components",
      31             :       "A list of 3 variables which represent the 3 components of the electric field");
      32         170 :   return params;
      33           0 : }
      34             : 
      35          86 : BorisStepper::BorisStepper(const InputParameters & parameters)
      36             :   : ParticleStepperBase(parameters),
      37          86 :     _efield_vars(getParam<std::vector<VariableName>>("efield_components")),
      38         344 :     _bfield_vars(getParam<std::vector<VariableName>>("bfield_components"))
      39             : {
      40          86 :   if (_efield_vars.size() != 3)
      41           1 :     paramError("efield_components",
      42             :                "You must provide 3 components representing the electric field");
      43             : 
      44          85 :   if (_bfield_vars.size() != 3)
      45           1 :     paramError("bfield_components",
      46             :                "You must provide 3 components representing the magnetic field");
      47             : 
      48         336 :   for (int i = 0; i < 3; ++i)
      49             :   {
      50         252 :     _efield_samplers.push_back(SALAMANDER::VariableSampler(_fe_problem, _efield_vars[i], _tid));
      51         504 :     _bfield_samplers.push_back(SALAMANDER::VariableSampler(_fe_problem, _bfield_vars[i], _tid));
      52             :   }
      53          84 : }
      54             : 
      55             : void
      56         320 : BorisStepper::setupStep(Ray & ray, Point & v, const Real q_m_ratio, const Real distance) const
      57             : {
      58             : 
      59         320 :   Real dt = _dt;
      60             :   // if it is the particles first step we need to take a half step
      61             :   // otherwise we take a full step
      62         320 :   if (distance == 0)
      63          32 :     dt /= 2;
      64             :   // let's sample the fields at the ray location so we can update the velocity properly
      65         320 :   Point E = sampleField(_efield_samplers, ray);
      66         320 :   Point B = sampleField(_bfield_samplers, ray);
      67             :   // calculate v^-
      68         320 :   v = linearImpulse(v, E, q_m_ratio, dt / 2);
      69             :   // calculate v^+
      70         320 :   v = magneticImpulse(v, B, q_m_ratio, dt);
      71             :   // calculate the final velocity and setup the particle for the step
      72         320 :   v = linearImpulse(v, E, q_m_ratio, dt / 2);
      73         320 :   setMaxDistanceAndDirection(ray, v, _dt);
      74         320 : }
      75             : 
      76             : Point
      77         320 : BorisStepper::magneticImpulse(const Point & v,
      78             :                               const Point & B,
      79             :                               const Real q_m_ratio,
      80             :                               const Real dt) const
      81             : {
      82             :   auto l = q_m_ratio * B * dt / 2;
      83         320 :   auto s = 2 * l / (1 + l * l);
      84             : 
      85         320 :   auto v_prime = v + v.cross(l);
      86             : 
      87         320 :   return v + v_prime.cross(s);
      88             : }