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