VERB_code_2.3
|
Making model matrixes, solving model matrixes. More...
#include "MatrixSolver.h"
#include <math.h>
#include <malloc/malloc.h>
#include <string>
#include <ctime>
#include "../Parameters/Parameters.h"
#include "../VariousFunctions/variousFunctions.h"
#include "../Exceptions/error.h"
#include "../Lapack/Include/cpplapack.h"
#include "../GMRES/itlin.h"
#include <iostream>
#include "../GMRES2/GMRES2.h"
Go to the source code of this file.
Classes | |
class | Preconditioner_jacobi |
Used for calculating matrix inversion with GMRES2. More... | |
class | Preconditioner_SOR |
Used for calculating matrix inversion with GMRES2. More... | |
Macros | |
#define | EE 1.e-19 |
Error value. | |
#define | I(l, k) (l)*n+(k) |
Specified Index [l][k]. | |
Functions | |
void | matvec (int n, double *x, double *b) |
Used for gmres. More... | |
void | preconr (int n, double *x, double *b) |
double | gmres_norm2 (int n, double *v) |
void | gmres (int n, double *y, MATVEC *matvec, PRECON *preconr, PRECON *preconl, double *b, struct ITLIN_OPT *opt, struct ITLIN_INFO *info) |
bool | MakeMatrix (double *A, double *B1, double *C, double *R, double *x, double tau, double n, double f_lower, double f_upper, int nx, double dt, double *Dxx, double taulc, double alc, int g_flag, string lower_border_condition_type, string upper_border_condition_type, string approximationMethod) |
Make matrix for 1d diffusion. More... | |
bool | SolveMatrix (double *f, double *A, double *B1, double *C, double *R, double f_lower, double f_upper, int nx, double dt) |
Solver for 1d general equation. More... | |
void | AddBoundary (DiagMatrix &Matrix_A, string type, int in, int id1, double dh) |
bool | MakeModelMatrix_3D (CalculationMatrix &matr_A, CalculationMatrix &matr_B, CalculationMatrix &matr_C, Matrix3D< double > &L, Matrix3D< double > &pc, Matrix3D< double > &alpha, int L_size, int pc_size, int alpha_size, Matrix2D< double > &L_lowerBoundaryCondition, Matrix2D< double > &L_upperBoundaryCondition, Matrix2D< double > &pc_lowerBoundaryCondition, Matrix2D< double > &pc_upperBoundaryCondition, Matrix2D< double > &alpha_lowerBoundaryCondition, Matrix2D< double > &alpha_upperBoundaryCondition, string L_lowerBoundaryCondition_calculationType, string L_upperBoundaryCondition_calculationType, string pc_lowerBoundaryCondition_calculationType, string pc_upperBoundaryCondition_calculationType, string alpha_lowerBoundaryCondition_calculationType, string alpha_upperBoundaryCondition_calculationType, Matrix3D< double > &DLL, Matrix3D< double > &Dpcpc, Matrix3D< double > &DpcpcLpp, Matrix3D< double > &Daa, Matrix3D< double > &DaaLpp, Matrix3D< double > &Dpca, Matrix3D< double > &DpcaLpp, Matrix3D< double > &Jacobian, double dt, double Lpp, double tau, double tauLpp) |
void | jacobi (int m_size, double *z, double *w) |
void | SOR (int m_size, double *z, double *w) |
void | for_norm_A (int n, double *x, double *b) |
void | gmres_wrapout (CalculationMatrix &A, Matrix1D< double > &B, Matrix1D< double > &X, Parameters_structure::PSD::GMRES_parameters_structure GMRES_parameters) |
void | gmres2_wrapout (CalculationMatrix &A, Matrix1D< double > &B, Matrix1D< double > &X, Parameters_structure::PSD::GMRES_parameters_structure GMRES_parameters) |
Used for calculating matrix inversion with GMRES2. | |
void | Lapack (DiagMatrix &A, Matrix1D< double > &B, Matrix1D< double > &X) |
void | over_relaxation_diag (DiagMatrix &A, Matrix1D< double > &B, Matrix1D< double > &X, int max_steps, double EE) |
void | GetDerivativeVector (string derivativeType, int &dL, int &dPc, int &dAlpha) |
void | SecondDerivativeApproximation_3D (CalculationMatrix &matr_A, int il, int im, int ia, string FirstDerivative, string SecondDerivative, Matrix3D< double > &L, Matrix3D< double > &pc, Matrix3D< double > &alpha, Matrix3D< double > &Dxx, Matrix3D< double > &Jacobian, double multiplicator) |
double | max2 (double a, double b) |
Function returns maximum between a and b. | |
void | mult_vect2 (double *af, double *vf, double *wf, int nf) |
Multiplication between vector and matrix. | |
void | gauss (DiagMatrix &A, Matrix1D< double > &B, Matrix1D< double > &X, int n) |
Gauss method - added to move out of if statements. | |
void | gauss_solve (double *a, double *b, int n) |
Gauss inversion. | |
bool | tridag (double a[], double b[], double c[], double r[], double u[], long n) |
Solve the AU=R system of equations, where A - tridiagonal matrix nxn with diagonals a[], b[], c[]. More... | |
bool | MakeModelMatrix_3D_KC (CalculationMatrix &matr_A, CalculationMatrix &matr_B, CalculationMatrix &matr_C, Matrix3D< double > &L, Matrix3D< double > &pc, Matrix3D< double > &alpha, int L_size, int pc_size, int alpha_size, Matrix2D< double > &L_lowerBoundaryCondition, Matrix2D< double > &L_upperBoundaryCondition, Matrix2D< double > &pc_lowerBoundaryCondition, Matrix2D< double > &pc_upperBoundaryCondition, Matrix2D< double > &alpha_lowerBoundaryCondition, Matrix2D< double > &alpha_upperBoundaryCondition, string L_lowerBoundaryCondition_calculationType, string L_upperBoundaryCondition_calculationType, string pc_lowerBoundaryCondition_calculationType, string pc_upperBoundaryCondition_calculationType, string alpha_lowerBoundaryCondition_calculationType, string alpha_upperBoundaryCondition_calculationType, Matrix3D< double > &DLL, Matrix3D< double > &Dpcpc, Matrix3D< double > &DpcpcLpp, Matrix3D< double > &Daa, Matrix3D< double > &DaaLpp, Matrix3D< double > &Dpca, Matrix3D< double > &DpcaLpp, Matrix3D< double > &Jacobian, double dt, double Lpp, double tau, double tauLpp) |
void | SecondDerivativeApproximation_3D_KC_Diagonal (CalculationMatrix &matr_A, int il, int im, int ia, string FirstDerivative, string SecondDerivative, Matrix3D< double > &L, Matrix3D< double > &pc, Matrix3D< double > &alpha, Matrix3D< double > &Dxx, Matrix3D< double > &Jacobian, double multiplicator) |
void | SecondDerivativeApproximation_3D_KC_Mixed (CalculationMatrix &matr_A, int il, int im, int ia, string FirstDerivative, string SecondDerivative, Matrix3D< double > &L, Matrix3D< double > &pc, Matrix3D< double > &alpha, Matrix3D< double > &Dxx, Matrix3D< double > &Jacobian, double multiplicator) |
Variables | |
static CalculationMatrix | A_copy |
Making model matrixes, solving model matrixes.
Matrix form of linear equations: A*X[t+1] = B*X[t], where A - model matrix, B - RHS, X[t] - known values of function (PSD), X[t+1] - unknown values of function.
In that file there are procedures for making model matrix for 1d-diffusion, 2d-diffusion, some ideas of 3d-diffusion and mixed terms. Solver for tridiagonal matrix, solver by gauss method and iteration method (upper relaxation).
This file is under development, and has a lot of commented code, unfinished etc code.
There is a checked version for 1d diffusion (or split method of 2d, 3d diffusions) - 1d_universal_solver.
Definition in file MatrixSolver.cpp.
void AddBoundary | ( | DiagMatrix & | Matrix_A, |
string | type, | ||
int | in, | ||
int | id1, | ||
double | dh | ||
) |
Supportive sub-function to add boundary conditions to model matrix
&Matrix_A | - model matrix to be updated |
type | - type of boundary condition |
in | - second index of matrix |
id1 | - upper boundary index |
dh | - derivative value? |
Definition at line 317 of file MatrixSolver.cpp.
void for_norm_A | ( | int | n, |
double * | x, | ||
double * | b | ||
) |
Simple preconditioner
Definition at line 731 of file MatrixSolver.cpp.
void GetDerivativeVector | ( | string | derivativeType, |
int & | dL, | ||
int & | dPc, | ||
int & | dAlpha | ||
) |
Get change in indexes according to the derivatives direction.
Used in approximation. If it gets derivativeType == DT_ALPHA_left, for example, that means we have alpha-left-derivative, which is ( f(alpha) - f(alpha-1) ) / ( delta alpha ). So it returns dAlpha = -1, as a derivative direction vector.
Definition at line 1185 of file MatrixSolver.cpp.
void gmres_wrapout | ( | CalculationMatrix & | A, |
Matrix1D< double > & | B, | ||
Matrix1D< double > & | X, | ||
Parameters_structure::PSD::GMRES_parameters_structure | GMRES_parameters | ||
) |
GMRES inversion. Potentially with bugs.
A * X = B - equation
Definition at line 745 of file MatrixSolver.cpp.
void jacobi | ( | int | m_size, |
double * | z, | ||
double * | w | ||
) |
Jacobi (or diagonal) preconditioner:
Definition at line 689 of file MatrixSolver.cpp.
void Lapack | ( | DiagMatrix & | A, |
Matrix1D< double > & | B, | ||
Matrix1D< double > & | X | ||
) |
bool MakeMatrix | ( | double * | A, |
double * | B1, | ||
double * | C, | ||
double * | R, | ||
double * | x, | ||
double | tau, | ||
double | n, | ||
double | f_lower, | ||
double | f_upper, | ||
int | nx, | ||
double | dt, | ||
double * | Dxx, | ||
double | taulc, | ||
double | alc, | ||
int | g_flag, | ||
string | lower_border_condition_type, | ||
string | upper_border_condition_type, | ||
string | approximationMethod | ||
) |
Make matrix for 1d diffusion.
Make matrix for 1D diffusion.
AM_Split_LR
\( Coef1 = \frac{x^n}{\hbar} \)
\( coef\_05m = 0.5 * (Dxx * x^{-n} * G(x)) + Dxx_{j-1} * G(x_{j-1}) * {x_{j-1}}^{-n} \)
\( coef\_05p = 0.5 * (Dxx * x^{-n} * G(x)) + Dxx_{j+1} * G(x_{j+1}) * {x_{j+1}}^{-n} \) where VF::G() is used if g_flag is set
\( A = \frac{Coef1 * coef\_05m}{h} \)
\( C = \frac{Coef1 * coef\_05p}{h+1} \)
AM_Split_C
\( Coef1 = \frac{x^n}{\hbar} \)
\( Dm\_05 = 0.5 * (Dxx_j + Dxx_{j-1}) \) \( Dp\_05 = 0.5 * (Dxx_j + Dxx_{j+1}) \)
\( A = Coef1 * \frac{Dm\_05 * G(\frac{x-h}{2}) * (\frac{x-h}{2})^{-n} }{h} \)
\( C = Coef1 * \frac{Dp\_05 * G(\frac{x-(h+1)}{2}) * (\frac{x-(h+1)}{2})^{-n} }{h_{j+1}} \)
Definition at line 72 of file MatrixSolver.cpp.
bool MakeModelMatrix_3D | ( | CalculationMatrix & | matr_A, |
CalculationMatrix & | matr_B, | ||
CalculationMatrix & | matr_C, | ||
Matrix3D< double > & | L, | ||
Matrix3D< double > & | pc, | ||
Matrix3D< double > & | alpha, | ||
int | L_size, | ||
int | pc_size, | ||
int | alpha_size, | ||
Matrix2D< double > & | L_lowerBoundaryCondition, | ||
Matrix2D< double > & | L_upperBoundaryCondition, | ||
Matrix2D< double > & | pc_lowerBoundaryCondition, | ||
Matrix2D< double > & | pc_upperBoundaryCondition, | ||
Matrix2D< double > & | alpha_lowerBoundaryCondition, | ||
Matrix2D< double > & | alpha_upperBoundaryCondition, | ||
string | L_lowerBoundaryCondition_calculationType, | ||
string | L_upperBoundaryCondition_calculationType, | ||
string | pc_lowerBoundaryCondition_calculationType, | ||
string | pc_upperBoundaryCondition_calculationType, | ||
string | alpha_lowerBoundaryCondition_calculationType, | ||
string | alpha_upperBoundaryCondition_calculationType, | ||
Matrix3D< double > & | DLL, | ||
Matrix3D< double > & | Dpcpc, | ||
Matrix3D< double > & | DpcpcLpp, | ||
Matrix3D< double > & | Daa, | ||
Matrix3D< double > & | DaaLpp, | ||
Matrix3D< double > & | Dpca, | ||
Matrix3D< double > & | DpcaLpp, | ||
Matrix3D< double > & | Jacobian, | ||
double | dt, | ||
double | Lpp, | ||
double | tau, | ||
double | tauLpp | ||
) |
Create matrix form of the Fokker-Planck equation: matr_A * PSD(t+1) = matr_B * PSD(t) + matr_C
Calculation matr_A, matr_B, and matr_C, i.e. numerical approximation of the derivatives Takes boundary conditions and diffusion coefficients as an input. L, pc, alpha should be orthogonal for 3D!!!
&matr_A | - matrix A |
&matr_B | - matrix B |
&matr_C | - matrix C |
&L | - matrix for L |
&pc | - matrix for pc |
&alpha | - matrix for alpha |
L_size | - size of L |
pc_size | - size of pc |
alpha_size | - size of alpha |
&L_lowerBoundaryCondition | - lower boundary condition for L |
&L_upperBoundaryCondition | - upper boundary condition for L |
&pc_lowerBoundaryCondition | - lower boundary condition for pc |
&pc_upperBoundaryCondition | - upper boundary condition for pc |
&alpha_lowerBoundaryCondition | - lower boundary condition for alpha |
&alpha_upperBoundaryCondition | - upper boundary condition for alpha |
L_lowerBoundaryCondition_calculationType | - calculation type of L's lower boundary condition |
L_upperBoundaryCondition_calculationType | - calculation type of L's upper boundary condition |
pc_lowerBoundaryCondition_calculationType | - calculation type of pc's lower boundary condition |
pc_upperBoundaryCondition_calculationType | - calculation type of pc's upper boundary condition |
alpha_lowerBoundaryCondition_calculationType | - calculation type of alpha's lower boundary condition |
alpha_upperBoundaryCondition_calculationType | - calculation type of alpha's upper boundary condition |
&DLL | - matrix of diffusion coefficients for L |
&Dpcpc | - diffusion coefficients group for Dpcpc |
&DpcpcLpp | - diffusion coefficients group for DpcpcLpp |
&Daa | - diffusion coefficients group for Daa |
&DaaLpp | - diffusion coefficients group for DaaLpp |
&Dpca | - diffusion coefficients group for Dpca |
&DpcaLpp | - diffusion coefficients group for DpcaLpp |
&Jacobian | - Jacobian values |
dt | - diffusion time(?) |
Lpp | - plasma pause location |
tau | - lifetime outside of plasmasphere |
tauLpp | - lifetime inside of plasmasphere |
Definition at line 371 of file MatrixSolver.cpp.
bool MakeModelMatrix_3D_KC | ( | CalculationMatrix & | matr_A, |
CalculationMatrix & | matr_B, | ||
CalculationMatrix & | matr_C, | ||
Matrix3D< double > & | L, | ||
Matrix3D< double > & | pc, | ||
Matrix3D< double > & | alpha, | ||
int | L_size, | ||
int | pc_size, | ||
int | alpha_size, | ||
Matrix2D< double > & | L_lowerBoundaryCondition, | ||
Matrix2D< double > & | L_upperBoundaryCondition, | ||
Matrix2D< double > & | pc_lowerBoundaryCondition, | ||
Matrix2D< double > & | pc_upperBoundaryCondition, | ||
Matrix2D< double > & | alpha_lowerBoundaryCondition, | ||
Matrix2D< double > & | alpha_upperBoundaryCondition, | ||
string | L_lowerBoundaryCondition_calculationType, | ||
string | L_upperBoundaryCondition_calculationType, | ||
string | pc_lowerBoundaryCondition_calculationType, | ||
string | pc_upperBoundaryCondition_calculationType, | ||
string | alpha_lowerBoundaryCondition_calculationType, | ||
string | alpha_upperBoundaryCondition_calculationType, | ||
Matrix3D< double > & | DLL, | ||
Matrix3D< double > & | Dpcpc, | ||
Matrix3D< double > & | DpcpcLpp, | ||
Matrix3D< double > & | Daa, | ||
Matrix3D< double > & | DaaLpp, | ||
Matrix3D< double > & | Dpca, | ||
Matrix3D< double > & | DpcaLpp, | ||
Matrix3D< double > & | Jacobian, | ||
double | dt, | ||
double | Lpp, | ||
double | tau, | ||
double | tauLpp | ||
) |
kckim test Create matrix form of the Fokker-Planck equation: matr_A * PSD(t+1) = matr_B * PSD(t) + matr_C
Calculation matr_A, matr_B, and matr_C, i.e. numerical approximation of the derivatives Takes boundary conditions and diffusion coefficients as an input. L, pc, alpha should be orthogonal for 3D!!!
Definition at line 1467 of file MatrixSolver.cpp.
void matvec | ( | int | m_size, |
double * | x, | ||
double * | b | ||
) |
Used for gmres.
Vector to matrix multiplication for GMRES
Definition at line 647 of file MatrixSolver.cpp.
void over_relaxation_diag | ( | DiagMatrix & | A, |
Matrix1D< double > & | B, | ||
Matrix1D< double > & | X, | ||
int | max_steps, | ||
double | EE | ||
) |
Over relaxation iteration method. Uses 0th diagonal row for A.
&A | - diagonal matrix A |
&B | - array B |
&X | - output matrix |
max_steps | - maximum steps of iteration |
EE | - error value |
Definition at line 1132 of file MatrixSolver.cpp.
void SecondDerivativeApproximation_3D | ( | CalculationMatrix & | matr_A, |
int | il, | ||
int | im, | ||
int | ia, | ||
string | FirstDerivative, | ||
string | SecondDerivative, | ||
Matrix3D< double > & | L, | ||
Matrix3D< double > & | pc, | ||
Matrix3D< double > & | alpha, | ||
Matrix3D< double > & | Dxx, | ||
Matrix3D< double > & | Jacobian, | ||
double | multiplicator | ||
) |
Second derivative approximation, returns coefficients to be putted into the model matrix. \(L_{\alpha \beta}(y) = (D_{\alpha \beta} \cdot y_{\bar{x}_\alpha})_{x_{\beta}}\)
Samarskiy, page 261
Returns coefficients to be put into model matrix for an approximation of a second derivative.
Definition at line 1212 of file MatrixSolver.cpp.
void SecondDerivativeApproximation_3D_KC_Diagonal | ( | CalculationMatrix & | matr_A, |
int | il, | ||
int | im, | ||
int | ia, | ||
string | FirstDerivative, | ||
string | SecondDerivative, | ||
Matrix3D< double > & | L, | ||
Matrix3D< double > & | pc, | ||
Matrix3D< double > & | alpha, | ||
Matrix3D< double > & | Dxx, | ||
Matrix3D< double > & | Jacobian, | ||
double | multiplicator | ||
) |
kckim test Second derivative approximation, returns coefficients to be put into the model matrix. \(L_{\alpha \beta}(y) = (D_{\alpha \beta} \cdot y_{\bar{x}_\alpha})_{x_{\beta}}\)
Samarskiy, page 261
Returns coefficients to be put into model matrix for an approximation of a second derivative.
Definition at line 1825 of file MatrixSolver.cpp.
void SecondDerivativeApproximation_3D_KC_Mixed | ( | CalculationMatrix & | matr_A, |
int | il, | ||
int | im, | ||
int | ia, | ||
string | FirstDerivative, | ||
string | SecondDerivative, | ||
Matrix3D< double > & | L, | ||
Matrix3D< double > & | pc, | ||
Matrix3D< double > & | alpha, | ||
Matrix3D< double > & | Dxx, | ||
Matrix3D< double > & | Jacobian, | ||
double | multiplicator | ||
) |
short form mixed term approximation takes approximations using alpha and pc left and right side data for derivatives
&matr_A | - calculation matrix |
il | - iterator for L |
im | - iterator for pc |
ia | - iterator for alpha |
FirstDerivative | - key term for first derivative approximation |
SecondDerivative | - key term for second derivative approximation |
&L | - matrix for L |
&pc | - matrix for pc |
&alpha | - matrix for alpha |
&Dxx | - matrix of derivatives |
&Jacobian | - jacobian |
multiplicator | - positive or negative multiplicator |
Definition at line 1904 of file MatrixSolver.cpp.
bool SolveMatrix | ( | double * | f, |
double * | A, | ||
double * | B1, | ||
double * | C, | ||
double * | R, | ||
double | f_lower, | ||
double | f_upper, | ||
int | nx, | ||
double | dt | ||
) |
Solver for 1d general equation.
Solver for 1d-diffusion matrix (tridiagonal).
Definition at line 271 of file MatrixSolver.cpp.
void SOR | ( | int | m_size, |
double * | z, | ||
double * | w | ||
) |
SOR preconditioner:
P = (D/rel + L) w = P^-1 * z; P * w = z (D/rel + L) * w = z
Definition at line 709 of file MatrixSolver.cpp.
bool tridag | ( | double | a[], |
double | b[], | ||
double | c[], | ||
double | r[], | ||
double | u[], | ||
long | n | ||
) |
Solve the AU=R system of equations, where A - tridiagonal matrix nxn with diagonals a[], b[], c[].
Solver for a system of equations with 3-diagonal matrix.
Au = r Where A = diag(a, b, c),
a[] | - array, diagonal '-1' of the matrix |
b[] | - array, diagonal '0' of the matrix |
c[] | - array, diagonal '+1' of the matrix |
r[] | - array, r-vector |
u[] | - array, result |
n | - size of the matrix |
Definition at line 1424 of file MatrixSolver.cpp.