A Priori Error Estimates and Computational Studies for a Fermi Pencil-Beam Equation

Abstract

We derive a priori error estimates for the standard Galerkin and streamline diffusion finite element methods for the Fermi pencil-beam equation obtained from a fully three-dimensional Fokker–Planck equation in space $\mathbf{x}=(x,y,z)$ and velocity $\tilde{\mathbf{v}}=(\mu ,\eta ,\xi )$ variables. For a constant transport cross-section, there is a closed form analytic solution available for the Fermi equation with a data as product of Dirac functions. Our objective is to study the case of nonconstant, nonincreasing transport cross-section. Therefore we start with a theoretical, that is, a priori, error analysis for a Fermi model with modified initial data in L₂. Then we construct semi-streamline-diffusion and characteristic streamline-diffusion schemes and consider an adaptive algorithm for local mesh refinements. To derive the stability estimates, for simplicity, we rely on the assumption of nonincreasing transport cross-section. Different numerical examples, in two space dimensions are justifying the theoretical results. Implementations show significant reduction of the computational error by using such adaptive procedure.

Keywords:

• Fermi and Fokker–Planck pencil-beam equations
• adaptive finite element method
• duality argument
• a priori error estimates
• efficiency
• reliability

1. Introduction

This work is a further development of studies in (Borgers and Larsen 1996; Asadzadeh 1997; Asadzadeh 2000; Asadzadeh and Sopasakis 2002; Asadzadeh and Larsen 2008) where adaptive finite element method was proposed (not performed) for a reduction of computational cost in numerical approximation for pencil-beam equations. However, focusing on theoretical convergence and stability aspects, except some special cases with limited amount of implementation in (Asadzadeh and Larsen 2008) and (Asadzadeh and Sopasakis 2002), the detailed numerical tests were postponed to future works. Here, first we construct and analyze fully discrete schemes using both standard Galerkin and flux correcting streamline diffusion (SD) finite element methods for the Fermi pencil-beam equation, with a modified L₂ data, in three dimensions. The Fermi model corresponds to considering the positive direction of x-axis as the penetration direction of the beam particles. This corresponds to a positive constant first component, μ, in velocity. The process of deriving three-dimensional Fermi (where the constant μ is 1) from the Boltzmann and Fokker–Planck equations is outlined through Equations (1(1) $$\mu \frac{\partial u}{\partial x}+\eta \frac{\partial u}{\partial y}+\xi \frac{\partial u}{\partial z}={\int }_{S^2}{\sigma }_s(\mathbf{v}·\mathbf{v}\prime )\left[u\right(\mathbf{x},\mathbf{v}\prime )-u(\mathbf{x},\mathbf{v}\left)\right]d^2\mathbf{v}\prime ,\hspace{1em}0<x<L,$$(1) .1)–(1.7). For detailed asymptotic derivations, we refer to Borgers and Larsen (1996), Larsen et al. (1985), and Pomraning (1992). In the applications, the quantity of interest is the deposit of energy from the particles into the two-dimensional transverse spatial domain ${\Omega }_{\perp }:=\left\{x_{\perp }\right|x_{\perp }:=(y,z)\}$, while they are moving with velocities $\tilde{\mathbf{v}}:=\{(1,\eta ,\xi )\}$. This corresponds to a two-dimensional model problem in ${\Omega }_{\perp }:=\left\{x_{\perp }\right|x_{\perp }:=(y,z)\}$ with ${\Omega }_{\mathbf{v}}:=\{(\eta ,\xi )\}$. We have derived our error estimates, in this geometry, while in numerical implementations, for the cost and visualization reasons, we have considered examples in lower dimensions.

More specifically, our study concerns a “pencil beam” of neutral or charged particles that are normally incident on a slab of finite thickness at the spatial origin (0, 0, 0) and in the direction of the positive x-axis. The governing equation for the pencil-beam problem is the Fermi equation which is obtained by two equivalent approaches (see Borgers and Larsen 1996): either as an asymptotic limit of the linear Boltzmann equation as the transport cross-section ${\sigma }_{tr}\to 0$ and the total cross-section ${\sigma }_t\to \infty$, or as an asymptotic limit in Taylor expansion of angular flux with respect to the velocity where the terms with derivatives of order three or higher are ignored. This procedure relies on an approach that follows the Fokker–Planck development.

The Boltzmann transport equation modeling the energy independent pencil-beam process can be written as a two-point boundary value problem, viz. (1) $$\mu \frac{\partial u}{\partial x}+\eta \frac{\partial u}{\partial y}+\xi \frac{\partial u}{\partial z}={\int }_{S^2}{\sigma }_s(\mathbf{v}·\mathbf{v}\prime )\left[u\right(\mathbf{x},\mathbf{v}\prime )-u(\mathbf{x},\mathbf{v}\left)\right]d^2\mathbf{v}\prime ,\hspace{1em}0<x<L,$$(1) where $\mathbf{x}=(x,y,z)$ and $\mathbf{v}=(\mu ,\eta ,\xi )$ are the space and velocity vectors, respectively. The model problem concerns sharply forward peaked beam of particles entering the spatial domain at x = 0: (2) $$u(0,y,z,\mathbf{v})=\delta \left(y\right)\delta \left(z\right)\frac{\delta (1-\mu )}{2\pi },\hspace{1em}0<\mu \le 1,$$(2) which are demising at x = L (we may assume, without loss of generality, that L = 1): (3) $$u(L,y,z,\mathbf{v})=0,\hspace{1em}\text{or}\hspace{1em}u(1,y,z,\mathbf{v})=0,\hspace{1em}-1\le \mu <0.$$(3)

In the realm of the Boltzmann transport equation (1)(1) $$\mu \frac{\partial u}{\partial x}+\eta \frac{\partial u}{\partial y}+\xi \frac{\partial u}{\partial z}={\int }_{S^2}{\sigma }_s(\mathbf{v}·\mathbf{v}\prime )\left[u\right(\mathbf{x},\mathbf{v}\prime )-u(\mathbf{x},\mathbf{v}\left)\right]d^2\mathbf{v}\prime ,\hspace{1em}0<x<L,$$(1) , an overview of the transport theory of charged particles can be found in (Luo and Brahme 1993). In this setting a few first coefficients in a Legendre polynomial expansion for σ_s and its integral σ_t are parameters corresponding to some physical quantities of vital importance. For instance, the slab width in the unit of mean free path: ${\sigma }_t^{-1}$ is the reciprocal of the total cross-section $${\sigma }_t=2\pi {\int }_{-1}^1{\sigma }_s(\omega )d\omega .$$

In the absorption-less case, the differential scattering cross-section is given by (4) $${\sigma }_s(\omega )={\sigma }_t\sum _{k=0}^{\infty }\frac{2k+1}{4\pi }{c}_k{P}_k(\omega ),\hspace{1em}{c}_0=1,\hspace{1em}{c}_1=2\pi {\int }_{-1}^1\omega {\sigma }_s(\omega )d\omega \equiv \overline{\omega },$$(4) with $P_k(\omega )$ being the Legendre polynomial of degree k, and $\overline{\omega }$ is the cosine of the mean scattering angle. The Fokker–Planck approximation to problem (1) is based on using spherical harmonics expansions and yields the following, degenerate-type partial differential equation (5) $$\mu \frac{\partial u}{\partial x}+\eta \frac{\partial u}{\partial y}+\xi \frac{\partial u}{\partial z}=\frac{{\sigma }_{tr}}{2}{\Delta }_{\mathbf{V}}u(\mathbf{x},\mathbf{v}),\hspace{1em}0<x<1,$$(5) associated with the same boundary data as (2) and (3), and with ${\Delta }_{\mathbf{V}}$ denoting the Laplace-Beltrami operator (6) $${\Delta }_{\mathbf{V}}:=\left[\frac{\partial }{\partial \mu }\right(1-{\mu }^2)\frac{\partial }{\partial \mu }+\frac{1}{1-{\mu }^2}\frac{{\partial }^2}{\partial {\varphi }^2}].$$(6)

Here, $\varphi$ is the angular variable appearing in the polar representation $\eta =\sqrt{1-{\mu }^2}\hspace{0.17em}\text{cos}\hspace{0.17em}\varphi$ and $\xi =\sqrt{1-{\mu }^2}\hspace{0.17em}\text{sin}\hspace{0.17em}\varphi$. Further, σ_tr is the transport cross-section defined by $${\sigma }_{tr}={\sigma }_t(c_0-c_1).$$

A thorough exposition of the Fokker–Planck operator as an asymptotic limit is given by Pomraning (1992). Due to successive asymptotic limits used in deriving the Fokker–Planck approximation, it is not obvious that this approximation is sufficiently accurate to be considered as a model for the pencil beams. However, for sufficiently small transport cross-section ${\sigma }_{tr}\ll 1$, Fermi proposed the following form of, projected, Fokker–Planck model: (7) $$\frac{\partial u}{\partial x}+\eta \frac{\partial u}{\partial y}+\xi \frac{\partial u}{\partial z}=\frac{{\sigma }_{tr}}{2}(\frac{{\partial }^2}{\partial {\eta }^2}+\frac{{\partial }^2}{\partial {\xi }^2})u(\mathbf{x},\mathbf{v}),\hspace{1em}0<x<1,$$(7)

with (8) $$u(0,y,z,\eta ,\xi )=\delta \left(y\right)\delta \left(z\right)\delta (\eta )\delta (\xi ).$$(8)

Fermi’s approach is different from the asymptotic ones and uses physical reasoning based on modeling cosmic rays. Note that the Fokker–Planck operator on the right hand side of (5), that is, (6), is the Laplacian on the unit sphere. The tangent plane to the unit sphere S² at the point ${\mu }_0:=(1,0,0)$ is an $\mathcal{O}({\eta }^2+{\xi }^2)$ approximation to the S² at the vicinity of μ₀. Extending $(\eta ,\xi )$ to ${\mathbb{R}}^2$, the Fourier transformation with respect to $y,z,\eta$, and ξ, assuming constant σ_tr, yields the following exact solution for the angular flux (9) $$u(x,y,z,\eta ,\xi )=\frac{3}{{\pi }^2{\sigma }_{tr}^2{x}^4}\text{exp}[-\frac{2}{{\sigma }_{tr}}(\frac{{\eta }^2+{\xi }^2}{x}-3\frac{y\eta +z\xi }{x^2}+3\frac{y^2+z^2}{x^3}\left)\right].$$(9)

The closed-form solution (9) was first derived by Fermi as referred in Rossi and Greisen (1941). Eyges (1948) has extended this exact solution to the case of an x-depending ${\sigma }_{tr}={\sigma }_{tr}(x)$. However, for the general case of ${\sigma }_{tr}(\mathbf{x})={\sigma }_{tr}(x,y,z)$, the closed-form analytic solution is not known. To obtain the scalar flux, we integrate (9) over $(\eta ,\xi )\in {\mathbb{R}}^2$: (10) $$\tilde{U}(x,y,z)={\int }_{{\mathbb{R}}^2}u(x,y,z,\eta ,\xi )d\eta d\xi =\frac{3}{2\pi {\sigma }_{tr}{x}^3}\text{exp}[-\frac{3}{2{\sigma }_{tr}}(\frac{y^2+z^2}{x^3}\left)\right].$$(10)

Equation (10)(10) $$\tilde{U}(x,y,z)={\int }_{{\mathbb{R}}^2}u(x,y,z,\eta ,\xi )d\eta d\xi =\frac{3}{2\pi {\sigma }_{tr}{x}^3}\text{exp}[-\frac{3}{2{\sigma }_{tr}}(\frac{y^2+z^2}{x^3}\left)\right].$$(10) satisfies the transverse diffusion equation (11) $$\frac{\partial \tilde{U}}{\partial x}=\frac{{\sigma }_{tr}{x}^2}{2}(\frac{{\partial }^2\tilde{U}}{\partial y^2}+\frac{{\partial }^2\tilde{U}}{\partial z^2}),$$(11) with (12) $$\tilde{U}(0,y,z)=\delta \left(y\right)\delta \left(z\right).$$(12)

Restricted to bounded phase-space domain, Fermi equation (7)(7) $$\frac{\partial u}{\partial x}+\eta \frac{\partial u}{\partial y}+\xi \frac{\partial u}{\partial z}=\frac{{\sigma }_{tr}}{2}(\frac{{\partial }^2}{\partial {\eta }^2}+\frac{{\partial }^2}{\partial {\xi }^2})u(\mathbf{x},\mathbf{v}),\hspace{1em}0<x<1,$$(7) can be written as the following “initial” boundary value problem (13) $$\{\begin{array}{cc}{u}_x+\mathbf{v}·{\nabla }_{\perp }u=\frac{{\sigma }_{tr}}{2}{\Delta }_{\mathbf{V}}u\hspace{1em}\hspace{1em}\hfill & \text{in}\hspace{1em}\Omega :={\Omega }_{\mathbf{x}}\times {\Omega }_{\mathbf{v}},\hfill \\ {\nabla }_{\mathbf{v}}u(x,x_{\perp },\mathbf{v})=0\hfill & \text{for}\hspace{1em}(x,x_{\perp },\mathbf{v})\in {\Omega }_{\mathbf{x}}\times \partial {\Omega }_{\mathbf{v}},\hfill \\ u(0,x_{\perp },\mathbf{v})=u_0(x_{\perp },\mathbf{v})\hfill & \text{for}\hspace{1em}(x_{\perp },\mathbf{v})\in {\Omega }_{x_{\perp }}\times {\Omega }_{\mathbf{v}}=:{\Omega }_{\perp },\hfill \\ u(x,x_{\perp },\mathbf{v})=0\hfill & \text{on}\hspace{1em}{\Gamma }_{\tilde{\beta }}^{-}\setminus \left\{\right(0,x_{\perp },\mathbf{v}\left)\right\},\hfill \end{array}$$(13) where $$\mathbf{v}=(\eta ,\xi ),{\nabla }_{\perp }=(\partial /\partial y,\partial /\partial z),{\Omega }_{\mathbf{x}}:=I_x\times {\Omega }_{\perp },{\Omega }_{\perp }\subset {\mathbb{R}}^2,{\Omega }_{\mathbf{v}}\subset {\mathbb{R}}^2$$ and (14) $${\Gamma }_{\tilde{\beta }}^{-}:=\left\{\right(x,x_{\perp },\mathbf{v})\in \partial \Omega ,\mathbf{n}·\tilde{\beta }<0\}$$(14) is the inflow boundary with respect to the characteristic line $\tilde{\beta }:=(1,\mathbf{v},0,0)$ and $\mathbf{n}$ is the outward unit normal to the boundary $\partial \Omega$. Note that, to derive energy estimates, the associated boundary data (viewed as a replacement for the initial data) at x = 0 is, in a sense, approximating the product of the Dirac’s delta functions on the right hand side of (8). Assuming that we can use separation of variables, we may write the data function u₀ as product of two functions $f(x_{\perp })$ and $g(\mathbf{v})$, $$u_0(x_{\perp },\mathbf{v})=f\left(x_{\perp }\right)g\left(\mathbf{v}\right).$$

The regularity of these functions has substantial impact in deriving theoretical stabilities and is essential in robustness of implemented results.

Finally, throughout the paper C will denote a generic constant, independent of the mesh parameters, unless it is explicitly expressed.

2. The phase-space standard Galerkin procedure

Below we introduce a framework that concerns a standard Galerkin discretization based on a quasi-uniform triangulation of the phase-space domain ${\Omega }_{\perp }:={\Omega }_{x_{\perp }}\times {\Omega }_{\mathbf{v}}:=I_{\perp }\times {\Omega }_{\mathbf{v}}$, where $I_{\perp }:=I_y\times I_z$. This is an extension of our studies in two dimensions in a flatland model (Asadzadeh and Larsen 2008). Previous numerical approaches are mostly devoted to the study of the one-dimensional problem see, for example, Larsen et al. (1985) and Prinja and Pomraning (1996).

Here we consider triangulation of the rectangular domains $I_{\perp }$ and ${\Omega }_{\mathbf{v}}:=I_{\eta }\times I_{\xi }$ into triangles ${\tau }_{\perp }$ and ${\tau }_{\mathbf{v}}$, and with the corresponding mesh parameters $h_{\perp }$ and $h_{\mathbf{v}}$, respectively. Then a general polynomial approximation of degree $\le r$ can be formulated in ${\mathbb{P}}_r(\tau ):={\mathbb{P}}_r({\tau }_{\perp })\otimes {\mathbb{P}}_r({\tau }_{\mathbf{v}})$. These polynomial spaces are more specified in the implementation section. We will assume a minimal angle condition on the triangles ${\tau }_{\perp }$ and ${\tau }_{\mathbf{v}}$ (see, e.g., Brenner and Scott 1994). Treating the beam’s entering direction x similar to a time variable, we let $\mathbf{n}:=\mathbf{n}(y,z,\mathbf{v})$ be the outward unit normal to the boundary of the phase-space domain ${\Omega }_{x_{\perp }}\times {\Omega }_{\mathbf{v}}$ at $(y,z,\mathbf{v})\in \partial {\Omega }_{\perp }$ where $\partial {\Omega }_{\perp }:=(\partial {\Omega }_{x_{\perp }}\times {\Omega }_{\mathbf{v}})\cup ({\Omega }_{x_{\perp }}\times \partial {\Omega }_{\mathbf{v}})$. Now set $\beta :=(\mathbf{v},0,0)$ and define the inflow (outflow) boundary as (15) $${\Gamma }_{\beta }^{-(+)}:=\left\{\right(x_{\perp },\mathbf{v})\in \Gamma :=\partial {\Omega }_{\perp }:\mathbf{n}·\beta <0(>0\left)\right\}.$$(15)

We shall also need an abstract finite element space as a subspace of a function space of Sobolev type, viz. (16) $${\mathcal{V}}_{h,\beta }\subset H_{\beta }^1:=\{w\in H^1(I_{\perp }\times {\Omega }_{\mathbf{v}}):w=0\text{on}{\Gamma }_{\beta }^{-}\}.$$(16)

Now for all $w\in H^1(I_{\perp }\times {\Omega }_{\mathbf{v}})\cap H^r(I_{\perp }\times {\Omega }_{\mathbf{v}})$ a classical standard estimate reads as (17) $$\underset{\chi \in {\mathcal{V}}_{h,\beta }}{\inf }\left|\right|w-\chi |{|}_j\le C{h}^{\alpha -j}|\left|w\right|{|}_{\alpha },\hspace{1em}j=1,2,\hspace{1em}1\le \alpha \le r\hspace{1em}\text{and}h=\max (h_{\perp },h_{\mathbf{v}}).$$(17)

To proceed let $\tilde{u}$ be an auxiliary interpolant of the solution u for Equation (13)(13) $$\{\begin{array}{cc}{u}_x+\mathbf{v}·{\nabla }_{\perp }u=\frac{{\sigma }_{tr}}{2}{\Delta }_{\mathbf{V}}u\hspace{1em}\hspace{1em}\hfill & \text{in}\hspace{1em}\Omega :={\Omega }_{\mathbf{x}}\times {\Omega }_{\mathbf{v}},\hfill \\ {\nabla }_{\mathbf{v}}u(x,x_{\perp },\mathbf{v})=0\hfill & \text{for}\hspace{1em}(x,x_{\perp },\mathbf{v})\in {\Omega }_{\mathbf{x}}\times \partial {\Omega }_{\mathbf{v}},\hfill \\ u(0,x_{\perp },\mathbf{v})=u_0(x_{\perp },\mathbf{v})\hfill & \text{for}\hspace{1em}(x_{\perp },\mathbf{v})\in {\Omega }_{x_{\perp }}\times {\Omega }_{\mathbf{v}}=:{\Omega }_{\perp },\hfill \\ u(x,x_{\perp },\mathbf{v})=0\hfill & \text{on}\hspace{1em}{\Gamma }_{\tilde{\beta }}^{-}\setminus \left\{\right(0,x_{\perp },\mathbf{v}\left)\right\},\hfill \end{array}$$(13) defined by (18) $$\mathcal{A}{(u-\tilde{u},\chi )}_{\perp }=0,\forall \chi \in {\mathcal{V}}_{h,\beta },$$(18) where (19) $$\mathcal{A}{(u,w)}_{\perp }={(u_x,w)}_{{\Omega }_{\perp }}+{(\mathbf{v}·{\nabla }_{\perp }u,w)}_{{\Omega }_{\perp }},$$(19) and ${(·,·)}_{\perp }:={(·,·)}_{{\Omega }_{\perp }}={(·,·)}_{I_{\perp }\times {\Omega }_{\mathbf{v}}}$.

With these notations the weak formulation for the problem (13) can be written as follows: for each $x\in (0,L]$, find $u(x,·)\in H_{\beta }^1$ such that, (20) $$\{\begin{array}{cc}\mathcal{A}{(u,\chi )}_{\perp }+\frac{1}{2}{({\sigma }_{tr}{\nabla }_{\mathbf{v}}u,{\nabla }_{\mathbf{v}}\chi )}_{\perp }=0\hspace{1em}\hspace{1em}\hfill & \forall \chi \in H_{\beta }^1,\hfill \\ u(0,x_{\perp },\mathbf{v})=u_0(x_{\perp },\mathbf{v})\hfill & \text{for}(x_{\perp },\mathbf{v})\in {\Gamma }_{\beta }^{+},\hfill \\ u(x,x_{\perp },\mathbf{v})=0\hfill & \text{on}{\Gamma }_{\beta }^{-}\setminus \left\{\right(0,x_{\perp },\mathbf{v}\left)\right\}.\hfill \end{array}$$(20)

Our objective is to solve the following finite element approximation for the problem (20): for each $x\in (0,L]$, find $u_h(x,·)\in {\mathcal{V}}_{h,\beta }$ such that, (21) $$\{\begin{array}{cc}\mathcal{A}{(u_h,\chi )}_{\perp }+\frac{1}{2}{({\sigma }_{tr}{\nabla }_{\mathbf{v}}{u}_h,{\nabla }_{\mathbf{v}}\chi )}_{\perp }=0\hspace{1em}\hspace{1em}\hfill & \forall \chi \in {\mathcal{V}}_{h,\beta },\hfill \\ u_h(0,x_{\perp },\mathbf{v})=u_{0,h}(x_{\perp },\mathbf{v})\hfill & \text{for}(x_{\perp },\mathbf{v})\in {\Gamma }_{\beta }^{+},\hfill \\ u_h(x,x_{\perp },\mathbf{v})=0\hfill & \text{on}{\Gamma }_{\beta }^{-}\setminus \left\{\right(0,x_{\perp },\mathbf{v}\left)\right\},\hfill \end{array}$$(21) where $u_{0,h}(x_{\perp },\mathbf{v})=\tilde{u}(0,x_{\perp },\mathbf{v})$.

2.1. A fully discrete scheme

For a partition of the interval $[0,L]$ into the subintervals $I_m:=(x_{m-1},x_m),m=1,2,\dots ,M$ with $k_m:=\left|I_m\right|:=x_m-x_{m-1}$, a finite element approximation U with continuous linear functions ${\psi }_m(x)$ on I_m can be written as (22) $$u_h(x,x_{\perp },\mathbf{v})=U_{m-1}(x_{\perp },\mathbf{v}){\psi }_{m-1}\left(x\right)+U_m(x_{\perp },\mathbf{v}){\psi }_m\left(x\right),$$(22) where $$x_{\perp }:=(y,z)$$ and (23) $${\psi }_{m-1}\left(x\right)=\frac{x_m-x}{k_m},\hspace{1em}{\psi }_m\left(x\right)=\frac{x-x_{m-1}}{k_m}.$$(23)

Hence, the setting (22) and (23) may be considered for an iterative, for example, backward Euler, scheme with continuous piecewise linear or discontinuous (with jump discontinuities at grid points x_m) piecewise linear functions for whole $I_x=[0,L]$.

To proceed we consider a normalized, rectangular domain ${\Omega }_{\mathbf{V}}$ for the velocity variable $\mathbf{v}$, as $(\eta ,\xi )\in [-1,1]\times [-1,1]$ and assume a uniform, “central adaptive” discretization mesh, viz. (24) $${\Omega }_{\mathbf{v}}^N:=\left\{{\mathbf{v}}_{i,j}\subset {\Omega }_{\mathbf{v}}|{\mathbf{v}}_{i,j}=({\eta }_i,{\xi }_j):=(\text{sin}\hspace{0.17em}\frac{i\pi }{2n},\hspace{0.17em}\text{sin}\hspace{0.17em}\frac{j\pi }{2n}),i,j=0,\pm 1,\dots \pm n\right\},$$(24) where $N={(2n+1)}^2$. Further, we assume that U has compact support in ${\Omega }_{\mathbf{V}}$. By a standard approach one can show that, for each $m=1,2,\dots ,M$, a finite element or finite difference solution $U_m^N$ obtained using the discretization (24) of the velocity domain ${\Omega }_{\mathbf{v}}$, satisfies the $L_2({\Omega }_{\mathbf{v}})$ error estimate (25) $$\left|\right|U_m-U_m^N|{|}_{L_2\left({\Omega }_{\mathbf{v}}\right)}\le \frac{C}{N^2}|\left|D_{\mathbf{V}}^2{U}_m\right(x_{\perp },·\left)\right|{|}_{L_2\left({\Omega }_{\mathbf{v}}\right)},$$(25) where $$D_{\mathbf{V}}^{2}w={(w_{\xi \xi }^2+2w_{\xi \eta }^2+w_{\eta \eta }^2)}^{1/2}.$$

Now we introduce a final, finite element, discretization using continuous piecewise linear basis functions ${\phi }_j(x_{\perp })$, on a partition ${\mathcal{T}}_h$ of the spatial domain ${\Omega }_{x_{\perp }}$, on a quasi-uniform triangulation with the mesh parameter h and obtain the fully discrete solution $U_m^{N,h}$. We introduce discontinuities on the direction of entering beam (on the x-direction). We also introduce jumps appearing in passing a collision site; say x_m, as the difference between the values at $x_m^{-}$ and $x_m^{+}$: (26) $$\left[U_m\right]:=U_m^{+}-U_m^{-},\hspace{1em}{U}_m^{\pm }:=\underset{s\to 0}{\lim }U(x_m\pm s,x_{\perp },\mathbf{v}).$$(26)

Due to the hyperbolic nature of the problem in $x_{\perp }$, for the solutions in the Sobolev space $H^{k+1}(·,x_{\perp },\mathbf{v})$ (see Adams 1975) for the exact definitions of the Sobolev norms and spaces) the final finite element approximation yields an $L_2({\Omega }_{x_{\perp }})$ error estimate, viz. (27) $$\left|\right|U_m^N-U_m^{N,k}|{|}_{L_2\left({\Omega }_{x_{\perp }}\right)}\le C{h}^{k+1/2}|\left|U_m^N\right(x_{\perp },·\left)\right|{|}_{H^{k+1}\left({\Omega }_{x_{\perp }}\right)}.$$(27)

To be specific, for each m and each ${\mathbf{v}}_{i,j}\in {\Omega }_{\mathbf{v}}$_, we obtain a spatially continuous version of the equations system (13) where, for u, we insert (28) $$U_m(x_{\perp },{\mathbf{v}}_{i,j})=\sum _{k=1}^K{U}_{m,k}\left({\mathbf{v}}_{i,j}\right){\phi }_k\left(x_{\perp }\right).$$(28)

Thus for each ${\mathbf{v}}_{i,j}\in {\Omega }_{\mathbf{v}}$ a variational formulation for a space-time like discretization in $(x,x_{\perp })$ of (22) reads as follows: find $U\in {\mathcal{V}}_{h,\beta }$ such that (29) $$\begin{array}{c}{\int }_{I_m}{\int }_{{\Omega }_{x_{\perp }}}{U}_x(x,x_{\perp },{\mathbf{v}}_{i,j}){\phi }_k\left(x_{\perp }\right)d{x}_{\perp }dx+{\int }_{I_m}{\int }_{{\Omega }_{x_{\perp }}}{\mathbf{v}}_{i,j}·{\nabla }_{\perp }U(x,x_{\perp },{\mathbf{v}}_{i,j}){\phi }_k\left(x_{\perp }\right)d{x}_{\perp }dx\hfill \\ ={\int }_{I_m}{\int }_{{\Omega }_{x_{\perp }}}\frac{{\sigma }_{tr}}{2}{\Delta }_{\mathbf{V}}U(x,x_{\perp },{\mathbf{v}}_{i,j}){\phi }_k\left(x_{\perp }\right)d{x}_{\perp }dx,\hspace{1em}\forall {\phi }_k\in {\mathcal{V}}_{h,\beta },\hfill \end{array}$$(29) where (30) $${\mathcal{V}}_{h,\beta }:=\{w\in {\mathcal{V}}_{\beta }|w_{{|}_{\tau }}\in {\mathbb{P}}_1(\tau ),w\text{is}\text{continuous}\}.$$(30)

This yields (31) $$\begin{array}{c}{\int }_{{\Omega }_{x_{\perp }}}(U_m(x_{\perp },{\mathbf{V}}_{i,j})-U_{m-1}(x_{\perp },{\mathbf{v}}_{i,j})\left){\phi }_k\right(x_{\perp })d{x}_{\perp }+\frac{k_m}{2}{\int }_{{\Omega }_{x_{\perp }}}{\mathbf{v}}_{i,j}·{\nabla }_{\perp }{U}_m(x_{\perp },{\mathbf{V}}_{i,j}\left){\phi }_k\right(x_{\perp })d{x}_{\perp }\hfill \\ +\frac{k_m}{2}{\int }_{{\Omega }_{x_{\perp }}}{\mathbf{v}}_{i,j}·{\nabla }_{\perp }{U}_{m-1}(x_{\perp },{\mathbf{v}}_{i,j}){\phi }_k\left(x_{\perp }\right)d{x}_{\perp }\hfill \\ =\frac{{\sigma }_{tr}}{2}\frac{k_m}{2}{\int }_{{\Omega }_{x_{\perp }}}{\Delta }_{\mathbf{v}}{U}_m(x_{\perp },{\mathbf{v}}_{i,j}){\phi }_k\left(x_{\perp }\right)d{x}_{\perp }+\frac{{\sigma }_{tr}}{2}\frac{k_m}{2}{\int }_{{\Omega }_{x_{\perp }}}{\Delta }_{\mathbf{v}}{U}_{m-1}(x_{\perp },{\mathbf{v}}_{i,j}){\phi }_k\left(x_{\perp }\right)d{x}_{\perp }.\hfill \end{array}$$(31)

Such an equation would lead to a linear system of equations which in compact form can be written as the following matrix equation: (32) $$\begin{array}{c}M{U}_m\left({\mathbf{v}}_{i,j}\right)-M{U}_{m-1}\left({\mathbf{v}}_{i,j}\right)+\frac{k_m}{2}{C}_{{\mathbf{v}}_{i,j}}{U}_m\left({\mathbf{v}}_{i,j}\right)+\frac{k_m}{2}{C}_{{\mathbf{v}}_{i,j}}{U}_{m-1}\left({\mathbf{v}}_{i,j}\right)\hfill \\ =\frac{{\sigma }_{tr}}{2}\frac{k_m}{2}{\Delta }_{\mathbf{v}}M{U}_m\left({\mathbf{v}}_{i,j}\right)+\frac{{\sigma }_{tr}}{2}\frac{k_m}{2}{\Delta }_{\mathbf{v}}M{U}_{m-1}\left({\mathbf{v}}_{i,j}\right).\hfill \end{array}$$(32)

Now considering $\mathbf{v}$-continuous version of (32): (33) $$\begin{array}{c}M{U}_m\left(\mathbf{v}\right)-M{U}_{m-1}\left(\mathbf{v}\right)+\frac{k_m}{2}{C}_{\mathbf{V}}{U}_m\left(\mathbf{v}\right)+\frac{k_m}{2}{C}_{\mathbf{V}}{U}_{m-1}\left(\mathbf{v}\right)\hfill \\ =\frac{{\sigma }_{tr}}{2}\frac{k_m}{2}{\Delta }_{\mathbf{V}}M{U}_m\left(\mathbf{v}\right)+\frac{{\sigma }_{tr}}{2}\frac{k_m}{2}{\Delta }_{\mathbf{V}}M{U}_{m-1}\left(\mathbf{v}\right),\hfill \end{array}$$(33) we may write (34) $$U_m(x_{\perp },\mathbf{v})=\sum _{k=1}^K\sum _{j=1}^J{U}_{m,k,j}{\chi }_j\left(\mathbf{v}\right){\phi }_k\left(x_{\perp }\right).$$(34)

Then a further variational form is obtained by multiplying (33) by χ_j, $j=1,2,\dots ,J$ and integrating over ${\Omega }_{\mathbf{V}}$: (35) $$\begin{array}{c}(M_{x_{\perp }}\otimes M_{\mathbf{v}})U_m+\frac{k_m}{2}(C_{x_{\perp }}\otimes {\tilde{M}}_{\mathbf{v}})U_m+\frac{{\sigma }_{tr}}{2}\frac{k_m}{2}(S_{\mathbf{v}}\otimes M_{x_{\perp }})U_m\hfill \\ \hspace{1em}\hspace{1em}=(M_{x_{\perp }}\otimes M_{\mathbf{v}})U_{m-1}-\frac{k_m}{2}(C_{x_{\perp }}\otimes {\tilde{M}}_{\mathbf{v}})U_{m-1}-\frac{{\sigma }_{tr}}{2}\frac{k_m}{2}(S_{\mathbf{v}}\otimes M_{x_{\perp }})U_{m-1},\hfill \end{array}$$(35) where ⊗ represents tensor products with the obvious notations for the coefficient matrices $M_{x_{\perp }},M_{\mathbf{v}}$ being the mass-matrices in spatial and velocity variables, $C_{x_{\perp }}$ is the convection matrix in space, ${\tilde{M}}_{\mathbf{v}}:=\mathbf{v}\otimes M_{\mathbf{v}}$ corresponds to the spatial convection terms with the coefficient $\mathbf{v}$: $\mathbf{v}·{\nabla }_{\perp }$, and finally $S_{\mathbf{v}}$ is the stiffness matrix in $\mathbf{v}$.

Now, given an initial beam configuration, U₀ = u₀, our objective is to use an iteration algorithm as the finite element version above or the corresponding equivalent backward Euler (or Crank-Nicholson) approach for discretization in the $x_{\perp }$ variable, and obtain successive U_m-values at the subsequent discrete $x_{\perp }$-levels. To this end the delicate issues of an initial data, viz. (8), as a product of Dirac delta functions, as well as the desired dose to the target that imposes the model to be transferred to a case having an inverse problem nature are challenging practicalities.

2.1.1. Standard stability estimates

We use the notion of the scalar products over a domain $\mathbf{D}$ and its boundary $\partial \mathbf{D}$ as ${(·,·)}_{\mathbf{D}}$ and ${〈·,·〉}_{\partial \mathbf{D}}$, respectively. Here, $\mathbf{D}$ can be $\Omega :=I_x\times {\Omega }_{\mathbf{x}}\times {\Omega }_{\mathbf{v}}$, $I_x\times {\Omega }_{\mathbf{x}},{\Omega }_{\mathbf{x}}\times {\Omega }_{\mathbf{v}}$, or possibly other relevant domains in the problem. Below we state and prove a stability lemma which, in some adequate norms, guarantees the control of the solution for the continuous problem by the data. The lemma is easily extended to the case of approximate solution.

We derive the stability estimate using the triple norm (36) $$\left|\right|\left|w\right||{|}_{\beta }^2={\int }_0^L{\int }_{{\Gamma }_{\beta }^{+}}{w}^2(\mathbf{n}·\beta )d\Gamma dx+|\left|{\sigma }_{tr}^{1/2}{\nabla }_{\mathbf{v}}w\right|{|}_{L_2(\Omega )}^2.$$(36)

Lemma 1

For u satisfying (13), we have the stability estimates (37) $$\underset{x\in I_x}{\sup }\left|\right|u(x,·,·)|{|}_{L_2({\Omega }_{\perp }\times {\Omega }_{\mathbf{v}})}\le |\left|u_0\right(·,·\left)\right|{|}_{L_2({\Omega }_{\perp }\times {\Omega }_{\mathbf{v}})},$$(37) and (38) $$\left|\right|\left|u\right||{|}_{\beta }\le |\left|u_0\right(·,·\left)\right|{|}_{L_2({\Omega }_{\perp }\times {\Omega }_{\mathbf{v}})}.$$(38)

Proof

We let $\chi =u$ in (20) and use (18) and (19) to obtain (39) $$\frac{1}{2}\frac{d}{dx}\left|\right|u|{|}_{L_2(I_{\perp }\times {\Omega }_{\mathbf{v}})}^2+{(\mathbf{v}·{\nabla }_{\perp }u,u)}_{I_{\perp }\times {\Omega }_{\mathbf{v}}}+\frac{1}{2}|\left|{\sigma }_{tr}^{1/2}{\nabla }_{\mathbf{v}}u\right|{|}_{L_2(I_{\perp }\times {\Omega }_{\mathbf{v}})}^2=0,$$(39) where using Green’s formula and with $\beta =(\mathbf{v},0)$ we have (40) $$\begin{array}{c}{(\mathbf{v}·{\nabla }_{\perp }u,u)}_{I_{\perp }\times {\Omega }_{\mathbf{v}}}={\int }_{{\Omega }_{\mathbf{v}}}({\int }_{I_{\perp }}(\mathbf{v}·{\nabla }_{\perp }u)u)d{x}_{\perp }d\mathbf{v}\\ =\frac{1}{2}{\int }_{{\Omega }_{\mathbf{v}}}(\mathbf{n}·\mathbf{v})u^{2}d\mathbf{v}\\ =\frac{1}{2}{\int }_{{\Gamma }_{\beta }^{+}}{u}^2(\mathbf{n}·\beta )d\Gamma \ge 0.\end{array}$$(40)

Thus (41) $$\frac{d}{dx}\left|\right|u|{|}_{L_2(I_{\perp }\times {\Omega }_{\mathbf{v}})}^2\le 0,$$(41) which yields (37) after integration over $(0,x)$ and taking supremum over $x\in I_x$. Integrating (39) over $x\in (0,L)$ and using (40) together with the definition of the triple norm $\left|\right||·||{|}_{\beta }$ we get (42) $$\left|\right|u(L,·,·)|{|}_{L_2(I_{\perp }\times {\Omega }_{\mathbf{v}})}^2+|\left|\right|u\left|\right|{|}_{\beta }^2=\left|\right|u(0,·,·)|{|}_{L_2(I_{\perp }\times {\Omega }_{\mathbf{v}})}^2,$$(42) and the estimate (38) follows. □

Using the same argument as above we obtain the semi-discrete version of the stability Lemma1.

Corollary 2

The semi-discrete solution u_h with $h=\max (h_{\perp },h_{\mathbf{v}})$ and standard Galerkin approximation in phase-space $I_{\perp }\times {\Omega }_{\mathbf{v}}$ satisfies the semi-discrete stability estimates: (43) $$\underset{x\in I_x}{\sup }\left|\right|u_h(x,·,·)|{|}_{L_2({\Omega }_{\perp }\times {\Omega }_{\mathbf{v}})}\le |\left|u_{0,h}\right(·,·\left)\right|{|}_{L_2({\Omega }_{\perp }\times {\Omega }_{\mathbf{v}})},$$(43) (44) $$\left|\right|\left|u_h\right||{|}_{\beta }\le |\left|u_{0,h}\right(·,·\left)\right|{|}_{L_2({\Omega }_{\perp }\times {\Omega }_{\mathbf{v}})}.$$(44)

2.1.2 Convergence

Below we state and prove an a priori error estimate for the finite element approximation u_h satisfying (21). The a priori error estimate will be stated in the triple norm defined by (36).

Lemma 3

[An a priori error estimate in the triple norm] Assume that u and u_h satisfy the continuous and discrete problems (20) and (21), respectively. Let $u\in H^r(\Omega )=H^r({\Omega }_{\mathbf{x}}\times {\Omega }_{\mathbf{v}})$, $r\ge 2$, then there is a constant C independent of v, u, and h such that (45) $$\left|\right||u-u_h||{|}_{\tilde{\beta }}\le C{h}^{r-1/2}|\left|u\right|{|}_{H^r(\Omega )}.$$(45)

Proof

. Taking the first equations in (20) and (21) and using (19) we end up with (46) $$\begin{array}{c}{({(u_h-\tilde{u})}_x,\chi )}_{I_{\perp }\times {\Omega }_{\mathbf{v}}}+{(\mathbf{v}·{\nabla }_{\perp }(u_h-\tilde{u}),\chi )}_{I_{\perp }\times {\Omega }_{\mathbf{v}}}+\frac{1}{2}{\left({\sigma }_{tr}{\nabla }_{\mathbf{v}}\right(u_h-\tilde{u}),{\nabla }_{\mathbf{v}}\chi )}_{I_{\perp }\times {\Omega }_{\mathbf{v}}}\hfill \\ =\frac{1}{2}{\left({\sigma }_{tr}{\nabla }_{\mathbf{v}}\right(u-\tilde{u}),{\nabla }_{\mathbf{v}}\chi )}_{I_{\perp }\times {\Omega }_{\mathbf{v}}}.\hfill \end{array}$$(46)

Let now $\chi =u_h-\tilde{u}$, then by the same argument as in the proof of Lemma 1 we get (47) $$\begin{array}{c}\frac{d}{dx}\left|\right|u_h-\tilde{u}|{|}_{L_2(I_{\perp }\times {\Omega }_{\mathbf{v}})}^2+{\int }_{{\Gamma }_{\beta }^{+}}(\mathbf{n}·\beta ){(u_h-\tilde{u})}^{2}d\Gamma +|\left|{\sigma }_{tr}^{1/2}{\nabla }_{\mathbf{v}}\right(u_h-\tilde{u}\left)\right|{|}_{L_2(I_{\perp }\times {\Omega }_{\mathbf{v}})}^2\hfill \\ \hspace{1em}\le \frac{1}{2}\left|\right|{\sigma }_{tr}^{1/2}{\nabla }_{\mathbf{v}}(u_h-\tilde{u})|{|}_{L_2(I_{\perp }\times {\Omega }_{\mathbf{v}})}^2+\frac{1}{2}|\left|{\sigma }_{tr}^{1/2}{\nabla }_{\mathbf{v}}\right(u-\tilde{u}\left)\right|{|}_{L_2(I_{\perp }\times {\Omega }_{\mathbf{v}})}^2,\hfill \end{array}$$(47) which yields (48) $$\begin{array}{c}\frac{d}{dx}\left|\right|u_h-\tilde{u}|{|}_{L_2(I_{\perp }\times {\Omega }_{\mathbf{v}})}^2+{\int }_{{\Gamma }_{\beta }^{+}}(\mathbf{n}·\beta ){(u_h-\tilde{u})}^{2}d\Gamma +\frac{1}{2}|\left|{\sigma }_{tr}^{1/2}{\nabla }_{\mathbf{v}}\right(u_h-\tilde{u}\left)\right|{|}_{L_2(I_{\perp }\times {\Omega }_{\mathbf{v}})}^2\hfill \\ \le \frac{1}{2}\left|\right|{\sigma }_{tr}^{1/2}{\nabla }_{\mathbf{v}}(u-\tilde{u})|{|}_{L_2(I_{\perp }\times {\Omega }_{\mathbf{v}})}^2.\hfill \end{array}$$(48)

Hence, integrating over $x\in (0,L)$ we obtain (49) $$\begin{array}{c}\left|\right|(u_h-\tilde{u})(L,·,·)|{|}_{L_2(I_{\perp }\times {\Omega }_{\mathbf{v}})}^2+{\int }_{{\Gamma }_{\beta }^{+}\setminus {\Gamma }_L}(\tilde{\mathbf{n}}·\tilde{\beta }){(u_h-\tilde{u})}^{2}d\Gamma +\frac{1}{2}|\left|{\sigma }_{tr}^{1/2}{\nabla }_{\mathbf{v}}\right(u_h-\tilde{u}\left)\right|{|}_{L_2(I_x\times I_{\perp }\times {\Omega }_{\mathbf{v}})}^2\hfill \\ \hspace{1em}\hspace{1em}\le \frac{1}{2}\left|\right|{\sigma }_{tr}^{1/2}{\nabla }_{\mathbf{v}}(u-\tilde{u})|{|}_{L_2(I_x\times I_{\perp }\times {\Omega }_{\mathbf{v}})}^2+|\left|\right(u_h-\tilde{u}\left)\right(0,·,·\left)\right|{|}_{L_2(I_{\perp }\times {\Omega }_{\mathbf{v}})}^2,\hfill \end{array}$$(49) where ${\Gamma }_{\beta }^{+}\setminus {\Gamma }_L:=\left\{\right\{L\}\times I_{\perp }\times {\Omega }_{\mathbf{v}}\}$. Now recalling that $u_h(0,·,·)=\tilde{u}(0,·,·)=u_{0,h}$ and the definition of $\left|\right||·||{|}_{\tilde{\beta }}$ we end up with (50) $$\left|\right||u_h-\tilde{u}||{|}_{\tilde{\beta }}\le |\left|{\sigma }_{tr}^{1/2}{\nabla }_{\mathbf{v}}\right(u-\tilde{u}\left)\right|{|}_{L_2(I_x\times I_{\perp }\times {\Omega }_{\mathbf{v}})}^2.$$(50)

Finally using the identity $u_h-u=(u_h-\tilde{u})+(\tilde{u}-u)$ and the interpolation estimate below we obtain the desired result. □

Proposition 4

(See Ciarlet 1941). Let $h^2\le {\sigma }_{tr}(\mathbf{x})\le h$, then there is an interpolation constant $\tilde{C}$ such that (51) $$\left|\right||u-\tilde{u}||{|}_{\tilde{\beta }}\le \tilde{C}{h}^{r-1/2}|\left|u\right|{|}_r.$$(51)

Proof

We rely on classical interpolation error estimates (see Brenner and Scott 1994 and Ciarlet 1941): Let $u\in H^r(\Omega )$, then there exists interpolation constants C₁ and C₂ such that for the nodal interpolant ${\pi }_h\in {\mathcal{V}}_{h,\tilde{\beta }}$ of u we have the interpolation error estimates (52) $$\left|\right|u-{\pi }_{h}u|{|}_s\le C_1{h}^{r-s}|\left|u\right|{|}_r,\hspace{1em}s=0,1,$$(52) (53) $$|u-{\pi }_{h}u{|}_{\tilde{\beta }}\le C_2{h}^{r-1/2}|\left|u\right|{|}_r,$$(53) where $$|\phi {|}_{\tilde{\beta }}:={\left({\int }_{{\Gamma }_{\tilde{\beta }}}{\phi }^2\right(\mathbf{n}·\tilde{\beta })d\Gamma )}^{1/2}.$$

Using the definition of the triple-norm we have that (54) $$\begin{array}{c}\left|\right||u-{\pi }_{h}u||{|}_{\tilde{\beta }}^2=|u-{\pi }_{h}u{|}_{\tilde{\beta }}^2+\left|\right|{\sigma }_{tr}^{1/2}{\nabla }_{\mathbf{v}}(u-{\pi }_{h}u)|{|}^2\\ \le |u-{\pi }_{h}u{|}_{\tilde{\beta }}^2+|\left|{\sigma }_{tr}^{1/2}\right|{|}_{\infty }^2\left|\right|u-{\pi }_{h}u|{|}_{H^1(\Omega )}^2\\ \le C_2^2{h}^{2r-1}\left|\right|u|{|}_r^2+C_1^2\underset{\mathbf{x}}{\sup }|{\sigma }_{tr}\left|h^{2r-1}\right|\left|u\right|{|}_r^2\\ =(C_2^2+C_1^2\underset{\mathbf{x}}{\sup }|{\sigma }_{tr}\left|\right)h^{2r-1}\left|\right|u|{|}_r^2,\end{array}$$(54) where in the last inequality we have used (52) and (53). Now choosing the constant $\tilde{C}={(C_2^2+C_1^2\underset{\mathbf{x}}{\sup }|{\sigma }_{tr}|)}^{1/2}$ we get the desired result.

This proposition yields the L₂ error estimate, viz.

Theorem 5

(L₂ error estimate). For $u\in H^r(\Omega )$ and $u_h\in {\mathcal{V}}_{h,\tilde{\beta }}$ satisfying (20) and (21), respectively, and with $h^2\le {\sigma }_{tr}\le h$, we have that there is a constant $C=C(\Omega ,f)$ such that (55) $$\left|\right|u-u_h|{|}_{L_2(\Omega )}\le C{h}^{r-3/2}|\left|u\right|{|}_r.$$(55)

Proof

Using the Poincaré inequality (56) $$\left|\right|u-u_h|{|}_{L_2(\Omega )}\le C|\left|{\nabla }_{\mathbf{v}}\right(u-u_h\left)\right|{|}_{L_2(\Omega )}\le \frac{C}{\min {\sigma }_{tr}^{1/2}}\left|\right|{\sigma }_{tr}^{1/2}{\nabla }_{\mathbf{v}}(u-u_h)|{|}_{L_2(\Omega )}.$$(56)

Further using Lemma 3 (57) $$\left|\right|{\sigma }_{tr}^{1/2}{\nabla }_{\mathbf{v}}(u-u_h)\left|\right|\le \left|\right||u-u_h||{|}_{\tilde{\beta }}\le C{h}^{r-1/2}|\left|u\right|{|}_r.$$(57)

Combining (56) and (57) and recalling that σ_tr is in the interval $[h^2,h]$ we end up with (58) $$\left|\right|u-u_h|{|}_{L_2(\Omega )}\le C{h}^{r-3/2}|\left|u\right|{|}_r,$$(58) and the proof is complete. □

3 Petrov-Galerkin approaches

Roughly speaking, in the Petrov-Galerkin method one adds a streaming term to the test function. The reason of such approach is described, motivated, and analyzed in the classical SD methods. Here, our objective is to briefly introduce a few cases of Petrov-Galerkin approaches in some lower dimensional geometry and implement them in both direct and adaptive settings. Some specific forms of the Petrov-Galerkin methods are studied in Johnson (1992) where the method of exact transport + projection is introduced. Also both the semi-streamline diffusion (SSD) as well as the characteristic streamline diffusion (CSD) methods, which in their simpler forms are implemented here, are studied in Asadzadeh (2002).

3.1 An SSD scheme

Here the main difference with the standard approach is that we employ modified test functions of the form $w+\delta \mathbf{v}·{\nabla }_{\perp }w$ with $\delta \ge {\sigma }_{tr}$. Further, we assume that w satisfies the vanishing inflow boundary condition of (13). Hence, multiplying the differential equation in (13) by $w+\delta (\mathbf{v}·{\nabla }_{\perp }w)$ and integrating over ${\Omega }_{\perp }={\Omega }_{{\mathbf{x}}_{\perp }}\times {\Omega }_{\mathbf{v}}$ we have a variational formulation, viz. (59) $$\begin{array}{c}{(u_x+\mathbf{v}·{\nabla }_{\perp }u-\frac{1}{2}{\sigma }_{tr}{\Delta }_{\mathbf{v}}u,w+\delta (\mathbf{v}·{\nabla }_{\perp }w\left)\right)}_{\perp }={(u_x,w)}_{\perp }+\delta {(u_x,\mathbf{v}·{\nabla }_{\perp }w)}_{\perp }\\ +{(\mathbf{v}·{\nabla }_{\perp }u,w)}_{\perp }+\delta {(\mathbf{v}·{\nabla }_{\perp }u,\mathbf{v}·{\nabla }_{\perp }w)}_{\perp }\\ +\frac{1}{2}{({\sigma }_{tr}{\nabla }_{\mathbf{v}}u,{\nabla }_{\mathbf{v}}w)}_{\perp }+\frac{\delta }{2}{({\sigma }_{tr}{\nabla }_{\mathbf{v}}u,{\nabla }_{\mathbf{v}}(\mathbf{v}·{\nabla }_{\perp }w\left)\right)}_{\perp }=0.\end{array}$$(59)

3.1.1 The SSD stability estimate

We let in (59) w = u and obtain the following identity: (60) $$\begin{array}{c}\frac{1}{2}\frac{d}{dx}\left|\right|u|{|}_{\perp }^2+\delta {(u_x,\mathbf{v}·{\nabla }_{\perp }u)}_{\perp }+\frac{1}{2}{\int }_{{\Gamma }_{\beta }^{+}}(\mathbf{n}·\beta )u^{2}d\Gamma +\delta ||\mathbf{v}·{\nabla }_{\perp }u|{|}_{\perp }^2\hfill \\ +\frac{1}{2}\left|\right|{\sigma }_{tr}^{1/2}{\nabla }_{\mathbf{v}}u|{|}_{\perp }^2+\frac{\delta }{2}{({\sigma }_{tr}{\nabla }_{\mathbf{v}}u,{\nabla }_{\mathbf{v}}(\mathbf{v}·{\nabla }_{\perp }u\left)\right)}_{\perp }=0.\hfill \end{array}$$(60)

Now it is easy to verify that the last term above can be written as (61) $$\delta {({\sigma }_{tr}{\nabla }_{\mathbf{v}}u,{\nabla }_{\mathbf{v}}(\mathbf{v}·{\nabla }_{\perp }u\left)\right)}_{\perp }=\delta {\int }_{{\Omega }_{\perp }}{\sigma }_{tr}({\nabla }_{\mathbf{v}}u·{\nabla }_{\perp }u+\frac{1}{2}\mathbf{v}·{\nabla }_{\perp }{\left(\right|{\nabla }_{\mathbf{v}}u\left|\right)}^2)d{x}_{\perp }d\mathbf{v}.$$(61)

Due to symmetry the second term in the integral above vanishes. Hence we end up with (62) $$\begin{array}{c}\frac{1}{2}\frac{d}{dx}\left|\right|u|{|}_{\perp }^2+\delta {(u_x,\mathbf{v}·{\nabla }_{\perp }u)}_{\perp }+\frac{1}{2}{\int }_{{\Gamma }_{\beta }^{+}}(\mathbf{n}·\beta )u^{2}d\Gamma +\delta ||\mathbf{v}·{\nabla }_{\perp }u|{|}_{\perp }^2\hfill \\ +\frac{1}{2}\left|\right|{\sigma }_{tr}^{1/2}{\nabla }_{\mathbf{v}}u|{|}_{\perp }^2+\frac{\delta }{2}({\sigma }_{tr}{\nabla }_{\mathbf{v}}u,{\nabla }_{\perp }u){}_{\perp }=0.\hfill \end{array}$$(62)

Next, we multiply the differential equation (13)(13) $$\{\begin{array}{cc}{u}_x+\mathbf{v}·{\nabla }_{\perp }u=\frac{{\sigma }_{tr}}{2}{\Delta }_{\mathbf{V}}u\hspace{1em}\hspace{1em}\hfill & \text{in}\hspace{1em}\Omega :={\Omega }_{\mathbf{x}}\times {\Omega }_{\mathbf{v}},\hfill \\ {\nabla }_{\mathbf{v}}u(x,x_{\perp },\mathbf{v})=0\hfill & \text{for}\hspace{1em}(x,x_{\perp },\mathbf{v})\in {\Omega }_{\mathbf{x}}\times \partial {\Omega }_{\mathbf{v}},\hfill \\ u(0,x_{\perp },\mathbf{v})=u_0(x_{\perp },\mathbf{v})\hfill & \text{for}\hspace{1em}(x_{\perp },\mathbf{v})\in {\Omega }_{x_{\perp }}\times {\Omega }_{\mathbf{v}}=:{\Omega }_{\perp },\hfill \\ u(x,x_{\perp },\mathbf{v})=0\hfill & \text{on}\hspace{1em}{\Gamma }_{\tilde{\beta }}^{-}\setminus \left\{\right(0,x_{\perp },\mathbf{v}\left)\right\},\hfill \end{array}$$(13) by $\delta u_x$ and integrate over $I_{\perp }\times {\Omega }_{\mathbf{v}}$ to get (63) $$\delta \left|\right|u_x|{|}^2+{(\delta u_x,\mathbf{v}·{\nabla }_{\perp }u)}_{\perp }+\frac{\delta }{2}{({\sigma }_{tr}{\nabla }_{\mathbf{v}}u,{\nabla }_{\mathbf{v}}{u}_x)}_{\perp }=0.$$(63)

The last inner product on the left hand side of (63) can be written as (64) $${({\sigma }_{tr}{\nabla }_{\mathbf{v}}u,{\nabla }_{\mathbf{v}}{u}_x)}_{\perp }=\frac{1}{2}\frac{d}{dx}{\int }_{I_{\perp }\times {\Omega }_{\mathbf{v}}}{\sigma }_{tr}|{\nabla }_{\mathbf{v}}u{|}^{2}d{x}_{\perp }d\mathbf{v}-\frac{1}{2}{\int }_{I_{\perp }\times {\Omega }_{\mathbf{v}}}\frac{\partial {\sigma }_{tr}}{\partial x}(\left|{\nabla }_{\mathbf{v}}u{|}^2\right)d{x}_{\perp }d\mathbf{v}.$$(64)

Now inserting (64) in (63) and adding the result to (62) we end up with (65) $$\begin{array}{c}\frac{1}{2}\frac{d}{dx}\left|\right|u|{|}_{\perp }^2+\delta ||u_x+\mathbf{v}·{\nabla }_{\perp }u|{|}_{\perp }^2+\frac{1}{2}{\int }_{{\Gamma }_{\beta }^{+}}(\mathbf{n}·\beta )u^{2}d\Gamma +\frac{1}{2}\left|\right|{\sigma }_{tr}^{1/2}{\nabla }_{\mathbf{v}}u|{|}_{\perp }^2\hfill \\ +\frac{\delta }{2}{({\sigma }_{tr}{\nabla }_{\mathbf{v}}u,{\nabla }_{\perp }u)}_{\perp }+\frac{\delta }{4}\frac{d}{dx}{\int }_{I_{\perp }\times {\Omega }_{\mathbf{v}}}{\sigma }_{tr}|{\nabla }_{\mathbf{v}}u{|}^{2}d{x}_{\perp }d\mathbf{v}\hfill \\ -\frac{\delta }{4}{\int }_{I_{\perp }\times {\Omega }_{\mathbf{v}}}\frac{\partial {\sigma }_{tr}}{\partial x}\left(\right|{\nabla }_{\mathbf{v}}u{|}^2)d{x}_{\perp }d\mathbf{v}=0.\hfill \end{array}$$(65)

Further, we use the Cauchy-Schwarz inequality to get (66) $$\left|{({\sigma }_{tr}{\nabla }_{\mathbf{v}}u,{\nabla }_{\perp }u)}_{\perp }\right|\le \frac{1}{2}\left|\right|{\sigma }_{tr}^{1/2}{\nabla }_{\perp }u|{|}_{\perp }^2+\frac{1}{2}|\left|{\sigma }_{tr}^{1/2}{\nabla }_{\mathbf{v}}u\right|{|}_{\perp }^2.$$(66)

Finally with an additional symmetry assumption on $x_{\perp }$ and $\mathbf{v}$ convection terms as (this is motivated by forward peakedness assumption in angle and energy which is used in deriving the Fokker–Planck/Fermi equations) (67) $$\left|\right|{\sigma }_{tr}^{1/2}{\nabla }_{\perp }u|{|}_{\perp }\sim |\left|{\sigma }_{tr}^{1/2}{\nabla }_{\mathbf{v}}u\right|{|}_{\perp },$$(67) and assuming that σ_tr is nonincreasing in the beam’s penetration direction, that is, $\frac{\partial {\sigma }_{tr}}{\partial x}\le 0$, we may write (65) as (68) $$\begin{array}{c}\frac{1}{2}\frac{d}{dx}\left(\right|\left|u\right|{|}_{\perp }^2+\frac{1}{2}\delta {\int }_{I_{\perp }\times {\Omega }_{\mathbf{v}}}{\sigma }_{tr}|{\nabla }_{\mathbf{v}}u{|}^{2}d{x}_{\perp }d\mathbf{v})+\frac{1}{2}{\int }_{{\Gamma }_{\beta }^{+}}(\mathbf{n}·\beta )u^{2}d\Gamma +\delta \left|\right|u_x+\mathbf{v}·{\nabla }_{\perp }u|{|}_{\perp }^2\hfill \\ +\frac{1}{2}(1-\delta )\left|\right|{\sigma }_{tr}^{1/2}{\nabla }_{\mathbf{v}}u|{|}_{\perp }^2\le 0.\hfill \end{array}$$(68)

As a consequence for sufficiently small δ (actually $\delta \approx {\sigma }_{tr}^{1/2}\ll 1$) we have, for example, (69) $$\frac{d}{dx}\left(\right|\left|u\right|{|}_{\perp }^2+\frac{1}{2}\delta {\int }_{I_{\perp }\times {\Omega }_{\mathbf{v}}}{\sigma }_{tr}|{\nabla }_{\mathbf{v}}u{|}^{2}d{x}_{\perp }d\mathbf{v})<0,$$(69) and hence $\left|\right|u|{|}_{\perp }^2+\frac{\delta }{2}{\int }_{I_{\perp }\times {\Omega }_{\mathbf{v}}}{\sigma }_{tr}|{\nabla }_{\mathbf{v}}u{|}^{2}d{x}_{\perp }d\mathbf{v}$ is strictly decreasing in x. Consequently, for each $x\prime \in [0,L]$ we have that (70) $$\left|\right|u(x\prime ,·,·)|{|}_{\perp }^2+\frac{\delta }{2}|\left|{\sigma }_{tr}^{1/2}{\nabla }_{\mathbf{v}}u\right(x\prime ,·,·\left)\right|{|}_{\perp }^2\le \left|\right|u(0,·,·)|{|}_{L_2(I_{\perp }\times {\Omega }_{\mathbf{v}})}^2+\frac{\delta }{2}|\left|{\sigma }_{tr}^{1/2}{\nabla }_{\mathbf{v}}u\right(0,·,·\left)\right|{|}_{\perp }^2.$$(70)

Thus, summing up we have proved the following stability estimates.

Proposition 6

Under the assumption (67) the following $L_2(I_{\perp }\times {\Omega }_{\mathbf{v}})$ stability holds true (71) $$\left|\right|u(L,·,·)|{|}_{\perp }^2+\frac{\delta }{2}|\left|{\sigma }_{tr}^{1/2}{\nabla }_{\mathbf{v}}u\right(L,·,·\left)\right|{|}_{\perp }^2\le \left|\right|u_0|{|}_{\perp }^2+\frac{\delta }{2}|\left|{\sigma }_{tr}^{1/2}{\nabla }_{\mathbf{v}}{u}_0\right|{|}_{\perp }^2.$$(71)

Moreover, we have also the second stability estimate (72) $$\left|\right|\left|u\right||{|}_{\tilde{\beta }}^2+\delta ||u_x+\mathbf{v}·{\nabla }_{\perp }u|{|}_{L_2(\Omega )}^2\le C\left(\right|\left|u_0\right|{|}_{L_2(I_{\perp }\times {\Omega }_{\mathbf{v}})}^2+\delta \left|\right|{\sigma }_{tr}^{1/2}{\nabla }_{\mathbf{v}}{u}_0\left|{|}_{L_2(I_{\perp }\times {\Omega }_{\mathbf{v}})}^2\right),$$(72) which is a consequence of (71) and (68).

4 Model problems in lower dimensions

We consider now a forward peaked narrow radiation beam entering into the symmetric domain $I_y\times I_{\eta }=[-y_0,y_0]\times [-{\eta }_0,{\eta }_0]$; $(y_0,{\eta }_0)\in {\mathbb{R}}_{+}^2$ at (0, 0) and penetrating in the direction of the positive x-axis. Then the computational domain $\Omega$ of our study is a three-dimensional slab with $(x,y,\eta )\in \Omega =I_x\times I_y\times I_{\eta }$ where $I_x=[0,L]$. In this way, the problem (13) will be transformed into the following lower dimensional model problem: (73) $$\{\begin{array}{cc}{u}_x+\eta u_y=\frac{1}{2}{\sigma }_{tr}{u}_{\eta \eta }\hspace{1em}& (x,y,\eta )\in \Omega ,\\ u_{\eta }(x,y,\pm {\eta }_0)=0\hspace{1em}& (x,y)\in I_x\times I_y,\\ u(0,y,\eta )=f(y,\eta )\hspace{1em}& (y,\eta )\in I_y\times I_{\eta },\\ u(x,y,\eta )=0\hspace{1em}& \text{on}{\Gamma }_{\beta }^{-}\setminus \left\{\right(0,y,\eta \left)\right\}.\end{array}$$(73) For this problem, we implement two different versions of the SD method: the SSD and the CSD.

4.1 The SSD method

In this version, we derive a discrete scheme for computing the approximate solution u_h of the exact solution u using the SD method for discretizing the $(y,\eta )$-variables (corresponding to multiplying the equation by test functions of the form $w+\delta \eta w_y$) combined with the backward Euler method for the x-variable. We start by introducing the bilinear forms $a(·,·)$ and $b(·,·)$ for the problem (73) as (74) $$\begin{array}{c}a(u,w)={(\eta u_y,w)}_{\perp }+\delta {(\eta u_y,\eta w_y)}_{\perp }+\frac{1}{2}{({\sigma }_{tr}{u}_{\eta },w_{\eta })}_{\perp }\hfill \\ +\frac{1}{2}\delta {({\sigma }_{tr}{u}_{\eta },w_y+\eta w_{y\eta })}_{\perp }-\frac{1}{2}\delta {\int }_{I_y}{\sigma }_{tr}\eta u_{\eta }{w}_y{|}_{\eta =-{\eta }_0}^{\eta ={\eta }_0}dy,\hfill \\ b(u,w)={(u,w)}_{\perp }+\delta {(u,\eta w_y)}_{\perp },\hfill \end{array}$$(74) where ${(·,·)}_{\perp }:={(·,·)}_{I_y\times I_{\eta }}$. Then the continuous problem reads as for each $x\in (0,L]$, find $u(x,·)\in H_{\beta }^1$ such that $$b(u_x,w)+a(u,w)=0,\hspace{1em}\forall w\in H_{\beta }^1,$$ where $$H_{\beta }^1:=\{w\in H^1(I_y\times I_{\eta });w=0\text{on}{\Gamma }_{\beta }^{-}\},$$

and (75) $${\Gamma }_{\beta }^{-}:=\left\{\right(y,\eta )\in \Gamma :=\partial (I_y\times I_{\eta }),\text{with}\mathbf{n}·\beta <0\},$$(75) with $\beta =(\eta ,0)$. Then the SSD method for the continuous problem (73) reads as follows: for each $x\in (0,L]$, find $u_h(x,·)\in {\mathcal{V}}_{h,\beta }$ such that (76) $$b(u_{h,x},w)+a(u_h,w)=0,\hspace{1em}\forall w\in {\mathcal{V}}_{h,\beta },$$(76) where ${\mathcal{V}}_{h,\beta }\subset H_{\beta }^1$ consists of continuous piecewise polynomials. Next, we write the global discrete solution by separation of variables as (77) $$u_h(x,y,\eta )=\sum _{j=1}^N{U}_j\left(x\right){\varphi }_j(y,\eta ),$$(77) where N is the number of degrees of freedom. Letting $w={\varphi }_i$ for $i=1,2,\dots ,N$ and inserting (77) into (76) we get the following discrete system of equations: (78) $$\sum _{j=1}^{N}U{\prime }_j\left(x\right)b({\varphi }_j,{\varphi }_i)+\sum _{j=1}^N{U}_j\left(x\right)a({\varphi }_j,{\varphi }_i)=0,\hspace{1em}i=1,2,\dots ,N.$$(78)

Equation (78)(78) $$\sum _{j=1}^{N}U{\prime }_j\left(x\right)b({\varphi }_j,{\varphi }_i)+\sum _{j=1}^N{U}_j\left(x\right)a({\varphi }_j,{\varphi }_i)=0,\hspace{1em}i=1,2,\dots ,N.$$(78) in matrix form can be written as (79) $$BU\prime \left(x\right)+AU\left(x\right)=0,$$(79) with $U={[U_1,\dots ,U_N]}^T,B=(b_{ij}),b_{ij}=b({\varphi }_j,{\varphi }_i)$_, and $A=(a_{ij}),a_{ij}=a({\varphi }_j,{\varphi }_i),i,j=1,2,\dots ,N$. We apply now the backward Euler method for further discretization of Equation (79)(79) $$BU\prime \left(x\right)+AU\left(x\right)=0,$$(79) in variable x, and with the step size k_m, to obtain an iterative form viz (80) $$B(U^{m+1}-U^m)+k_{m}A{U}^{m+1}=0.$$(80)

The equation above can be rewritten as a system of equations for finding the solution $U^{m+1}$ (at “time” level $x=x_{m+1}$) on iteration m + 1 from the known solution U^m from the previous iteration step m: (81) $$[B+k_{m}A]U^{m+1}=B{U}^m.$$(81)

4.2 CSD method

In this part, we construct an oriented phase-space mesh to obtain the CSD method. Before formulating this method, we need to construct a new subdivision of $\Omega =I_x\times I_y\times I_{\eta }$. To this end and for $m=1,2,\dots ,M$, we define a subdivision of ${\Omega }_m:=I_m\times I_y\times I_{\eta },I_m:=[x_{m-1},x_m]$, into elements $${\widehat{\tau }}_m=\left\{\right(x,y+(x-x_m)\eta ,\eta ):(y,\eta )\in \tau \in {\mathcal{T}}_h,x\in I_m\},$$

where ${\mathcal{T}}_h$ is a previous triangulation of $I_{\perp }$. Then we introduce, slabwise, the function spaces $${\widehat{\mathcal{V}}}_m=\{\widehat{w}\in C({\Omega }_m):\widehat{w}(x,y,\eta )=w(y+(x-x_m)\eta ,\eta ),w\in {\mathcal{V}}_{h,\beta }\}.$$

In other words ${\widehat{\mathcal{V}}}_m$ consists of continuous functions $\widehat{w}(x,y,\eta )$ on Ω_m such that $\widehat{w}$ is constant along characteristics $(\widehat{y},\widehat{\eta })=(y+x\eta ,\eta )$ parallel to the sides of the elements ${\widehat{\tau }}_m$, meaning that the derivative in the characteristic direction vanishes: ${\widehat{w}}_x+\eta {\widehat{w}}_y=0$. The SD method can now be reduced to the following formulation (where only the σ_tr-term survives): find ${\widehat{u}}_h$ such that, for each $m=1,2,\dots ,M,$ ${\widehat{u}}_h{|}_{{\Omega }_m}\in {\widehat{\mathcal{V}}}_m$ and (82) $$\begin{array}{c}\frac{1}{2}{\int }_{{\Omega }_m}{\sigma }_{tr}{\widehat{u}}_{h,\eta }{w}_{\eta }\mathit{\text{dxdyd}}\eta +{\int }_{I_{\perp }}{\widehat{u}}_{h,+}(x_{m-1},y,\eta )w_{+}(x_{m-1},y,\eta )dyd\eta \hfill \\ ={\int }_{I_{\perp }}{\widehat{u}}_{h,-}(x_{m-1},y,\eta )w_{+}(x_{m-1},y,\eta )dyd\eta ,\hspace{1em}\forall w\in {\widehat{\mathcal{V}}}_m.\hfill \end{array}$$(82)

Here, for definition of ${\widehat{u}}_{h,+},{\widehat{u}}_{h,-},w_{+}$ we refer to (26).

5 Adaptive algorithm

In this section, we formulate an adaptive algorithm, which is used in computations of the numerical examples studied in Section 6. This algorithm improves the accuracy of the computed solution u_h of the model problem (73). In the sequel for simplicity we denote $I_y\times I_{\eta }$ also by ${\Omega }_{\perp }$ (this, however, should not be mixed with the notation in the theoretical Sections 1–3).

The mesh refinement recommendation: We refine the mesh in neighborhoods of those points in $I_y\times I_{\eta }$ where the error ${\epsilon }_k=|u-u_h^k|$ attains its maximal values. More specifically, we refine the mesh in such subdomains of $I_y\times I_{\eta }$ where (83) $${\epsilon }_k\ge \tilde{\gamma }\underset{{\Omega }_{\perp }}{\max }{\epsilon }_k.$$(83)

Here $\tilde{\gamma }\in (0,1)$ is a number which should be chosen computationally and $u_h^k$ denotes the computed solution on the kth refinement of the mesh.

The steps in adaptive algorithm

Step 0. Choose a fixed $\tilde{\gamma }\in (0,1)$ and an initial coarse mesh $I_m^0\times {\tau }^0$ in $I_x\times I_y\times I_{\eta }$, and obtain the numerical solution $u_h^0,(n=0)$.

Step 1. Assume that considering (83), we have made $n-1,n-1>0$ successive iterations and computed numerical solution $u_h^{n-1}$ on ${\tau }^{n-1}$ using any of the finite element methods introduced in Section 4.

Step 2. Refine those elements in the mesh ${\tau }^{n-1}$ for which (84) $${\epsilon }_{n-1}\ge \tilde{\gamma }\underset{{\Omega }_{\perp }}{\max }{\epsilon }_{n-1},$$(84) and call this new mesh τⁿ.

Step 3. Compute $u_h^n$ on the new mesh τⁿ and stop mesh refinements when $\left|\right|u_h^n-u_h^{n-1}|{|}_{L_2({\Omega }_{\perp })}<\text{tol}$, where tol is a total tolerance chosen by the user. Otherwise go to Step 2, relabel n as n + 1 and continue.

6 Numerical examples

In this section, we present numerical examples which show the performance of an adaptive finite element method for the solution of the model problem (73). Here, all computations are performed in Matlab COMSOL Multiphysics using module LiveLink for MATLAB. As initial data we have chosen even functions replacing the Dirac delta functions. Hence we choose symmetric phase-space domain of computation as $${\Omega }_{\perp }=I_y\times I_{\eta }:=\left\{(y,\eta )\in (-1,1)\times (-1,1)\right\}.$$

Our tests are performed with a fixed diffusion coefficient ${\sigma }_{tr}=0.002$. Further, due to smallness of the parameters δ and σ_tr, the terms that involve the product $\delta {\sigma }_{tr}$ are assumed to be negligible. In the backward Euler scheme, used for discretization in x-variable, we solve the system of equations (67)(67) $$\left|\right|{\sigma }_{tr}^{1/2}{\nabla }_{\perp }u|{|}_{\perp }\sim |\left|{\sigma }_{tr}^{1/2}{\nabla }_{\mathbf{v}}u\right|{|}_{\perp },$$(67) which ends up with a discrete (computed) solution $U^{m+1}$ of (73) at the time iteration m + 1 and with the time step k_m which has been chosen to be $k_m=0.01$.

Previous computational studies, for example Asadzadeh and Larsen (2008), have shown oscillatory behavior of the solution u_h when the SSD method was used, and layer behavior when the standard Galerkin method was applied to solve the model problem (73). Formation of strong boundary layers appears in our Galerkin approaches for the electron beams. This is the case even for a solution with Maxwellian initial data at different depths, see, for example, Naqos (2005).

In this work, we significantly improve results of Asadzadeh and Larsen (2008), and in particular Naqos (2005), by using the adaptive algorithm of Section 5 on the locally adaptively refined meshes. More specifically, the oscillations as well as the boundary layers in non-refined meshes, were removed using the adaptive mesh refinement and also reducing the dominance of the convection term.

All our computational results are for the depth x = 1 (the broadening phenomenon, which is not reported in here, was checked for x = 2).

For the model problem (73), with the initial data $u(0,y,\eta )=\delta (y)\delta (\eta )$, the analytic solution is given by (85) $$u(x,y,\eta )=\frac{\sqrt{3}}{\pi {\sigma }_{tr}{x}^2}\text{exp}\hspace{0.17em}\left[-\frac{2}{{\sigma }_{tr}}\left(\frac{3y^2}{x^3}-\frac{3y\eta }{x^2}+\frac{{\eta }^2}{x}\right)\right].$$(85)

We have performed the following computational tests:

• Test 1. Solution of the model problem (73) with a “Dirac type” initial condition (86) $$u(0,y,\eta )=f(y,\eta )=1/(y^2+{\eta }^2+\alpha ),\hspace{1em}(y,\eta )\in {\Omega }_{\perp },$$(86)

  for different values of the parameter $\alpha \in (0,1)$.

• Test 2. Solution of the model problem (73) with “Maxwellian type” initial condition (87) $$u(0,y,\eta )=f(y,\eta )={\hspace{0.17em}\text{exp}\hspace{0.17em}}^{-(y^2+{\eta }^2+\alpha )},\hspace{1em}(y,\eta )\in {\Omega }_{\perp },$$(87) for different values of $\alpha \in (0,1)$.

• Test 3. Solution of the model problem (73) with a “hyperbolic type” initial condition (88) $$u(0,y,\eta )=f(y,\eta )=\frac{1}{\sqrt{y^2+{\eta }^2+\alpha }},\hspace{1em}(y,\eta )\in {\Omega }_{\perp },$$(88) for $\alpha =0.19$.

The initial data above are all approximations for $\delta (y)\delta (\eta )$ and the computed solutions are compared with (85), which are mainly of the form (85).

The analytic solutions with the initial data as in Tests 1–3 can be obtained through Green’s function approach which is a rather involved procedure. Since we have a stable system these solutions should not be very different from the expression for u in (85).

We used piecewise linear polynomial approximation in Test 1. This, however, was not giving satisfactory results in the other two tests. Therefore, we choose higher order polynomials (of degree p = 2 for Test 2 and p = 3 for Test 3) which performed with corresponding higher convergence rates in, for example, comparing $e_{n-1}/e_n$.

6.1. Test 1

In this test, we compute numerical simulations for the problem (73) with a “Dirac type” initial condition (86) and for different values of the parameter $\alpha \in (0,1)$ in (86), where we use adaptive algorithm of Section 4 on the locally adaptively refined meshes. These meshes were refined according to the error indicator (84) in the adaptive algorithm. For computation of the finite element solution we employ SSD method of Section 4.1. We performed two set of numerical experiments:

• Test 1(a). We take $\tilde{\gamma }=0.5$ in (84). This choice of the parameter allows to refine the mesh τ not only at the center of the domain ${\Omega }_{\perp }$, but also at the boundaries of ${\Omega }_{\perp }$.

• Test 1(b). We take $\tilde{\gamma }=0.7$ in (84). Such choice of the parameter allows to refine the mesh τ only at the middle of the domain ${\Omega }_{\perp }$.

Our computational tests have shown that the values for $\alpha \in (0.05,0.1)$ give smaller computational errors $e_n=\left|\right|u-u_h^n|{|}_{L_2({\Omega }_{\perp })}$ than the other α-values. The results of the computations for $\alpha =0.1$ are presented in Tables 1 and 2 for Tests 1(a) and 1(b), respectively. Using these tables and Figures 1 and 2 we observe that we have obtained significant reduction of the computational error $e_n=\left|\right|u-u_h^n|{|}_{L_2({\Omega }_{\perp })}$ on the adaptively refined meshes. Using Tables 1 and 2 we observe that the reduction of the computational error is faster and more significant in the case (a) than in the case (b). Thus, choosing the parameter $\tilde{\gamma }=0.5$ in (84) gives a better computational result and smaller error e_n than $\tilde{\gamma }=0.7$. This allows us to conclude that the refinement of the mesh τ not only at the center of the domain ${\Omega }_{\perp }$, but also at the boundaries of ${\Omega }_{\perp }$ give significantly smaller computational error $e_n=\left|\right|u-u_h^n|{|}_{L_2({\Omega }_{\perp })}$.

Figure 1. Test 1(a). (a)–(d) Locally adaptively refined meshes of Table 1. (e) Computed solution on the four times adaptively refined mesh (d).

Figure 2. Test 1(b). (a)–(d) Locally adaptively refined meshes of Table 2. (e) Computed solution on the four times adaptively refined mesh (d).

Table 1. Test 1(a). Computed errors $e_n=\left|\right|u-u_h^n|{|}_{L_2({\Omega }_{\perp })}$ and $e_{n-1}/e_n$ on the adaptively refined meshes. Here, the solution $u_h^n$ is computed using semi-streamline diffusion method of Section 4.1 with $\tilde{\gamma }=0.5$ in the adaptive algorithm and $\alpha =0.1$ in (86).

No. of refinement, n	No. of elements	No. of vertices	DOF	$e_n=\left|\right|u-u_h^n|{|}_{L_2}$	$e_{n-1}/e_n$
0	272	157	157	8.364e-03
1	1176	591	591	2.134e-03	3.92
2	4704	2268	2268	5.368e-04	3.97
3	17616	8878	8878	1.345e-04	3.99
4	69864	35231	35231	3.363e-05	4.00

Table 2. Test 1(b). Computed errors $e_n=\left|\right|u-u_h^n|{|}_{L_2({\Omega }_{\perp })}$ and $e_{n-1}/e_n$ on the adaptively refined meshes. Here, the solution $u_h^n$ is computed using semi-streamline diffusion method of Section 4.1 with $\tilde{\gamma }=0.7$ in the adaptive algorithm and $\alpha =0.1$ in (86).

No. of refinement, n	No. of elements	No. of vertices	DOF	$e_n=\left|\right|u-u_h^n|{|}_{L_2}$	$e_{n-1}/e_n$
0	272	157	157	8.364e-03
1	1088	585	585	8.278e-03	1.01
2	4352	2257	2257	2.105e-03	3.93
3	17408	8865	8865	5.290e-04	3.98
4	69632	35137	35137	1.325e-04	3.99

We present the final solution $u_h^4$ computed on the four times adaptively refined mesh on the Figure 1(e) for Test 1(a) and on the Figure 2(e) for Test 1(b). We note that in both cases we have obtained smoother computed solution $u_h^4$ without any oscillatory behavior. This is a significant improvement of the result of Asadzadeh and Larsen (2008) where mainly oscillatory solution could be obtained.

6.2. Test 2

In this test, we perform numerical simulations for the problem (73) with Maxwellian initial condition (87) and for different values of the parameter $\alpha \in (0,1)$. Again we use the error indicator (84) in the adaptive algorithm for local refinement of meshes and perform two set of tests as in the case of Test 1 and with the same values on the parameter $\tilde{\gamma }$.

For finite element discretization we use CSD method as in the Test 1. To be able to control the formation of the layer which appears at the central point $(y,\eta )=(0,0)$ we use different values of $\alpha \in (0,1)$ inside the function (87). Our computational tests show that the value of the parameter $\alpha =0.19$ is optimal.

We present results of our computations for $\alpha =0.19$ in Tables 3 and 4 for Tests 2(a) and 2(b), respectively. Using these tables and Figures 3 and 4 once again, we observe significant reduction of the computational error $e_n=\left|\right|u-u_h^n|{|}_{L_2({\Omega }_{\perp })}$ on the adaptively refined meshes. Using Tables 3 and 4 again, we observe more significant reduction of the computational error in the case (a) than in the case (b). Thus, choosing the parameter $\tilde{\gamma }=0.5$ in (84) yields better computational results.

Figure 3. Test 2(a). (a)–(d) Locally adaptively refined meshes of Table 3. (e) Computed solution on the four times adaptively refined mesh (d).

Figure 4. Test 2(b). (a)–(d) Locally adaptively refined meshes of Table 4. (e) Computed solution on the four times adaptively refined mesh (d).

Table 3. Test 2(a). Computed values of errors $e_n=\left|\right|u-u_h^n|{|}_{L_2({\Omega }_{\perp })}$ and $e_{n-1}/e_n$ on the adaptively refined meshes. Here, the solution $u_h^n$ is computed using characteristic streamline diffusion method with $\tilde{\gamma }=0.5$ in the adaptive algorithm.

No. of refinement, n	No. of elements	No. of vertices	DOF	$e_n=\left|\right|u-u_h^n|{|}_{L_2}$	$e_{n-1}/e_n$
0	272	157	585	2.582e-04
1	1288	592	2267	3.242e-05	7.97
2	4552	2267	8911	4.062e-06	7.98
3	17628	8911	35432	5.085e-07	7.99
4	69812	35432	141612	6.362e-08	7.99

Table 4. Test 2(b). Computed values of errors $e_n=\left|\right|u-u_h^n|{|}_{L_2({\Omega }_{\perp })}$ and $e_{n-1}/e_n$ on the adaptively refined meshes. Here, the solution $u_h^n$ is computed using characteristic streamline diffusion method with $\tilde{\gamma }=0.7$ in the adaptive algorithm.

No. of refinement, n	No. of elements	No. of vertices	DOF	$e_n=\left|\right|u-u_h^n|{|}_{L_2}$	$e_{n-1}/e_n$
0	272	157	585	2.114e-02
1	1088	585	2257	2.706e-03	7.81
2	4352	2257	8865	3.427e-04	7.90
3	17408	8865	35137	4.307e-05	7.96
4	69632	35137	139905	5.398e-06	7.98

Final solution $u_h^4$ computed on the four times adaptively refined mesh is shown on Figure 3(e) for Test 2(a) and on Figure 4 for Test 2(b). Again we observe that, with the above numerical values for the parameters, α and $\tilde{\gamma }$ we have avoided the formation of layers and in both tests we have obtained smooth computed solution $u_h^4$.

6.3. Test 3

In this test, we perform numerical simulations of the problem (73) with hyperbolic initial condition (88) on the locally adaptively refined meshes. Taking into account results of our previous Tests 1 and 2 we take fixed value of $\alpha =0.19$ in (88). For finite element discretization we used the SSD method of Section 4. We again perform two set of tests with different values of $\tilde{\gamma }$ in (84): in Test 3(a) we choose $\tilde{\gamma }=0.5$, and in Test 3(b) we assign this parameter to be $\tilde{\gamma }=0.7$.

We present results of our computations in Tables 5 and 6. Using these tables and Figures 5 and 6 we observe significant reduction of the computational error $e_n=\left|\right|u-u_h^n|{|}_{L_2({\Omega }_{\perp })}$ on the adaptively refined meshes. Final solutions $u_h^4$ computed on the four times adaptively refined meshes are shown on Figure 5(e) for Test 3(a) and on Figure 6(e) for Test 3(b), respectively.

Figure 5. Test 3(a). (a)–(d) Locally adaptively refined meshes of Table 5. (e) Computed solution on the four times adaptively refined mesh (d).

Figure 6. Test 3(b). (a)–(d) Locally adaptively refined meshes of Table 6. (e) Computed solution on the four times adaptively refined mesh (d).

Table 5. Test 3(a). Computed values of errors $e_n=\left|\right|u-u_h^n|{|}_{L_2({\Omega }_{\perp })}$ and $e_{n-1}/e_n$ on the adaptively refined meshes. Here, the solution $u_h^n$ is computed using semi-streamline diffusion method with $\tilde{\gamma }=0.5$ in the adaptive algorithm.

No. of refinement, n	No. of elements	No. of vertices	DOF	$e_n=\left|\right|u-u_h^n|{|}_{L_2}$	$e_{n-1}/e_n$
0	272	157	1285	1.565e-05
1	1271	597	5115	9.732e-07	16.08
2	5084	2267	20937	6.052e-08	16.08
3	20336	9075	79825	3.771e-09	16.05

Table 6. Test 3(b). Computed values of errors $e_n=\left|\right|u-u_h^n|{|}_{L_2({\Omega }_{\perp })}$ and $e_{n-1}/e_n$ on the adaptively refined meshes. Here, the solution $u_h^n$ is computed using semi-streamline diffusion method with $\tilde{\gamma }=0.7$ in the adaptive algorithm.

No. of refinement, n	No. of elements	No. of vertices	DOF	$e_n=\left|\right|u-u_h^n|{|}_{L_2}$	$e_{n-1}/e_n$
0	272	157	1285	1.484e-05
1	1088	585	5017	9.289e-07	15.98
2	4352	2257	19825	5.799e-08	16.02
3	17408	8865	78817	3.620e-09	16.02

7. Conclusion

In this work, FEM is applied to compute approximate solution of a, degenerate type, convection dominated convection-diffusion problem. We studied different finite element discretizations for the solutions of pencil-beam models based on Fermi and Fokker–Planck equations. The objective was to derive stability estimates and prove optimal convergence rates (due to the maximal available regularity of the exact solution). We specified some “goal-oriented” numerical schemes derived using a variety of Galerkin methods such as standard Galerkin, SSD, characteristic Galerkin, and CSD methods.

Our focus has been in two of these approximation schemes: (i) the SSD and (ii) the CSD methods. The SSD scheme automatically adds extra diffusion to the system, whereas in the CSD the shock capturing property is more pronounced. Both methods are efficient in improving stability and also layer resolution.

For these two settings, we formulated an adaptive algorithm. The Fermi equation, with the initial data of the form $\delta (y)\delta (\eta )$ (however with constant σ_tr) has a closed-form analytic solution. This analytic solution is used as an “adequate approximation” to the analytic solutions of our test problems. This is due to the following facts that we have proved stability of the considered schemes and in our numerical tests the initial data are approximating $\delta (y)\delta (\eta )$, the initial data for the Fermi equation. Also the closed-form solution is used to make local mesh refinements.

Numerically, we tested our adaptive algorithm for different types of initial data in (73) in three tests with different mesh refinement parameters $\tilde{\gamma }$ in the mesh refinement criterion (84). The initial data in the tests, although all approximating the product of Dirac functions, are chosen with different kind of singularities.

The goal of our numerical experiments was to remove oscillatory behavior of the computed solution in Asadzadeh and Larsen (2008) as well as removing of the formation of the artificial and boundary layers that appeared, for example, in Naqos (2005).

Based on the results of this study, we conclude that the SD approach decreases the dominance of the convection. Further the, local refinement based-adaptivity strategy can remove both the oscillatory behavior of the computed solutions and layer formations appearing for problems with nonsmooth initial data.

Additional information

Funding

The research of first and second authors is supported by Swedish Research Council (VR).