VERB4D
Documentation for the VERB4D_SOLVER

Includes classes, namespaces, and files

For use in conjunction with Matlab in order to solve for Phase Space Density, Diffusion, Convection etc.

Diffusion + convection code, VERB4D.

The code solves the Fokker-Planck equation with convection terms:

\[ \frac{\partial f}{\partial t} = v_{\Phi} \frac{\partial f}{\partial \Pi} + v_{R} \frac{\partial f}{\partial R} + \frac{1}{G} \frac{\partial}{\partial L} G D_{LL} \frac{\partial f}{\partial L} + \frac{1}{G} \frac{\partial}{\partial V} G (D_{VV} \frac{\partial f}{\partial V} + D_{VK} \frac{\partial f}{\partial K}) + \frac{1}{G} \frac{\partial}{\partial K} G (D_{KV} \frac{\partial f}{\partial V} + D_{KK} \frac{\partial f}{\partial K}) + Losses * f + Sources \]

where f is a Phase Space Density (PSD); Φ is MLT angle, R is radial distance from Earth center, $ V≡μ∙(K+0.5)^2 $ is a combination of adiabatic invariants μ and $ K≡ {J}{\sqrt{8μ}} $; L is the third adiabatic invariant; D are bounce-averaged diffusion coefficients; $ v_Φ $, $ v_R $ are drift velocities; G is the Jacobian of the transformation from (μ,J,L) to (V,K,L)

Books, required to understand the code:

VARIABLES for understanding code

Convection is solved using Convection_1D_ULTIMATE_QUICKEST6.h which implements Leonard BP (1988) Universal Limiter for transient interpolation modeling of the advective transport equations: the ULTIMATE conservative difference scheme, NASA technical Memorandum 100916 ICOMP-88-11 http://www.hadian.ir/teaching/CompHydr/3.pdf

Diffusion is solved using Diffusion_1D() and Diffusion_2D() for the standard cases.

Three other methods are also used for finding diffusion in 2D: these can be found in Diffusion_ADI1.h, Diffusion_ADI2.h, and Diffusion_ADI3.h

The parameters.ini file created in matlab will specify the inversion method. If "Lapack" is specified then Diffusion_2D() will be used. Otherwise 1 of the 3 ADI methods will be used to inversion. All 3 ADI methods can use multiple threads while Lapack can only use 1. The current chosen method of inversion is Diffusion_ADI3() if "ADI" is specified. It can be changed in VERB4D_SOLVER.cpp

From current testing all of the ADI methods produce bad results - namely negative PSD values. All examples should thus be run with Lapack until the ADI methods are fixed. Unfortunately Lapack is not compatable with multiple threads so running the examples with Lapack will be much slower than ADI and not scalable. Newer versions of Lapack may allow multithreading however.

The main solver function can be found in VERB4D_Solver.cpp It consists of a series of PSD calculations done for every time step that is specified.

This will be used in order to take the input data from matlab store them into Matrix1D, Matrix2D, Matrix3D, and Matrix4D using ReadInitialData.h

The matrix classes listed above are used to store the data in 1,2,3 or 4 dimension. One example is a 4D matrix containing P,R,V,K

There is also the standardized MatrixND for any size matrices which is used in conjunction with UpdatableMatrix which can be updated with new data.

Most of the computation such as numerically approximating derivatives are done with MatrixSolver.h using CalculationMatrix. These calculation matrices are made up of DiagMatrix which is a mapping of an int (which diagonal number) to a 1d matrix of values (for that diagonal).

There is also a Parameters class which holds paramters and their value as defined in the file they came from. These parameters get the values specified in the parameters.ini text file and set variables to determine inversion method, number of threads, using .mat files, etc.

This code was built with MATLAB capabilities including reading and writing .mat files. In order to switch all files to .mat for increased speed and percision change the use_matlab variable in the MATLAB Conv_Dif.m file for whichever example you are running to true. If the machine that this code is being run on does not have matlab installed change the MATLAB_CAPABLE variable in Matrix.h to false