25 #include "../VariousFunctions/variousFunctions.h"
30 #include "../Exceptions/error.h"
50 if (!arr.initialized) {
65 if (grid.
L.
size != 1) {
66 throw error_msg(
"INITIAL_PSD_2D_METHOD_FOR_3D_GRID",
"Initial PSD error: impossible to use ORBIT_FLUX_2D if L.size > 1");
72 if (grid.
L.
size != 1) {
73 throw error_msg(
"INITIAL_PSD_2D_METHOD_FOR_3D_GRID",
"Initial PSD error: impossible to use FLUX_2D_MAXWELLIAN if L.size > 1");
75 cout<<
"kim................................"<<endl;
76 Load_initial_f_maxwell(grid.
L, grid.
pc, grid.
alpha);
82 if (grid.
L.
size == 1) {
83 throw error_msg(
"INITIAL_PSD_3D_METHOD_FOR_2D_GRID",
"Initial PSD error: impossible to use STEADY_STATE if L.size == 1");
93 if (grid.
L.
size == 1) {
94 throw error_msg(
"INITIAL_PSD_3D_METHOD_FOR_2D_GRID",
"Initial PSD error: impossible to use STEADY_STATE_L7 if L.size == 1");
103 if (grid.
L.
size == 1) {
104 throw error_msg(
"INITIAL_PSD_3D_METHOD_FOR_2D_GRID",
"Initial PSD error: impossible to use STEADY_STATE_TWO_ZONE if L.size == 1");
109 else if (parameters.
initial_PSD_Type ==
"IPSDT_STEADY_STATE_FROM_BOUNDARY") {
111 if (grid.
L.
size == 1) {
112 throw error_msg(
"INITIAL_PSD_3D_METHOD_FOR_2D_GRID",
"Initial PSD error: impossible to use STEADY_STATE_FROM_BOUNDARY if L.size == 1");
115 throw error_msg(
"INITIAL_PSD_ERROR",
"Initial PSD error: upper boundary is not initialized, can not be used for steady state calculation.");
119 Load_initial_f_from_outer_L(grid.
L, grid.
pc, grid.
alpha,
131 for (il = 0; il < grid.
L.
size; il++)
132 for (im = 0; im < grid.
pc.
size; im++)
145 throw error_msg(
"INITIAL_PSD_TYPE_UNKNOWN",
"Unknown initial PSD type");
153 for (il = 0; il < grid.
L.
size; il++)
154 for (im = 0; im < grid.
pc.
size; im++)
155 for (ia = 0; ia < grid.
alpha.
size; ia++) {
158 double a = grid.
alpha.
arr[il][im][ia];
159 if (a < lc || a >
VC::pi - lc)
165 arr[il][im][ia] = arr[il][im][ia] * pow(sinh(a/lc),2)/pow(sinh(1.0),2);
175 output_without_grid_file <<
"VARIABLES = \"Phase_Space_Density\" " << endl;
176 output_without_grid_file.setf(ios_base::scientific, ios_base::floatfield);
323 double dt,
double Lpp,
324 double tau,
double tauLpp,
double Kp)
331 double tau_p,tau_c,tauLpp0;
334 for (il = 0; il < L.
size; il++) {
335 for (im = 0; im < pc.
size; im++) {
336 for (ia = 0; ia < alpha.
size; ia++) {
354 if (L.
arr[il][im][ia] >= Lpp) {
356 arr[il][im][ia] = arr[il][im][ia] - arr[il][im][ia] * (dt / tau);
366 }
else if (parameters.
usetauLpp ==
"precipitaion") {
368 }
else if (parameters.
usetauLpp ==
"combined") {
371 tauLpp0 = 1.0/(1.0/tau_p+1.0/tau_c);
375 arr[il][im][ia] = arr[il][im][ia] - arr[il][im][ia] * (dt / tauLpp0);
384 arr[il][im][ia] = arr[il][im][ia] * SL.
arr[il][im][ia];
421 string pc_lowerBoundaryCondition_calculationType,
422 string pc_upperBoundaryCondition_calculationType,
423 string alpha_lowerBoundaryCondition_calculationType,
424 string alpha_upperBoundaryCondition_calculationType)
437 double Dpca_, Dpca_a_L, Dpca_a_R, Dpca_m_L, Dpca_m_R;
438 double J_, J_a_L, J_a_R, J_m_L, J_m_R;
440 double Laml, Lmal, Lamr, Lmar;
441 double PSD_m_r, PSD_m_l, PSD_a_r, PSD_a_l, Lam, Lma;
452 if (PSD_parameters.approximationMethod ==
"AM_Split_LR") {
455 for (il = 0; il < L.
size; il++) {
456 for (im = 0; im < pc.
size; im++) {
457 for (ia = 1; ia < alpha.
size-1; ia++) {
459 PSD_Der_a_L[il][im][ia] =
VF::G(alpha.
arr[il][im][ia]) * Dpca.
arr[il][im][ia]*(arr[il][im][ia] - arr[il][im][ia-1])/(alpha.
arr[il][im][ia] - alpha.
arr[il][im][ia-1]);
461 PSD_Der_a_R[il][im][ia] =
VF::G(alpha.
arr[il][im][ia]) * Dpca.
arr[il][im][ia]*(arr[il][im][ia+1] - arr[il][im][ia])/(alpha.
arr[il][im][ia+1] - alpha.
arr[il][im][ia]);
464 for (im = 1; im < pc.
size-1; im++) {
465 for (ia = 0; ia < alpha.
size; ia++) {
467 PSD_Der_m_L[il][im][ia] = pc.
arr[il][im][ia] * pc.
arr[il][im][ia] * Dpca.
arr[il][im][ia]*(arr[il][im][ia] - arr[il][im-1][ia])/(pc.
arr[il][im][ia] - pc.
arr[il][im-1][ia]);
469 PSD_Der_m_R[il][im][ia] = pc.
arr[il][im][ia] * pc.
arr[il][im][ia] * Dpca.
arr[il][im][ia]*(arr[il][im+1][ia] - arr[il][im][ia])/(pc.
arr[il][im+1][ia] - pc.
arr[il][im][ia]);
475 for (il = 0; il < L.
size; il++) {
476 for (im = 1; im < pc.
size-1; im++) {
477 for (ia = 1; ia < alpha.
size-1; ia++) {
479 Dpca_ = Dpca.
arr[il][im][ia];
480 Dpca_a_L = Dpca.
arr[il][im][ia-1];
481 Dpca_a_R = Dpca.
arr[il][im][ia+1];
482 Dpca_m_L = Dpca.
arr[il][im-1][ia];
483 Dpca_m_R = Dpca.
arr[il][im+1][ia];
484 J_ = Jacobian[il][im][ia];
485 J_a_L = Jacobian[il][im][ia-1];
486 J_a_R = Jacobian[il][im][ia+1];
487 J_m_L = Jacobian[il][im-1][ia];
488 J_m_R = Jacobian[il][im+1][ia];
496 Lmar = 1.0/(pc.
arr[il][im][ia] * pc.
arr[il][im][ia]) * (PSD_Der_a_R[il][im][ia] - PSD_Der_a_R[il][im-1][ia])/(pc.
arr[il][im][ia] - pc.
arr[il][im-1][ia]);
497 Lmal = 1.0/(pc.
arr[il][im][ia] * pc.
arr[il][im][ia]) * (PSD_Der_a_L[il][im+1][ia] - PSD_Der_a_L[il][im][ia])/(pc.
arr[il][im+1][ia] - pc.
arr[il][im][ia]);
504 Lamr = 1.0/(
VF::G(alpha.
arr[il][im][ia])) * (PSD_Der_m_R[il][im][ia] - PSD_Der_m_R[il][im][ia-1])/(alpha.
arr[il][im][ia] - alpha.
arr[il][im][ia-1]);
505 Laml = 1.0/(
VF::G(alpha.
arr[il][im][ia])) * (PSD_Der_m_L[il][im][ia+1] - PSD_Der_m_L[il][im][ia])/(alpha.
arr[il][im][ia+1] - alpha.
arr[il][im][ia]);
509 arr[il][im][ia] = arr[il][im][ia] + dt * 0.5 * (Lmar + Lamr + Lmal + Laml);
537 else if (PSD_parameters.approximationMethod ==
"AM_Split_C") {
538 for (il = 0; il < L.
size; il++) {
539 for (im = 1; im < pc.
size-1; im++) {
540 for (ia = 1; ia < alpha.
size-1; ia++) {
542 Dpca_ = Dpca.
arr[il][im][ia];
543 Dpca_a_L = Dpca.
arr[il][im][ia-1];
544 Dpca_a_R = Dpca.
arr[il][im][ia+1];
545 Dpca_m_L = Dpca.
arr[il][im-1][ia];
546 Dpca_m_R = Dpca.
arr[il][im+1][ia];
547 J_ = Jacobian[il][im][ia];
548 J_a_L = Jacobian[il][im][ia-1];
549 J_a_R = Jacobian[il][im][ia+1];
550 J_m_L = Jacobian[il][im-1][ia];
551 J_m_R = Jacobian[il][im+1][ia];
554 PSD_m_r = Dpca.
arr[il][im+1][ia]*Jacobian[il][im+1][ia]*(arr[il][im+1][ia+1] - arr[il][im+1][ia-1])/(alpha.
arr[il][im+1][ia+1] - alpha.
arr[il][im+1][ia-1]);
555 PSD_m_l = Dpca.
arr[il][im-1][ia]*Jacobian[il][im-1][ia]*(arr[il][im-1][ia+1] - arr[il][im-1][ia-1])/(alpha.
arr[il][im-1][ia+1] - alpha.
arr[il][im-1][ia-1]);
556 Lma = 1.0/Jacobian[il][im][ia] * (PSD_m_r - PSD_m_l) / (pc.
arr[il][im+1][ia] - pc.
arr[il][im-1][ia]);
559 PSD_a_r = Dpca.
arr[il][im][ia+1]*Jacobian[il][im][ia+1]*(arr[il][im+1][ia+1] - arr[il][im-1][ia+1])/(pc.
arr[il][im+1][ia+1] - pc.
arr[il][im-1][ia+1]);
560 PSD_a_l = Dpca.
arr[il][im][ia-1]*Jacobian[il][im][ia-1]*(arr[il][im+1][ia-1] - arr[il][im-1][ia-1])/(pc.
arr[il][im+1][ia-1] - pc.
arr[il][im-1][ia-1]);
561 Lam = 1.0/Jacobian[il][im][ia] * (PSD_a_r - PSD_a_l) / (alpha.
arr[il][im][ia+1] - alpha.
arr[il][im][ia-1]);
565 arr[il][im][ia] = arr[il][im][ia] + dt * (Lma + Lam);
571 throw error_msg(
"APPROXIMATION",
"Unknown approximation method for mixed terms");
604 string alpha_lowerBoundaryCondition_calculationType,
605 string alpha_upperBoundaryCondition_calculationType)
619 matr_A.Initialize(alpha.
size, 1, 1, 1);
620 matr_B.Initialize(alpha.
size, 1, 1, 0),
621 matr_C.Initialize(alpha.
size, 1, 1, 0);
639 static Matrix3D<double> Daa_1d(1, 1, alpha.
size), DaaLpp_1d(1, 1, alpha.
size), Jacobian_1d(1, 1, alpha.
size), alpha_1d(1, 1, alpha.
size), L_1d(1, 1, alpha.
size), pc_1d(1, 1, alpha.
size);
641 static Matrix2D<double> alpha_lowerBoundaryCondition_1d(1, 1), alpha_upperBoundaryCondition_1d(1, 1);
645 DiagMatrix::iterator it;
646 for (il = 0; il < L.
size; il++) {
647 for (im = 0; im < pc.
size; im++) {
650 for (ia = 0; ia < alpha.
size; ia++) {
651 slicePSD1D[ia] = arr[il][im][ia];
653 Daa_1d[0][0][ia] = Daa.
arr[il][im][ia];
654 DaaLpp_1d[0][0][ia] = DaaLpp.
arr[il][im][ia];
656 L_1d[0][0][ia] = L.
arr[il][im][ia];
657 pc_1d[0][0][ia] = pc.
arr[il][im][ia];
658 alpha_1d[0][0][ia] = alpha.
arr[il][im][ia];
660 Jacobian_1d[0][0][ia] = Jacobian[il][im][ia];
663 alpha_upperBoundaryCondition_1d[0][0] = alpha_upperBoundaryCondition[il][im];
664 alpha_lowerBoundaryCondition_1d[0][0] = alpha_lowerBoundaryCondition[il][im];
668 matr_A, matr_B, matr_C,
669 L_1d, pc_1d, alpha_1d,
671 Zero_Matrix_2d, Zero_Matrix_2d,
672 Zero_Matrix_2d, Zero_Matrix_2d,
673 alpha_lowerBoundaryCondition_1d, alpha_upperBoundaryCondition_1d,
674 "BCT_CONSTANT_VALUE",
"BCT_CONSTANT_VALUE",
675 "BCT_CONSTANT_VALUE",
"BCT_CONSTANT_VALUE",
676 alpha_lowerBoundaryCondition_calculationType, alpha_upperBoundaryCondition_calculationType,
678 Zero_Matrix, Zero_Matrix,
680 Zero_Matrix, Zero_Matrix,
681 Jacobian_1d, dt, Lpp);
693 for (in = 0; in < alpha.
size; in++) {
696 RHS[in] = matr_C[0][in];
697 for (it = matr_B.begin(); it != matr_B.end(); it++) {
699 if (in + it->first >= 0 && in + it->first < matr_B.total_size) {
700 RHS[in] += it->second[in] * slicePSD1D[in + it->first];
720 for (ia = 0; ia < alpha.
size; ia++) {
721 arr[il][im][ia] = slicePSD1D[ia];
752 string pc_lowerBoundaryCondition_calculationType,
753 string pc_upperBoundaryCondition_calculationType)
767 matr_A.Initialize(1, pc.
size, 1, 1);
768 matr_B.Initialize(1, pc.
size, 1, 0),
769 matr_C.Initialize(1, pc.
size, 1, 0);
786 static Matrix3D<double> Dpcpc_1d(1, pc.
size, 1), DpcpcLpp_1d(1, pc.
size, 1), Jacobian_1d(1, pc.
size, 1), L_1d(1, pc.
size, 1), pc_1d(1, pc.
size, 1), alpha_1d(1, pc.
size, 1);
788 static Matrix2D<double> pc_lowerBoundaryCondition_1d(1, 1), pc_upperBoundaryCondition_1d(1, 1);
792 DiagMatrix::iterator it;
793 for (il = 0; il < L.
size; il++) {
794 for (ia = 0; ia < alpha.
size; ia++) {
797 for (im = 0; im < pc.
size; im++) {
799 slicePSD1D[im] = arr[il][im][ia];
801 Dpcpc_1d[0][im][0] = Dpcpc.
arr[il][im][ia];
802 DpcpcLpp_1d[0][im][0] = DpcpcLpp.
arr[il][im][ia];
804 L_1d[0][im][0] = L.
arr[il][im][ia];
805 pc_1d[0][im][0] = pc.
arr[il][im][ia];
806 alpha_1d[0][im][0] = alpha.
arr[il][im][ia];
808 Jacobian_1d[0][im][0] = Jacobian[il][im][ia];
811 pc_upperBoundaryCondition_1d[0][0] = pc_upperBoundaryCondition[il][ia];
812 pc_lowerBoundaryCondition_1d[0][0] = pc_lowerBoundaryCondition[il][ia];
816 matr_A, matr_B, matr_C,
817 L_1d, pc_1d, alpha_1d,
819 Zero_Matrix_2d, Zero_Matrix_2d,
820 pc_lowerBoundaryCondition_1d, pc_upperBoundaryCondition_1d,
821 Zero_Matrix_2d, Zero_Matrix_2d,
822 "BCT_CONSTANT_VALUE",
"BCT_CONSTANT_VALUE",
823 pc_lowerBoundaryCondition_calculationType, pc_upperBoundaryCondition_calculationType,
824 "BCT_CONSTANT_VALUE",
"BCT_CONSTANT_VALUE",
826 Dpcpc_1d, DpcpcLpp_1d,
827 Zero_Matrix, Zero_Matrix,
828 Zero_Matrix, Zero_Matrix,
829 Jacobian_1d, dt, Lpp);
841 DiagMatrix::iterator it;
842 for (in = 0; in < pc.
size; in++) {
845 RHS[in] = matr_C[0][in];
846 for (it = matr_B.begin(); it != matr_B.end(); it++) {
848 if (in + it->first >= 0 && in + it->first < matr_B.total_size) {
849 RHS[in] += it->second[in] * slicePSD1D[in + it->first];
869 for (im = 0; im < pc.
size; im++) {
870 arr[il][im][ia] = slicePSD1D[im];
900 string L_lowerBoundaryCondition_calculationType,
901 string L_upperBoundaryCondition_calculationType,
902 double tau,
double tauLpp)
916 matr_A.Initialize(1, 1, L.
size, 1);
917 matr_B.Initialize(1, 1, L.
size, 0),
918 matr_C.Initialize(1, 1, L.
size, 0);
937 static Matrix2D<double> L_lowerBoundaryCondition_1d(1, 1), L_upperBoundaryCondition_1d(1, 1);
941 DiagMatrix::iterator it;
942 for (ia = 0; ia < alpha.
size; ia++) {
943 for (im = 0; im < pc.
size; im++) {
946 for (il = 0; il < L.
size; il++) {
948 slicePSD1D[il] = arr[il][im][ia];
950 DLL_1d[il][0][0] = DLL.
arr[il][im][ia];
952 L_1d[il][0][0] = L.
arr[il][im][ia];
956 Jacobian_1d[il][0][0] = Jacobian[il][im][ia];
959 L_upperBoundaryCondition_1d[0][0] = L_upperBoundaryCondition[im][ia];
960 L_lowerBoundaryCondition_1d[0][0] = L_lowerBoundaryCondition[im][ia];
964 matr_A, matr_B, matr_C,
965 L_1d, Zero_Matrix, Zero_Matrix,
967 L_lowerBoundaryCondition_1d, L_upperBoundaryCondition_1d,
968 Zero_Matrix_2d, Zero_Matrix_2d,
969 Zero_Matrix_2d, Zero_Matrix_2d,
970 L_lowerBoundaryCondition_calculationType, L_upperBoundaryCondition_calculationType,
971 "BCT_CONSTANT_VALUE",
"BCT_CONSTANT_VALUE",
972 "BCT_CONSTANT_VALUE",
"BCT_CONSTANT_VALUE",
974 Zero_Matrix, Zero_Matrix,
975 Zero_Matrix, Zero_Matrix,
976 Zero_Matrix, Zero_Matrix,
977 Jacobian_1d, dt, Lpp,
992 DiagMatrix::iterator it;
993 for (in = 0; in < L.
size; in++) {
996 RHS[in] = matr_C[0][in];
997 for (it = matr_B.begin(); it != matr_B.end(); it++) {
999 if (in + it->first >= 0 && in + it->first < matr_B.total_size) {
1000 RHS[in] += it->second[in] * slicePSD1D[in + it->first];
1020 for (il = 0; il < L.
size; il++) {
1021 arr[il][im][ia] = slicePSD1D[il];
1062 string pc_lowerBoundaryCondition_calculationType,
1063 string pc_upperBoundaryCondition_calculationType,
1064 string alpha_lowerBoundaryCondition_calculationType,
1065 string alpha_upperBoundaryCondition_calculationType)
1076 matr_A.Initialize(alpha.
size, pc.
size, 1, 1);
1077 matr_B.Initialize(alpha.
size, pc.
size, 1, 0),
1078 matr_C.Initialize(alpha.
size, pc.
size, 1, 0);
1090 static Matrix1D<double> slicePSD1D(matr_A.total_size), RHS(matr_A.total_size);
1093 static Matrix3D<double> Dpcpc_slice(1, pc.
size, alpha.
size), DpcpcLpp_slice(1, pc.
size, alpha.
size), Daa_slice(1, pc.
size, alpha.
size), DaaLpp_slice(1, pc.
size, alpha.
size), Dpca_slice(1, pc.
size, alpha.
size), DpcaLpp_slice(1, pc.
size, alpha.
size);
1094 static Matrix3D<double> L_slice(1, pc.
size, alpha.
size), pc_slice(1, pc.
size, alpha.
size), alpha_slice(1, pc.
size, alpha.
size), Jacobian_slice(1, pc.
size, alpha.
size);
1096 static Matrix2D<double> pc_lowerBoundaryCondition_slice(1, alpha.
size), pc_upperBoundaryCondition_slice(1, alpha.
size);
1097 static Matrix2D<double> alpha_lowerBoundaryCondition_slice(1, pc.
size), alpha_upperBoundaryCondition_slice(1, pc.
size);
1101 DiagMatrix::iterator it;
1102 for (il = 0; il < L.
size; il++) {
1104 for (im = 0; im < pc.
size; im++) {
1105 for (ia = 0; ia < alpha.
size; ia++) {
1107 in = matr_A.index1d(ia, im);
1108 slicePSD1D[in] = arr[il][im][ia];
1110 Dpcpc_slice[0][im][ia] = Dpcpc.
arr[il][im][ia];
1111 DpcpcLpp_slice[0][im][ia] = DpcpcLpp.
arr[il][im][ia];
1112 Daa_slice[0][im][ia] = Daa.
arr[il][im][ia];
1113 DaaLpp_slice[0][im][ia] = DaaLpp.
arr[il][im][ia];
1114 Dpca_slice[0][im][ia] = Dpca.
arr[il][im][ia];
1115 DpcaLpp_slice[0][im][ia] = DpcaLpp.
arr[il][im][ia];
1117 L_slice[0][im][ia] = L.
arr[il][im][ia];
1118 pc_slice[0][im][ia] = pc.
arr[il][im][ia];
1119 alpha_slice[0][im][ia] = alpha.
arr[il][im][ia];
1121 Jacobian_slice[0][im][ia] = Jacobian[il][im][ia];
1125 for (ia = 0; ia < alpha.
size; ia++) {
1126 pc_upperBoundaryCondition_slice[0][ia] = pc_upperBoundaryCondition[il][ia];
1127 pc_lowerBoundaryCondition_slice[0][ia] = pc_lowerBoundaryCondition[il][ia];
1129 for (im = 0; im < pc.
size; im++) {
1130 alpha_upperBoundaryCondition_slice[0][im] = alpha_upperBoundaryCondition[il][im];
1131 alpha_lowerBoundaryCondition_slice[0][im] = alpha_lowerBoundaryCondition[il][im];
1136 matr_A, matr_B, matr_C,
1137 L_slice, pc_slice, alpha_slice,
1139 Zero_Matrix_2d, Zero_Matrix_2d,
1140 pc_lowerBoundaryCondition_slice, pc_upperBoundaryCondition_slice,
1141 alpha_lowerBoundaryCondition_slice, alpha_upperBoundaryCondition_slice,
1142 "BCT_CONSTANT_VALUE",
"BCT_CONSTANT_VALUE",
1143 pc_lowerBoundaryCondition_calculationType, pc_upperBoundaryCondition_calculationType,
1144 alpha_lowerBoundaryCondition_calculationType, alpha_upperBoundaryCondition_calculationType,
1146 Dpcpc_slice, DpcpcLpp_slice,
1147 Daa_slice, DaaLpp_slice,
1148 Dpca_slice, DpcaLpp_slice,
1149 Jacobian_slice, dt, Lpp);
1158 DiagMatrix::iterator it;
1159 for (im = 0; im < pc.
size; im++) {
1160 for (ia = 0; ia < alpha.
size; ia++) {
1161 in = matr_B.index1d(ia, im);
1164 RHS[in] = matr_C[0][in];
1165 for (it = matr_B.begin(); it != matr_B.end(); it++) {
1167 if (in + it->first >= 0 && in + it->first < matr_B.total_size) {
1168 RHS[in] += it->second[in] * slicePSD1D[in + it->first];
1181 if (PSD_parameters.solutionMethod ==
"SM_Gauss") {
1183 gauss(matr_A, RHS, slicePSD1D, n);
1189 }
else if (PSD_parameters.solutionMethod ==
"SM_Relaxation") {
1192 Output::echo(5,
"Inverting A matrix (relaxation)... ");
1197 }
else if (PSD_parameters.solutionMethod ==
"SM_Lapack") {
1202 Lapack(matr_A, RHS, slicePSD1D);
1205 }
else if (PSD_parameters.solutionMethod ==
"SM_GMRES") {
1211 gmres_wrapout(matr_A, RHS, slicePSD1D, PSD_parameters.GMRES_parameters);
1214 }
else if (PSD_parameters.solutionMethod ==
"SM_GMRES2") {
1220 gmres2_wrapout(matr_A, RHS, slicePSD1D, PSD_parameters.GMRES_parameters);
1224 }
else if (PSD_parameters.solutionMethod ==
"SM_Tridiag") {
1226 throw error_msg(
"SOLUTIONERR",
"Tridiagonal method can't be used");
1228 throw error_msg(
"SOLUTIONERR",
"Unknown solution method");
1232 for (im = 0; im < pc.
size; im++) {
1233 for (ia = 0; ia < alpha.
size; ia++) {
1234 in = matr_A.index1d(ia, im);
1235 arr[il][im][ia] = slicePSD1D[in];
1278 string pc_lowerBoundaryCondition_calculationType,
1279 string pc_upperBoundaryCondition_calculationType,
1280 string alpha_lowerBoundaryCondition_calculationType,
1281 string alpha_upperBoundaryCondition_calculationType)
1292 matr_A.Initialize(alpha.
size, pc.
size, 1, 1);
1293 matr_B.Initialize(alpha.
size, pc.
size, 1, 0),
1294 matr_C.Initialize(alpha.
size, pc.
size, 1, 0);
1309 static Matrix1D<double> slicePSD1D(matr_A.total_size), RHS(matr_A.total_size);
1312 static Matrix3D<double> Dpcpc_slice(1, pc.
size, alpha.
size), DpcpcLpp_slice(1, pc.
size, alpha.
size), Daa_slice(1, pc.
size, alpha.
size), DaaLpp_slice(1, pc.
size, alpha.
size), Dpca_slice(1, pc.
size, alpha.
size), DpcaLpp_slice(1, pc.
size, alpha.
size);
1313 static Matrix3D<double> L_slice(1, pc.
size, alpha.
size), pc_slice(1, pc.
size, alpha.
size), alpha_slice(1, pc.
size, alpha.
size), Jacobian_slice(1, pc.
size, alpha.
size);
1315 static Matrix2D<double> pc_lowerBoundaryCondition_slice(1, alpha.
size), pc_upperBoundaryCondition_slice(1, alpha.
size);
1316 static Matrix2D<double> alpha_lowerBoundaryCondition_slice(1, pc.
size), alpha_upperBoundaryCondition_slice(1, pc.
size);
1320 DiagMatrix::iterator it;
1321 for (il = 0; il < L.
size; il++) {
1323 for (im = 0; im < pc.
size; im++) {
1324 for (ia = 0; ia < alpha.
size; ia++) {
1326 in = matr_A.index1d(ia, im);
1327 slicePSD1D[in] = arr[il][im][ia];
1329 Dpcpc_slice[0][im][ia] = Dpcpc.
arr[il][im][ia];
1330 DpcpcLpp_slice[0][im][ia] = DpcpcLpp.
arr[il][im][ia];
1331 Daa_slice[0][im][ia] = Daa.
arr[il][im][ia];
1332 DaaLpp_slice[0][im][ia] = DaaLpp.
arr[il][im][ia];
1333 Dpca_slice[0][im][ia] = Dpca.
arr[il][im][ia];
1334 DpcaLpp_slice[0][im][ia] = DpcaLpp.
arr[il][im][ia];
1336 L_slice[0][im][ia] = L.
arr[il][im][ia];
1337 pc_slice[0][im][ia] = pc.
arr[il][im][ia];
1338 alpha_slice[0][im][ia] = alpha.
arr[il][im][ia];
1340 Jacobian_slice[0][im][ia] = Jacobian[il][im][ia];
1344 for (ia = 0; ia < alpha.
size; ia++) {
1345 pc_upperBoundaryCondition_slice[0][ia] = pc_upperBoundaryCondition[il][ia];
1346 pc_lowerBoundaryCondition_slice[0][ia] = pc_lowerBoundaryCondition[il][ia];
1348 for (im = 0; im < pc.
size; im++) {
1349 alpha_upperBoundaryCondition_slice[0][im] = alpha_upperBoundaryCondition[il][im];
1350 alpha_lowerBoundaryCondition_slice[0][im] = alpha_lowerBoundaryCondition[il][im];
1355 matr_A, matr_B, matr_C,
1356 L_slice, pc_slice, alpha_slice,
1358 Zero_Matrix_2d, Zero_Matrix_2d,
1359 pc_lowerBoundaryCondition_slice, pc_upperBoundaryCondition_slice,
1360 alpha_lowerBoundaryCondition_slice, alpha_upperBoundaryCondition_slice,
1361 "BCT_CONSTANT_VALUE",
"BCT_CONSTANT_VALUE",
1362 pc_lowerBoundaryCondition_calculationType, pc_upperBoundaryCondition_calculationType,
1363 alpha_lowerBoundaryCondition_calculationType, alpha_upperBoundaryCondition_calculationType,
1365 Dpcpc_slice, DpcpcLpp_slice,
1366 Daa_slice, DaaLpp_slice,
1367 Dpca_slice, DpcaLpp_slice,
1368 Jacobian_slice, dt, Lpp);
1377 DiagMatrix::iterator it;
1378 for (im = 0; im < pc.
size; im++) {
1379 for (ia = 0; ia < alpha.
size; ia++) {
1380 in = matr_B.index1d(ia, im);
1383 RHS[in] = matr_C[0][in];
1384 for (it = matr_B.begin(); it != matr_B.end(); it++) {
1386 if (in + it->first >= 0 && in + it->first < matr_B.total_size) {
1387 RHS[in] += it->second[in] * slicePSD1D[in + it->first];
1402 int modelMatrixLineNumber;
1404 double sum, error, error_max;
1405 static bool recalculate_A =
true;
1406 if (PSD_parameters.solutionMethod ==
"SM_Gauss") {
1408 if (recalculate_A) {
1420 for (it = matr_A.begin(); it != matr_A.end(); it++) {
1421 for (modelMatrixLineNumber = 0; modelMatrixLineNumber < n; modelMatrixLineNumber++) {
1424 if (modelMatrixLineNumber + it->first >= 0 && modelMatrixLineNumber + it->first < n) {
1428 A[modelMatrixLineNumber][modelMatrixLineNumber + it->first] = it->second[modelMatrixLineNumber];
1438 recalculate_A =
true;
1443 slicePSD1D[n-1] = B[n-1] / A[n-1][n-1];
1444 for (k = n - 2; k >= 0; k--) {
1446 for (l = k + 1; l < n; l++)
1447 sum += A[k][l] * slicePSD1D[l];
1448 slicePSD1D[k] = B[k] - sum;
1456 for (in = 0; in < n; in++) {
1458 for (it = matr_A.begin(); it != matr_A.end(); it++) {
1460 if (in + it->first >= 0 && in + it->first < n) {
1461 error -= it->second[in] * slicePSD1D[in + it->first];
1464 error_max =
VF::max( error_max, error );
1469 }
else if (PSD_parameters.solutionMethod ==
"SM_Relaxation") {
1472 Output::echo(5,
"Inverting A matrix (relaxation)... ");
1477 }
else if (PSD_parameters.solutionMethod ==
"SM_Lapack") {
1482 Lapack(matr_A, RHS, slicePSD1D);
1485 }
else if (PSD_parameters.solutionMethod ==
"SM_GMRES") {
1491 gmres_wrapout(matr_A, RHS, slicePSD1D, PSD_parameters.GMRES_parameters);
1494 }
else if (PSD_parameters.solutionMethod ==
"SM_GMRES2") {
1500 gmres2_wrapout(matr_A, RHS, slicePSD1D, PSD_parameters.GMRES_parameters);
1504 }
else if (PSD_parameters.solutionMethod ==
"SM_Tridiag") {
1506 throw error_msg(
"SOLUTIONERR",
"Tridiagonal method can't be used");
1508 throw error_msg(
"SOLUTIONERR",
"Unknown solution method");
1512 for (im = 0; im < pc.
size; im++) {
1513 for (ia = 0; ia < alpha.
size; ia++) {
1514 in = matr_A.index1d(ia, im);
1515 arr[il][im][ia] = slicePSD1D[in];
1535 if (arr[il][im][ia] != arr[il][im][ia] || arr[il][im][ia] > maxNum || arr[il][im][ia] < -maxNum) {
1536 throw error_msg(
"INTERPOLATION_ERROR",
"Interpolation error at [%d][%d][%d]", il, im, ia);
1551 if ( arr[im] != arr[im] || arr[im] > maxNum || arr[im] < -maxNum) {
1552 throw error_msg(
"INTERPOLATION_ERROR",
"Interpolation error at [%d][%d][%d]", il, im, ia);
1561 double lowerBoundary1d=0, upperBoundary1d=0;
1565 if (interpolationParameters_structure.
type ==
"IT_NONE") {
1567 }
else if (oldGrid.
type == newGrid.
type) {
1570 for (il = 0; il < oldGrid.
L.
size; il++) {
1571 for (ia = 0; ia < oldGrid.
alpha.
size; ia++) {
1572 for (im = 0; im < oldGrid.
pc.
size; im++) {
1574 arr[il][im][ia] = otherPSD.
arr[il][im][ia];
1591 for (il = 0; il < oldGrid.
L.
size; il++) {
1592 for (ia = 0; ia < oldGrid.
alpha.
size; ia++) {
1593 for (im = 0; im < oldGrid.
pc.
size; im++) {
1595 if (newGrid.
L.
arr[il][im][ia] != oldGrid.
L.
arr[il][im][ia] || newGrid.
alpha.
arr[il][im][ia] != oldGrid.
alpha.
arr[il][im][ia]) {
1597 throw error_msg(
"INTERPOLATION_ERROR",
"Interpolation error: grids must have the same L and Alpha");
1605 if (interpolationParameters_structure.
useLog ==
"Log_10") {
1608 matrix_array_1d_old[im] = log10(otherPSD.
arr[il][im][ia]);
1616 lowerBoundary1d = log10(
VF::max(newGrid_pc_lowerBoundaryCondition[il][ia],
VC::zero_f));
1617 upperBoundary1d = log10(
VF::max(newGrid_pc_upperBoundaryCondition[il][ia],
VC::zero_f));
1619 else if (interpolationParameters_structure.
useLog ==
"Log_2") {
1622 matrix_array_1d_old[im] = log(otherPSD.
arr[il][im][ia])/log(2.0);
1625 matrix_array_1d_old[im] = log(
VC::zero_f)/log(2.0);
1632 lowerBoundary1d = log(
VF::max(newGrid_pc_lowerBoundaryCondition[il][ia],
VC::zero_f)) / log(2.0);
1633 upperBoundary1d = log(
VF::max(newGrid_pc_upperBoundaryCondition[il][ia],
VC::zero_f)) / log(2.0);
1635 else if (interpolationParameters_structure.
useLog ==
"Log_E") {
1638 matrix_array_1d_old[im] = log(otherPSD.
arr[il][im][ia]);
1648 lowerBoundary1d = log(
VF::max(newGrid_pc_lowerBoundaryCondition[il][ia],
VC::zero_f));
1649 upperBoundary1d = log(
VF::max(newGrid_pc_upperBoundaryCondition[il][ia],
VC::zero_f));
1651 else if (interpolationParameters_structure.
useLog ==
"NoLog") {
1652 matrix_array_1d_old[im] = otherPSD.
arr[il][im][ia];
1653 lowerBoundary1d =
VF::max(newGrid_pc_lowerBoundaryCondition[il][ia],
VC::zero_f);
1654 upperBoundary1d =
VF::max(newGrid_pc_upperBoundaryCondition[il][ia],
VC::zero_f);
1657 throw error_msg(
"INTERPOLATION_ERROR",
"Usage of log(PSD) is not define properly: %s", interpolationParameters_structure.
useLog.c_str());
1666 grid_1d_old[im] = oldGrid.
pc.
arr[il][im][ia];
1667 grid_1d_new[im] = newGrid.
pc.
arr[il][im][ia];
1672 if (interpolationParameters_structure.
type ==
"IT_LINEAR") {
1674 matrix_array_1d_new.
Interpolate(matrix_array_1d_old, grid_1d_old, grid_1d_new);
1675 }
else if (interpolationParameters_structure.
type ==
"IT_POLYNOMIAL") {
1677 matrix_array_1d_new.
Polilinear(matrix_array_1d_old, grid_1d_old, grid_1d_new, lowerBoundary1d, upperBoundary1d);
1678 }
else if (interpolationParameters_structure.
type ==
"IT_SPLINE") {
1680 matrix_array_1d_new.
Spline(matrix_array_1d_old, grid_1d_old, grid_1d_new, lowerBoundary1d, upperBoundary1d, interpolationParameters_structure.
linearSplineCoef, interpolationParameters_structure.
maxSecondDerivative);
1681 }
else if (interpolationParameters_structure.
type ==
"IT_COPY") {
1683 matrix_array_1d_new = matrix_array_1d_old;
1685 throw error_msg(
"INTERPOLATION_TYPE_ERROR",
"Unknown interpolation type");
1688 for (im = 0; im < newGrid.
pc.
size; im++) {
1695 if (interpolationParameters_structure.
useLog ==
"Log_10") arr[il][im][ia] = pow(10, matrix_array_1d_new[im]);
1696 else if (interpolationParameters_structure.
useLog ==
"Log_2") arr[il][im][ia] = pow(2, matrix_array_1d_new[im]);
1697 else if (interpolationParameters_structure.
useLog ==
"Log_E") arr[il][im][ia] = pow(
VC::exp, matrix_array_1d_new[im]);
1698 else arr[il][im][ia] = matrix_array_1d_new[im];
1722 output_without_grid_file <<
"ZONE T=\"" <<
internal << time <<
"\" I=" << arr.size_z <<
", J=" << arr.size_y <<
", K=" << arr.size_x << endl;
1723 for (il = 0; il < arr.size_x; il++) {
1724 for (im = 0; im < arr.size_y; im++) {
1725 for (ia = 0; ia < arr.size_z; ia++) {
1726 output_without_grid_file << arr[il][im][ia] << endl;
1745 double Load_L=-1, Load_epc=-1, Load_alpha=-1;
1749 ifstream input_f(filename);
1750 if (input_f != NULL && !input_f.eof()) {
1752 getline(input_f, inBuf);
1753 getline(input_f, inBuf);
1754 for (il = 0; il < L.
size; il++) {
1755 for (im = 0; im < epc.
size; im++) {
1756 for (ia = 0; ia < alpha.
size; ia++) {
1757 if (withGrid) input_f >> Load_L >> Load_epc >> Load_alpha;
1758 input_f >> arr[il][im][ia];
1760 if (fabs(Load_L - L.
arr[il][im][ia]) > err || fabs(Load_epc - epc.
arr[il][im][ia]) > err || fabs(Load_alpha - alpha.
arr[il][im][ia]) > err) {
1761 throw error_msg(
"INITIAL_PSD_LOAD_GRID_ERR",
"Loading %s: grid mismatch.\nLoaded: L=%e, epc=%e, alpha=%e\nGrid: L=%e, epc=%e, alpha=%e\n", filename, Load_L, Load_epc, Load_alpha, L.
arr[il][im][ia], epc.
arr[il][im][ia], alpha.
arr[il][im][ia]);
1767 throw error_msg(
"INITIAL_PSD_ERROR",
"Error reading PSD input file %s.\n", filename);
1792 double p1, Lmax, Lmin, scale_factor, Lgrid;
1793 int Lsize_7, addGrid;
1801 int ia_90=alpha.
size-1;
1802 while (alpha.
arr[L.
size-1][0][ia_90] >
VC::pi/2 && ia_90 > 0) ia_90--;
1806 if (L.
arr[L.
size-1][0][ia_90] < 7) {
1812 Lmax = L.
arr[L.
size-1][0][ia_90];
1813 Lmin = L.
arr[0][0][ia_90];
1816 Lmin = L.
arr[0][0][0];
1819 Lgrid = (Lmax-Lmin)/(L.
size-1);
1820 Lsize_7 = (7.0-Lmin)/Lgrid + 1;
1821 addGrid = Lsize_7 - L.
size;
1827 for (il = 0; il < L.
size; il++) {
1828 L_1d[il] = L.
arr[il][0][ia_90];
1830 for (il = L.
size; il < Lsize_7; il++) {
1831 L_1d[il] = Lmin+Lgrid*il;
1838 for (il = 0; il < L.
size; il++) L_1d[il] = L.
arr[il][0][ia_90];
1842 steady_state(fL, tau, Kp, Lsize_7, L_1d, fb_out, fb_in);
1846 scale_factor = fb_out/fL[L.
size-1];
1850 if (L_1d[L_1d.
size_x-1] > 7) {
1851 int il_7 = L.
size-1;
1852 while (L_1d[il_7] > 7 && il_7 > 0) il_7--;
1853 Output::echo(5,
"Steady state normalized by fL(L[%d] = %f) = %f.\n", il_7, L_1d[il_7], fL(il_7));
1860 for (
int i = 0; i < fL.
size_x; i++)
if (fL[i] < min_f) fL[i] = min_f;
1873 for (im = 0; im < pc.
size; im++) {
1874 for (il = 0; il < L.
size; il++) {
1875 for (ia = 0; ia < alpha.
size; ia++) {
1886 while (alpha.
arr[il][im][ia_90] >
VC::pi/2 && ia_90 > 0) ia_90--;
1897 arr[il][im][ia] = arr[il][im][ia] * pow(sinh(alpha.
arr[il][im][ia]/
VF::alc(L.
arr[il][im][ia])),2)/pow(sinh(1.0),2);
1928 double Lmax, Lmin, Lgrid, Lpp;
1929 int Lsize_7, addGrid;
1930 double tau_p, tau_c, alpha1, p1;
1942 int ia_90 = alpha.
size - 1;
1943 while (alpha.
arr[L.
size - 1][0][ia_90] >
VC::pi / 2 && ia_90 > 0) ia_90--;
1947 if (L.
arr[L.
size - 1][0][ia_90] < 7.0) {
1950 Lmax = L.
arr[L.
size - 1][0][ia_90];
1951 Lmin = L.
arr[0][0][ia_90];
1952 Lgrid = (Lmax - Lmin) / (L.
size - 1);
1953 Lsize_7 = (7.0 - Lmin) / Lgrid + 1;
1954 addGrid = Lsize_7 - L.
size;
1961 for (il = 0; il < L.
size; il++) {
1962 for (im = 0; im < pc.
size; im++) {
1963 L_1d_i[il] = L.
arr[il][0][ia_90];
1964 pc_i[il][im] = pc.
arr[il][im][ia_90];
1969 for (il = L.
size; il < Lsize_7; il++) {
1970 for (im = 0; im < pc.
size; im++) {
1971 L_1d_i[il] = Lmin + Lgrid * il;
1985 for (il = 0; il < L.
size; il++) {
1986 for (im = 0; im < pc.
size; im++) {
1987 L_1d_i[il] = L.
arr[il][0][ia_90];
1988 pc_i[il][im] = pc.
arr[il][im][ia_90];
1995 Lpp = (5.6 - 0.46 * Kp);
2001 for (im = 0; im < pc.
size; im++) {
2002 for (il = 0; il < Lsize_7; il++) {
2005 alpha1 = alpha.
arr[L.
size - 1][0][alpha.
size - 1];
2007 if(PSD_parameters.usetauLpp_SteadyState ==
"coulomb") {
2009 }
else if(PSD_parameters.usetauLpp_SteadyState ==
"precipitation") {
2011 }
else if (PSD_parameters.usetauLpp_SteadyState ==
"combined") {
2014 tau_0[il] = 1.0 / (1.0/tau_p+1.0/tau_c);
2029 for (il = 0; il < Lsize_7; il++) {
2030 if (fL_i[il] < min_f) {
2031 fLm[il][im] = min_f;
2033 fLm[il][im] = fL_i[il];
2039 for (im = 0; im < pc.
size; im++) {
2040 scale_factor[im] = fb_out / fLm[L.
size - 1][im];
2043 for (im = 0; im < pc.
size; im++) {
2044 for (il = 0; il < L.
size; il++) {
2045 for (ia = 0; ia < alpha.
size; ia++) {
2054 int ia_90 = alpha.
size - 1;
2055 while (alpha.
arr[il][im][ia_90] >
VC::pi / 2 && ia_90 > 0)
2059 L.
arr[il][im][ia_90]));
2060 arr[il][im][ia] = sin(alpha.
arr[il][im][ia]) * scale_factor[im]
2062 VF::Kfunc(p1), J_L7_function) * pow(1.0 / p1, 2);
2067 L.
arr[il][im][ia])) {
2068 arr[il][im][ia] = arr[il][im][ia] * pow(
2070 L.
arr[il][im][ia])), 2) / pow(sinh(1.0), 2);
2116 for (i = 0; i <= nx-1; i++) {
2118 Dxx[i] = pow(10, 0.506*Kp-9.325)*pow(L[i],10) +
VF::Dfe(alpha,Ke[i],L[i],Kp);
2123 for (i = 1; i < nx-1; i++) {
2124 alfa[i] = 1.0/((L[i+1] - L[i-1])/2)*pow(L[i],2);
2126 C_1[i] = ((Dxx[i] + Dxx[i-1])/2)/(L[i] - L[i-1])/(pow(((L[i] + L[i-1])/2), 2));
2127 C_2[i] = ((Dxx[i+1] + Dxx[i])/2)/(L[i+1] - L[i])/(pow(((L[i+1] + L[i])/2), 2));
2133 for (i = 1; i < nx-1; i++) {
2134 A[i] = alfa[i] * C_1[i];
2135 B[i] = -alfa[i] * (C_1[i] + C_2[i]) - C_3[i];
2136 C[i] = alfa[i] * C_2[i];
2139 for (i = 0; i < nx; i++) D[i] = 0;;
2162 tridag(&A[0], &B[0], &C[0], &D[0], &f[0], nx);
2192 for (il = 0; il < L.
size; il++) L_1d[il] = L.
arr[il][0][alpha.
size-1];
2194 for (
int i = 0; i < fL.
size_x; i++)
if (fL[i] < min_f) fL[i] = min_f;
2201 for (im = 0; im < pc.
size; im++) {
2202 for (il = 0; il < L.
size; il++) {
2203 for (ia = 0; ia < alpha.
size; ia++) {
2205 arr[il][im][ia] = sin(alpha.
arr[il][im][ia]) * fL[il] * L_UpperBoundaryCondition[im][ia] * pow(1.0/p1,2);
2242 ifstream input_f(filename);
2243 if (input_f != NULL && !input_f.eof()) {
2244 for (k = 0; k <= lmax; k++) {
2245 input_f >> orbit[k] >> dayno[k] >> eventtime[k];
2246 for (j = 0; j <= jmax; j++) {
2247 input_f >> i >> energy[j] >> jperp[k][j];
2248 energy[j] = energy[j]/1000;
2252 throw error_msg(
"INITIAL_PSD_ERROR",
"Error reading initial conditions input file %s.\n", filename);
2260 for (i = 0; i <= jmax; i++) jcopy[i] = jperp[2][i];
2279 j_0m.Spline(jcopy, energy, epc1D, jcopy[0], log10(
VC::zero_f), 0, 1);
2280 for (i = 0; i < pc.
size; i++) {
2290 for (i = 0; i < jmax+2; i++) {
2291 jcopy[i] = log10(jcopy[i]);
2293 j_0m.Spline(jcopy, energy, epc1D, jcopy[0], log10(
VC::zero_f), 0, 1);
2296 for (i = 0; i < pc.
size; i++) {
2297 j_0m[i] = pow(10, j_0m[i]);
2301 for (i = 0; i < pc.
size; i++) {
2303 f_0m[i] = j_0m[i]/pow(pc.
arr[0][i][0], 2);
2307 j_0m.writeToFile(
"./Debug/j_0m.plt", epc1D);
2312 for (il = 0; il < L.
size; il++) {
2313 for (im = 0; im < pc.
size; im++) {
2314 for (ia = 0; ia < alpha.
size; ia++) {
2316 arr[il][im][ia] = f_0m[im] * sin( alpha.
arr[il][im][ia] );
2318 arr[il][im][ia] = arr[il][im][ia] * pow(sinh(alpha.
arr[il][im][ia]/
VF::alc(L.
arr[il][im][ia])),2)/pow(sinh(1.0),2);
2337 double Ke, alf, flux;
2339 for (il = 0; il < L.
size; il++) {
2340 for (im = 0; im < pc.
size; im++) {
2341 for (ia = 0; ia < alpha.
size; ia++) {
2344 alf = alpha.
arr[il][im][ia];
2345 flux =
exp(-(Ke-0.2)/0.1);
2347 arr[il][im][ia] = flux / pow(pc.
arr[0][im][0], 2);
2383 for (i = 0; i <= nx-1; i++) {
2385 Dxx[i] = pow(10, 0.506*Kp-9.325)*pow(L[i],10);
2388 for (i = 1; i < nx-1; i++) {
2389 alfa[i] = 1.0/((L[i+1] - L[i-1])/2)*pow(L[i],2);
2391 C_1[i] = ((Dxx[i] + Dxx[i-1])/2)/(L[i] - L[i-1])/(pow(((L[i] + L[i-1])/2), 2));
2392 C_2[i] = ((Dxx[i+1] + Dxx[i])/2)/(L[i+1] - L[i])/(pow(((L[i+1] + L[i])/2), 2));
2398 for (i = 1; i < nx-1; i++) {
2399 A[i] = alfa[i] * C_1[i];
2400 B[i] = -alfa[i] * (C_1[i] + C_2[i]) - C_3[i];
2401 C[i] = alfa[i] * C_2[i];
2404 for (i = 0; i < nx; i++) D[i] = 0;;
2427 tridag(&A[0], &B[0], &C[0], &D[0], &f[0], nx);
bool initialized
Flag, equal true if initialized.
GridElement alpha
grid elements to be used
bool useLmp
magnetopause position flag
Matrix3D< double > arr
array of grid points
void Spline(Matrix1D< T > &old_function, Matrix1D< T > &old_grid, Matrix1D< T > &new_grid, double lb, double ub, double lin_spline_coef=0, double max_second_der=0)
void Diffusion_L(double dt, double Lpp, DiffusionCoefficient &DLL, GridElement &L, GridElement &pc, GridElement &alpha, Matrix3D< double > Jacobian, Matrix2D< double > lowerBoundaryCondition, Matrix2D< double > upperBoundaryCondition, string lowerBoundaryCondition_calculationType, string upperBoundaryCondition_calculationType, double tau, double tauLpp)
void Interpolate(PSD &otherPSD, Parameters_structure::Interpolation interpolationParameters, Grid &oldGrid, Grid &newGrid, Matrix2D< double > newGrid_pc_lowerBoundaryCondition, Matrix2D< double > newGrid_pc_upperBoundaryCondition)
Interpolation function.
Interpolation parameters structure
void Polilinear(Matrix1D< T > &old_function, Matrix1D< T > &old_grid, Matrix1D< T > &new_grid, double lb, double ub)
void SourcesAndLosses(Parameters_structure parameters, GridElement &L, GridElement &pc, GridElement &alpha, AdditionalSourcesAndLosses &SL, double dt, double Lpp, double tau, double tauLpp, double Kp)
double initial_PSD_inner_psd
Inner PSD boundary value for steady state, default 0.
Array of values of coordinate axes.
void Load_initial_f_file(GridElement &L, GridElement &pc, GridElement &alpha, const char *filename, bool with_grid)
void Lapack(DiagMatrix &A, Matrix1D< double > &B, Matrix1D< double > &X)
Calculates any additional sources and losses due to magnetopause.
string output_PSD_fileName4D
File name for psd output.
double max(double v1, double v2)
Return maximum.
void steady_state(Matrix1D< double > &f, double tau, double Kp, int nx, Matrix1D< double > &L, double f_bnd_out, double f_bnd_in)
void steady_state_two_zone(Matrix1D< double > &f, Matrix1D< double > &tau, double Kp, double alpha, Matrix1D< double > &Ke, int nx, Matrix1D< double > &L, double f_bnd_out, double f_bnd_in)
void Diffusion_alpha(double dt, double Lpp, DiffusionCoefficient &Daa, DiffusionCoefficient &DaaLpp, GridElement &L, GridElement &pc, GridElement &alpha, Matrix3D< double > Jacobian, Matrix2D< double > alpha_lowerBoundaryCondition, Matrix2D< double > alpha_upperBoundaryCondition, string alpha_lowerBoundaryCondition_calculationType, string alpha_upperBoundaryCondition_calculationType)
Holds diffusion coefficient matrix and routines to load and calculate it.
string usetauLpp
key word to define if and how to use tau/tauLpp
double initial_PSD_Kp0
Kp for steady state?
double Precipitation(double Kp, double L, double energy)
Empirical electron lifetime.
void Diffusion_pc_alpha(double dt, double Lpp, DiffusionCoefficient &Dpcpc, DiffusionCoefficient &DpcpcLpp, DiffusionCoefficient &Daa, DiffusionCoefficient &DaaLpp, DiffusionCoefficient &Dpca, DiffusionCoefficient &DpcaLpp, GridElement &L, GridElement &pc, GridElement &alpha, Matrix3D< double > Jacobian, Matrix2D< double > pc_lowerBoundaryCondition, Matrix2D< double > pc_upperBoundaryCondition, Matrix2D< double > alpha_lowerBoundaryCondition, Matrix2D< double > alpha_upperBoundaryCondition, string pc_lowerBoundaryCondition_calculationType, string pc_upperBoundaryCondition_calculationType, string alpha_lowerBoundaryCondition_calculationType, string alpha_upperBoundaryCondition_calculationType)
void AllocateMemory(int size_x, int size_y)
void Create_Initial_PSD(Parameters_structure::PSD parameters, Grid &grid, BoundaryCondition L_UpperBoundaryCondition)
Create inital PSD (Steady State)
void Output_without_grid(double time)
void gmres_wrapout(CalculationMatrix &A, Matrix1D< double > &B, Matrix1D< double > &X, Parameters_structure::PSD::GMRES_parameters_structure GMRES_parameters)
Matrix2D< double > arr
2D matrix of boundary conditions
double initial_PSD_tauSteadyState
Tau for steady state, if we are calculation initial values as steady state solution.
Parameters_structure::GridElement GridElement_parameters
parameters for grid element
double Outer_Boundary_Flux_at_L7(double E7, string J_L7_function)
Call various flux model of the outer boundary.
static const double exp
exponential
string initial_PSD_Type
Tells us where to get initial PSD values. Check StrToVal(string input, InitialPSDTypes &place) for kn...
void echo(int msglvl, const char *format,...)
Basic function! Write the message to the file.
void AllocateMemory(int x_size)
Matrix3D< double > arr
array of sources and losses values
Makes operations with PSD (Phase Space Density (PSD).) (like, diffusion).
void Load_initial_f_two_zone(GridElement &L, GridElement &pc, GridElement &alpha, double tau, double tauLpp, double Kp, double min_psd=1.e-99, string J_L7_function="J_L7", double fb_out=1, double fb_in=0)
two zone initial profile (kckim)
string type
Interpolation type: spline, linear etc.
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[]...
struct Parameters_structure::SL_structure SL
instance of the struct
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.
void Load_initial_f_maxwell(GridElement &L, GridElement &pc, GridElement &alpha)
double maxSecondDerivative
void Load_initial_f_from_outer_L(GridElement &L, GridElement &pc, GridElement &alpha, double tau, double Kp, Matrix2D< double > L_UpperBoundaryCondition, double min_psd=1.e-99, double fb_out=1, double fb_in=0)
Calculate initial PSD from steady state using boundary conditions. Load initial PSD from outer L...
void DiffusionMixTermExplicit(double dt, double Lpp, DiffusionCoefficient &Dpca, DiffusionCoefficient &DpcaLpp, GridElement &L, GridElement &pc, GridElement &alpha, Matrix3D< double > Jacobian, Matrix2D< double > pc_lowerBoundaryCondition, Matrix2D< double > pc_upperBoundaryCondition, Matrix2D< double > alpha_lowerBoundaryCondition, Matrix2D< double > alpha_upperBoundaryCondition, string pc_lowerBoundaryCondition_calculationType, string pc_upperBoundaryCondition_calculationType, string alpha_lowerBoundaryCondition_calculationType, string alpha_upperBoundaryCondition_calculationType)
double B(double lambda, double L)
void Diffusion_pc(double dt, double Lpp, DiffusionCoefficient &Dpcpc, DiffusionCoefficient &DpcpcLpp, GridElement &L, GridElement &pc, GridElement &alpha, Matrix3D< double > Jacobian, Matrix2D< double > pc_lowerBoundaryCondition, Matrix2D< double > pc_upperBoundaryCondition, string pc_lowerBoundaryCondition_calculationType, string pc_upperBoundaryCondition_calculationType)
Main parameters structure that holds smaller structures for individual parameters.
double Coulomb(double L, double energy)
Coulomb scattering electron lifetime.
int size
size of grid element
void Load_initial_f(GridElement &L, GridElement &pc, GridElement &alpha, double tau, double Kp, double min_f=1.e-99, string J_L7_function="J_L7", double fb_out=1, double fb_in=0, bool using_ip_90=true)
double initial_PSD_some_constant_value
Some psd value. Used as initial PSD value everywhere if initial_PSD_Type = IPSDT_CONSTANT or as minim...
Matrix3D< double > arr
array of diffusion coefficients
double Dfe(double alpha, double Ke, double L, double Kp)
Electrostatic radial diffusion coefficient.
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.
double initial_PSD_tauLppSteadyState
Tau for steady state, if we are calculation initial values as steady state solution.
double B(double Lparam)
Dipole magnetic field.
void writeToFile(string filename)
void Diffusion_pc_alpha_KC(double dt, double Lpp, DiffusionCoefficient &Dpcpc, DiffusionCoefficient &DpcpcLpp, DiffusionCoefficient &Daa, DiffusionCoefficient &DaaLpp, DiffusionCoefficient &Dpca, DiffusionCoefficient &DpcaLpp, GridElement &L, GridElement &pc, GridElement &alpha, Matrix3D< double > Jacobian, Matrix2D< double > pc_lowerBoundaryCondition, Matrix2D< double > pc_upperBoundaryCondition, Matrix2D< double > alpha_lowerBoundaryCondition, Matrix2D< double > alpha_upperBoundaryCondition, string pc_lowerBoundaryCondition_calculationType, string pc_upperBoundaryCondition_calculationType, string alpha_lowerBoundaryCondition_calculationType, string alpha_upperBoundaryCondition_calculationType)
Matrix3D< double > arr
array of PSD values
void Interpolate(Matrix1D< T > &old_function, Matrix1D< T > &old_grid, Matrix1D< T > &new_grid)
Phase Space Density class.
GridElement L
grid elements to be used
double alc(double L)
Loss cone calculations.
void checkInf_3D(Matrix3D< double > arr, int il, int im, int ia, double maxNum=1e99)
Holds upper and lower boundary conditions.
string initial_PSD_fileName
File name, if we need to load initial values from file.
string output_PSD_folderName
Folder name for psd output.
Parameters_structure parameters
Parameters structure, with all parameters from the parameters.ini file. The default parameters define...
string useLog
Flag, if we should interpolate logarithm of the values.
GridElement epc
grid elements to be used
Error message - stack of single_errors.
void gauss_solve(double *a, double *b, int n)
Gauss inversion.
string initial_PSD_J_L7_function
Name of the function for flux at L=7. Can be 'J_L7' or 'J_L7_corrected'.
double initial_PSD_outer_psd
Outer PSD boundary value for steady state, default 1.
Computational grid composed of 3 different GridElement.
PSD parameters structure.
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)
void checkInf_1D(Matrix1D< double > arr, int il, int im, int ia, double maxNum=1e99)
static const double zero_f
Minimum value for PSD.
double Kfunc(double pc)
Computation of Kinetic energy from given momentum .
void Load_initial_f_2d(GridElement &L, GridElement &pc, GridElement &alpha, const char *filename)
GridElement pc
grid elements to be used
double G(double alpha)
G-function.