Universal Boltzmann Equation Solver Ver. 1.0

1. Overview

   This is a FORTRAN library for "Universal Boltzmann Equation Solver" (UBER), which
   solves the general form of Fokker-Planck equation and Boltzmann equation, 
   diffusive or non-diffusive, that appear in modeling planetary radiation belts. 
   Users can freely specify (1) the coordinate system, (2) boundary geometry and
   boundary conditions, and (3) the equation terms and coefficients. The solver
   works for problems in one to three spatial dimensions.

   The solver is based upon the mathematical theory of stochastic differential
   equations. By its nature, the solver scheme is intrinsically Monte Carlo, and the
   solutions thus contain stochastic uncertainty, though the user may dictate an 
   arbitrarily small relative tolerance of the stochastic uncertainty at the cost of
   longer Monte Carlo iterations. Ref [1] gives a thorough description of the theory
   and applications of the solver.

2. Compilation and configuration

   The package of UBER contains the following files and directories:

      README         This file
      LICENSE        Software license
      CHANGELOG      Log of changes
      Makefile       Makefile of the UBER library
      Makefile.in    Makefile include that defines directories and compilers
      src            Source code
      scripts        MATLAB scripts for processing input and output data files
      examples       Examples of UBER usage

   UBER depends on the following well-maintained libraries:

      DCMT              (github.com/MersenneTwister-Lab/dcmt)
      REFERENCE BLAS    (www.netlib.org/blas/)
      LAPACK            (www.netlib.org/lapack/)

   The current release has been tested using dcmt version 0.6.1, refblas version
   3.8.0 and lapack version 3.9.0. The user must edit Makefile.in to specify the
   directories of these libraries and their include files.

   To compile UBER, in the directory of this README, edit Makefile.in to specify
   your compilers and the relevant library directories, and then do "make". Library
   archive file "libuber.a" will be created in this directory, and the compilation 
   is done. To use the library, "use uber" in your source code, and link your code 
   with libuber.a, libdcmt.a, librefblas.a and liblapack.a when compiling the source.
   In addition, "make oclean" cleans the intermediate object files generated in the
   compiling process, and "make clean" removes everything that "make" would create.

   UBER is parallelized by OpenMP. For its proper execution, set these environmental
   variables in bash

      export OMP_STACKSIZE=256M (or higher)
      export OMP_NUM_THREADS= the number of threads you want to use

   in csh or tcsh

      setenv OMP_STACKSIZE=256M (or higher)
      setenv OMP_NUM_THREADS= the number of threads you want to use

   Make sure you add these environmental variables to your .bashrc, .cshrc, or
   .tcshrc script.

3. The library structure and functions

   The library contains a physics branch and a mathematics branch. The physics
   branch is a collection of physical constants and functions that often appear in
   radiation belt studies. They are contained and explained in the source file
   ./src/utilities/functions.F90.

   The mathematics branch fully defines and solves the partial differential equation
   (PDE) problem. It is consisted of three parts: (1) a domain part that states the
   computational domain, the initial condition, the boundary geometry and the 
   boundary conditions; (2) an equation part that specifies the coordinate system 
   (in terms of the gradient of logarithm of the Jacobian determinant), the equation
   terms and the equation coefficients; and (3) the solver, with the solver data 
   types and parameters collectively defined and explained in the source file
   ./src/params/mod_typedef_params.F90.

   To specify the domain and the equation, the user must edit the self-explanatory
   source file user_input.F90 following the templates provided in the examples 
   directory. This source file is to be compiled together with your code making use
   of UBER. See the Makefile in the examples directory for an illustration.

   There are two solver subroutines in the library: a lower-level subroutine, 
   SOLVER_ENGINE, that directly takes a spacetime point and gives the solution there
   (see example0); and a higher-level subroutine, SOLVER, that reads lists of
   spacetimes from a binary input file and writes solutions into a binary output 
   file (see example1-3). The input and output files could be processed by the 
   MATLAB scripts provided in the ./scripts directory. Solver parameters, some of
   which must be adjusted per the concrete problem, are read in at run time from the
   text file UBER_params.in. This will be further explained in the next section.

4. Examples 

   In each ./examples/example? directory, do "make" to compile the example
   executable "main.x". Then typing "./main.x" runs the example. 

   4.0 Example0

   This example solves a diffusion problem in a three-dimensional unit sphere with
   internal sources and a partially reflective boundary. The PDE of the problem
   is

         df      d2    d2    d2
         -- = D(--- + --- + ---)f + Sf + v(x,y,z)                             (1)
         dt     dx2   dy2   dz2

   where the differential operators shall all be understood as partial 
   differentials, and the coefficients D = 0.1, S = 2.5, and

                            10^(-6)
         v(x,y,z) = -------------------------                                 (2)
                    1 + sqrt(x^2 + y^2 + z^2)
   
   Initial condition is

                            x^2 + y^2 + z^2
         f(0,x,y,z) = exp(- ---------------)                                  (3)
                                 0.02

   and the boundary condition is

           df   1  |
         (D-- + -f)|   = 0                                                    (4)
           dr   2  |r=1

   which says a half of the diffusive flux is reflected whereas the other half is
   lost on the sphere surface. r is the radius of the sphere.

   Point solutions are obtained at t = 0.05 at three locations along the sphere
   radii: (0.44, 0, 0), (0, 0.64, 0), and (0, 0, 0.84). The subroutine
   SOLVER_ENGINE is called in its verbose mode as indicated by the existence of the
   character argument "V", which would generate a detailed statistical diagnosis 
   file for the solution. For more background information of this example, see Ref
   [1] Section 4.1.

   Values of the first three parameters in UBER_params.in should be adjusted for 
   each problem to yield optimal performance. The first two parameters dictate the 
   range of temporal stepsizes the solver may take in simulating stochastic 
   trajectories, and the third parameter indicates the root-mean-square (RMS) 
   spatial stepsize of the simulated stochastic trajectories the solver should 
   achieve. The RMS spatial stepsize shall be small enough compared to the size of 
   the domain and any scale length of the equation coefficients to allow sufficient
   resolution of the problem. In this example, the size of the domain is O(1) and
   all equation coefficients are slowly varying, therefore an RMS spatial stepsize
   0.01 is sufficient.

   For a stochastic process, the relation between the RMS spatial stepsize (dS) and
   the temporal stepsize (dt) is

         dS = sqrt[2tr(D)dt]                                                  (5)
                       =

   where tr(D) is trace of the diffusion tensor D which in this example is a rank-3
            =                                   =
   diagonal matrix with each diagonal component 0.1. From Eq. (5) one can calculate
   that the temporal stepsize corresponding to dS = 0.01 is about 10^(-4). The range
   of temporal stepsizes should enclose this value, and be properly widened so that
   the adaptive algorithms of the solver can take full advantage. 

   4.1 Example1

   This example solves exactly the same problem as Example0 but in a spherical
   coordinate system, in which the problem is reduced to one-dimensional. The
   Jacobian determinant is r^2, and the gradient of its logarithm is (2/r, 0, 0) in
   spherical coordinates, which is provided in the source file user_input.F90 to
   indicate the coordinate system.

   With a singularity of the equation at the pole of the spherical coordinates, a
   boundary condition must be imposed at r = 0 as

         df|
         --|   = 0                                                            (6)
         dr|r=0

   In its numerical implementation, to avoid the pole, this boundary condition is 
   actually imposed at r = epsBnd, which is a very small number effecting the 
   numerical thickness of all boundaries in UBER. The value of the parameter epsBnd
   is related to the given RMS spatial stepsize (which should be much smaller than 
   the dimension of the domain), and is printed to screen together with other 
   parameters at the beginning of the code execution.

   This example uses the higher-level subroutine SOLVER which takes solution
   spacetimes from an input file. To generate the input file, run the MATLAB script
   example1_input.m. After the main program execution, run the script 
   example1_plot.m to plot the UBER solutions against reference ones obtained from 
   finite difference method.

   4.2 Example2

   This example solves a one-dimensional Fokker-Planck equation with additional 
   jumping (non-diffusive) transport. The equation is given as

          d          d2         /                   /
         --f(t,x) = ---f(t,x) - | P(x->y)f(t,x)dy + | P(y->x)f(t,y)dy         (7)
         dt         dx2         /                   /

   where the integration kernel P(x->y) is the jumping rate from x to y. The
   computational domain is x from 0 to 1, with reflective boundaries on both ends.
   To simplify the situation, the jumping rate is chosen as a delta function

         P(x->y) = 10*delta(x - y + 1/2)                                      (8)

   which says that a particle initially at x (x<1/2) has a transition rate of 10 to
   jump to y = x + 1/2. With this simplification, the integrals in Eq. (7) can be
   worked out explicitly, which yield 

         /              ( 10, for x < 1/2
         | P(x->y)dy = <                                                      (9)
         /              ( 0,  for x>= 1/2

   and

         /                    ( 0,            for x < 1/2
         | P(y->x)f(t,y)dy = <                                                (10)
         /                    ( 10f(t,x-1/2), for x>= 1/2

   For the general situations where P(x->y) is not a delta function, functions
   carrying out these integrals shall be called at assigning these equation terms.
   For more discussion on UBER solving integro-differential equations, see Ref [1]
   Section 5.

   Eq. (7) is solved with an initial Gaussian distribution centered at x = 1/2. From
   the plot of solutions generated by example2_plot.m, one can see two signatures of
   the jumping transport: (1) immediately after start, the distribution is non-zero
   near x = 1; and (2) as time goes by, the distribution approaches an equilibrium
   that is asymmetric about x = 1/2.

   This example takes about 10 minutes to run using 86 CPUs with 2.1 GHz floating
   point operation frequency.

   4.3 Example3

   This example illustrates how gridded data are used in setting boundary geometry
   and equation coefficients in UBER, and how periodic boundary conditions are
   implemented. The equation is a Boltzmann equation arising from Earth's inner
   radiation belt dynamics, which has the form

         df    df   1  d       df
         -- + w-- = - --(G Dxx --) + S                                        (11)
         dt    dy   G dx       dx

   where f is electron phase space density, x is electron equatorial pitch angle, y
   is electron drift phase, w is electron drift frequency, Dxx is pitch-angle
   diffusion coefficient, S is the source rate caused by cosmic ray albedo neutron
   decay, and G is the Jacobian determinant of the (x,y) coordinates. Periodic
   boundary condition is imposed for y such that

         f(t,x,y=0) = f(t,x,y=2pi)                                            (12) 

   Ref [1] Section 4.4 provides more background information about this example.
   
   The auxilliary data files ex3_bounce_loss_cone.dat and ex3_dipole_lshell.dat,
   located in ./examples/data/, store the local bounce loss cone angle and dipole
   L-shell as functions of geomagnetic longitude along the electron drift shell.
   They were generated using the uber_create_data_file.m script (in ./scripts/) to
   conform with the UBER input binary file format.

   The script example3_plot.m plots two snapshots of f versus y and x at t = 2 and 4
   hr. One should see distinct west-east gradient of f below x = 60 deg that changes
   very little in time, and increased f above x = 60 deg that is a result of the
   source.

   This example takes about 45 minutes to run using 86 CPUs with 2.1 GHz floating
   point operation frequency.

5. Comment and bug reports

   If you find any comment or bug, please kindly send an email to Liheng Zheng
   ([email protected]).

6. How to cite

   [1] Zheng, L., L. Chen, A. A. Chan, P. Wang, Z. Xia and X. Liu (2021),
       UBER v1.0: a universal kinetic equation solver for radiation belts,
       Geosci. Model Dev., 14, 5825-5842, doi:10.5194/gmd-14-5825-2021.