VERB_code_2.2  2
 All Classes Namespaces Files Functions Variables Typedefs Macros Pages
MatrixSolver.h
Go to the documentation of this file.
1 /**
2  * \file MatrixSolver.h
3  *
4  * Making model matrixes, solving model matrixes.
5  *
6  * Matrix form of linear equations: A*X[t+1] = B*X[t],
7  * where A - model matrix, B - RHS, X[t] - known values of function (PSD), X[t+1] - unknown values of function.
8  *
9  * In that file there are procedures for making model matrix for 1d-diffusion, 2d-diffusion, some ideas of 3d-diffusion and mixed terms.
10  * Solver for tridiagonal matrix, solver by gauss method and iteration method (upper relaxation).
11  *
12  * This file is under development, and has a lot of commented code, unfinished etc code.
13  *
14  * There is a checked version for 1d diffusion (or split method of 2d, 3d diffusions) - 1d_universal_solver.
15  *
16  * \author Developed under supervision of the PI Yuri Shprits
17  */
18 
19 #ifndef MatrixSolver_H
20 #define MatrixSolver_H
21 
22 #include "../Grid/Grid.h"
23 #include "../DiffusionCoefficient/DiffusionCoefficient.h"
24 
25 #include "../GMRES/itlin.h"
26 
27 
28 
29 /// Make matrix for 1d diffusion.
30 bool MakeMatrix(double *A,
31  double *B1,
32  double *C,
33  double *R,
34  double *x, // one dimentional grid
35  double tau, // lifetime
36  double n, // power of x
37  double f_lower, // lower boundary conditions
38  double f_upper, // upper boundary conditions
39  int nx, // number of grid points
40  double dt, // time step
41  double *Dxx, // diffusion coefficient
42  double taulc, // lifetime in loss cone
43  double alc, // loss cone size
44  int g_flag, // diffusion flag
45  string lower_border_condition_type,// type of lower boundary conditions: 0 - for value, 1 - for derivative
46  string upper_border_condition_type,// type of upper boundary conditions: 0 - for value, 1 - for derivative
47  string approximationMethod = "AM_Split_C");
48 
49 
50 /// Solver for 1d general equation
51 bool SolveMatrix(double *f, // phase space density function
52  double *A,
53  double *B1,
54  double *C,
55  double *R,
56  double f_lower, // lower boundary conditions
57  double f_upper, // upper boundary conditions
58  int nx, // number of grid points
59  double dt); // time step
60 
61 
62 // Matrix solvers
63 void over_relaxation_diag(DiagMatrix &A, Matrix1D<double> &B, Matrix1D<double> &X, int max_steps = 1e5, double EE = 1e-3);
65 void gauss_solve(double *a, double *b, int n);
66 
67 /** Solver for a system of equations with 3-diagonal matrix.
68  *
69  * Au = r
70  * Where A = diag(a, b, c),
71  *
72  * \param a[] - array, diagonal '-1' of the matrix
73  * \param b[] - array, diagonal '0' of the matrix
74  * \param c[] - array, diagonal '+1' of the matrix
75  * \param r[] - array, r-vector
76  * \param u[] - array, result
77  * \param n - size of the matrix
78  */
79 bool tridag(double a[], double b[], double c[], double r[], double u[], long n);
80 
81 //void gmres_wrapout(CalculationMatrix &A, Matrix1D<double> &B, Matrix1D<double> &X, int maxiter = 100, int i_max = 10, double tol = 1e-10);
83 
84 // A test
86 
88  CalculationMatrix &matr_A, CalculationMatrix &matr_B, CalculationMatrix &matr_C,
90  int L_size, int pc_size, int alpha_size,
91  Matrix2D<double> &L_lowerBoundaryCondition,
92  Matrix2D<double> &L_upperBoundaryCondition,
93  Matrix2D<double> &pc_lowerBoundaryCondition,
94  Matrix2D<double> &pc_upperBoundaryCondition,
95  Matrix2D<double> &alpha_lowerBoundaryCondition,
96  Matrix2D<double> &alpha_upperBoundaryCondition,
97  string L_lowerBoundaryCondition_calculationType,
98  string L_upperBoundaryCondition_calculationType,
99  string pc_lowerBoundaryCondition_calculationType,
100  string pc_upperBoundaryCondition_calculationType,
101  string alpha_lowerBoundaryCondition_calculationType,
102  string alpha_upperBoundaryCondition_calculationType,
103  Matrix3D<double> &DLL,
104  Matrix3D<double> &Dpcpc, Matrix3D<double> &DpcpcLpp,
105  Matrix3D<double> &Daa, Matrix3D<double> &DaaLpp,
106  Matrix3D<double> &Dpca, Matrix3D<double> &DpcaLpp,
107  Matrix3D<double> &Jacobian, double dt, double Lpp,
108  double tau = 1e99, double tauLpp = 1e99);
109 
111  CalculationMatrix &matr_A, CalculationMatrix &matr_B, CalculationMatrix &matr_C,
113  int L_size, int pc_size, int alpha_size,
114  Matrix2D<double> &L_lowerBoundaryCondition,
115  Matrix2D<double> &L_upperBoundaryCondition,
116  Matrix2D<double> &pc_lowerBoundaryCondition,
117  Matrix2D<double> &pc_upperBoundaryCondition,
118  Matrix2D<double> &alpha_lowerBoundaryCondition,
119  Matrix2D<double> &alpha_upperBoundaryCondition,
120  string L_lowerBoundaryCondition_calculationType,
121  string L_upperBoundaryCondition_calculationType,
122  string pc_lowerBoundaryCondition_calculationType,
123  string pc_upperBoundaryCondition_calculationType,
124  string alpha_lowerBoundaryCondition_calculationType,
125  string alpha_upperBoundaryCondition_calculationType,
126  Matrix3D<double> &DLL,
127  Matrix3D<double> &Dpcpc, Matrix3D<double> &DpcpcLpp,
128  Matrix3D<double> &Daa, Matrix3D<double> &DaaLpp,
129  Matrix3D<double> &Dpca, Matrix3D<double> &DpcaLpp,
130  Matrix3D<double> &Jacobian, double dt, double Lpp,
131  double tau = 1e99, double tauLpp = 1e99);
132 
133 
135  int il, int im, int ia,
136  string FirstDerivative, string SecondDerivative,
138  Matrix3D<double> &Dxx, Matrix3D<double> &Jacobian,
139  double multiplicator);
140 
142  int il, int im, int ia,
143  string FirstDerivative, string SecondDerivative,
145  Matrix3D<double> &Dxx, Matrix3D<double> &Jacobian,
146  double multiplicator);
147 
149  int il, int im, int ia,
150  string FirstDerivative, string SecondDerivative,
152  Matrix3D<double> &Dxx, Matrix3D<double> &Jacobian,
153  double multiplicator);
154 
155 
156 #endif
157