23 #if defined(_WIN32) || defined(_WIN64)
24 #define _CRTDBG_MAP_ALLOC
30 #include <malloc/malloc.h>
41 #include "../Parameters/Parameters.h"
44 #include "../VariousFunctions/variousFunctions.h"
45 #include "../Exceptions/error.h"
48 #include "../Lapack/Include/cpplapack.h"
52 #include "../GMRES/itlin.h"
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);
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);
66 #include "../GMRES2/GMRES2.h"
87 string lower_border_condition_type,
88 string upper_border_condition_type,
89 string approximationMethod)
93 double Coef1, Dm_05, Dp_05, tau_cur;
97 double *
B =
new double[nx];
100 double *h =
new double[nx];
101 double *h_bar =
new double[nx];
105 double coef_05m, coef_05p;
108 for (j = 1; j <= nx-1; j++) {
109 h[j] = x[j] - x[j-1];
111 for (j = 0; j <= nx-2; j++) {
112 h_bar[j] = (h[j] + h[j+1])/2;
115 for (j = 1; j <= nx-2; j++) {
130 if (approximationMethod ==
"AM_Split_LR") {
132 Coef1 = pow(x[j], n)/h_bar[j];
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));
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));
141 A[j] = Coef1 * coef_05m / h[j];
142 C[j] = Coef1 * coef_05p / h[j+1];
144 }
else if (approximationMethod ==
"AM_Split_LR_oldway") {
146 Coef1 = pow(x[j], n);
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));
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));
155 A[j] = Coef1 * coef_05m / h[j] / h[j+1];
156 C[j] = Coef1 * coef_05p / h[j+1] / h[j];
171 else if (approximationMethod ==
"AM_Split_C") {
173 Coef1 = pow(x[j], n)/h_bar[j];
176 Coef1 = Coef1/
VF::G(x[j]);
178 Dm_05 = 0.5*(Dxx[j] + Dxx[j-1]);
179 Dp_05 = 0.5*(Dxx[j] + Dxx[j+1]);
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]);
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]);
189 }
else if (approximationMethod ==
"AM_Split_C_oldway") {
191 Coef1 = pow(x[j], n)/h_bar[j];
193 Coef1 = Coef1/
VF::G(x[j]);
195 Dm_05 = 0.5*(Dxx[j] + Dxx[j-1]);
196 Dp_05 = 0.5*(Dxx[j] + Dxx[j+1]);
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]);
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]);
206 throw error_msg(
"APPRERR",
"Unknown approximation type: %d", approximationMethod.c_str());
213 if ((x[j] < alc || x[j] >
VC::pi - alc) && (g_flag == 1)) {
219 B1[j] = B[j] - (1.0/dt + 1/tau_cur);
223 for (i = 1; i <= nx-2; i++) {
231 if (lower_border_condition_type ==
"BCT_CONSTANT_VALUE") {
232 A[i] = 0; B1[i] = 1; C[i] = 0;
233 }
else if (lower_border_condition_type ==
"BCT_CONSTANT_DERIVATIVE") {
234 A[i] = 0; B1[i] = -1; C[i] = 1;
236 throw error_msg(
"1D_DIFF_BOUNDARY",
"1d_universal_solver: unknown lower boundary type: %d", lower_border_condition_type.c_str());
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;
249 throw error_msg(
"1D_DIFF_BOUNDARY",
"1d_universal_solver: unknown lower boundary type: %d", upper_border_condition_type.c_str());
282 double *R1 =
new double[nx];
287 for (i = 1; i <= nx-2; i++) {
301 tridag(A, B1, C, R1, f, nx);
319 if (type ==
"BCT_CONSTANT_VALUE") {
320 Matrix_A[id0][in] = 1;
321 Matrix_A[id1][in] = 0;
322 }
else if (type ==
"BCT_CONSTANT_DERIVATIVE") {
323 Matrix_A[id0][in] = -1/dh;
324 Matrix_A[id1][in] = 1/dh;
326 throw error_msg(
"2D_DIFF_BOUNDARY",
"1d_universal_solver: unknown boundary type: %d", type.c_str());
374 int L_size,
int pc_size,
int alpha_size,
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,
392 double tau,
double tauLpp) {
398 int il, im, ia, id, in;
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;
429 for (il = 0; il < L_size; il++) {
430 for (im = 0; im < pc_size; im++) {
431 for (ia = 0; ia < alpha_size; ia++) {
433 in = matr_A.
index1d(ia, im, il);
438 if (il == 0 && L_size >= 3) {
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);
447 }
else if (il == L_size-1 && L_size >= 3) {
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);
454 }
else if (im == 0 && pc_size >= 3) {
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);
461 }
else if (im == pc_size-1 && pc_size >= 3) {
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) {
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);
474 }
else if (ia == alpha_size-1 && alpha_size >= 3) {
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);
486 matr_A[0][in] += 1.0/dt;
488 matr_B[0][in] += 1.0/dt;
494 "DT_L_left",
"DT_L_right",
495 L, pc, alpha, DLL, Jacobian, -0.5);
497 "DT_L_right",
"DT_L_left",
498 L, pc, alpha, DLL, Jacobian, -0.5);
503 if (L[il][im][ia] >= Lpp) {
505 "DT_PC_left",
"DT_PC_right",
506 L, pc, alpha, Dpcpc, Jacobian, -0.5);
508 "DT_PC_right",
"DT_PC_left",
509 L, pc, alpha, Dpcpc, Jacobian, -0.5);
512 "DT_PC_left",
"DT_PC_right",
513 L, pc, alpha, DpcpcLpp, Jacobian, -0.5);
515 "DT_PC_right",
"DT_PC_left",
516 L, pc, alpha, DpcpcLpp, Jacobian, -0.5);
520 if (alpha_size >=3) {
529 if ( (alpha[il][im][ia] <=
VF::alc(L[il][im][0]) || alpha[il][im][ia] >=
VC::pi -
VF::alc(L[il][im][0]))) {
534 matr_A[0][in] += 1.0/taulc;
538 if (L[il][im][ia] >= Lpp) {
540 "DT_ALPHA_left",
"DT_ALPHA_right",
541 L, pc, alpha, Daa, Jacobian, -0.5);
543 "DT_ALPHA_right",
"DT_ALPHA_left",
544 L, pc, alpha, Daa, Jacobian, -0.5);
547 "DT_ALPHA_left",
"DT_ALPHA_right",
548 L, pc, alpha, DaaLpp, Jacobian, -0.5);
550 "DT_ALPHA_right",
"DT_ALPHA_left",
551 L, pc, alpha, DaaLpp, Jacobian, -0.5);
555 if (alpha_size >= 3 && pc_size >= 3) {
558 if (L[il][im][ia] >= Lpp) {
562 "DT_ALPHA_left",
"DT_PC_right",
563 L, pc, alpha, Dpca, Jacobian, -0.25);
565 "DT_ALPHA_right",
"DT_PC_left",
566 L, pc, alpha, Dpca, Jacobian, -0.25);
570 "DT_PC_left",
"DT_ALPHA_right",
571 L, pc, alpha, Dpca, Jacobian, -0.25);
573 "DT_PC_right",
"DT_ALPHA_left",
574 L, pc, alpha, Dpca, Jacobian, -0.25);
578 "DT_ALPHA_left",
"DT_PC_left",
579 L, pc, alpha, Dpca, Jacobian, -0.25);
581 "DT_ALPHA_right",
"DT_PC_right",
582 L, pc, alpha, Dpca, Jacobian, -0.25);
586 "DT_PC_left",
"DT_ALPHA_left",
587 L, pc, alpha, Dpca, Jacobian, -0.25);
589 "DT_PC_right",
"DT_ALPHA_right",
590 L, pc, alpha, Dpca, Jacobian, -0.25);
596 "DT_ALPHA_left",
"DT_PC_right",
597 L, pc, alpha, DpcaLpp, Jacobian, -0.25);
599 "DT_ALPHA_right",
"DT_PC_left",
600 L, pc, alpha, DpcaLpp, Jacobian, -0.25);
602 "DT_ALPHA_left",
"DT_PC_left",
603 L, pc, alpha, DpcaLpp, Jacobian, -0.25);
605 "DT_ALPHA_right",
"DT_PC_right",
606 L, pc, alpha, DpcaLpp, Jacobian, -0.25);
611 "DT_PC_left",
"DT_ALPHA_right",
612 L, pc, alpha, DpcaLpp, Jacobian, -0.25);
614 "DT_PC_right",
"DT_ALPHA_left",
615 L, pc, alpha, DpcaLpp, Jacobian, -0.25);
617 "DT_PC_left",
"DT_ALPHA_left",
618 L, pc, alpha, DpcaLpp, Jacobian, -0.25);
620 "DT_PC_right",
"DT_ALPHA_right",
621 L, pc, alpha, DpcaLpp, Jacobian, -0.25);
647 void matvec(
int m_size,
double *x,
double *b)
649 DiagMatrix::iterator it;
651 for (in = 0; in < m_size; in++) {
655 for (it = A_copy.begin(); (it != A_copy.end()); it++) {
661 if (ik >= 0 && ik < m_size) {
667 b[in] += it->second[in] * x[ik];
689 void jacobi(
int m_size,
double *z,
double *w)
692 for (in = 0; in < m_size; in++) {
693 if (A_copy[0][in] != 0) {
694 w[in] = z[in] / A_copy[0][in];
709 void SOR(
int m_size,
double *z,
double *w)
712 DiagMatrix::iterator it;
715 for (in = 0; in < m_size; in++) {
716 diag = A_copy[0][in] / rel;
717 w[in] = z[in] / diag;
719 for (it = A_copy.begin(); (it->first < 0 && it != A_copy.end()); it++) {
720 if (in + it->first >= 0 && it->first < 0) {
722 w[in] -= w[in + it->first] * it->second[in] / diag;
734 for (j=0;j<n;j++) b[j] = x[j];
746 int m_size = A[0].size_x;
752 DiagMatrix::iterator it;
755 std::clock_t start_time, current_time;
783 for (in = 0; in < m_size; 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;
814 struct ITLIN_OPT *opt =
new struct ITLIN_OPT;
815 struct ITLIN_INFO *info =
new struct ITLIN_INFO;
817 TERM_CHECK termcheck = CheckOnRestart;
822 opt->termcheck = termcheck;
824 opt->errorlevel = None;
825 opt->monitorlevel = None;
826 opt->datalevel = None;
827 opt->errorfile = stdout;
828 opt->monitorfile = stdout;
829 opt->datafile = NULL;
832 opt->iterfile = NULL;
834 opt->miscfile = NULL;
836 start_time = std::clock();
840 gmres(m_size, &X[0], &
matvec, NULL, NULL, &B_copy[0], opt, info);
842 gmres(m_size, &X[0], &
matvec, &
jacobi, NULL, &B_copy[0], opt, info);
844 gmres(m_size, &X[0], &
matvec, &
SOR, NULL, &B_copy[0], opt, info);
850 current_time = std::clock();
851 Output::echo(5,
"Done (%.2f sec)\n", (( current_time - start_time ) / (
double)CLOCKS_PER_SEC ));
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);
864 error = 0; err_max = 0;
865 for (in = 0; in < m_size; 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];
873 err_max =
max(err_max, fabs(error));
887 const CPPL::dgbmatrix *A;
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) {
904 res(in) = RHS(in) / (*A)(in, in);
917 const CPPL::dgbmatrix *A;
927 CPPL::dcovector
solve(
const CPPL::dcovector &RHS)
const {
928 int in, it, m_size = RHS.l;
929 CPPL::dcovector res(m_size);
932 for (in = 0; in < m_size; in++) {
933 diag = (*A)(in, in) / rel;
934 res(in) = RHS(in) / diag;
937 for (it = (*A).kl; it <= (*A).ku; it++) {
939 res(in) -= res(it) * (*A)(in, it) / diag;
956 CPPL::dcovector X_lp(X.
size_x);
958 for (i = 0; i < X.
size_x; i++) X_lp(i) = X(i);
960 CPPL::dcovector B_lp(B.
size_x);
962 for (i = 0; i < B.
size_x; i++) B_lp(i) =
B(i);
964 long m_size = A[0].size_x;
965 DiagMatrix::iterator it;
967 long kl = -it->first;
971 CPPL::dgbmatrix A_lp(m_size, m_size, kl, ku);
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;
985 for (modelMatrixLineNumber = 0; modelMatrixLineNumber < m_size; modelMatrixLineNumber++) {
986 for (it = A.begin(); it != A.end(); it++) {
988 if (modelMatrixLineNumber + it->first >= 0 && modelMatrixLineNumber + it->first < m_size) {
990 A_lp(modelMatrixLineNumber, modelMatrixLineNumber + it->first) = it->second[modelMatrixLineNumber];
1001 GMRES(A_lp, X_lp, B_lp,
1002 M, H_lp, GMRES_parameters.
SOL_i_max, iter,
1006 GMRES(A_lp, X_lp, B_lp,
1007 M, H_lp, GMRES_parameters.
SOL_i_max, iter,
1020 for (modelMatrixLineNumber = 0; modelMatrixLineNumber < m_size; modelMatrixLineNumber++) {
1021 X[modelMatrixLineNumber] = X_lp(modelMatrixLineNumber);
1035 std::clock_t start_time, current_time;
1038 DiagMatrix::iterator it;
1041 int modelMatrixLineNumber;
1044 long m_size = A[0].size_x;
1046 long kl = -it->first;
1049 long ku = it->first;
1053 CPPL::dgbmatrix lap_AB(m_size, m_size, kl, ku);
1054 CPPL::dgbmatrix lap_AB_res(m_size, m_size, kl, ku);
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;
1068 for (modelMatrixLineNumber = 0; modelMatrixLineNumber < m_size; modelMatrixLineNumber++) {
1069 for (it = A.begin(); it != A.end(); it++) {
1074 if (modelMatrixLineNumber + it->first >= 0 && modelMatrixLineNumber + it->first < m_size) {
1076 lap_AB(modelMatrixLineNumber, modelMatrixLineNumber + it->first) = it->second[modelMatrixLineNumber];
1083 lap_AB_res = lap_AB;
1086 CPPL::dcovector lap_B(m_size);
1087 CPPL::dcovector lap_B_res(m_size);
1088 for(i = 0; i < lap_B.l; i++){
1093 start_time = std::clock();
1097 lap_AB.dgbsv(lap_B);
1099 current_time = std::clock();
1100 Output::echo(5,
"done (%.2f sec)\n", (( current_time - start_time ) / (
double)CLOCKS_PER_SEC ));
1103 lap_B_res = lap_AB_res*lap_B - lap_B_res;
1107 for (i = 0; i < lap_B_res.l; i++) {
1108 max = (max>lap_B_res(i))?max:lap_B_res(i);
1114 for (modelMatrixLineNumber = 0; modelMatrixLineNumber < m_size; modelMatrixLineNumber++) {
1115 X[modelMatrixLineNumber] = lap_B(modelMatrixLineNumber);
1135 double err_max = 1e99;
1136 int n = A[0].size_x;
1142 DiagMatrix::iterator it;
1145 while(err_max > EE && nsteps < max_steps) {
1147 for (i = 0; i < n; i++) {
1149 for (it = A.begin(); it != A.end(); it++) {
1152 if (j >= 0 && j < n && j != i) {
1155 sum = sum + it->second[i]/A[0][i] * X[j];
1158 X[i] = (1 - omega) * X[i] - omega * ( sum - B[i]/A[0][i]);
1164 for (j = 0; j < n; 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];
1172 err_max =
max(err_max, fabs(error));
1175 Output::echo(0,
"[%d]Error_norm = %e \n", nsteps, err_max);
1186 if (derivativeType ==
"DT_ALPHA_left") {
1188 }
else if (derivativeType ==
"DT_ALPHA_right") {
1190 }
else if (derivativeType ==
"DT_PC_left") {
1192 }
else if (derivativeType ==
"DT_PC_right") {
1194 }
else if (derivativeType ==
"DT_L_left") {
1196 }
else if (derivativeType ==
"DT_L_right") {
1199 throw error_msg(
"DERIVATIVE_TYPE",
"Unknown derivative approximation: %s\n", derivativeType.c_str());
1213 int il,
int im,
int ia,
1214 string FirstDerivative,
string SecondDerivative,
1217 double multiplicator) {
1219 int dL1 = 0, dPc1 = 0, dAlpha1 = 0;
1220 int dL2 = 0, dPc2 = 0, dAlpha2 = 0;
1231 int in = matr_A.
index1d(ia, im, il);
1233 double dh1, dh11, dh12;
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]);
1250 double common_part = multiplicator / Jacobian[il][im][ia] / dh1 ;
1253 id = matr_A.
index1d(ia, im, il) - in;
1255 matr_A[id][in] += + common_part / dh11 * Dxx[il][im][ia] * Jacobian[il][im][ia];
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];
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];
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];
1277 #define I(l,k) (l)*n+(k)
1288 for (k = 0; k < nf; k++) {
1290 for(l=0; l < nf; l++)
1291 wf[k]+=af[( (k)*n+(l) )]*vf[l];
1302 int modelMatrixLineNumber;
1304 DiagMatrix::iterator it;
1305 double sum, error, error_max;
1306 static bool recalculate_m_A =
true;
1308 if (recalculate_m_A) {
1321 for (it = A.begin(); it != A.end(); it++) {
1322 for (modelMatrixLineNumber = 0; modelMatrixLineNumber < n; modelMatrixLineNumber++) {
1325 if (modelMatrixLineNumber + it->first >= 0 && modelMatrixLineNumber + it->first < n) {
1329 m_A[modelMatrixLineNumber][modelMatrixLineNumber + it->first] = it->second[modelMatrixLineNumber];
1339 recalculate_m_A =
true;
1344 X[n-1] = m_B[n-1] / m_A[n-1][n-1];
1345 for (k = n-2; k >= 0; k--) {
1347 for (l = k+1; l < n; l++)
1348 sum += m_A[k][l] * X[l];
1349 X[k] = m_B[k] - sum;
1357 for (in = 0; in < n; in++) {
1359 for (it = A.begin(); it != A.end(); it++) {
1361 if (in + it->first >= 0 && in + it->first < n) {
1362 error -= it->second[in] * X[in + it->first];
1365 error_max =
VF::max( error_max, error );
1379 dum = (
double*) malloc(n*
sizeof(
double));
1381 for (k = 0; k < (n-1); k++) {
1384 a_max = fabs ( a[
I(k,k)] );
1387 for(l = k + 1; l < n; l++)
1389 if( fabs( a[
I(l,k)] ) > a_max ) {
1390 a_max = fabs( a[
I(l,k)] );
1394 printf(
"GAUSS: zero element string=%d size=%d \n", k, n);
1395 printf(
" a(k,k)=%e \n",a[
I(k,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) );
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)];
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)];
1424 bool tridag(
double a[],
double b[],
double c[],
double r[],
double u[],
long n) {
1427 gam =
new double[n];
1430 throw error_msg(
"TRIDAG_MATRIX_ERROR",
"tridag: error, b[0] = 0");
1435 u[0]=r[0]/(bet=b[0]);
1437 for (j=1;j<=n-1;j++) {
1439 bet=b[j]-a[j]*gam[j];
1442 throw error_msg(
"TRIDAG_MATRIX_ERROR",
"tridag: error, bet = 0");
1447 u[j]=(r[j]-a[j]*u[j-1])/bet;
1450 for (j=(n-2);j>=0;j--)
1451 u[j] -= gam[j+1]*u[j+1];
1470 int L_size,
int pc_size,
int alpha_size,
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,
1488 double tau,
double tauLpp) {
1494 int il, im, ia, id, in;
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;
1525 for (il = 0; il < L_size; il++) {
1526 for (im = 0; im < pc_size; im++) {
1527 for (ia = 0; ia < alpha_size; ia++) {
1529 in = matr_A.
index1d(ia, im, il);
1534 if (il == 0 && L_size >= 3) {
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);
1543 }
else if (il == L_size-1 && L_size >= 3) {
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);
1550 }
else if (im == 0 && pc_size >= 3) {
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);
1557 }
else if (im == pc_size-1 && pc_size >= 3) {
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) {
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);
1570 }
else if (ia == alpha_size-1 && alpha_size >= 3) {
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);
1582 matr_A[0][in] += 1.0/dt;
1584 matr_B[0][in] += 1.0/dt;
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);
1599 if (L[il][im][ia] >= Lpp) {
1609 "DT_PC_right",
"DT_PC_left",
1610 L, pc, alpha, Dpcpc, Jacobian, -1.0);
1625 "DT_PC_right",
"DT_PC_left",
1626 L, pc, alpha, DpcpcLpp, Jacobian, -1.0);
1631 if (alpha_size >=3) {
1648 if (L[il][im][ia] >= Lpp) {
1658 "DT_ALPHA_right",
"DT_ALPHA_left",
1659 L, pc, alpha, Daa, Jacobian, -1.0);
1672 "DT_ALPHA_right",
"DT_ALPHA_left",
1673 L, pc, alpha, DaaLpp, Jacobian, -1.0);
1678 if (alpha_size >= 3 && pc_size >= 3) {
1681 if (L[il][im][ia] >= Lpp) {
1704 "DT_ALPHA_right",
"DT_PC_right",
1705 L, pc, alpha, Dpca, Jacobian, -1.0);
1707 "DT_ALPHA_left",
"DT_PC_right",
1708 L, pc, alpha, Dpca, Jacobian, +1.0);
1710 "DT_ALPHA_right",
"DT_PC_left",
1711 L, pc, alpha, Dpca, Jacobian, +1.0);
1713 "DT_ALPHA_left",
"DT_PC_left",
1714 L, pc, alpha, Dpca, Jacobian, -1.0);
1720 "DT_PC_right",
"DT_ALPHA_right",
1721 L, pc, alpha, Dpca, Jacobian, -1.0);
1723 "DT_PC_left",
"DT_ALPHA_right",
1724 L, pc, alpha, Dpca, Jacobian, +1.0);
1726 "DT_PC_right",
"DT_ALPHA_left",
1727 L, pc, alpha, Dpca, Jacobian, +1.0);
1729 "DT_PC_left",
"DT_ALPHA_left",
1730 L, pc, alpha, Dpca, Jacobian, -1.0);
1771 "DT_ALPHA_right",
"DT_PC_right",
1772 L, pc, alpha, DpcaLpp, Jacobian, -1.0);
1774 "DT_ALPHA_right",
"DT_PC_left",
1775 L, pc, alpha, DpcaLpp, Jacobian, +1.0);
1777 "DT_ALPHA_left",
"DT_PC_right",
1778 L, pc, alpha, DpcaLpp, Jacobian, +1.0);
1780 "DT_ALPHA_left",
"DT_PC_left",
1781 L, pc, alpha, DpcaLpp, Jacobian, -1.0);
1787 "DT_PC_right",
"DT_ALPHA_right",
1788 L, pc, alpha, DpcaLpp, Jacobian, -1.0);
1790 "DT_PC_right",
"DT_ALPHA_left",
1791 L, pc, alpha, DpcaLpp, Jacobian, +1.0);
1793 "DT_PC_left",
"DT_ALPHA_right",
1794 L, pc, alpha, DpcaLpp, Jacobian, +1.0);
1796 "DT_PC_left",
"DT_ALPHA_left",
1797 L, pc, alpha, DpcaLpp, Jacobian, -1.0);
1826 int il,
int im,
int ia,
1827 string FirstDerivative,
string SecondDerivative,
1830 double multiplicator) {
1832 int dL1 = 0, dPc1 = 0, dAlpha1 = 0;
1833 int dL2 = 0, dPc2 = 0, dAlpha2 = 0;
1843 int in = matr_A.
index1d(ia, im, il);
1849 double Dxx_0, Dxx_1, Dxx_2, Dxx_11, Dxx_21;
1850 double J_0, J_1, J_2;
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]));
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];
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];
1868 common = multiplicator / J_0 / (dh_1+dh_2) * 2;
1870 id = matr_A.
index1d(ia, im, il) - in;
1872 matr_A[id][in] += -common / dh_1 * (Dxx_0+Dxx_1)/2 * (J_0+J_1)/2;
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;
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;
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;
1905 int il,
int im,
int ia,
1906 string FirstDerivative,
string SecondDerivative,
1909 double multiplicator) {
1911 int dL1 = 0, dPc1 = 0, dAlpha1 = 0;
1912 int dL2 = 0, dPc2 = 0, dAlpha2 = 0;
1922 int in = matr_A.
index1d(ia, im, il);
1927 double dh_1, dh_2, dh_11, dh_21;
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;
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]));
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];
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];
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];
1961 common = multiplicator / J_0 / (dh_1+dh_11) / (dh_2+dh_21);
1963 id = matr_A.
index1d(ia+dAlpha1, im+dPc1, il+dL1) - in;
1965 matr_A[id][in] += +common * (Dxx_0+Dxx_2)/2 * (J_0+J_2)/2;
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;
bool initialized
Flag, equal true if initialized.
string preconditioner_type
Type of preconditioner to use.
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)
void AllocateMemory(int size_x, int size_y)
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.
void echo(int msglvl, const char *format,...)
Basic function! Write the message to the file.
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)
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.
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...
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
bool use_normalization
To normalize (or not) matrix A from Ax=B system.
double alc(double L)
Loss cone calculations.
#define I(l, k)
Specified Index [l][k].
GMRES parameters structure.
Error message - stack of single_errors.
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.
double G(double alpha)
G-function.