VERB_code_2.2  2
 All Classes Namespaces Files Functions Variables Typedefs Macros Pages
MatrixSolver.cpp
Go to the documentation of this file.
1 /**
2  * \file MatrixSolver.cpp
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 //#define CPPL_DEBUG
20 //#define CPPL_VERBOSE
21 
22 // Memory leaks
23 #if defined(_WIN32) || defined(_WIN64)
24 #define _CRTDBG_MAP_ALLOC
25 #include <crtdbg.h>
26 #endif
27 
28 #include "MatrixSolver.h"
29 #include <math.h>
30 #include <malloc.h>
31 #include <string>
32 #include <ctime>
33 
34 // Lapack debugging
35 //#define CPPL_VERBOSE
36 //#define CPPL_DEBUG
37 
38 using namespace std;
39 
40 // parameters class
41 #include "../Parameters/Parameters.h"
42 
43 //#include "PSD.h"
44 #include "../VariousFunctions/variousFunctions.h"
45 #include "../Exceptions/error.h"
46 
47 // Lapack
48 #include "../Lapack/Include/cpplapack.h"
49 
50 // gmres
51 extern "C" {
52  #include "../GMRES/itlin.h"
53 }
54 
55 void matvec(int n, double *x, double *b);
56 void preconr(int n, double *x, double *b);
57 extern "C" double gmres_norm2(int n, double *v);
58 
59 extern "C" void gmres(int n, double *y, MATVEC *matvec,
60  PRECON *preconr, PRECON *preconl, double *b,
61  struct ITLIN_OPT *opt, struct ITLIN_INFO *info);
62 
63 #include <iostream>
64 
65 #include "../GMRES2/GMRES2.h"
66 
67 
68 /**
69  * Make matrix for 1D diffusion.
70  */
71 bool MakeMatrix(double *A,
72  double *B1,
73  double *C,
74  double *R,
75  double *x, // one dimentional grid
76  double tau, // lifetime
77  double n, // power of x
78  double f_lower, // lower boundary conditions
79  double f_upper, // upper boundary conditions
80  int nx, // number of grid points
81  double dt, // time step
82  double *Dxx, // diffusion coefficient
83  double taulc, // lifetime in loss cone
84  double alc, // loss cone size
85  int g_flag, // diffusion flag
86  string lower_border_condition_type,// type of lower boundary conditions: 0 - for value, 1 - for derivative
87  string upper_border_condition_type,// type of upper boundary conditions: 0 - for value, 1 - for derivative
88  string approximationMethod)
89 {
90 
91  // local variables
92  double Coef1, Dm_05, Dp_05, tau_cur;
93 
94  // initialisation of temp matrix_arrays
95 // double *A = new double[nx];
96  double *B = new double[nx];
97 // double *C = new double[nx];
98 // double *B1 = new double[nx];
99  double *h = new double[nx];
100  double *h_bar = new double[nx];
101 
102  // iterators i, j
103  int i, j;
104  double coef_05m, coef_05p;
105 
106  // making matrix
107  for (j = 1; j <= nx-1; j++) {
108  h[j] = x[j] - x[j-1];
109  }
110  for (j = 0; j <= nx-2; j++) {
111  h_bar[j] = (h[j] + h[j+1])/2;
112  }
113 
114  for (j = 1; j <= nx-2; j++) {
115  // making model matrix, according to the approximation method
116  if (approximationMethod == "AM_Split_LR") {
117  //Ver1
118  Coef1 = pow(x[j], n)/h_bar[j];
119  if (g_flag == 1) {
120  Coef1 = Coef1/VF::G(x[j]);
121  coef_05m = 0.5*(Dxx[j]*VF::G(x[j])*pow(x[j], -n) + Dxx[j-1]*VF::G(x[j-1])*pow(x[j-1], -n));
122  coef_05p = 0.5*(Dxx[j]*VF::G(x[j])*pow(x[j], -n) + Dxx[j+1]*VF::G(x[j+1])*pow(x[j+1], -n));
123  } else {
124  coef_05m = 0.5*(Dxx[j]*pow(x[j], -n) + Dxx[j-1]*pow(x[j-1], -n));
125  coef_05p = 0.5*(Dxx[j]*pow(x[j], -n) + Dxx[j+1]*pow(x[j+1], -n));
126  }
127  A[j] = Coef1 * coef_05m / h[j];
128  C[j] = Coef1 * coef_05p / h[j+1];
129  // end ver1
130  } else if (approximationMethod == "AM_Split_LR_oldway") {
131  //Ver1
132  Coef1 = pow(x[j], n);
133  if (g_flag == 1) {
134  Coef1 = Coef1/VF::G(x[j]);
135  coef_05m = 0.5*(Dxx[j]*VF::G(x[j])*pow(x[j], -n) + Dxx[j-1]*VF::G(x[j-1])*pow(x[j-1], -n));
136  coef_05p = 0.5*(Dxx[j]*VF::G(x[j])*pow(x[j], -n) + Dxx[j+1]*VF::G(x[j+1])*pow(x[j+1], -n));
137  } else {
138  coef_05m = 0.5*(Dxx[j]*pow(x[j], -n) + Dxx[j-1]*pow(x[j-1], -n));
139  coef_05p = 0.5*(Dxx[j]*pow(x[j], -n) + Dxx[j+1]*pow(x[j+1], -n));
140  }
141  A[j] = Coef1 * coef_05m / h[j] / h[j+1];
142  C[j] = Coef1 * coef_05p / h[j+1] / h[j];
143  // end ver1
144  } else if (approximationMethod == "AM_Split_C") {
145  //Ver2
146  Coef1 = pow(x[j], n)/h_bar[j];
147  if (g_flag == 1) {
148  Coef1 = Coef1/VF::G(x[j]);
149  }
150  Dm_05 = 0.5*(Dxx[j] + Dxx[j-1]);
151  Dp_05 = 0.5*(Dxx[j] + Dxx[j+1]);
152  if (g_flag == 0) {
153  A[j] = Coef1 * (Dm_05 * pow(x[j] - h[j]/2, -n) / h[j]);
154  C[j] = Coef1 * (Dp_05 * pow(x[j] + h[j+1]/2, -n) / h[j+1]);
155  } else {
156  A[j] = Coef1 * (Dm_05 * VF::G(x[j] - h[j]/2) * pow(x[j] - h[j]/2, -n) / h[j]);
157  C[j] = Coef1 * (Dp_05 * VF::G(x[j] + h[j+1]/2) * pow(x[j] + h[j+1]/2, -n) / h[j+1]);
158  }
159  // end ver2
160  } else if (approximationMethod == "AM_Split_C_oldway") {
161  //Ver2
162  Coef1 = pow(x[j], n)/h_bar[j];
163  if (g_flag == 1) {
164  Coef1 = Coef1/VF::G(x[j]);
165  }
166  Dm_05 = 0.5*(Dxx[j] + Dxx[j-1]);
167  Dp_05 = 0.5*(Dxx[j] + Dxx[j+1]);
168  if (g_flag == 0) {
169  A[j] = Coef1 * (Dm_05 * pow(x[j] - h[j]/2, -n) / h[j]);
170  C[j] = Coef1 * (Dp_05 * pow(x[j] + h[j]/2, -n) / h[j+1]);
171  } else {
172  A[j] = Coef1 * (Dm_05 * VF::G(x[j] - h[j]/2) * pow(x[j] - h[j]/2, -n) / h[j]);
173  C[j] = Coef1 * (Dp_05 * VF::G(x[j] + h[j]/2) * pow(x[j] + h[j]/2, -n) / h[j+1]);
174  }
175  // end ver2
176  } else {
177  throw error_msg("APPRERR", "Unknown approximation type: %d", approximationMethod.c_str());
178  }
179 
180  B[j] = -A[j] - C[j];
181 
182  //if (x[j] < (alc * pi / 180) && (g_flag == 1)) {
183  // Loss cone close to 0 and 180 deg pitch-angle
184  if ((x[j] < alc || x[j] > VC::pi - alc) && (g_flag == 1)) {
185  tau_cur = taulc;
186  } else {
187  tau_cur = tau;
188  }
189 
190  B1[j] = B[j] - (1.0/dt + 1/tau_cur);
191  }
192 
193  for (i = 1; i <= nx-2; i++) {
194  R[i] = - 1.0 / dt;
195  }
196 
197  // Using boundary conditions.
198  // lower boundary
199  i = 0;
200  R[i] = f_lower;
201  if (lower_border_condition_type == "BCT_CONSTANT_VALUE") { // for condition on value
202  A[i] = 0; B1[i] = 1; C[i] = 0;
203  } else if (lower_border_condition_type == "BCT_CONSTANT_DERIVATIVE") { // for condition on derivative
204  A[i] = 0; B1[i] = -1; C[i] = 1;
205  } else {
206  throw error_msg("1D_DIFF_BOUNDARY", "1d_universal_solver: unknown lower boundary type: %d", lower_border_condition_type.c_str());
207  //log1::fcout << "1d_universal_solver: unknown lower boundary type: " << lower_border_condition_type << "\n";
208  //return -1;
209  }
210 
211  // upper boundary
212  i = nx-1;
213  R[i] = f_upper;
214  if (upper_border_condition_type == "BCT_CONSTANT_VALUE") {
215  A[i] = 0; B1[i] = 1; C[i] = 0;
216  } else if (upper_border_condition_type == "BCT_CONSTANT_DERIVATIVE") {
217  A[i] = -1; B1[i] = 1; C[i] = 0;
218  } else {
219  throw error_msg("1D_DIFF_BOUNDARY", "1d_universal_solver: unknown lower boundary type: %d", upper_border_condition_type.c_str());
220  //log1::fcout << "1d_universal_solver: unknown upper boundary type: " << upper_border_condition_type << "\n";
221  //return -1;
222  }
223 
224  // solving matrix
225 // tridag(A, B1, C, R, f, nx);
226 
227  // deleting matrix_arrays
228 // delete A;
229  delete B;
230 // delete C;
231 // delete B1;
232  delete h;
233  delete h_bar;
234 
235  return true;
236 }
237 
238 /**
239  * Solver for 1d-diffusion matrix (tridiagonal).
240  */
241 bool SolveMatrix(double *f, // phase space density function
242  double *A,
243  double *B1,
244  double *C,
245  double *R,
246  double f_lower, // lower boundary conditions
247  double f_upper, // upper boundary conditions
248  int nx, // number of grid points
249  double dt) // time step
250 {
251  // initialisation of temp matrix_arrays
252  double *R1 = new double[nx];
253 
254  // iterators i, j
255  int i;
256 
257  for (i = 1; i <= nx-2; i++) {
258  R1[i] = R[i] * f[i];
259  }
260 
261  // Using borders conditions.
262  // lower boundary
263  i = 0;
264  R1[i] = R[i];
265 
266  // upper boundary
267  i = nx-1;
268  R1[i] = R[i];
269 
270  // solving matrix
271  tridag(A, B1, C, R1, f, nx);
272 
273  // deleting matrix_arrays
274  delete R1;
275 
276  return true;
277 }
278 
279 
280 /** Supportive sub-function to add boundary conditions to model matrix
281  */
282 void AddBoundary(DiagMatrix &Matrix_A, string type, int in, int id1, double dh) {
283  int id0 = 0;
284  if (type == "BCT_CONSTANT_VALUE") { // for condition on value
285  Matrix_A[id0][in] = 1;
286  Matrix_A[id1][in] = 0;
287  } else if (type == "BCT_CONSTANT_DERIVATIVE") { // for condition on derivative
288  Matrix_A[id0][in] = -1/dh;
289  Matrix_A[id1][in] = 1/dh;
290  } else {
291  throw error_msg("2D_DIFF_BOUNDARY", "1d_universal_solver: unknown boundary type: %d", type.c_str());
292  }
293 }
294 
295 /**
296  * Create matrix form of the Fokker-Planck equation: matr_A * PSD(t+1) = matr_B * PSD(t) + matr_C
297  *
298  * Calculation matr_A, matr_B, and matr_C, i.e. numerical approximation of the derivatives
299  * Takes boundary conditions and diffusion coefficients as an input.
300  * L, pc, alpha should be orthogonal for 3D!!!
301  */
303  CalculationMatrix &matr_A, CalculationMatrix &matr_B, CalculationMatrix &matr_C,
305  int L_size, int pc_size, int alpha_size,
306  Matrix2D<double> &L_lowerBoundaryCondition,
307  Matrix2D<double> &L_upperBoundaryCondition,
308  Matrix2D<double> &pc_lowerBoundaryCondition,
309  Matrix2D<double> &pc_upperBoundaryCondition,
310  Matrix2D<double> &alpha_lowerBoundaryCondition,
311  Matrix2D<double> &alpha_upperBoundaryCondition,
312  string L_lowerBoundaryCondition_calculationType,
313  string L_upperBoundaryCondition_calculationType,
314  string pc_lowerBoundaryCondition_calculationType,
315  string pc_upperBoundaryCondition_calculationType,
316  string alpha_lowerBoundaryCondition_calculationType,
317  string alpha_upperBoundaryCondition_calculationType,
318  Matrix3D<double> &DLL,
319  Matrix3D<double> &Dpcpc, Matrix3D<double> &DpcpcLpp,
320  Matrix3D<double> &Daa, Matrix3D<double> &DaaLpp,
321  Matrix3D<double> &Dpca, Matrix3D<double> &DpcaLpp,
322  Matrix3D<double> &Jacobian, double dt, double Lpp,
323  double tau, double tauLpp) {
324 
325  // variables, for future use
326  // il - L-index, im - pc-index, ia - alpha-index,
327  // id - index of a diagonal of the model matrix (counting from main diagonal, main is zero),
328  // in - model matrix line number
329  int il, im, ia, id, in;
330 
331  // create A-matrix, B-matrix, C-matrix
332  // Output::echo("Recalculating model matrix... ");
333 
334  // make everything zero
335  // TODO: make zero not everything, but existing diagonals
336  /* !!! matr_C[0] = 0;
337  for (il = -1; il <= 1; il++) {
338  for (im = -1; im <= 1; im++) {
339  for (ia = -1; ia <= 1; ia++) {
340  // calculating a diagonal number (id)
341  id = matr_A.index1d(ia, im, il);
342  // make the diagonal zero
343  matr_A[id] = 0;
344  matr_B[id] = 0;
345  }
346  }
347  }*/
348 
349  // Make diagonals to be equal to zero
350  DiagMatrix::iterator it;
351  for (it = matr_A.begin(); it != matr_A.end(); it++) it->second = 0;
352  for (it = matr_B.begin(); it != matr_B.end(); it++) it->second = 0;
353  for (it = matr_C.begin(); it != matr_C.end(); it++) it->second = 0;
354 
355  // create a new model matrix
356  // (f^{t+1} - f^{t})/dt = L1(f^{t+1}) + L2(f^{t+1}) + L3(f^{t+1})[main equation, losses are calculated separately]
357  // L1, L2, L3 - diffusion operators
358  // + need to add there the equations for boundary conditions also
359  double dh;
360  for (il = 0; il < L_size; il++) {
361  for (im = 0; im < pc_size; im++) {
362  for (ia = 0; ia < alpha_size; ia++) {
363  // calculating current line number (in)
364  in = matr_A.index1d(ia, im, il);
365 
366  // Describe all boundaries first. The 'inner' area will be approximated later.
367  // If L_size > 1 we should use L-boundaries, right?
368  // If not, then just calculate diffusion
369  if (il == 0 && L_size >= 3) {
370  // Boundary conditions
371  // Points at the corners are used only at one of the boundaries
372  // L-bottom boundary
373  matr_C[0][in] = L_lowerBoundaryCondition[im][ia];
374  id = matr_A.index1d(ia, im, il+1) - in;
375  dh = L[il+1][im][ia] - L[il][im][ia];
376  AddBoundary(matr_A, L_lowerBoundaryCondition_calculationType, in, id, dh);
377 
378  } else if (il == L_size-1 && L_size >= 3) {
379  // L-top boundary
380  matr_C[0][in] = L_upperBoundaryCondition[im][ia];
381  id = matr_A.index1d(ia, im, il-1) - in;
382  dh = L[il][im][ia] - L[il-1][im][ia];
383  AddBoundary(matr_A, L_upperBoundaryCondition_calculationType, in, id, dh);
384 
385  } else if (im == 0 && pc_size >= 3) {
386  // pc-low boundary
387  matr_C[0][in] = pc_lowerBoundaryCondition[il][ia];
388  id = matr_A.index1d(ia, im+1, il) - in;
389  dh = pc[il][im+1][ia] - pc[il][im][ia];
390  AddBoundary(matr_A, pc_lowerBoundaryCondition_calculationType, in, id, dh);
391 
392  } else if (im == pc_size-1 && pc_size >= 3) {
393  // pc-high boundary
394  matr_C[0][in] = pc_upperBoundaryCondition[il][ia];
395  id = matr_A.index1d(ia, im-1, il) - in;
396  dh = pc[il][im][ia] - pc[il][im-1][ia];
397  AddBoundary(matr_A, pc_upperBoundaryCondition_calculationType, in, id, dh);
398  } else if (ia == 0 && alpha_size >= 3) {
399  // alpha-loss cone boundary
400  matr_C[0][in] = alpha_lowerBoundaryCondition[il][im];
401  id = matr_A.index1d(ia+1, im, il) - in;
402  dh = alpha[il][im][ia+1] - alpha[il][im][ia];
403  AddBoundary(matr_A, alpha_lowerBoundaryCondition_calculationType, in, id, dh);
404 
405  } else if (ia == alpha_size-1 && alpha_size >= 3) {
406  // alpha-equatorial boundary
407  matr_C[0][in] = alpha_upperBoundaryCondition[il][im];
408  id = matr_A.index1d(ia-1, im, il) - in;
409  dh = alpha[il][im][ia] - alpha[il][im][ia-1];
410  AddBoundary(matr_A, alpha_upperBoundaryCondition_calculationType, in, id, dh);
411 
412  } else {
413 
414  // now we are sure we are not on a boundary, can do the Fokker-Planck equation approximation in the inner area
415 
416  // f^{t+1}/dt
417  matr_A[0][in] += 1.0/dt;
418  // f^{t}/dt
419  matr_B[0][in] += 1.0/dt;
420 
421  // calculate radial diffusion, if L_size >= 3
422  if (L_size >=3) {
423  // L term, add approximated derivatives into the calculation matrix
424  SecondDerivativeApproximation_3D(matr_A, il, im, ia,
425  "DT_L_left", "DT_L_right",
426  L, pc, alpha, DLL, Jacobian, -0.5);
427  SecondDerivativeApproximation_3D(matr_A, il, im, ia,
428  "DT_L_right", "DT_L_left",
429  L, pc, alpha, DLL, Jacobian, -0.5);
430  }
431 
432  if (pc_size >=3) {
433  // Energy term
434  if (L[il][im][ia] >= Lpp) {
435  SecondDerivativeApproximation_3D(matr_A, il, im, ia,
436  "DT_PC_left", "DT_PC_right",
437  L, pc, alpha, Dpcpc, Jacobian, -0.5);
438  SecondDerivativeApproximation_3D(matr_A, il, im, ia,
439  "DT_PC_right", "DT_PC_left",
440  L, pc, alpha, Dpcpc, Jacobian, -0.5);
441  } else {
442  SecondDerivativeApproximation_3D(matr_A, il, im, ia,
443  "DT_PC_left", "DT_PC_right",
444  L, pc, alpha, DpcpcLpp, Jacobian, -0.5);
445  SecondDerivativeApproximation_3D(matr_A, il, im, ia,
446  "DT_PC_right", "DT_PC_left",
447  L, pc, alpha, DpcpcLpp, Jacobian, -0.5);
448  }
449  }
450 
451  if (alpha_size >=3) {
452  // Pitch angle term(s)
453 
454  // Loss cone.
455  // Warning! Can cause the loss cone to be added during each diffusion, but
456  // now it's ok, cause alpha_size >= 3 only for pitch-angle diffusion
457  // f^{n+1}/taulc
458 
459 
460  if ( (alpha[il][im][ia] <= VF::alc(L[il][im][0]) || alpha[il][im][ia] >= VC::pi - VF::alc(L[il][im][0]))) {
461  double taulc;
462  //double taulc = VF::bounce_time(VF::Kfunc(pc[il][im][ia]), L[il][im][ia])/60./60./24.;
463  // Quarter bounce time
464  taulc = 0.25 * VF::bounce_time_new(L[il][im][ia], pc[il][im][ia], alpha[il][im][ia]);
465  matr_A[0][in] += 1.0/taulc; // Implicit
466  //matr_B[0][in] -= 1.0/taulc; // Explicit
467  }
468 
469  if (L[il][im][ia] >= Lpp) {
470  SecondDerivativeApproximation_3D(matr_A, il, im, ia,
471  "DT_ALPHA_left", "DT_ALPHA_right",
472  L, pc, alpha, Daa, Jacobian, -0.5);
473  SecondDerivativeApproximation_3D(matr_A, il, im, ia,
474  "DT_ALPHA_right", "DT_ALPHA_left",
475  L, pc, alpha, Daa, Jacobian, -0.5);
476  } else {
477  SecondDerivativeApproximation_3D(matr_A, il, im, ia,
478  "DT_ALPHA_left", "DT_ALPHA_right",
479  L, pc, alpha, DaaLpp, Jacobian, -0.5);
480  SecondDerivativeApproximation_3D(matr_A, il, im, ia,
481  "DT_ALPHA_right", "DT_ALPHA_left",
482  L, pc, alpha, DaaLpp, Jacobian, -0.5);
483  }
484  }
485 
486  if (alpha_size >= 3 && pc_size >= 3) {
487  // Mixed pitch angle and energy term
488  // There are several possibilities for approximation
489  if (L[il][im][ia] >= Lpp) {
490 
491  // d/dp D d/da f
492  SecondDerivativeApproximation_3D(matr_A, il, im, ia,
493  "DT_ALPHA_left", "DT_PC_right",
494  L, pc, alpha, Dpca, Jacobian, -0.25);
495  SecondDerivativeApproximation_3D(matr_A, il, im, ia,
496  "DT_ALPHA_right", "DT_PC_left",
497  L, pc, alpha, Dpca, Jacobian, -0.25);
498 
499  // d/da D d/dp f
500  SecondDerivativeApproximation_3D(matr_A, il, im, ia,
501  "DT_PC_left", "DT_ALPHA_right",
502  L, pc, alpha, Dpca, Jacobian, -0.25);
503  SecondDerivativeApproximation_3D(matr_A, il, im, ia,
504  "DT_PC_right", "DT_ALPHA_left",
505  L, pc, alpha, Dpca, Jacobian, -0.25);
506 
507  // d/dp D d/da f
508  SecondDerivativeApproximation_3D(matr_A, il, im, ia,
509  "DT_ALPHA_left", "DT_PC_left",
510  L, pc, alpha, Dpca, Jacobian, -0.25);
511  SecondDerivativeApproximation_3D(matr_A, il, im, ia,
512  "DT_ALPHA_right", "DT_PC_right",
513  L, pc, alpha, Dpca, Jacobian, -0.25);
514 
515  // d/da D d/dp f
516  SecondDerivativeApproximation_3D(matr_A, il, im, ia,
517  "DT_PC_left", "DT_ALPHA_left",
518  L, pc, alpha, Dpca, Jacobian, -0.25);
519  SecondDerivativeApproximation_3D(matr_A, il, im, ia,
520  "DT_PC_right", "DT_ALPHA_right",
521  L, pc, alpha, Dpca, Jacobian, -0.25);
522 
523  } else {
524 
525  // d/dp D d/da f
526  SecondDerivativeApproximation_3D(matr_A, il, im, ia,
527  "DT_ALPHA_left", "DT_PC_right",
528  L, pc, alpha, DpcaLpp, Jacobian, -0.25);
529  SecondDerivativeApproximation_3D(matr_A, il, im, ia,
530  "DT_ALPHA_right", "DT_PC_left",
531  L, pc, alpha, DpcaLpp, Jacobian, -0.25);
532  SecondDerivativeApproximation_3D(matr_A, il, im, ia,
533  "DT_ALPHA_left", "DT_PC_left",
534  L, pc, alpha, DpcaLpp, Jacobian, -0.25);
535  SecondDerivativeApproximation_3D(matr_A, il, im, ia,
536  "DT_ALPHA_right", "DT_PC_right",
537  L, pc, alpha, DpcaLpp, Jacobian, -0.25);
538 
539 
540  // d/da D d/dp f
541  SecondDerivativeApproximation_3D(matr_A, il, im, ia,
542  "DT_PC_left", "DT_ALPHA_right",
543  L, pc, alpha, DpcaLpp, Jacobian, -0.25);
544  SecondDerivativeApproximation_3D(matr_A, il, im, ia,
545  "DT_PC_right", "DT_ALPHA_left",
546  L, pc, alpha, DpcaLpp, Jacobian, -0.25);
547  SecondDerivativeApproximation_3D(matr_A, il, im, ia,
548  "DT_PC_left", "DT_ALPHA_left",
549  L, pc, alpha, DpcaLpp, Jacobian, -0.25);
550  SecondDerivativeApproximation_3D(matr_A, il, im, ia,
551  "DT_PC_right", "DT_ALPHA_right",
552  L, pc, alpha, DpcaLpp, Jacobian, -0.25);
553 
554  }
555  }
556  }
557  }
558  }
559  }
560  // Output::echo("recalculated.\n");
561 
562  // save the time of matrix change
563  matr_A.change_ind = clock();
564  matr_B.change_ind = clock();
565  matr_C.change_ind = clock();
566 
567  return true;
568 }
569 
570 // GMRES stuff starts from here
571 
572 //static double diag = 6.01;
574 
575 /**
576  * Vector to matrix multiplication for GMRES
577  */
578 void matvec(int m_size, double *x, double *b)
579 {
580  DiagMatrix::iterator it;
581  int in, ik;
582  for (in = 0; in < m_size; in++) {
583  // cout << in << ": ";
584  b[in] = 0;
585  // for (it = A_copy.begin(); ( (ik = in + it->first) < m_size && it != A_copy.end()); it++) { // 1 - 43 sec
586  for (it = A_copy.begin(); (it != A_copy.end()); it++) { // 2 - 46 sec
587  // it->first is diagonal number
588  // it->second is the pointer to diagonal array
589  // First I check, if the element at line (modelMatrixLineNumber) and diagonal (it->first) is inside the matrix.
590  // Cause it can be outside, if line number is 1 and diagonal number is -5, for example.
591  ik = in + it->first; // column number
592  if (ik >= 0 && ik < m_size) {
593  // if (in + it->first >= 0 && in + it->first < m_size) { // 2 - 46 sec
594  // if (in >= -it->first) {
595  // if (ik >= 0 ) { // 1 - 43 sec
596  // multiply
597  // cout << it->second[in] << " * " << x[in + it->first] << " + ";
598  b[in] += it->second[in] * x[ik];
599  // b[in] += A_copy[ik - in][in] * x[ik];
600  }
601  }
602  // cout << " = " << b[in] << endl;
603  }
604 }
605 
606 /*
607  * A pointer to the right preconditioner user routine, which computes
608  * w=B_r*z, where B_r should be an approximation of the inverse of the
609  * matrix A. If a null pointer is supplied, then a dummy preconditioner
610  * routine will be used which realizes the preconditioner matrix
611  * B_r=identity.
612  * int n input Number of vector components.
613  * double *z input Preconditioned iterate, of size n .
614  * double *w output Vector which holds the matrix-vector product B_r*z,
615  * i.e. the original iterate. */
616 
617 /**
618  * Jacobi (or diagonal) preconditioner:
619  */
620 void jacobi(int m_size, double *z, double *w)
621 {
622  int in;
623  for (in = 0; in < m_size; in++) {
624  if (A_copy[0][in] != 0) {
625  w[in] = z[in] / A_copy[0][in];
626  } else {
627  w[in] = z[in];
628  }
629  }
630 }
631 
632 /**
633  * SOR preconditioner:
634  *
635  * P = (D/rel + L)
636  * w = P^-1 * z;
637  * P * w = z
638  * (D/rel + L) * w = z
639  */
640 void SOR(int m_size, double *z, double *w)
641 {
642  int in;
643  DiagMatrix::iterator it;
644  double diag;
645  double rel = 1; // relaxation parameter - between >0 and <2
646  for (in = 0; in < m_size; in++) {
647  diag = A_copy[0][in] / rel;
648  w[in] = z[in] / diag;
649  // Lower triangular part of A
650  for (it = A_copy.begin(); (it->first < 0 && it != A_copy.end()); it++) {
651  if (in + it->first >= 0 && it->first < 0) { // double (redundant) check
652  // cout << (in + it->first) << ", " << it->second[in] << ", " << w[in + it->first] << endl;
653  w[in] -= w[in + it->first] * it->second[in] / diag;
654  }
655  }
656  }
657 }
658 
659 /**
660  * Simple preconditioner
661  */
662 void for_norm_A(int n, double *x, double *b)
663 {
664  int j;
665  for (j=0;j<n;j++) b[j] = x[j];
666 }
667 
668 
669 /**
670  * GMRES inversion.
671  * Potentially with bugs.
672  *
673  * A * X = B - equation
674  */
675 //void gmres_wrapout(CalculationMatrix &A, Matrix1D<double> &B, Matrix1D<double> &X, int maxiter, int i_max, double max_err) {
677  int m_size = A[0].size_x;
678  // int maxiter = 1000;
679 
680  int in;
681  //double sum = 0;
682  //double b_norm;
683  DiagMatrix::iterator it;
684 
685  // variables to store time
686  std::clock_t start_time, current_time;
687 
688  double err_max = 0.;
689  double error;
690 /*
691  for (in = 0; in < m_size; in++) {
692  // multiplying vector to matrix and calculating the difference
693  error = B[in];
694  for (it = A.begin(); it != A.end(); it++) {
695  if (in + it->first >= 0 && in + it->first < m_size) {
696  error -= it->second[in] * X[in + it->first];
697  }
698  }
699  err_max = max(err_max, fabs(error));
700  }
701  Output::echo(0, "Error_norm = %e \n", err_max);
702  err_max = 0;
703 */
704  A_copy = A;
705  static Matrix1D<double> B_copy;
706  B_copy = B;
707 /*
708  A_copy.writeToFile("./Debug/A_copy.dat");
709  B_copy.writeToFile("./Debug/B_copy.dat");
710 */
711  // normalization A
712  if (GMRES_parameters.use_normalization) {
713  double A0;
714  for (in = 0; in < m_size; in++) {
715  A0 = A_copy[0][in];
716  B_copy[in] = B_copy[in] / A0;
717  for (it = A_copy.begin(); it != A_copy.end(); it++) {
718  it->second[in] = it->second[in]/A0;
719  }
720  }
721  }
722 
723  // normalization of B
724  /* for (in = 0; in < m_size; in++) {
725  sum += B_copy[in] * B_copy[in];
726  }
727  b_norm = sqrt(sum);
728  cout << b_norm << endl;
729 
730  for (in = 0; in < m_size; in++) {
731  B_copy[in] = B_copy[in]/b_norm;
732  for (it = A_copy.begin(); it != A_copy.end(); it++) {
733  it->second[in] = it->second[in] / b_norm;
734  }
735  }
736 
737  A_copy.writeToFile("./Debug/A_copy.dat");
738  B_copy.writeToFile("./Debug/B_copy.dat");
739 
740  A.writeToFile("./Debug/A_check.dat"); */
741 
742  // struct ITLIN_OPT *opt = malloc(sizeof(struct ITLIN_OPT));
743  // struct ITLIN_INFO *info = malloc(sizeof(struct ITLIN_INFO));
744 
745  struct ITLIN_OPT *opt = new struct ITLIN_OPT;
746  struct ITLIN_INFO *info = new struct ITLIN_INFO;
747 
748  TERM_CHECK termcheck = CheckOnRestart;//CheckEachIter;
749 
750  /* initialize solver options, see gmres.c for help */
751  opt->tol = GMRES_parameters.SOL_max_iter_err; //max_err;
752  opt->i_max = GMRES_parameters.SOL_i_max;//i_max;
753  opt->termcheck = termcheck;
754  opt->maxiter = GMRES_parameters.SOL_maxiter;//maxiter;
755  opt->errorlevel = None;//Verbose;
756  opt->monitorlevel = None;//Verbose;
757  opt->datalevel = None;//Verbose;
758  opt->errorfile = stdout;
759  opt->monitorfile = stdout;
760  opt->datafile = NULL; //fopen("test_gmres.data","w");
761  /*if (!opt->datafile)
762  fprintf(stdout, "\n open of test_gmres.data failed!\n");*/
763  opt->iterfile = NULL;//fopen("test_gmres_iter.data","w");
764  opt->resfile = NULL;//fopen("test_gmres_res.data","w");
765  opt->miscfile = NULL;//fopen("test_gmres_misc.data","w");
766 
767  start_time = std::clock();
768  Output::echo(5, "Inverting... \n");
769 
770  if (GMRES_parameters.preconditioner_type == "none")
771  gmres(m_size, &X[0], &matvec, NULL, NULL, &B_copy[0], opt, info);
772  else if (GMRES_parameters.preconditioner_type == "jacobi")
773  gmres(m_size, &X[0], &matvec, &jacobi, NULL, &B_copy[0], opt, info);
774  else if (GMRES_parameters.preconditioner_type == "SOR")
775  gmres(m_size, &X[0], &matvec, &SOR, NULL, &B_copy[0], opt, info);
776  else if (GMRES_parameters.preconditioner_type == "for_norm_A")
777  gmres(m_size, &X[0], &matvec, &for_norm_A, NULL, &B_copy[0], opt, info);
778  else
779  throw error_msg("GMRES", "Unknown preconditioner type: %s", GMRES_parameters.preconditioner_type.c_str());
780 
781  current_time = std::clock();
782  Output::echo(5, "Done (%.2f sec)\n", (( current_time - start_time ) / (double)CLOCKS_PER_SEC ));
783 
784 
785  //fclose(opt->datafile);
786  //fclose(opt->iterfile);
787  //fclose(opt->miscfile);
788  fprintf(stdout,"\n Return code: %6i\n",info->rcode);
789  fprintf(stdout," Iterations: %6i\n",info->iter);
790  fprintf(stdout," Matvec calls: %6i\n",info->nomatvec);
791  fprintf(stdout," Right Precon calls: %6i\n",info->noprecr);
792  fprintf(stdout," Left Precon calls: %6i\n",info->noprecl);
793  fprintf(stdout," Estimated residual reduction: %e\n",info->precision);
794 
795  error = 0; err_max = 0;
796  for (in = 0; in < m_size; in++) {
797  // multiplying vector to matrix and calculating the difference
798  error = B[in];
799  for (it = A.begin(); it != A.end(); it++) {
800  if (in + it->first >= 0 && in + it->first < m_size) {
801  error -= it->second[in] * X[in + it->first];
802  }
803  }
804  err_max = max(err_max, fabs(error));
805  }
806  Output::echo(2, "Error_norm = %e \n", err_max);
807  err_max = 0;
808 
809  delete opt;
810  delete info;
811 
812 }
813 
815  const CPPL::dgbmatrix *A;
816 public:
817 
818  Preconditioner_jacobi(CPPL::dgbmatrix *A){
819  //this->A.copy(A);
820  this->A = A;
821  }
822 
823  // APPROXIMATELY solves M * X = RHS for X.
824  //const Matrix1D<double>& solve(const Matrix1D<double> &RHS) const;
825  CPPL::dcovector solve(const CPPL::dcovector &RHS) const {
826  int in, m_size = RHS.l;
827  CPPL::dcovector res(m_size);
828  for (in = 0; in < m_size; in++) {
829  if ((*A)(in, in) != 0) {
830  // divide to diagonal element
831  res(in) = RHS(in) / (*A)(in, in);
832  } else {
833  res(in) = RHS(in);
834  }
835  }
836  return res;
837  }
838 };
839 
841  const CPPL::dgbmatrix *A;
842 public:
843 
844  Preconditioner_SOR(CPPL::dgbmatrix *A){
845  //this->A.copy(A);
846  this->A = A;
847  }
848 
849  // APPROXIMATELY solves M * X = RHS for X.
850  CPPL::dcovector solve(const CPPL::dcovector &RHS) const {
851  int in, it, m_size = RHS.l;
852  CPPL::dcovector res(m_size);
853  double diag;
854  double rel = 1; // relaxation parameter - between >0 and <2
855  for (in = 0; in < m_size; in++) {
856  diag = (*A)(in, in) / rel;
857  res(in) = RHS(in) / diag;
858  // Lower triangular part of A
859  //for (it = A_copy.begin(); (it->first < 0 && it != A_copy.end()); it++) {
860  for (it = (*A).kl; it <= (*A).ku; it++) {
861  //if (in + it->first >= 0 && it->first < 0) { // ???
862  res(in) -= res(it) * (*A)(in, it) / diag;
863  //}
864  }
865  }
866  return res;
867  }
868 };
869 
870 
872  int i;
873 
874  // Matrix2D<double> H(GMRES_parameters.SOL_i_max, GMRES_parameters.SOL_i_max);
875  CPPL::dgematrix H_lp(GMRES_parameters.SOL_i_max+1, GMRES_parameters.SOL_i_max+1);
876 
877  CPPL::dcovector X_lp(X.size_x);
878  // copy X
879  for (i = 0; i < X.size_x; i++) X_lp(i) = X(i);
880 
881  CPPL::dcovector B_lp(B.size_x);
882  // copy B
883  for (i = 0; i < B.size_x; i++) B_lp(i) = B(i);
884 
885  long m_size = A[0].size_x;
886  DiagMatrix::iterator it;
887  it = A.begin();
888  long kl = -it->first; // first diagonal
889  it = A.end();
890  it--;
891  long ku = it->first; // last diagonal
892  CPPL::dgbmatrix A_lp(m_size, m_size, kl, ku); // it's our matrix A
893 
894  // zero A
895  int modelMatrixLineNumber;
896  for (modelMatrixLineNumber = 0; modelMatrixLineNumber < m_size; modelMatrixLineNumber++) {
897  for (int modelMatrixDiagNumber = -kl; modelMatrixDiagNumber <= ku; modelMatrixDiagNumber++) {
898  if (modelMatrixDiagNumber + modelMatrixLineNumber >= 0 && modelMatrixDiagNumber + modelMatrixLineNumber < m_size) {
899  A_lp(modelMatrixLineNumber, modelMatrixDiagNumber + modelMatrixLineNumber) = 0;
900  }
901  }
902  }
903 
904  // Fast-Copy values into the lapack library matrices
905  Output::echo(5, "Copying values...\n");
906  for (modelMatrixLineNumber = 0; modelMatrixLineNumber < m_size; modelMatrixLineNumber++) {
907  for (it = A.begin(); it != A.end(); it++) {
908  // check, if the element at line (modelMatrixLineNumber) and diagonal (it->first) is inside the matrix.
909  if (modelMatrixLineNumber + it->first >= 0 && modelMatrixLineNumber + it->first < m_size) {
910  // converting matrix, stored as diagonals, into lapack matrix (also diagonal)
911  A_lp(modelMatrixLineNumber, modelMatrixLineNumber + it->first) = it->second[modelMatrixLineNumber];
912  }
913  }
914  }
915 
916  Output::echo(3, "Inverting...\n");
917  int iter = GMRES_parameters.SOL_maxiter;
918  double tol = GMRES_parameters.SOL_max_iter_err;
919 
920  if (GMRES_parameters.preconditioner_type == "jacobi") {
921  Preconditioner_jacobi M(&A_lp);
922  GMRES(A_lp, X_lp, B_lp,
923  M, H_lp, GMRES_parameters.SOL_i_max, iter,
924  tol);
925  } else if (GMRES_parameters.preconditioner_type == "SOR") {
926  Preconditioner_SOR M(&A_lp);
927  GMRES(A_lp, X_lp, B_lp,
928  M, H_lp, GMRES_parameters.SOL_i_max, iter,
929  tol);
930  } else {
931  throw error_msg("GMRES", "Unknown preconditioner type: %s", GMRES_parameters.preconditioner_type.c_str());
932  }
933  Output::echo(3, "Done: iter %d, err %e\n", iter, tol);
934 
935  /*GMRES(A, X, B,
936  M, H, GMRES_parameters.SOL_i_max, GMRES_parameters.SOL_maxiter,
937  GMRES_parameters.SOL_max_iter_err);*/
938 
939  // save result into X
940  Output::echo(5, "Coping back... ");
941  for (modelMatrixLineNumber = 0; modelMatrixLineNumber < m_size; modelMatrixLineNumber++) {
942  X[modelMatrixLineNumber] = X_lp(modelMatrixLineNumber);
943  }
944 
945 }
946 
947 
948 /**
949  * Lapack inversion.
950  *
951  * A * X = B - equation
952  */
954 
955  // variables to store time
956  std::clock_t start_time, current_time;
957 
958  // iterator for diagonals of the diagonal matrix
959  DiagMatrix::iterator it;
960 
961  // model matrix line number
962  int modelMatrixLineNumber;
963 
964  long i;
965  long m_size = A[0].size_x;
966  it = A.begin();
967  long kl = -it->first; // first diagonal
968  it = A.end();
969  it--;
970  long ku = it->first; // last diagonal
971 
972  // Lapack library diagonal matrices
973  Output::echo(5, "Creating matrixes...\n");
974  CPPL::dgbmatrix lap_AB(m_size, m_size, kl, ku); // it's our matrix A
975  CPPL::dgbmatrix lap_AB_res(m_size, m_size, kl, ku);
976 
977  // !!! Make zero
978  for (modelMatrixLineNumber = 0; modelMatrixLineNumber < m_size; modelMatrixLineNumber++) {
979  for (int modelMatrixDiagNumber = -kl; modelMatrixDiagNumber <= ku; modelMatrixDiagNumber++) {
980  if (modelMatrixDiagNumber + modelMatrixLineNumber >= 0 && modelMatrixDiagNumber + modelMatrixLineNumber < m_size) {
981  lap_AB(modelMatrixLineNumber, modelMatrixDiagNumber + modelMatrixLineNumber) = 0;
982  }
983  }
984  }
985 
986 
987  // Copy values into the lapack library matrices
988  Output::echo(5, "Copying values...\n");
989  for (modelMatrixLineNumber = 0; modelMatrixLineNumber < m_size; modelMatrixLineNumber++) {
990  for (it = A.begin(); it != A.end(); it++) {
991  // it->first is diagonal number
992  // it->second is the pointer to diagonal array
993  // First I check, if the element at line (modelMatrixLineNumber) and diagonal (it->first) is inside the matrix.
994  // Cause it can be outside, if line number is 1 and diagonal number is -5, for example.
995  if (modelMatrixLineNumber + it->first >= 0 && modelMatrixLineNumber + it->first < m_size) {
996  // converting matrix, stored as diagonals, into lapack matrix (also diagonal)
997  lap_AB(modelMatrixLineNumber, modelMatrixLineNumber + it->first) = it->second[modelMatrixLineNumber];
998  }
999  }
1000  }
1001 
1002  // solution with lapack
1003  Output::echo(5, "matr --> res\n");
1004  lap_AB_res = lap_AB;
1005 
1006  // lapack vectors
1007  CPPL::dcovector lap_B(m_size);
1008  CPPL::dcovector lap_B_res(m_size);
1009  for(i = 0; i < lap_B.l; i++){
1010  lap_B(i) = B[i];
1011  }
1012  lap_B_res = lap_B;
1013 
1014  start_time = std::clock();
1015  Output::echo(5, "inverting... ");
1016 
1017  // lapack solution
1018  lap_AB.dgbsv(lap_B);
1019 
1020  current_time = std::clock();
1021  Output::echo(5, "done (%.2f sec)\n", (( current_time - start_time ) / (double)CLOCKS_PER_SEC ));
1022 
1023  Output::echo(5, "Checking... ");
1024  lap_B_res = lap_AB_res*lap_B - lap_B_res; // should be zero
1025 
1026  // calculate max error
1027  double max = 0;
1028  for (i = 0; i < lap_B_res.l; i++) {
1029  max = (max>lap_B_res(i))?max:lap_B_res(i);
1030  }
1031  Output::echo(5, " Max error: %e.\n", max);
1032 
1033  // save result into X
1034  Output::echo(5, "Coping back... ");
1035  for (modelMatrixLineNumber = 0; modelMatrixLineNumber < m_size; modelMatrixLineNumber++) {
1036  X[modelMatrixLineNumber] = lap_B(modelMatrixLineNumber);
1037  }
1038  Output::echo(5, "done. \n");
1039 
1040 }
1041 
1042 
1043 /**
1044  * Over relaxation iteration method.
1045  */
1046 void over_relaxation_diag(DiagMatrix &A, Matrix1D<double> &B, Matrix1D<double> &X, int max_steps, double EE) {
1047  int i, j;
1048  double omega, sum;
1049  double err_max = 1e99;
1050  int n = A[0].size_x;
1051 
1052  omega = 1.75;
1053  //omega = 1.1;
1054 
1055  DiagMatrix::iterator it;
1056  int nsteps = 0;
1057 
1058  while(err_max > EE && nsteps < max_steps) {
1059  nsteps++;
1060  for (i = 0; i < n; i++) {
1061  sum = 0;
1062  for (it = A.begin(); it != A.end(); it++) {
1063  j = i + it->first;
1064  if (j >= 0 && j < n && j != i) {
1065  sum = sum + it->second[i]/A[0][i] * X[j];
1066  }
1067  }
1068  X[i] = (1 - omega) * X[i] - omega * ( sum - B[i]/A[0][i]);
1069  }
1070 
1071  err_max = 0.;
1072  DiagMatrix::iterator it;
1073  double error;
1074  for (j = 0; j < n; j++) {
1075  // multiplying vector to matrix and calculating the difference
1076  error = B[j];
1077  for (it = A.begin(); it != A.end(); it++) {
1078  if (j + it->first >= 0 && j + it->first < n) {
1079  error -= it->second[j] * X[j + it->first];
1080  }
1081  }
1082  err_max = max(err_max, fabs(error));
1083  }
1084  }
1085  Output::echo(0, "[%d]Error_norm = %e \n", nsteps, err_max);
1086 }
1087 
1088 
1089 /**
1090  * Get change in indexes according to the derivatives direction.
1091  *
1092  * Used in approximation.
1093  * 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.
1094  */
1095 void GetDerivativeVector(string derivativeType, int &dL, int &dPc, int &dAlpha) {
1096  if (derivativeType == "DT_ALPHA_left") {
1097  dAlpha = -1;
1098  } else if (derivativeType == "DT_ALPHA_right") {
1099  dAlpha = 1;
1100  } else if (derivativeType == "DT_PC_left") {
1101  dPc = -1;
1102  } else if (derivativeType == "DT_PC_right") {
1103  dPc = 1;
1104  } else if (derivativeType == "DT_L_left") {
1105  dL = -1;
1106  } else if (derivativeType == "DT_L_right") {
1107  dL = 1;
1108  } else {
1109  throw error_msg("DERIVATIVE_TYPE", "Unknown derivative approximation: %s\n", derivativeType.c_str());
1110  }
1111 }
1112 
1113 
1114 /**
1115  * Second derivative approximation, returns coefficients to be putted into the model matrix.
1116  * \f$L_{\alpha \beta}(y) = (D_{\alpha \beta} \cdot y_{\bar{x}_\alpha})_{x_{\beta}}\f$
1117  *
1118  * Samarskiy, page 261
1119  *
1120  * Returns coefficients to be put into model matrix for an approximation of a second derivative.
1121  */
1123  int il, int im, int ia,
1124  string FirstDerivative, string SecondDerivative,
1126  Matrix3D<double> &Dxx, Matrix3D<double> &Jacobian,
1127  double multiplicator) {
1128 
1129  int dL1 = 0, dPc1 = 0, dAlpha1 = 0;
1130  int dL2 = 0, dPc2 = 0, dAlpha2 = 0;
1131  // calculating derivative directions according derivative types
1132  // each second derivative has two directions, like
1133  // (alpha, alpha) - is a second derivative in alpha direction
1134  // (pc, pc) - second in pc direction
1135  // (pc, alpha) - mixed term
1136  // here we are getting that directions for derivative
1137  GetDerivativeVector(FirstDerivative, dL1, dPc1, dAlpha1);
1138  GetDerivativeVector(SecondDerivative, dL2, dPc2, dAlpha2);
1139 
1140  // matr_A line number
1141  int in = matr_A.index1d(ia, im, il);
1142 
1143  double dh1, dh11, dh12;
1144  // The following trick work only for orthogonal grid
1145  // first and second dh in derivatove calculation
1146 
1147  // TODO: Rewrite in some normal way
1148  dh1 = (L[il][im][ia] + pc[il][im][ia] + alpha[il][im][ia]) - (L[il+dL1][im+dPc1][ia+dAlpha1] + pc[il+dL1][im+dPc1][ia+dAlpha1] + alpha[il+dL1][im+dPc1][ia+dAlpha1]);
1149  dh11 = (L[il][im][ia] + pc[il][im][ia] + alpha[il][im][ia]) - (L[il+dL2][im+dPc2][ia+dAlpha2] + pc[il+dL2][im+dPc2][ia+dAlpha2] + alpha[il+dL2][im+dPc2][ia+dAlpha2]);
1150  dh12 = (L[il+dL1][im+dPc1][ia+dAlpha1] + pc[il+dL1][im+dPc1][ia+dAlpha1] + alpha[il+dL1][im+dPc1][ia+dAlpha1]) - (L[il+dL1+dL2][im+dPc1+dPc2][ia+dAlpha1+dAlpha2] + pc[il+dL1+dL2][im+dPc1+dPc2][ia+dAlpha1+dAlpha2] + alpha[il+dL1+dL2][im+dPc1+dPc2][ia+dAlpha1+dAlpha2]);
1151 
1152  // for 2D
1153  //dh1 = (pc[il][im][ia] + alpha[il][im][ia]) - (pc[il][im+dPc1][ia+dAlpha1] + alpha[il][im+dPc1][ia+dAlpha1]);
1154  //dh11 = (pc[il][im][ia] + alpha[il][im][ia]) - (pc[il][im+dPc2][ia+dAlpha2] + alpha[il][im+dPc2][ia+dAlpha2]);
1155  //dh12 = (pc[il][im+dPc1][ia+dAlpha1] + alpha[il][im+dPc1][ia+dAlpha1]) - (pc[il][im+dPc1+dPc2][ia+dAlpha1+dAlpha2] + alpha[il][im+dPc1+dPc2][ia+dAlpha1+dAlpha2]);
1156 
1157  int id;
1158 
1159  // The same for all four coefficients, what's stending befire derivatives pretty much
1160  double common_part = multiplicator / Jacobian[il][im][ia] / dh1 ;
1161 
1162  // getting model matrix diagonal number according to derivative
1163  id = matr_A.index1d(ia, im, il) - in;
1164  // add corresponding coefficient for corresponding matrix cell
1165  matr_A[id][in] += + common_part / dh11 * Dxx[il][im][ia] * Jacobian[il][im][ia];
1166  // test - delete: matr_A[id][in] += + common_part / dh11 * Dxx[il][im][ia] * Jacobian[il][im][ia];
1167 
1168  id = matr_A.index1d(ia+dAlpha2, im+dPc2, il+dL2) - in;
1169  matr_A[id][in] += - common_part / dh11 * Dxx[il][im][ia] * Jacobian[il][im][ia];
1170  // test - delete: matr_A[id][in] += - common_part / dh12 * Dxx[il][im][ia] * Jacobian[il][im][ia];
1171  //matr_A[id][in] += -multiplicator / Jacobian[il][im][ia] / dh1 / dh11 * Dxx[il+dL2][im+dPc2][ia+dAlpha2] * Jacobian[il+dL2][im+dPc2][ia+dAlpha2];
1172 
1173  id = matr_A.index1d(ia+dAlpha1, im+dPc1, il+dL1) - in;
1174  matr_A[id][in] += - common_part / dh12 * Dxx[il+dL1][im+dPc1][ia+dAlpha1] * Jacobian[il+dL1][im+dPc1][ia+dAlpha1];
1175  // test - delete: matr_A[id][in] += - common_part / dh11 * Dxx[il+dL1][im+dPc1][ia+dAlpha1] * Jacobian[il+dL1][im+dPc1][ia+dAlpha1];
1176 
1177  id = matr_A.index1d(ia+dAlpha1+dAlpha2, im+dPc1+dPc2, il+dL1+dL2) - in;
1178  matr_A[id][in] += + common_part / dh12 * Dxx[il+dL1][im+dPc1][ia+dAlpha1] * Jacobian[il+dL1][im+dPc1][ia+dAlpha1];
1179  // test - delete: matr_A[id][in] += + common_part / dh12 * Dxx[il+dL1][im+dPc1][ia+dAlpha1] * Jacobian[il+dL1][im+dPc1][ia+dAlpha1];
1180  //matr_A[id][in] += +multiplicator / Jacobian[il][im][ia] / dh1 / dh12 * Dxx[il+dL1+dL2][im+dPc1+dPc2][ia+dAlpha1+dAlpha2] * Jacobian[il+dL1+dL2][im+dPc1+dPc2][ia+dAlpha1+dAlpha2];
1181 
1182 }
1183 
1184 
1185 
1186 #define EE 1.e-19
1187 #define I(l,k) (l)*n+(k)
1188 
1189 /// Function returns maximum between a and b.
1190 double max2(double a, double b) {
1191  return (a>b)?a:b;
1192 }
1193 /// Multiplocation between vector and matrix.
1194 void mult_vect2(double *af, double *vf, double *wf, int nf) {
1195  int k,l,n;
1196  n = nf;
1197  for (k = 0; k < nf; k++) {
1198  wf[k] = 0.;
1199  for(l=0; l < nf; l++)
1200  wf[k]+=af[( (k)*n+(l) )]*vf[l];
1201  }
1202 }
1203 
1204 
1205 /// Gauss inversion.
1206 void gauss_solve(double *a, double *b, int n) {
1207  int j, k, l, l_max;
1208  double a_max;
1209  double *dum;
1210 
1211  dum = (double*) malloc(n*sizeof(double));
1212 
1213  for (k = 0; k < (n-1); k++) {
1214  /* find max diag el */
1215 
1216  a_max = fabs ( a[I(k,k)] );
1217  l_max = k;
1218 
1219  for(l = k + 1; l < n; l++)
1220 
1221  if( fabs( a[I(l,k)] ) > a_max ) {
1222  a_max = fabs( a[I(l,k)] );
1223  l_max = l;
1224  }
1225  if( a_max < EE ) {
1226  printf("GAUSS: zero element string=%d size=%d \n", k, n);
1227  printf(" a(k,k)=%e \n",a[I(k,k)]);
1228  return;
1229  }
1230  if( l_max != k) {
1231  memcpy( dum, &(a[I(k,k)]), (n-k)*sizeof(double) );
1232  memcpy( &(a[I(k,k)]), &(a[I(l_max,k)]), (n-k)*sizeof(double) );
1233  memcpy( &(a[I(l_max,k)]), dum, (n-k)*sizeof(double) );
1234  a_max = b[k];
1235  b[k] = b[l_max];
1236  b[l_max] = a_max;
1237  }
1238 
1239  /* k-elemination */
1240  b[k] = b[k] / a[I(k,k)];
1241  for (l = k + 1; l < n; l++)
1242  a[( (k)*n+(l) )] = a[( (k)*n+(l) )] / a[I(k,k)];
1243  a[I(k,k)] = 1.;
1244  for (l = k + 1; l < n; l++) {
1245  for(j = k + 1; j < n; j++)
1246  a[I(l,j)] = a[I(l,j)] - a[I(k,j)] * a[I(l,k)];
1247  b[l] = b[l] - b[k] * a[I(l,k)];
1248  a[I(l,k)] = 0.;
1249  }
1250  }
1251 
1252  free(dum);
1253 }
1254 
1255 /// Solve the AU=R system of equations, where A - tridiagonal matrix nxn with diagonals a[], b[], c[].
1256 bool tridag(double a[], double b[], double c[], double r[], double u[], long n) {
1257  long j;
1258  double bet, *gam;
1259  gam = new double[n];
1260 
1261  if (b[0] == 0.0) {
1262  throw error_msg("TRIDAG_MATRIX_ERROR", "tridag: error, b[0] = 0");
1263  //log1::fcout << "error 1 in tridag\n";
1264  //return -1;
1265  }
1266 
1267  u[0]=r[0]/(bet=b[0]);
1268 
1269  for (j=1;j<=n-1;j++) {
1270  gam[j]=c[j-1]/bet;
1271  bet=b[j]-a[j]*gam[j];
1272 
1273  if (bet == 0.0) {
1274  throw error_msg("TRIDAG_MATRIX_ERROR", "tridag: error, bet = 0");
1275  //log1::fcout << "tridag(: Error 2 in tridag\n";
1276  //return -1;
1277  }
1278 
1279  u[j]=(r[j]-a[j]*u[j-1])/bet;
1280  }
1281 
1282  for (j=(n-2);j>=0;j--)
1283  u[j] -= gam[j+1]*u[j+1];
1284 
1285  delete gam;
1286 
1287  return true;
1288 }
1289 
1290 
1291 /** kckim test **/
1292 /**
1293  * Create matrix form of the Fokker-Planck equation: matr_A * PSD(t+1) = matr_B * PSD(t) + matr_C
1294  *
1295  * Calculation matr_A, matr_B, and matr_C, i.e. numerical approximation of the derivatives
1296  * Takes boundary conditions and diffusion coefficients as an input.
1297  * L, pc, alpha should be orthogonal for 3D!!!
1298  */
1300  CalculationMatrix &matr_A, CalculationMatrix &matr_B, CalculationMatrix &matr_C,
1302  int L_size, int pc_size, int alpha_size,
1303  Matrix2D<double> &L_lowerBoundaryCondition,
1304  Matrix2D<double> &L_upperBoundaryCondition,
1305  Matrix2D<double> &pc_lowerBoundaryCondition,
1306  Matrix2D<double> &pc_upperBoundaryCondition,
1307  Matrix2D<double> &alpha_lowerBoundaryCondition,
1308  Matrix2D<double> &alpha_upperBoundaryCondition,
1309  string L_lowerBoundaryCondition_calculationType,
1310  string L_upperBoundaryCondition_calculationType,
1311  string pc_lowerBoundaryCondition_calculationType,
1312  string pc_upperBoundaryCondition_calculationType,
1313  string alpha_lowerBoundaryCondition_calculationType,
1314  string alpha_upperBoundaryCondition_calculationType,
1315  Matrix3D<double> &DLL,
1316  Matrix3D<double> &Dpcpc, Matrix3D<double> &DpcpcLpp,
1317  Matrix3D<double> &Daa, Matrix3D<double> &DaaLpp,
1318  Matrix3D<double> &Dpca, Matrix3D<double> &DpcaLpp,
1319  Matrix3D<double> &Jacobian, double dt, double Lpp,
1320  double tau, double tauLpp) {
1321 
1322  // variables, for future use
1323  // il - L-index, im - pc-index, ia - alpha-index,
1324  // id - index of a diagonal of the model matrix (counting from main diagonal, main is zero),
1325  // in - model matrix line number
1326  int il, im, ia, id, in;
1327 
1328  // create A-matrix, B-matrix, C-matrix
1329  // Output::echo("Recalculating model matrix... ");
1330 
1331  // make everything zero
1332  // TODO: make zero not everything, but existing diagonals
1333  /* !!! matr_C[0] = 0;
1334  for (il = -1; il <= 1; il++) {
1335  for (im = -1; im <= 1; im++) {
1336  for (ia = -1; ia <= 1; ia++) {
1337  // calculating a diagonal number (id)
1338  id = matr_A.index1d(ia, im, il);
1339  // make the diagonal zero
1340  matr_A[id] = 0;
1341  matr_B[id] = 0;
1342  }
1343  }
1344  }*/
1345 
1346  // Make diagonals to be equal to zero
1347  DiagMatrix::iterator it;
1348  for (it = matr_A.begin(); it != matr_A.end(); it++) it->second = 0;
1349  for (it = matr_B.begin(); it != matr_B.end(); it++) it->second = 0;
1350  for (it = matr_C.begin(); it != matr_C.end(); it++) it->second = 0;
1351 
1352  // create a new model matrix
1353  // (f^{t+1} - f^{t})/dt = L1(f^{t+1}) + L2(f^{t+1}) + L3(f^{t+1})[main equation, losses are calculated separately]
1354  // L1, L2, L3 - diffusion operators
1355  // + need to add there the equations for boundary conditions also
1356  double dh;
1357  for (il = 0; il < L_size; il++) {
1358  for (im = 0; im < pc_size; im++) {
1359  for (ia = 0; ia < alpha_size; ia++) {
1360  // calculating current line number (in)
1361  in = matr_A.index1d(ia, im, il);
1362 
1363  // Describe all boundaries first. The 'inner' area will be approximated later.
1364  // If L_size > 1 we should use L-boundaries, right?
1365  // If not, then just calculate diffusion
1366  if (il == 0 && L_size >= 3) {
1367  // Boundary conditions
1368  // Points at the corners are used only at one of the boundaries
1369  // L-bottom boundary
1370  matr_C[0][in] = L_lowerBoundaryCondition[im][ia];
1371  id = matr_A.index1d(ia, im, il+1) - in;
1372  dh = L[il+1][im][ia] - L[il][im][ia];
1373  AddBoundary(matr_A, L_lowerBoundaryCondition_calculationType, in, id, dh);
1374 
1375  } else if (il == L_size-1 && L_size >= 3) {
1376  // L-top boundary
1377  matr_C[0][in] = L_upperBoundaryCondition[im][ia];
1378  id = matr_A.index1d(ia, im, il-1) - in;
1379  dh = L[il][im][ia] - L[il-1][im][ia];
1380  AddBoundary(matr_A, L_upperBoundaryCondition_calculationType, in, id, dh);
1381 
1382  } else if (im == 0 && pc_size >= 3) {
1383  // pc-low boundary
1384  matr_C[0][in] = pc_lowerBoundaryCondition[il][ia];
1385  id = matr_A.index1d(ia, im+1, il) - in;
1386  dh = pc[il][im+1][ia] - pc[il][im][ia];
1387  AddBoundary(matr_A, pc_lowerBoundaryCondition_calculationType, in, id, dh);
1388 
1389  } else if (im == pc_size-1 && pc_size >= 3) {
1390  // pc-high boundary
1391  matr_C[0][in] = pc_upperBoundaryCondition[il][ia];
1392  id = matr_A.index1d(ia, im-1, il) - in;
1393  dh = pc[il][im][ia] - pc[il][im-1][ia];
1394  AddBoundary(matr_A, pc_upperBoundaryCondition_calculationType, in, id, dh);
1395  } else if (ia == 0 && alpha_size >= 3) {
1396  // alpha-loss cone boundary
1397  matr_C[0][in] = alpha_lowerBoundaryCondition[il][im];
1398  id = matr_A.index1d(ia+1, im, il) - in;
1399  dh = alpha[il][im][ia+1] - alpha[il][im][ia];
1400  AddBoundary(matr_A, alpha_lowerBoundaryCondition_calculationType, in, id, dh);
1401 
1402  } else if (ia == alpha_size-1 && alpha_size >= 3) {
1403  // alpha-equatorial boundary
1404  matr_C[0][in] = alpha_upperBoundaryCondition[il][im];
1405  id = matr_A.index1d(ia-1, im, il) - in;
1406  dh = alpha[il][im][ia] - alpha[il][im][ia-1];
1407  AddBoundary(matr_A, alpha_upperBoundaryCondition_calculationType, in, id, dh);
1408  //cout<<in<<" "<<ia<<" "<<id<<" "<<dh<<" "<<alpha_upperBoundaryCondition[il][im]<<endl;
1409  } else {
1410 
1411  // now we are sure we are not on a boundary, can do the Fokker-Planck equation approximation in the inner area
1412 
1413  // f^{t+1}/dt
1414  matr_A[0][in] += 1.0/dt;
1415  // f^{t}/dt
1416  matr_B[0][in] += 1.0/dt;
1417 
1418  // calculate radial diffusion, if L_size >= 3
1419  if (L_size >=3) {
1420  // L term, add approximated derivatives into the calculation matrix
1422  "DT_L_left", "DT_L_right",
1423  L, pc, alpha, DLL, Jacobian, -0.5);
1425  "DT_L_right", "DT_L_left",
1426  L, pc, alpha, DLL, Jacobian, -0.5);
1427  }
1428 
1429  if (pc_size >=3) {
1430  // Energy term
1431  if (L[il][im][ia] >= Lpp) {
1432  /*
1433  SecondDerivativeApproximation_3D(matr_A, il, im, ia,
1434  "DT_PC_left", "DT_PC_right",
1435  L, pc, alpha, Dpcpc, Jacobian, -0.5);
1436  SecondDerivativeApproximation_3D(matr_A, il, im, ia,
1437  "DT_PC_right", "DT_PC_left",
1438  L, pc, alpha, Dpcpc, Jacobian, -0.5);
1439  */
1441  "DT_PC_right", "DT_PC_left",
1442  L, pc, alpha, Dpcpc, Jacobian, -1.0);
1443 
1444  } else {
1445 
1446  /*
1447  SecondDerivativeApproximation_3D(matr_A, il, im, ia,
1448  "DT_PC_left", "DT_PC_right",
1449  L, pc, alpha, DpcpcLpp, Jacobian, -0.5);
1450 
1451  SecondDerivativeApproximation_3D(matr_A, il, im, ia,
1452  "DT_PC_right", "DT_PC_left",
1453  L, pc, alpha, DpcpcLpp, Jacobian, -0.5);
1454  */
1455 
1457  "DT_PC_right", "DT_PC_left",
1458  L, pc, alpha, DpcpcLpp, Jacobian, -1.0);
1459 
1460  }
1461  }
1462 
1463  if (alpha_size >=3) {
1464  // Pitch angle term(s)
1465 
1466  // Loss cone.
1467  // Warning! Can cause the loss cone to be added during each diffusion, but
1468  // now it's ok, cause alpha_size >= 3 only for pitch-angle diffusion
1469  // f^{n+1}/taulc
1470 
1471 
1472  /*
1473  if ( (alpha[il][im][ia] <= VF::alc(L[il][im][0]) || alpha[il][im][ia] >= VC::pi - VF::alc(L[il][im][0]))) {
1474  double taulc = VF::bounce_time(VF::Kfunc(pc[il][im][ia]), L[il][im][ia])/60./60./24.;
1475  matr_A[0][in] += 1.0/taulc; // Implicit
1476  //matr_B[0][in] -= 1.0/taulc; // Explicit
1477  }*/
1478 
1479 
1480  if (L[il][im][ia] >= Lpp) {
1481  /*
1482  SecondDerivativeApproximation_3D(matr_A, il, im, ia,
1483  "DT_ALPHA_left", "DT_ALPHA_right",
1484  L, pc, alpha, Daa, Jacobian, -0.5);
1485  SecondDerivativeApproximation_3D(matr_A, il, im, ia,
1486  "DT_ALPHA_right", "DT_ALPHA_left",
1487  L, pc, alpha, Daa, Jacobian, -0.5);
1488  */
1490  "DT_ALPHA_right", "DT_ALPHA_left",
1491  L, pc, alpha, Daa, Jacobian, -1.0);
1492 
1493  } else {
1494  /*
1495  SecondDerivativeApproximation_3D(matr_A, il, im, ia,
1496  "DT_ALPHA_left", "DT_ALPHA_right",
1497  L, pc, alpha, DaaLpp, Jacobian, -0.5);
1498  SecondDerivativeApproximation_3D(matr_A, il, im, ia,
1499  "DT_ALPHA_right", "DT_ALPHA_left",
1500  L, pc, alpha, DaaLpp, Jacobian, -0.5);
1501  */
1502 
1504  "DT_ALPHA_right", "DT_ALPHA_left",
1505  L, pc, alpha, DaaLpp, Jacobian, -1.0);
1506 
1507  }
1508  }
1509 
1510  if (alpha_size >= 3 && pc_size >= 3) {
1511  // Mixed pitch angle and energy term
1512  // There are several possibilities for approximation
1513  if (L[il][im][ia] >= Lpp) {
1514  /*
1515  // d/dp D d/da f
1516  SecondDerivativeApproximation_3D_Mixed(matr_A, il, im, ia,
1517  "DT_ALPHA_right", "DT_PC_right",
1518  L, pc, alpha, Dpca, Jacobian, -1.0);
1519  SecondDerivativeApproximation_3D_Mixed(matr_A, il, im, ia,
1520  "DT_ALPHA_right", "DT_PC_left",
1521  L, pc, alpha, Dpca, Jacobian, +1.0);
1522 
1523 
1524  // d/da D d/dp f
1525  SecondDerivativeApproximation_3D_Mixed(matr_A, il, im, ia,
1526  "DT_PC_right", "DT_ALPHA_right",
1527  L, pc, alpha, Dpca, Jacobian, -1.0);
1528  SecondDerivativeApproximation_3D_Mixed(matr_A, il, im, ia,
1529  "DT_PC_right", "DT_ALPHA_left",
1530  L, pc, alpha, Dpca, Jacobian, +1.0);
1531  */
1532 
1533 
1534  // d/dp D d/da f
1535  SecondDerivativeApproximation_3D_KC_Mixed(matr_A, il, im, ia,
1536  "DT_ALPHA_right", "DT_PC_right",
1537  L, pc, alpha, Dpca, Jacobian, -1.0);
1538  SecondDerivativeApproximation_3D_KC_Mixed(matr_A, il, im, ia,
1539  "DT_ALPHA_left", "DT_PC_right",
1540  L, pc, alpha, Dpca, Jacobian, +1.0);
1541  SecondDerivativeApproximation_3D_KC_Mixed(matr_A, il, im, ia,
1542  "DT_ALPHA_right", "DT_PC_left",
1543  L, pc, alpha, Dpca, Jacobian, +1.0);
1544  SecondDerivativeApproximation_3D_KC_Mixed(matr_A, il, im, ia,
1545  "DT_ALPHA_left", "DT_PC_left",
1546  L, pc, alpha, Dpca, Jacobian, -1.0);
1547 
1548 
1549  // d/da D d/dp f
1550 
1551  SecondDerivativeApproximation_3D_KC_Mixed(matr_A, il, im, ia,
1552  "DT_PC_right", "DT_ALPHA_right",
1553  L, pc, alpha, Dpca, Jacobian, -1.0);
1554  SecondDerivativeApproximation_3D_KC_Mixed(matr_A, il, im, ia,
1555  "DT_PC_left", "DT_ALPHA_right",
1556  L, pc, alpha, Dpca, Jacobian, +1.0);
1557  SecondDerivativeApproximation_3D_KC_Mixed(matr_A, il, im, ia,
1558  "DT_PC_right", "DT_ALPHA_left",
1559  L, pc, alpha, Dpca, Jacobian, +1.0);
1560  SecondDerivativeApproximation_3D_KC_Mixed(matr_A, il, im, ia,
1561  "DT_PC_left", "DT_ALPHA_left",
1562  L, pc, alpha, Dpca, Jacobian, -1.0);
1563 
1564 
1565  } else {
1566 
1567  // full form
1568  // d/dp D d/da f
1569  /*
1570  SecondDerivativeApproximation_3D_Mixed(matr_A, il, im, ia,
1571  "DT_ALPHA_right", "DT_PC_right",
1572  L, pc, alpha, DpcaLpp, Jacobian, -1.0);
1573  SecondDerivativeApproximation_3D_Mixed(matr_A, il, im, ia,
1574  "DT_ALPHA_left", "DT_PC_right",
1575  L, pc, alpha, DpcaLpp, Jacobian, +1.0);
1576  SecondDerivativeApproximation_3D_Mixed(matr_A, il, im, ia,
1577  "DT_ALPHA_right", "DT_PC_left",
1578  L, pc, alpha, DpcaLpp, Jacobian, +1.0);
1579  SecondDerivativeApproximation_3D_Mixed(matr_A, il, im, ia,
1580  "DT_ALPHA_left", "DT_PC_left",
1581  L, pc, alpha, DpcaLpp, Jacobian, -1.0);
1582 
1583 
1584  // d/da D d/dp f
1585 
1586  SecondDerivativeApproximation_3D_Mixed(matr_A, il, im, ia,
1587  "DT_PC_right", "DT_ALPHA_right",
1588  L, pc, alpha, DpcaLpp, Jacobian, -1.0);
1589  SecondDerivativeApproximation_3D_Mixed(matr_A, il, im, ia,
1590  "DT_PC_left", "DT_ALPHA_right",
1591  L, pc, alpha, DpcaLpp, Jacobian, +1.0);
1592  SecondDerivativeApproximation_3D_Mixed(matr_A, il, im, ia,
1593  "DT_PC_right", "DT_ALPHA_left",
1594  L, pc, alpha, DpcaLpp, Jacobian, +1.0);
1595  SecondDerivativeApproximation_3D_Mixed(matr_A, il, im, ia,
1596  "DT_PC_left", "DT_ALPHA_left",
1597  L, pc, alpha, DpcaLpp, Jacobian, -1.0);
1598  */
1599 
1600  //short form
1601  // d/dp D d/da f
1602  SecondDerivativeApproximation_3D_KC_Mixed(matr_A, il, im, ia,
1603  "DT_ALPHA_right", "DT_PC_right",
1604  L, pc, alpha, DpcaLpp, Jacobian, -1.0);
1605  SecondDerivativeApproximation_3D_KC_Mixed(matr_A, il, im, ia,
1606  "DT_ALPHA_right", "DT_PC_left",
1607  L, pc, alpha, DpcaLpp, Jacobian, +1.0);
1608  SecondDerivativeApproximation_3D_KC_Mixed(matr_A, il, im, ia,
1609  "DT_ALPHA_left", "DT_PC_right",
1610  L, pc, alpha, DpcaLpp, Jacobian, +1.0);
1611  SecondDerivativeApproximation_3D_KC_Mixed(matr_A, il, im, ia,
1612  "DT_ALPHA_left", "DT_PC_left",
1613  L, pc, alpha, DpcaLpp, Jacobian, -1.0);
1614 
1615 
1616 
1617  // d/da D d/dp f
1618  SecondDerivativeApproximation_3D_KC_Mixed(matr_A, il, im, ia,
1619  "DT_PC_right", "DT_ALPHA_right",
1620  L, pc, alpha, DpcaLpp, Jacobian, -1.0);
1621  SecondDerivativeApproximation_3D_KC_Mixed(matr_A, il, im, ia,
1622  "DT_PC_right", "DT_ALPHA_left",
1623  L, pc, alpha, DpcaLpp, Jacobian, +1.0);
1624  SecondDerivativeApproximation_3D_KC_Mixed(matr_A, il, im, ia,
1625  "DT_PC_left", "DT_ALPHA_right",
1626  L, pc, alpha, DpcaLpp, Jacobian, +1.0);
1627  SecondDerivativeApproximation_3D_KC_Mixed(matr_A, il, im, ia,
1628  "DT_PC_left", "DT_ALPHA_left",
1629  L, pc, alpha, DpcaLpp, Jacobian, -1.0);
1630 
1631  }
1632  }
1633  }
1634  }
1635  }
1636  }
1637  // Output::echo("recalculated.\n");
1638 
1639  // save the time of matrix change
1640  matr_A.change_ind = clock();
1641  matr_B.change_ind = clock();
1642  matr_C.change_ind = clock();
1643 
1644  return true;
1645 }
1646 
1647 /** kckim test **/
1648 /**
1649  * Second derivative approximation, returns coefficients to be putted into the model matrix.
1650  * \f$L_{\alpha \beta}(y) = (D_{\alpha \beta} \cdot y_{\bar{x}_\alpha})_{x_{\beta}}\f$
1651  *
1652  * Samarskiy, page 261
1653  *
1654  * Returns coefficients to be put into model matrix for an approximation of a second derivative.
1655  */
1656 
1658  int il, int im, int ia,
1659  string FirstDerivative, string SecondDerivative,
1661  Matrix3D<double> &Dxx, Matrix3D<double> &Jacobian,
1662  double multiplicator) {
1663 
1664  int dL1 = 0, dPc1 = 0, dAlpha1 = 0;
1665  int dL2 = 0, dPc2 = 0, dAlpha2 = 0;
1666  // calculating derivative directions according derivative types
1667  // each second derivative has two directions, like
1668  // (alpha, alpha) - is a second derivative in alpha direction
1669  // (pc, pc) - second in pc direction
1670  // here we are getting that directions for derivative
1671  GetDerivativeVector(FirstDerivative, dL1, dPc1, dAlpha1);
1672  GetDerivativeVector(SecondDerivative, dL2, dPc2, dAlpha2);
1673 
1674  // matr_A line number
1675  int in = matr_A.index1d(ia, im, il);
1676 
1677  // The following trick work only for orthogonal grid
1678  // first and second dh in derivatove calculation
1679 
1680  double dh_1, dh_2; // _0: current point, _1: +1, _2: -1
1681  double Dxx_0, Dxx_1, Dxx_2, Dxx_11, Dxx_21;
1682  double J_0, J_1, J_2;
1683  // TODO: Rewrite in some normal way
1684 
1685  dh_1 = abs((L[il][im][ia] + pc[il][im][ia] + alpha[il][im][ia]) - (L[il+dL1][im+dPc1][ia+dAlpha1] + pc[il+dL1][im+dPc1][ia+dAlpha1] + alpha[il+dL1][im+dPc1][ia+dAlpha1]));
1686  dh_2 = abs((L[il][im][ia] + pc[il][im][ia] + alpha[il][im][ia]) - (L[il+dL2][im+dPc2][ia+dAlpha2] + pc[il+dL2][im+dPc2][ia+dAlpha2] + alpha[il+dL2][im+dPc2][ia+dAlpha2]));
1687 
1688  Dxx_0 = Dxx[il][im][ia];
1689  Dxx_1 = Dxx[il+dL1][im+dPc1][ia+dAlpha1];
1690  Dxx_2 = Dxx[il+dL2][im+dPc2][ia+dAlpha2];
1691  Dxx_21 = Dxx[il+dL1+dL2][im+dPc1+dPc2][ia+dAlpha1+dAlpha2];
1692 
1693  J_0 = Jacobian[il][im][ia];
1694  J_1 = Jacobian[il+dL1][im+dPc1][ia+dAlpha1];
1695  J_2 = Jacobian[il+dL2][im+dPc2][ia+dAlpha2];
1696 
1697  int id;
1698  double common;
1699 
1700  common = multiplicator / J_0 / (dh_1+dh_2) * 2;
1701  // getting model matrix diagonal number according to derivative
1702  id = matr_A.index1d(ia, im, il) - in;
1703  // add corresponding coefficient for corresponding matrix cell
1704  matr_A[id][in] += -common / dh_1 * (Dxx_0+Dxx_1)/2 * (J_0+J_1)/2; // f_0
1705 
1706  id = matr_A.index1d(ia+dAlpha1, im+dPc1, il+dL1) - in;
1707  matr_A[id][in] += +common / dh_1 * (Dxx_0+Dxx_1)/2 * (J_0+J_1)/2; // f_+1
1708 
1709  id = matr_A.index1d(ia+dAlpha2, im+dPc2, il+dL2) - in;
1710  matr_A[id][in] += +common / dh_2 * (Dxx_0+Dxx_2)/2 * (J_0+J_2)/2; // f_-1
1711 
1712  id = matr_A.index1d(ia+dAlpha1+dAlpha2, im+dPc1+dPc2, il+dL1+dL2) - in;
1713  matr_A[id][in] += -common / dh_2 * (Dxx_0+Dxx_2)/2 * (J_0+J_2)/2; // f_0 again because dA1+dA2 = 0, dPc1+dPc2 = 0
1714 
1715 }
1716 
1717 
1718 // short form mixed term
1719 
1721  int il, int im, int ia,
1722  string FirstDerivative, string SecondDerivative,
1724  Matrix3D<double> &Dxx, Matrix3D<double> &Jacobian,
1725  double multiplicator) {
1726 
1727  int dL1 = 0, dPc1 = 0, dAlpha1 = 0;
1728  int dL2 = 0, dPc2 = 0, dAlpha2 = 0;
1729  // calculating derivative directions according derivative types
1730  // each second derivative has two directions, like
1731  // (alpha, alpha) - is a second derivative in alpha direction
1732  // (pc, pc) - second in pc direction
1733  // here we are getting that directions for derivative
1734  GetDerivativeVector(FirstDerivative, dL1, dPc1, dAlpha1);
1735  GetDerivativeVector(SecondDerivative, dL2, dPc2, dAlpha2);
1736 
1737  // matr_A line number
1738  int in = matr_A.index1d(ia, im, il);
1739 
1740  // The following trick work only for orthogonal grid
1741  // first and second dh in derivatove calculation
1742 
1743  double dh_1, dh_2, dh_11, dh_21; // _0: current point, _1: +1, _2: -1
1744  double Dxx_0, Dxx_1, Dxx_2, Dxx_11, Dxx_21;
1745  double Dxx_c1, Dxx_c2, Dxx_c3, Dxx_c4;
1746  double J_0, J_1, J_2, J_11, J_21;
1747  double p_0;
1748  // TODO: Rewrite in some normal way
1749 
1750  dh_1 = abs((L[il][im][ia] + pc[il][im][ia] + alpha[il][im][ia]) - (L[il+dL1][im+dPc1][ia+dAlpha1] + pc[il+dL1][im+dPc1][ia+dAlpha1] + alpha[il+dL1][im+dPc1][ia+dAlpha1]));
1751  dh_2 = abs((L[il][im][ia] + pc[il][im][ia] + alpha[il][im][ia]) - (L[il+dL2][im+dPc2][ia+dAlpha2] + pc[il+dL2][im+dPc2][ia+dAlpha2] + alpha[il+dL2][im+dPc2][ia+dAlpha2]));
1752  dh_11 = abs((L[il][im][ia] + pc[il][im][ia] + alpha[il][im][ia]) - (L[il-dL1][im-dPc1][ia-dAlpha1] + pc[il-dL1][im-dPc1][ia-dAlpha1] + alpha[il-dL1][im-dPc1][ia-dAlpha1]));
1753  dh_21 = abs((L[il][im][ia] + pc[il][im][ia] + alpha[il][im][ia]) - (L[il-dL2][im-dPc2][ia-dAlpha2] + pc[il-dL2][im-dPc2][ia-dAlpha2] + alpha[il-dL2][im-dPc2][ia-dAlpha2]));
1754 
1755  Dxx_0 = Dxx[il][im][ia];
1756  Dxx_1 = Dxx[il+dL1][im+dPc1][ia+dAlpha1];
1757  Dxx_2 = Dxx[il+dL2][im+dPc2][ia+dAlpha2];
1758  Dxx_11 = Dxx[il-dL1][im-dPc1][ia-dAlpha1];
1759  Dxx_21 = Dxx[il-dL2][im-dPc2][ia-dAlpha2];
1760 
1761  // four points at the outer corners
1762  Dxx_c1 = Dxx[il+dL1+dL2][im+dPc1+dPc2][ia+dAlpha1+dAlpha2];
1763  Dxx_c2 = Dxx[il+dL1-dL2][im+dPc1-dPc2][ia+dAlpha1-dAlpha2];
1764  Dxx_c3 = Dxx[il-dL1+dL2][im-dPc1+dPc2][ia-dAlpha1+dAlpha2];
1765  Dxx_c4 = Dxx[il-dL1-dL2][im-dPc1-dPc2][ia-dAlpha1-dAlpha2];
1766 
1767 
1768  J_0 = Jacobian[il][im][ia];
1769  J_1 = Jacobian[il+dL1][im+dPc1][ia+dAlpha1];
1770  J_2 = Jacobian[il+dL2][im+dPc2][ia+dAlpha2];
1771  J_11 = Jacobian[il-dL1][im-dPc1][ia-dAlpha1];
1772  J_21 = Jacobian[il-dL2][im-dPc2][ia-dAlpha2];
1773 
1774  int id;
1775  double common;
1776 
1777  common = multiplicator / J_0 / (dh_1+dh_11) / (dh_2+dh_21);
1778  // getting model matrix diagonal number according to derivative
1779  id = matr_A.index1d(ia+dAlpha1, im+dPc1, il+dL1) - in;
1780  // add corresponding coefficient for corresponding matrix cell
1781  matr_A[id][in] += +common * (Dxx_0+Dxx_2)/2 * (J_0+J_2)/2; // f_0
1782  id = matr_A.index1d(ia+dAlpha1+dAlpha2, im+dPc1+dPc2, il+dL1+dL2) - in;
1783  matr_A[id][in] += +common * (Dxx_0+Dxx_2)/2 * (J_0+J_2)/2; // f_-1
1784 
1785 }
1786