Commit 39ae30c2 authored by Steffen Schotthöfer's avatar Steffen Schotthöfer
Browse files

Merge branch 'develop' of https://git.scc.kit.edu/rtsn/rtsn into develop


Former-commit-id: da2db94b
parents 9f26f65e 0882c593
......@@ -60,7 +60,7 @@ class CSDSolverTrafoFP : public SNSolver
void FVMUpdate( unsigned idx_energy ) override final;
void FluxUpdate() override final;
void IterPreprocessing( unsigned idx_pseudotime ) override final;
void virtual IterPostprocessing() override final;
void virtual IterPostprocessing( unsigned idx_pseudotime ) override final;
void SolverPreprocessing() override final;
};
......
......@@ -72,7 +72,7 @@ class CSDSolverTrafoFP2D : public SNSolver
void FVMUpdate( unsigned idx_energy ) override final;
void FluxUpdate() override final;
void IterPreprocessing( unsigned idx_pseudotime ) override final;
void virtual IterPostprocessing() override final;
void virtual IterPostprocessing( unsigned idx_pseudotime ) override final;
void SolverPreprocessing() override final;
};
......
......@@ -89,7 +89,7 @@ class CSDSolverTrafoFPSH2D : public SNSolver
void FVMUpdate( unsigned idx_energy ) override final;
void FluxUpdate() override final;
void IterPreprocessing( unsigned idx_pseudotime ) override final;
void virtual IterPostprocessing() override final;
void virtual IterPostprocessing( unsigned idx_pseudotime ) override final;
void SolverPreprocessing() override final;
};
......
......@@ -55,7 +55,7 @@ class MNSolver : public SolverBase
void FVMUpdate( unsigned idx_energy ) override;
void FluxUpdate() override;
void IterPreprocessing( unsigned idx_pseudotime ) override;
void IterPostprocessing();
void IterPostprocessing( unsigned idx_pseudotime );
/*! @brief : Construct flux by computing the Moment of the sum of FVM discretization at the interface of cell
* @param : idx_cell = current cell id
* @returns : sum over all neighbors of flux for all moments at interface of idx_cell, idx_neighbor */
......
......@@ -54,7 +54,7 @@ class PNSolver : public SolverBase
void FVMUpdate( unsigned idx_energy ) override;
void FluxUpdate() override;
void IterPreprocessing( unsigned idx_pseudotime ) override;
void IterPostprocessing();
void IterPostprocessing( unsigned idx_pseudotime );
// Helper
void ComputeRadFlux();
......
......@@ -30,7 +30,7 @@ class SNSolver : public SolverBase
void virtual FVMUpdate( unsigned idx_energy ) override;
void virtual FluxUpdate() override;
void virtual IterPreprocessing( unsigned idx_pseudotime ) override;
void virtual IterPostprocessing() override;
void virtual IterPostprocessing( unsigned idx_pseudotime ) override;
// Helper
void ComputeRadFlux();
......
......@@ -90,7 +90,7 @@ class SolverBase
@param idx_iter : current (peudo) time iteration */
virtual void IterPreprocessing( unsigned idx_iter ) = 0;
/*! @brief Performs postprocessing for the current solver iteration */
virtual void IterPostprocessing() = 0;
virtual void IterPostprocessing( unsigned idx_pseudotime ) = 0;
/*! @brief Constructs the flux update for the current iteration and stores it in psiNew*/
virtual void FluxUpdate() = 0;
/*! @brief Computes the finite Volume update step for the current iteration
......
......@@ -11,13 +11,14 @@ LOG_DIR = ../result/logs
MESH_FILE = meshes/1DMesh.su2
PROBLEM = WATERPHANTOM
SOLVER = CSD_SN_FOKKERPLANCK_TRAFO_SOLVER
CONTINUOUS_SLOWING_DOWN = YES
HYDROGEN_FILE = ENDL_H.txt
OXYGEN_FILE = ENDL_O.txt
KERNEL = ISOTROPIC_1D
CFL_NUMBER = 0.008
TIME_FINAL = 1.0
CLEAN_FLUX_MATRICES = NO
MAX_ENERGY_CSD = 5
VOLUME_OUTPUT = (MINIMAL, MEDICAL)
BC_DIRICHLET = ( dirichlet )
BC_NEUMANN = ( wall_low, wall_up )
......
......@@ -17,7 +17,7 @@ CSDSolverTrafoFP::CSDSolverTrafoFP( Config* settings ) : SNSolver( settings ) {
// Set angle and energies
_energies = Vector( _nEnergies, 0.0 ); // equidistant
_energyMin = 1e-4 * 0.511;
_energyMax = 10e0;
_energyMax = _settings->GetMaxEnergyCSD();
// write equidistant energy grid (false) or refined grid (true)
GenerateEnergyGrid( false );
......@@ -183,10 +183,24 @@ void CSDSolverTrafoFP::FVMUpdate( unsigned /*idx_energy*/ ) {
}
}
void CSDSolverTrafoFP::IterPostprocessing() {
void CSDSolverTrafoFP::IterPostprocessing( unsigned idx_pseudotime ) {
// --- Update Solution ---
_sol = _solNew;
unsigned n = idx_pseudotime;
for( unsigned j = 0; j < _nCells; ++j ) {
_fluxNew[j] = dot( _sol[j], _weights );
if( n > 0 ) {
_dose[j] += 0.5 * _dE * ( _fluxNew[j] * _s[_nEnergies - n - 1] + _flux[j] * _s[_nEnergies - n] ) /
_density[j]; // update dose with trapezoidal rule
}
else {
_dose[j] += _dE * _fluxNew[j] * _s[_nEnergies - n - 1] / _density[j];
}
_solverOutput[j] = _fluxNew[j];
_flux[j] = _fluxNew[j];
}
// --- Compute Flux for solution and Screen Output ---
ComputeRadFlux();
}
......@@ -340,16 +354,7 @@ void CSDSolverTrafoFP::WriteVolumeOutput( unsigned idx_pseudoTime ) {
case MEDICAL:
// Compute Dose
for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
if( idx_cell > 0 ) {
_outputFields[idx_group][0][idx_cell] +=
0.5 * _dE *
( _fluxNew[idx_cell] * _s[_nEnergies - idx_pseudoTime - 1] + _flux[idx_cell] * _s[_nEnergies - idx_pseudoTime] ) /
_density[idx_cell]; // update dose with trapezoidal rule
}
else {
_outputFields[idx_group][0][idx_cell] +=
_dE * _fluxNew[idx_cell] * _s[_nEnergies - idx_pseudoTime - 1] / _density[idx_cell];
}
_outputFields[idx_group][0][idx_cell] = _dose[idx_cell];
}
// Compute normalized dose
_outputFields[idx_group][1] = _outputFields[idx_group][0];
......
......@@ -247,16 +247,7 @@ void CSDSolverTrafoFP2D::WriteVolumeOutput( unsigned idx_pseudoTime ) {
case MEDICAL:
// Compute Dose
for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
if( idx_cell > 0 ) {
_outputFields[idx_group][0][idx_cell] +=
0.5 * _dE *
( _fluxNew[idx_cell] * _s[_nEnergies - idx_pseudoTime - 1] + _flux[idx_cell] * _s[_nEnergies - idx_pseudoTime] ) /
_density[idx_cell]; // update dose with trapezoidal rule
}
else {
_outputFields[idx_group][0][idx_cell] +=
_dE * _fluxNew[idx_cell] * _s[_nEnergies - idx_pseudoTime - 1] / _density[idx_cell];
}
_outputFields[idx_group][0][idx_cell] += _dose[idx_cell];
}
// Compute normalized dose
_outputFields[idx_group][1] = _outputFields[idx_group][0];
......@@ -358,13 +349,27 @@ void CSDSolverTrafoFP2D::IterPreprocessing( unsigned idx_pseudotime ) {
}
}
void CSDSolverTrafoFP2D::IterPostprocessing() {
void CSDSolverTrafoFP2D::IterPostprocessing( unsigned idx_pseudotime ) {
unsigned n = idx_pseudotime;
// --- Update Solution ---
for( unsigned j = 0; j < _nCells; ++j ) {
if( _boundaryCells[j] == BOUNDARY_TYPE::DIRICHLET ) continue;
_sol[j] = _solNew[j];
}
for( unsigned j = 0; j < _nCells; ++j ) {
_fluxNew[j] = dot( _sol[j], _weights );
if( n > 0 ) {
_dose[j] += 0.5 * _dE * ( _fluxNew[j] * _s[_nEnergies - n - 1] + _flux[j] * _s[_nEnergies - n] ) /
_density[j]; // update dose with trapezoidal rule
}
else {
_dose[j] += _dE * _fluxNew[j] * _s[_nEnergies - n - 1] / _density[j];
}
_solverOutput[j] = _fluxNew[j];
_flux[j] = _fluxNew[j];
}
// --- Compute Flux for solution and Screen Output ---
ComputeRadFlux();
}
......
......@@ -357,16 +357,7 @@ void CSDSolverTrafoFPSH2D::WriteVolumeOutput( unsigned idx_pseudoTime ) {
case MEDICAL:
// Compute Dose
for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
if( idx_cell > 0 ) {
_outputFields[idx_group][0][idx_cell] +=
0.5 * _dE *
( _fluxNew[idx_cell] * _s[_nEnergies - idx_pseudoTime - 1] + _flux[idx_cell] * _s[_nEnergies - idx_pseudoTime] ) /
_density[idx_cell]; // update dose with trapezoidal rule
}
else {
_outputFields[idx_group][0][idx_cell] +=
_dE * _fluxNew[idx_cell] * _s[_nEnergies - idx_pseudoTime - 1] / _density[idx_cell];
}
_outputFields[idx_group][0][idx_cell] += _dose[idx_cell];
}
// Compute normalized dose
_outputFields[idx_group][1] = _outputFields[idx_group][0];
......@@ -468,7 +459,8 @@ void CSDSolverTrafoFPSH2D::IterPreprocessing( unsigned idx_pseudotime ) {
}
}
void CSDSolverTrafoFPSH2D::IterPostprocessing() {
void CSDSolverTrafoFPSH2D::IterPostprocessing( unsigned idx_pseudotime ) {
unsigned n = idx_pseudotime;
for( unsigned j = 0; j < _nCells; ++j ) {
if( _boundaryCells[j] == BOUNDARY_TYPE::DIRICHLET ) continue;
_sol[j] = _solNew[j];
......@@ -476,6 +468,19 @@ void CSDSolverTrafoFPSH2D::IterPostprocessing() {
// --- Compute Flux for solution and Screen Output ---
ComputeRadFlux();
for( unsigned j = 0; j < _nCells; ++j ) {
_fluxNew[j] = dot( _sol[j], _weights );
if( n > 0 ) {
_dose[j] += 0.5 * _dE * ( _fluxNew[j] * _s[_nEnergies - n - 1] + _flux[j] * _s[_nEnergies - n] ) /
_density[j]; // update dose with trapezoidal rule
}
else {
_dose[j] += _dE * _fluxNew[j] * _s[_nEnergies - n - 1] / _density[j];
}
_solverOutput[j] = _fluxNew[j];
_flux[j] = _fluxNew[j];
}
}
void CSDSolverTrafoFPSH2D::SolverPreprocessing() {
......
......@@ -152,7 +152,7 @@ void MNSolver::IterPreprocessing( unsigned /*idx_pseudotime*/ ) {
}
}
void MNSolver::IterPostprocessing() {
void MNSolver::IterPostprocessing( unsigned idx_pseudotime ) {
// --- Update Solution ---
_sol = _solNew;
......
......@@ -57,7 +57,7 @@ void PNSolver::IterPreprocessing( unsigned /*idx_pseudotime*/ ) {
// Nothing to preprocess for PNSolver
}
void PNSolver::IterPostprocessing() {
void PNSolver::IterPostprocessing( unsigned idx_pseudotime ) {
// --- Update Solution ---
_sol = _solNew;
......
......@@ -27,7 +27,7 @@ void SNSolver::IterPreprocessing( unsigned /*idx_pseudotime*/ ) {
// Nothing to do for SNSolver
}
void SNSolver::IterPostprocessing() {
void SNSolver::IterPostprocessing( unsigned idx_pseudotime ) {
// --- Update Solution ---
_sol = _solNew;
......
......@@ -146,7 +146,7 @@ void SolverBase::Solve() {
FVMUpdate( iter );
// --- Iter Postprocessing ---
IterPostprocessing();
IterPostprocessing( iter );
// --- Solver Output ---
WriteVolumeOutput( iter );
......
......@@ -7,15 +7,15 @@ Authors
KiT-RT Software
*****************
+--------------------+------------------------+
|Main contributors | - Jannick Wolters |
| | - Jonas Kusch |
| | - Steffen Schotthöfer |
| | - Tianbai Xiao |
| | - Pia Stammer |
+--------------------+------------------------+
|Other contributors | |
+--------------------+------------------------+
+--------------------+----------------------------------------+
|Main contributors | - Jannick Wolters |
| | - Jonas Kusch |
| | - Steffen Schotthöfer |
| | - Tianbai Xiao (tianbai.xiao@kit.edu) |
| | - Pia Stammer |
+--------------------+----------------------------------------+
|Other contributors | |
+--------------------+----------------------------------------+
*********************
KiT-RT Dokumentation
......
===========================================================
KiT-RT: a kinetic transport solver for radiation therapy
KiT-RT: A kinetic transport solver for radiation therapy
===========================================================
The KiT-RT framework is a powerful open source platform for radiation transport.
Its main focus is on radiotherapy planning in cancer treatment.
To allow problem-specific method selection, the framework provides different deterministic solver types.
This not only facilitates treatment planning, but also provides tools to investigate various research questions in the field of radiative transfer.
This goal is supported by an easily extensible code structure that allows straightforward implementation of additional methods and techniques.
The KiT-RT framework is an open source project for radiation transport written in C++ programming language.
It is interested in the evolution of many-particle systems, e.g. photons, neutrons, electrons, etc.
Based on the finite volume method (FVM), it provides an efficient tool to solve the Boltzmann and related equations in multiple dimensions.
Special attention has been paid to the application of radiation therapy and treatment planning.
The framework provides rich deterministic solver types for different domain-specific problems.
A list of current supported models and equations is as follows.
- linear Boltzmann (:math:`S_N`) equation
- spherical harmonics (:math:`P_N`) equations
- entropy-closure moment (:math:`M_N`) equations
- continuous slowing down equation
Design philosophy
------------------------
The code hierarchy is designed as intuitive and neat as possible.
It's dedicated to providing a friendly interface for educational usage in kinetic theory and rich functionality for scientific research.
Benefiting from the brilliant expressiveness and low-overhead abstraction provided by the C++ programming language, we implement the easily extensible code structure,
which allow the users to focus on physics or easily extend the codes with their own methods.
What is new?
------------------------
Finite volume method is a proven approach for simulating conservation laws.
Compared with the existing open-source softwares, e.g. OpenFOAM, SU2 and Clawpack, Kit-RT holds the novelty through the following points:
- Comprehensive support for kinetic theory and phase-space equations
- Special focus on radiation therapy
- Lightweight design to ensure the flexibility for secondary development
How to get help?
------------------------
The software is being developed by members of the group `CSMM <https://www.scc.kit.edu/en/aboutus/rg-csmm.php>`_ at the Karlsruhe Institute of Technology (KIT).
For more information, please feel free to get in touch with us (link to authors page).
--------
Contents
--------
If you are interested in using KiT-RT or are trying to figure out how to use it, please feel free to get in touch with `us <authors.html>`_.
Do open an issue or pull request if you have questions, suggestions or solutions.
Table of contents
------------------------
.. toctree::
:maxdepth: 1
......
......@@ -57,9 +57,6 @@ The resulting executable will automatically be placed in the `code/bin` folder.
Run
**********
Local
===========
Execute the compiled binary from the `bin` folder and hand over a valid *TOML*-styled config file.
Example from inside the `code` directory:
......@@ -75,56 +72,6 @@ OMP_NUM_THREADS=N mpirun -np J ./KiT-RT ../input/example.cfg
with `N` equal to the number of shared memory threads and `J` equal to the number of distrubuted memory threads.
BwUniCluster
==============
As VTK is not available on the bwUniCluster, it needs to be installed first. This just needs to be done once. Example:
.. code-block:: bash
module load devel/cmake/3.16
module load compiler/gnu/9.2
wget --no-check-certificate --quiet https://www.vtk.org/files/release/8.2/VTK-8.2.0.tar.gz
tar xzf VTK-8.2.0.tar.gz
mkdir VTK-build
cd VTK-build
cmake -DCMAKE_BUILD_TYPE=Release -DBUILD_DOCUMENTATION=OFF -DBUILD_TESTING=OFF -
DCMAKE_INSTALL_PREFIX=~/VTK-install ../VTK-8.2.0
make -j
make install
cd -
rm -r VTK-8.2.0 VTK-build
Example for build and run on bwUniCluster:
Get the code
.. code-block:: bash
git clone https://git.scc.kit.edu/rtsn/rtsn.git KiT-RT
cd KiT-RT/
git submodule init
git submodule update
Append ``HINTS VTK_INSTALL_DIR` to the ``find_package( VTK ... )`` line in the CMakeLists.txt. E.g.:
.. code-block:: bash
find_package( VTK REQUIRED COMPONENTS vtkIOGeometry vtkFiltersCore HINTS ~/VTK-install )
Compile it
.. code-block:: bash
module load devel/cmake/3.16
module load compiler/gnu/9.2
module load mpi/openmpi/4.0
cd code/build/release/
cmake -DCMAKE_BUILD_TYPE=Release ../../
make -j
---------------------------------------------------------------
Tests
......@@ -138,7 +85,3 @@ After compiling the framework as described above just run:
The ``unit_tests`` executable will also be placed in in the build folder.
......@@ -2,128 +2,113 @@
Theory
================
The kinetic theory is dedicated to describe the dynamical behavior of a many-particle system through ensemble averaging.
Particles, e.g. molecules, photons, neutrons, electrons and plasmas, travel along the trajectories and undergo occasional collisions that change their directions and energies.
Such dynamics can be formulated via operator splitting approach in the Boltzmann equation, i.e.
.. math::
:label: boltzmann
\frac{\partial \psi}{\partial t}+v \cdot \nabla_x \psi = Q(\psi)
The Boltzmann equation
---------
where :math:`\psi` is the one-particle probability distribution function, :math:`v` is velocity and :math:`Q(\psi)` is the scattering term.
In the Kit-RT solver, we consider the liner Boltzmann equation
The particle transport phenomena enjoy rich academic research value and application prospects.
A many-particle system can exhibit different behaviors at characteristic different scales.
Down to the finest scale of a many-particle system, the Newton’s second law depicts particle motions via
.. math::
:label: linearbz
\partial_{t} f(v)+v \cdot \nabla_{x} f(v)=\sigma \int k\left(v, v^{\prime}\right)\left(f\left(v^{\prime}\right)-f(v)\right) d v^{\prime}-\tau f(v)
F = m a,
where the particles don't interact with one another but scatter with the background material.
For convenience, we reformulate the particle velocity into polar coordinates.
The physical world shows a diverse set of behaviors on different
characteristic scales. Consider the molecular motion of gases as an
example. Down to the finest scale of a many-particle system, the
Newton’s second law depicts particle motions via
which leads
.. math::
\mathbf{F} = m \mathbf{a}
As a first order system it reads
.. math::
\frac{d x}{dt} = v, \ \frac{d v}{dt} = \frac{F}{m}.
\frac{d \mathbf x}{dt} = \mathbf v, \ \frac{d \mathbf v}{dt} = \frac{\mathbf F}{m}
An intuitive numerical solution algorithm is to get the numerous particles on board and track the trajectories of them.
A typical example is the molecular dynamics (MD) method.
This is not going to be efficient since there are more than :math:`2\times 10^{25}` molecules per cubic meter in normal atmosphere,
and things get extremely complicated if the N-body interactions are counted all the time.
An intuitive numerical algorithm is to get the numerous particles on
board and track the trajectories of them. A typical example is the
`Molecular Dynamics`_. This is not going to be efficient since there are
more than ``2e25`` molecules per cubic meter in normal atmosphere, and
things get even more complicated when you count on the N-body
interactions all the time. Some methods have been proposed to simplify
the computation. As an example, the `Direct simulation Monte Carlo`_
employs certain molecular models and conduct the intermolecular
collisions in a stochastic manner. It significantly reduces the
computational cost, while the trade-off is the artificial fluctuations.
Many realizations must be simulated successively to average the
solutions and reduce the errors.
Simplifications can be conducted to accelerate the numerical computation.
As an example, the Monte Carlo method employs certain particle models and conduct the interactions in a stochastic manner.
It significantly reduces the computational cost, while the trade-off is the artificial fluctuations.
Many realizations must be simulated successively to average the solutions and reduce the errors.
An alternative strategy is made from ensemble averaging, where the
coarse-grained modeling is used to provide a bottom-up view. At the mean
free path and collision time scale of molecules, particles travel freely
during most of time with mild intermolecular collisions. Such dynamics
can be described with an operator splitting approach, i.e. the kinetic
transport equation
An alternative strategy can be made from ensemble averaging, where the
coarse-grained modeling is used to provide a bottom-up view.
At the mean free path and collision time scale of particles. Such dynamics can be described with kinetic theory.
The Boltzmann equation can be formulated via an operator splitting approach, i.e.
.. math::
\frac{\partial f}{\partial t}+ \mathbf v \cdot \nabla_\mathbf x f + \mathbf a \cdot \nabla_\mathbf v f = Q(f)
\partial_{t} f(v)+v \cdot \nabla_{x} f(v)=\int_{\mathcal R^3} \int_{\mathcal S^2} k\left(v, v^{\prime}\right) \left(f\left(v^{\prime}\right)f\left(v_*^{\prime}\right)-f(v)f(v_*)\right) d\Omega d v_*,
where the left and right hand sides model particle transports and
collisions correspondingly. Different collision models can be inserted
into such equation. If the particles only collide with a background
material one obtains linear Boltzmann collision operator
where the left and right hand sides model particle transports and collisions correspondingly.
The distribution function :math:`f` is the probability of finding a particle with certain location, and :math:`\{v, v_*\}` denotes the velocities of two classes of colliding particles.
The collision kernel :math:`k` models the strength of collisions at different velocities.
Different collision models can be inserted into the Boltzmann equation.
In the KiT-RT solver, we are interested in the linear Boltzmann equation, where the particles don't interact with one another but scatter with the background material.
Therefore, the Boltzmann can be simplified as the linear equation with respect to :math:`f`
.. math::
Q(f)=\int_{\mathbb R^3} \mathcal B(\mathbf v_*, \mathbf v) \left[ f(\mathbf v_*)-f(\mathbf v)\right] d\mathbf v_*
\partial_{t} f(v)+v \cdot \nabla_{x} f(v)=\int k\left(v, v^{\prime}\right)\left(f\left(v^{\prime}\right)-f(v)\right) d v^{\prime}-\tau f(v).
where the collision kernel ``\mathcal B`` models the strength of
collisions at different velocities. If the interactions among particles
are considered, the collision operator becomes nonlinear. For example,
the two-body collision results in nonlinear Boltzmann equation
For convenience, it is often reformulated with polar coordinates :math:`\{r, \phi, \theta \}`,
.. math::
Q(f)=\int_{\mathbb R^3} \int_{\mathcal S^2} \mathcal B(\cos \beta, |\mathbf{v}-\mathbf{v_*}|) \left[ f(\mathbf v')f(\mathbf v_*')-f(\mathbf v)f(\mathbf v_*)\right] d\mathbf \Omega d\mathbf v_*
&\left[\frac{1}{v(E)} \partial_{t} +\Omega \cdot \nabla+\Sigma_t (r, E, t)\right] \psi(r, \Omega, E, t) \\
&=\int_{0}^{\infty} d E^{\prime} \int_{\mathcal R^2} d \Omega^{\prime} \Sigma_{s}\left(r, \Omega^{\prime} \bullet \Omega, E^{\prime} \rightarrow E\right) \psi\left(r, \Omega^{\prime}, E^{\prime}, t\right) + Q(r, \Omega, E, t).
The particle distribution :math:`\psi(r, \Omega, E, t)` here is often named as angular flux, :math:`\{\Sigma_s, \Sigma_t \}` are the scattering and total cross sections correspondingly, and :math:`Q` denotes a source term.
The continuous slowing down approximation
---------
Our main goal is to compute the radiation dose
For the radiation therapy, the main goal is to compute the radiation dose accurately, which is defined as
.. math::
D(x)=\frac{1}{\rho(x)}\int_0^{\infty}\int_{\mathbb{S}^2}S(E,x)\psi(E,x,\Omega)\,d\Omega dE.
The angular flux $\psi$ can be approximated by the continuous slowing down (CSD) equation, which reads
where :math:`\rho(x)` is the background material density.
The angular flux :math:`\psi` can be approximated by a further approximation equation, i.e. the continuous slowing down (CSD) equation which reads
.. math::
-\partial_E\left(S(E,x)\psi(E,x,\Omega)\right)+\Omega\cdot\nabla_x\psi(E,x,\Omega)+\Sigma_t(E,x)\psi(E,x,\Omega) = \int_{\mathbb{S}^2}\Sigma_s(E,x,\Omega\cdot\Omega')\psi(E,x,\Omega')d\Omega'.
&-\partial_E\left(S(E,x)\psi(E,x,\Omega)\right)+\Omega\cdot\nabla_x\psi(E,x,\Omega)+\Sigma_t(E,x)\psi(E,x,\Omega) \\
&= \int_{\mathbb{S}^2}\Sigma_s(E,x,\Omega\cdot\Omega')\psi(E,x,\Omega')d\Omega'.
Here $E\in\mathbb{R}_+$ is energy, $x\in D\subset \mathbb{R}^3$ is the spatial domain and $\Omega\in\mathbb{S}^2$ is the direction of travel. The stopping power $S$ is given by
Here :math:`E\in\mathbb{R}_+` is energy, :math:`x\in D\subset \mathbb{R}^3` is the spatial domain and :math:`\Omega\in\mathbb{S}^2` is the direction of travel.
The stopping power :math:`S` is given by
.. math::
S(E,x) = \int_0^{\infty} E'\int_{-1}^1\Sigma(E,E',x,\mu)d\mu dE'.
Since there are no absorption effects, the total cross section is given by
Let us assume there are no absorption effects, and thus the total cross section is given by
.. math::
\Sigma_t(E,x) = \Sigma_{s,0}(E,x)=2\pi \int_{-1}^1\Sigma_s(E,x,\mu)d\mu.
With a given background material density $\rho(x)$ now make the following assumptions
With a given :math:`\rho(x)`, we now make the following assumptions
.. math::
S(E,x) = S^{H_2O}(E)\rho(x), \\
\Sigma_t(E,x) = \Sigma_t^{H_2O}(E)\rho(x), \\
\Sigma_s(E,x,\Omega\cdot\Omega') = \Sigma_s(E,\Omega\cdot\Omega')\rho(x).
&S(E,x) = S^{H_2O}(E)\rho(x), \\
&\Sigma_t(E,x) = \Sigma_t^{H_2O}(E)\rho(x), \\
&\Sigma_s(E,x,\Omega\cdot\Omega') = \Sigma_s(E,\Omega\cdot\Omega')\rho(x).
Leaving out the superscript $H_2O$, the CSD equation simplifies to
Leaving out the superscript :math:`H_2O`, the CSD equation can be simplified as
.. math::
:label: CSD2
-\partial_E\left(\rho(x)S(E)\psi(E,x,\Omega)\right)+\Omega\cdot\nabla_x\psi(E,x,\Omega)+\rho(x)\Sigma_t(E)\psi(E,x,\Omega) = \int_{\mathbb{S}^2}\rho(x)\Sigma_s(E,\Omega\cdot\Omega')\psi(E,x,\Omega')d\Omega'.
&-\partial_E\left(\rho(x)S(E)\psi(E,x,\Omega)\right)+\Omega\cdot\nabla_x\psi(E,x,\Omega)+\rho(x)\Sigma_t(E)\psi(E,x,\Omega) \\
&= \int_{\mathbb{S}^2}\rho(x)\Sigma_s(E,\Omega\cdot\Omega')\psi(E,x,\Omega')d\Omega'.
Now, we bring this system in a form which resembles the standard Boltzmann equation. Multiplying \eqref{eq:CSD2} with $S(E)$ gives
Now, we bring this system in a form which resembles the standard Boltzmann equation.
Multiplying :eq:`CSD2` with :math:`S(E)` gives
.. math::
:label: CSD3
\begin{align}
-S(E)\partial_E\left(S(E)\rho(x)\psi(E,x,\Omega)\right)+&\Omega\cdot\nabla_x S(E)\psi(E,x,\Omega)+\Sigma_t(E)S(E)\rho(x)\psi(E,x,\Omega)\\
&= \int_{\mathbb{S}^2}\Sigma_s(E,\Omega\cdot\Omega')S(E)\rho(x)\psi(E,x,\Omega')d\Omega'.
......@@ -134,11 +119,13 @@ Then, we substitute
.. math::
\widehat{\psi}(E,x,\Omega):= S(E)\rho(x)\psi(E,x,\Omega)
into \eqref{eq:CSD3}, which yields
into :eq:`CSD3`, which yields
.. math::
:label: CSD4