Main program for the Fokker-Planck example problem.
For example, the program is run as ./Fokker_Planck 2 4 Fokker_Planck_2D.input glex This sets dimension = 2 , number of quadrature points = 4, Fokker_Planck_2D.input as the input file to be read and glex as the moment ordering option. int main (int argc, char **argv) { Declare variable solver_params as an object of the inputs class. solver_params can then access the member function input_reader that reads the input file inputs solver_params Declare as extern variables defined elsewhere in inputs.c. num_steps : total number of time steps, print_interval : intervals at which output is to be printed, dt : time step, tau : time scale, mu : mean of initial normal distribution, std : standard deviation of initial normal distribution extern int num_steps, print_interval; extern double dt; extern double mu, std; extern matrix_int dtuple_set; Invoke the the input_reader function which is a member of the inputs class of which solver_params is declared to be a member. The argv[3] is the third variable in the command line options when running the example. solver_params.input_reader(argv[3]); Declare the dimension of the problem. atoi(argv[1]) assigns the first variable in the command line options to dimension. size_t dimension = atoi(argv[1]);
Declare the number of quadrature points to be used in the problem. atoi(argv[2]) assigns the second variable in the command line options to num_points. size_t num_points = atoi(argv[2]);
The total number of moments required for the chosen dimension and num_points. size_t itot = num_points + dimension*num_points;
Declare the solution variable quadwts as a member of the quadpts class. quadpts quadwts(itot); Generate initial conditions for the quadrature points by taking a univariate Gaussian. There are many ways of specifying initial distributions. For this problem in which we are looking at the stationary state, the initial condition is less important. A valid initial bivariate pDF is generated from a univariate Gaussian. Compute quadrature points and weights associated with it. Assign the same quadrature points to both variables in the bivariate case. The quadrature weights for the initial bivariate case are the quadrature weights of the univaruate Gaussian. moments moms_init(2*num_points); quadpts quads_init(2*num_points); @ endcode Moments of initial univariate Gaussian with known mean and standard deviation @code moms_init = moms_init.univariate_gaussian_raw(1, num_points, mu, std_dev); Quadrature points and weights from zeros of univariate Gauss Hermite polynomials quads_init= moms_init.get_quadpts_univariate_Gaussian(1, num_points); Now assign initial conditions to the quadwts vector as described above. First assign the weights and the weights*quadrature points for the first variable. for (size_t i = 0; i < 2*num_points; i++){ quadwts.set(i,quads_init.get(i))} Now assign the product of weights*quadrature points for the second variable. for (size_t i = 2*num_points; i < 3*num_points; i++){ quadwts.set(i,drand48()*quads_init.get(i-num_points)); Set containing indices for a valid moment set (non-trivial for a bivariate distribution). The argv[4] is the fourth variable entered in the command line. Only glex option is supported at present. matrix_int valid_dtuple; dtuple_set = valid_dtuple_set(dimension, num_points, argv[4]); In this example we choose to solve the quadrature points and weights. The function solve is called with the dimension, the number of quadrature points num_points, the time step dt, the number of steps num_steps, the interval to print print_interval. dimension, num_points can be changed in the command line. This example only works with dimension = 1. The other variables can be changed in the input file. Exit the program return 0;
}
|