VERB_code_2.3
MatrixSolver.h File Reference
#include "../Grid/Grid.h"
#include "../DiffusionCoefficient/DiffusionCoefficient.h"
#include "../GMRES/itlin.h"

Go to the source code of this file.

Functions

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="AM_Split_C")
 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 over_relaxation_diag (DiagMatrix &A, Matrix1D< double > &B, Matrix1D< double > &X, int max_steps=1e5, double EE=1e-3)
 
void Lapack (DiagMatrix &A, Matrix1D< double > &B, Matrix1D< double > &X)
 
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...
 
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.
 
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=1e99, double tauLpp=1e99)
 
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=1e99, double tauLpp=1e99)
 
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)
 
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)
 

Detailed Description

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.

Author
Developed under supervision of the PI Yuri Shprits

Definition in file MatrixSolver.h.

Function Documentation

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 Lapack ( DiagMatrix A,
Matrix1D< double > &  B,
Matrix1D< double > &  X 
)

Lapack inversion.

A * X = B - equation

Definition at line 1032 of file MatrixSolver.cpp.

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!!!

Parameters
&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 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.

Parameters
&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

Note
SecondDerivativeApproximation_3D_Mixed is deprecated in favor of this function
Parameters
&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.

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),

Parameters
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.