00001
00015
00016
00017 #include "PSD.h"
00018
00019 #include <math.h>
00020 #include "../VariousFunctions/variousFunctions.h"
00021 #include "../Exceptions/error.h"
00022
00023 #include <iostream>
00024 #include <string>
00025 #include <ctime>
00026
00027 using namespace std;
00028
00031 PSD::~PSD()
00032 {
00033
00034
00035 }
00036
00051 PSD::PSD(ParamStructure::PSD parameters, Grid &grid) : Matrix3D<double>(grid.L.size, grid.pc.size, grid.alpha.size)
00052 {
00053 this->Initialize(parameters, grid);
00054 }
00055
00062 PSD::PSD(ParamStructure::PSD parameters, Grid &grid, BoundaryCondition L_UpperBoundaryCondition) : Matrix3D<double>(grid.L.size, grid.pc.size, grid.alpha.size)
00063 {
00064
00065
00066 if (parameters.initial_PSD_Type == "IPSDT_STEADY_STATE_FROM_BOUNDARY"
00067 && (L_UpperBoundaryCondition.initialType != "BCT_FILE"
00068 && L_UpperBoundaryCondition.initialType != "BCT_FILE_GRID"
00069 && L_UpperBoundaryCondition.initialType != "BCT_FILE_CHANGES")
00070 ) throw error_msg("PSD_INITIALIZATION", "PSD initialization: top L boundary can not be used for steady state, cause was not loaded");
00071 this->Initialize(parameters, grid, L_UpperBoundaryCondition.xSlice(0));
00072 }
00073
00074
00083 void PSD::Initialize(ParamStructure::PSD parameters, Grid &grid, Matrix2D<double> L_UpperBoundaryCondition)
00084 {
00085
00086 PSD_parameters = parameters;
00087
00088
00089 this->LoadInitialValue(parameters, grid, L_UpperBoundaryCondition);
00090
00091
00092 output_without_grid_file = new ofstream((parameters.output_PSD_folderName + parameters.output_PSD_fileName4D).c_str());
00093 (*output_without_grid_file) << "VARIABLES = \"Phase_Space_Density\" " << endl;
00094 output_without_grid_file->setf(ios_base::scientific, ios_base::floatfield);
00095
00096 initialized = true;
00097 }
00098
00099
00112 void PSD::SourcesAndLosses(
00113 GridElement &L, GridElement &pc, GridElement &alpha,
00114 Matrix3D<double> &SL,
00115 double dt, double Lpp,
00116 double tau, double tauLpp)
00117 {
00118 int il, im, ia;
00119
00120
00121 for (il = 0; il < L.size; il++) {
00122 for (im = 0; im < pc.size; im++) {
00123 for (ia = 0; ia < alpha.size; ia++) {
00124
00125
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135 if (L[il][im][ia] >= Lpp) {
00136
00137 (*this)[il][im][ia] = (*this)[il][im][ia] - (*this)[il][im][ia] * (dt / tau);
00138 } else {
00139
00140 (*this)[il][im][ia] = (*this)[il][im][ia] - (*this)[il][im][ia] * (dt / tauLpp);
00141 }
00142 }
00143 }
00144 }
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157 }
00158
00159
00180 void PSD::DiffusionMixTermExplicit(double dt, double Lpp,
00181 DiffusionCoefficient &Dpca, DiffusionCoefficient &DpcaLpp,
00182 GridElement &L, GridElement &pc, GridElement &alpha, Matrix3D<double> Jacobian,
00183 Matrix2D<double> pc_lowerBoundaryCondition,
00184 Matrix2D<double> pc_upperBoundaryCondition,
00185 Matrix2D<double> alpha_lowerBoundaryCondition,
00186 Matrix2D<double> alpha_upperBoundaryCondition,
00187 string pc_lowerBoundaryCondition_calculationType,
00188 string pc_upperBoundaryCondition_calculationType,
00189 string alpha_lowerBoundaryCondition_calculationType,
00190 string alpha_upperBoundaryCondition_calculationType)
00191 {
00192
00193 int il, im, ia;
00194
00195
00196 Matrix3D<double> PSD_Der_a_L(L.size, pc.size, alpha.size);
00197 Matrix3D<double> PSD_Der_a_R(L.size, pc.size, alpha.size);
00198 Matrix3D<double> PSD_Der_m_L(L.size, pc.size, alpha.size);
00199 Matrix3D<double> PSD_Der_m_R(L.size, pc.size, alpha.size);
00200
00201
00202 double Dpca_, Dpca_a_L, Dpca_a_R, Dpca_m_L, Dpca_m_R;
00203 double J_, J_a_L, J_a_R, J_m_L, J_m_R;
00204 double Laml, Lmal, Lamr, Lmar;
00205 double PSD_m_r, PSD_m_l, PSD_a_r, PSD_a_l, Lam, Lma;
00206
00207
00208 if (PSD_parameters.approximationMethod == "AM_Split_LR") {
00209
00210
00211 for (il = 0; il < L.size; il++) {
00212 for (im = 0; im < pc.size; im++) {
00213 for (ia = 1; ia < alpha.size-1; ia++) {
00214 PSD_Der_a_L[il][im][ia] = VF::G(alpha[il][im][ia]) * Dpca[il][im][ia]*((*this)[il][im][ia] - (*this)[il][im][ia-1])/(alpha[il][im][ia] - alpha[il][im][ia-1]);
00215 PSD_Der_a_R[il][im][ia] = VF::G(alpha[il][im][ia]) * Dpca[il][im][ia]*((*this)[il][im][ia+1] - (*this)[il][im][ia])/(alpha[il][im][ia+1] - alpha[il][im][ia]);
00216 }
00217 }
00218 for (im = 1; im < pc.size-1; im++) {
00219 for (ia = 0; ia < alpha.size; ia++) {
00220 PSD_Der_m_L[il][im][ia] = pc[il][im][ia] * pc[il][im][ia] * Dpca[il][im][ia]*((*this)[il][im][ia] - (*this)[il][im-1][ia])/(pc[il][im][ia] - pc[il][im-1][ia]);
00221 PSD_Der_m_R[il][im][ia] = pc[il][im][ia] * pc[il][im][ia] * Dpca[il][im][ia]*((*this)[il][im+1][ia] - (*this)[il][im][ia])/(pc[il][im+1][ia] - pc[il][im][ia]);
00222 }
00223 }
00224 }
00225
00226
00227 for (il = 0; il < L.size; il++) {
00228 for (im = 1; im < pc.size-1; im++) {
00229 for (ia = 1; ia < alpha.size-1; ia++) {
00230
00231 Dpca_ = Dpca[il][im][ia];
00232 Dpca_a_L = Dpca[il][im][ia-1];
00233 Dpca_a_R = Dpca[il][im][ia+1];
00234 Dpca_m_L = Dpca[il][im-1][ia];
00235 Dpca_m_R = Dpca[il][im+1][ia];
00236 J_ = Jacobian[il][im][ia];
00237 J_a_L = Jacobian[il][im][ia-1];
00238 J_a_R = Jacobian[il][im][ia+1];
00239 J_m_L = Jacobian[il][im-1][ia];
00240 J_m_R = Jacobian[il][im+1][ia];
00241
00242
00243
00244 Lmar = 1.0/(pc[il][im][ia] * pc[il][im][ia]) * (PSD_Der_a_R[il][im][ia] - PSD_Der_a_R[il][im-1][ia])/(pc[il][im][ia] - pc[il][im-1][ia]);
00245 Lmal = 1.0/(pc[il][im][ia] * pc[il][im][ia]) * (PSD_Der_a_L[il][im+1][ia] - PSD_Der_a_L[il][im][ia])/(pc[il][im+1][ia] - pc[il][im][ia]);
00246
00247
00248 Lamr = 1.0/(VF::G(alpha[il][im][ia])) * (PSD_Der_m_R[il][im][ia] - PSD_Der_m_R[il][im][ia-1])/(alpha[il][im][ia] - alpha[il][im][ia-1]);
00249 Laml = 1.0/(VF::G(alpha[il][im][ia])) * (PSD_Der_m_L[il][im][ia+1] - PSD_Der_m_L[il][im][ia])/(alpha[il][im][ia+1] - alpha[il][im][ia]);
00250
00251 (*this)[il][im][ia] = (*this)[il][im][ia] + dt * 0.5 * (Lmar + Lamr + Lmal + Laml);
00252
00253
00254
00255
00256
00257
00258
00259
00260
00261
00262 }
00263 }
00264 }
00265 } else if (PSD_parameters.approximationMethod == "AM_Split_C") {
00266 for (il = 0; il < L.size; il++) {
00267 for (im = 1; im < pc.size-1; im++) {
00268 for (ia = 1; ia < alpha.size-1; ia++) {
00269
00270 Dpca_ = Dpca[il][im][ia];
00271 Dpca_a_L = Dpca[il][im][ia-1];
00272 Dpca_a_R = Dpca[il][im][ia+1];
00273 Dpca_m_L = Dpca[il][im-1][ia];
00274 Dpca_m_R = Dpca[il][im+1][ia];
00275 J_ = Jacobian[il][im][ia];
00276 J_a_L = Jacobian[il][im][ia-1];
00277 J_a_R = Jacobian[il][im][ia+1];
00278 J_m_L = Jacobian[il][im-1][ia];
00279 J_m_R = Jacobian[il][im+1][ia];
00280
00281
00282 PSD_m_r = Dpca[il][im+1][ia]*Jacobian[il][im+1][ia]*((*this)[il][im+1][ia+1] - (*this)[il][im+1][ia-1])/(alpha[il][im+1][ia+1] - alpha[il][im+1][ia-1]);
00283 PSD_m_l = Dpca[il][im-1][ia]*Jacobian[il][im-1][ia]*((*this)[il][im-1][ia+1] - (*this)[il][im-1][ia-1])/(alpha[il][im-1][ia+1] - alpha[il][im-1][ia-1]);
00284 Lma = 1.0/Jacobian[il][im][ia] * (PSD_m_r - PSD_m_l) / (pc[il][im+1][ia] - pc[il][im-1][ia]);
00285
00286
00287 PSD_a_r = Dpca[il][im][ia+1]*Jacobian[il][im][ia+1]*((*this)[il][im+1][ia+1] - (*this)[il][im-1][ia+1])/(pc[il][im+1][ia+1] - pc[il][im-1][ia+1]);
00288 PSD_a_l = Dpca[il][im][ia-1]*Jacobian[il][im][ia-1]*((*this)[il][im+1][ia-1] - (*this)[il][im-1][ia-1])/(pc[il][im+1][ia-1] - pc[il][im-1][ia-1]);
00289 Lam = 1.0/Jacobian[il][im][ia] * (PSD_a_r - PSD_a_l) / (alpha[il][im][ia+1] - alpha[il][im][ia-1]);
00290
00291
00292 (*this)[il][im][ia] = (*this)[il][im][ia] + dt * (Lma + Lam);
00293 }
00294 }
00295 }
00296
00297 } else {
00298 throw error_msg("APPROXIMATION", "Unknown approximation method for mixed terms");
00299 }
00300
00301 }
00302
00303
00323 void PSD::Diffusion_alpha(double dt, double Lpp,
00324 DiffusionCoefficient &Daa, DiffusionCoefficient &DaaLpp,
00325 GridElement &L, GridElement &pc, GridElement &alpha, Matrix3D<double> Jacobian,
00326 Matrix2D<double> alpha_lowerBoundaryCondition,
00327 Matrix2D<double> alpha_upperBoundaryCondition,
00328 string alpha_lowerBoundaryCondition_calculationType,
00329 string alpha_upperBoundaryCondition_calculationType)
00330 {
00331
00332
00333
00334 static CalculationMatrix
00335 matr_A(alpha.size, 1, 1, 1),
00336 matr_B(alpha.size, 1, 1, 0),
00337 matr_C(alpha.size, 1, 1, 0);
00338
00339
00340 static Matrix3D<double> Zero_Matrix(1, 1, alpha.size);
00341 static Matrix2D<double> Zero_Matrix_2d(1, 1);
00342 static Matrix1D<double> Zero_Matrix_1d(alpha.size);
00343 Zero_Matrix = 0;
00344 Zero_Matrix_2d = 0;
00345 Zero_Matrix_1d = 0;
00346
00347
00348
00349 static Matrix1D<double> slicePSD1D(alpha.size), RHS(alpha.size);
00350
00351
00352
00353 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);
00354
00355 static Matrix2D<double> alpha_lowerBoundaryCondition_1d(1, 1), alpha_upperBoundaryCondition_1d(1, 1);
00356
00357
00358 int il, im, ia, in;
00359 DiagMatrix::iterator it;
00360 for (il = 0; il < L.size; il++) {
00361 for (im = 0; im < pc.size; im++) {
00362
00363
00364 for (ia = 0; ia < alpha.size; ia++) {
00365 slicePSD1D[ia] = (*this)[il][im][ia];
00366
00367 Daa_1d[0][0][ia] = Daa[il][im][ia];
00368 DaaLpp_1d[0][0][ia] = DaaLpp[il][im][ia];
00369
00370 L_1d[0][0][ia] = L[il][im][ia];
00371 pc_1d[0][0][ia] = pc[il][im][ia];
00372 alpha_1d[0][0][ia] = alpha[il][im][ia];
00373
00374 Jacobian_1d[0][0][ia] = Jacobian[il][im][ia];
00375 }
00376
00377 alpha_upperBoundaryCondition_1d[0][0] = alpha_upperBoundaryCondition[il][im];
00378 alpha_lowerBoundaryCondition_1d[0][0] = alpha_lowerBoundaryCondition[il][im];
00379
00380
00381 MakeModelMatrix_3D(
00382 matr_A, matr_B, matr_C,
00383 L_1d, pc_1d, alpha_1d,
00384 1, 1, alpha.size,
00385 Zero_Matrix_2d, Zero_Matrix_2d,
00386 Zero_Matrix_2d, Zero_Matrix_2d,
00387 alpha_lowerBoundaryCondition_1d, alpha_upperBoundaryCondition_1d,
00388 "BCT_CONSTANT_VALUE", "BCT_CONSTANT_VALUE",
00389 "BCT_CONSTANT_VALUE", "BCT_CONSTANT_VALUE",
00390 alpha_lowerBoundaryCondition_calculationType, alpha_upperBoundaryCondition_calculationType,
00391 Zero_Matrix,
00392 Zero_Matrix, Zero_Matrix,
00393 Daa_1d, DaaLpp_1d,
00394 Zero_Matrix, Zero_Matrix,
00395 Jacobian_1d, dt, Lpp);
00396
00397
00398
00399
00400
00401
00402
00403
00404 for (in = 0; in < alpha.size; in++) {
00405
00406
00407 RHS[in] = matr_C[0][in];
00408 for (it = matr_B.begin(); it != matr_B.end(); it++) {
00409
00410 if (in + it->first >= 0 && in + it->first < matr_B.total_size) {
00411 RHS[in] += it->second[in] * slicePSD1D[in + it->first];
00412 }
00413 }
00414 }
00415
00416
00417
00418
00419
00420
00421
00422 tridag(
00423 &matr_A[-1][0],
00424 &matr_A[0][0],
00425 &matr_A[+1][0],
00426 &RHS[0],
00427 &slicePSD1D[0],
00428 alpha.size);
00429
00430
00431 for (ia = 0; ia < alpha.size; ia++) {
00432 (*this)[il][im][ia] = slicePSD1D[ia];
00433 }
00434 }
00435 }
00436 }
00437
00438
00458 void PSD::Diffusion_pc(double dt, double Lpp,
00459 DiffusionCoefficient &Dpcpc, DiffusionCoefficient &DpcpcLpp,
00460 GridElement &L, GridElement &pc, GridElement &alpha, Matrix3D<double> Jacobian,
00461 Matrix2D<double> pc_lowerBoundaryCondition,
00462 Matrix2D<double> pc_upperBoundaryCondition,
00463 string pc_lowerBoundaryCondition_calculationType,
00464 string pc_upperBoundaryCondition_calculationType)
00465 {
00466
00467
00468
00469 static CalculationMatrix
00470 matr_A(1, pc.size, 1, 1),
00471 matr_B(1, pc.size, 1, 0),
00472 matr_C(1, pc.size, 1, 0);
00473
00474
00475 static Matrix3D<double> Zero_Matrix(1, pc.size, 1);
00476 static Matrix2D<double> Zero_Matrix_2d(1, 1);
00477 static Matrix1D<double> Zero_Matrix_1d(alpha.size);
00478 Zero_Matrix = 0;
00479 Zero_Matrix_2d = 0;
00480 Zero_Matrix_1d = 0;
00481
00482
00483
00484 static Matrix1D<double> slicePSD1D(pc.size), RHS(pc.size);
00485
00486
00487
00488 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);
00489
00490 static Matrix2D<double> pc_lowerBoundaryCondition_1d(1, 1), pc_upperBoundaryCondition_1d(1, 1);
00491
00492
00493 int il, im, ia, in;
00494 DiagMatrix::iterator it;
00495 for (il = 0; il < L.size; il++) {
00496 for (ia = 0; ia < alpha.size; ia++) {
00497
00498
00499 for (im = 0; im < pc.size; im++) {
00500
00501 slicePSD1D[im] = (*this)[il][im][ia];
00502
00503 Dpcpc_1d[0][im][0] = Dpcpc[il][im][ia];
00504 DpcpcLpp_1d[0][im][0] = DpcpcLpp[il][im][ia];
00505
00506 L_1d[0][im][0] = L[il][im][ia];
00507 pc_1d[0][im][0] = pc[il][im][ia];
00508 alpha_1d[0][im][0] = alpha[il][im][ia];
00509
00510 Jacobian_1d[0][im][0] = Jacobian[il][im][ia];
00511 }
00512
00513 pc_upperBoundaryCondition_1d[0][0] = pc_upperBoundaryCondition[il][ia];
00514 pc_lowerBoundaryCondition_1d[0][0] = pc_lowerBoundaryCondition[il][ia];
00515
00516
00517 MakeModelMatrix_3D(
00518 matr_A, matr_B, matr_C,
00519 L_1d, pc_1d, alpha_1d,
00520 1, pc.size, 1,
00521 Zero_Matrix_2d, Zero_Matrix_2d,
00522 pc_lowerBoundaryCondition_1d, pc_upperBoundaryCondition_1d,
00523 Zero_Matrix_2d, Zero_Matrix_2d,
00524 "BCT_CONSTANT_VALUE", "BCT_CONSTANT_VALUE",
00525 pc_lowerBoundaryCondition_calculationType, pc_upperBoundaryCondition_calculationType,
00526 "BCT_CONSTANT_VALUE", "BCT_CONSTANT_VALUE",
00527 Zero_Matrix,
00528 Dpcpc_1d, DpcpcLpp_1d,
00529 Zero_Matrix, Zero_Matrix,
00530 Zero_Matrix, Zero_Matrix,
00531 Jacobian_1d, dt, Lpp);
00532
00533
00534
00535
00536
00537
00538
00539
00540 DiagMatrix::iterator it;
00541 for (in = 0; in < pc.size; in++) {
00542
00543
00544 RHS[in] = matr_C[0][in];
00545 for (it = matr_B.begin(); it != matr_B.end(); it++) {
00546
00547 if (in + it->first >= 0 && in + it->first < matr_B.total_size) {
00548 RHS[in] += it->second[in] * slicePSD1D[in + it->first];
00549 }
00550 }
00551 }
00552
00553
00554
00555
00556
00557
00558
00559 tridag(
00560 &matr_A[-1][0],
00561 &matr_A[0][0],
00562 &matr_A[+1][0],
00563 &RHS[0],
00564 &slicePSD1D[0],
00565 pc.size);
00566
00567
00568 for (im = 0; im < pc.size; im++) {
00569 (*this)[il][im][ia] = slicePSD1D[im];
00570 }
00571 }
00572 }
00573 }
00574
00575
00594 void PSD::Diffusion_L( double dt, double Lpp,
00595 DiffusionCoefficient &DLL,
00596 GridElement &L, GridElement &pc, GridElement &alpha, Matrix3D<double> Jacobian,
00597 Matrix2D<double> L_lowerBoundaryCondition,
00598 Matrix2D<double> L_upperBoundaryCondition,
00599 string L_lowerBoundaryCondition_calculationType,
00600 string L_upperBoundaryCondition_calculationType,
00601 double tau, double tauLpp)
00602 {
00603
00604
00605
00606 static CalculationMatrix
00607 matr_A(1, 1, L.size, 1),
00608 matr_B(1, 1, L.size, 0),
00609 matr_C(1, 1, L.size, 0);
00610
00611
00612 static Matrix3D<double> Zero_Matrix(L.size, 1, 1);
00613 static Matrix2D<double> Zero_Matrix_2d(1, 1);
00614 static Matrix1D<double> Zero_Matrix_1d(1);
00615 Zero_Matrix = 0;
00616 Zero_Matrix_2d = 0;
00617 Zero_Matrix_1d = 0;
00618
00619
00620
00621 static Matrix1D<double> slicePSD1D(L.size), RHS(L.size);
00622
00623
00624
00625 static Matrix3D<double> DLL_1d(L.size, 1, 1), L_1d(L.size, 1, 1), Jacobian_1d(L.size, 1, 1);
00626
00627 static Matrix2D<double> L_lowerBoundaryCondition_1d(1, 1), L_upperBoundaryCondition_1d(1, 1);
00628
00629
00630 int il, im, ia, in;
00631 DiagMatrix::iterator it;
00632 for (ia = 0; ia < alpha.size; ia++) {
00633 for (im = 0; im < pc.size; im++) {
00634
00635
00636 for (il = 0; il < L.size; il++) {
00637
00638 slicePSD1D[il] = (*this)[il][im][ia];
00639
00640 DLL_1d[il][0][0] = DLL[il][im][ia];
00641
00642 L_1d[il][0][0] = L[il][im][ia];
00643
00644
00645
00646 Jacobian_1d[il][0][0] = Jacobian[il][im][ia];
00647 }
00648
00649 L_upperBoundaryCondition_1d[0][0] = L_upperBoundaryCondition[im][ia];
00650 L_lowerBoundaryCondition_1d[0][0] = L_lowerBoundaryCondition[im][ia];
00651
00652
00653 MakeModelMatrix_3D(
00654 matr_A, matr_B, matr_C,
00655 L_1d, Zero_Matrix, Zero_Matrix,
00656 L.size, 1, 1,
00657 L_lowerBoundaryCondition_1d, L_upperBoundaryCondition_1d,
00658 Zero_Matrix_2d, Zero_Matrix_2d,
00659 Zero_Matrix_2d, Zero_Matrix_2d,
00660 L_lowerBoundaryCondition_calculationType, L_upperBoundaryCondition_calculationType,
00661 "BCT_CONSTANT_VALUE", "BCT_CONSTANT_VALUE",
00662 "BCT_CONSTANT_VALUE", "BCT_CONSTANT_VALUE",
00663 DLL_1d,
00664 Zero_Matrix, Zero_Matrix,
00665 Zero_Matrix, Zero_Matrix,
00666 Zero_Matrix, Zero_Matrix,
00667 Jacobian_1d, dt, Lpp,
00668 tau, tauLpp);
00669
00670
00671
00672
00673
00674
00675
00676
00677 DiagMatrix::iterator it;
00678 for (in = 0; in < L.size; in++) {
00679
00680
00681 RHS[in] = matr_C[0][in];
00682 for (it = matr_B.begin(); it != matr_B.end(); it++) {
00683
00684 if (in + it->first >= 0 && in + it->first < matr_B.total_size) {
00685 RHS[in] += it->second[in] * slicePSD1D[in + it->first];
00686 }
00687 }
00688 }
00689
00690
00691
00692
00693
00694
00695
00696 tridag(
00697 &matr_A[-1][0],
00698 &matr_A[0][0],
00699 &matr_A[+1][0],
00700 &RHS[0],
00701 &slicePSD1D[0],
00702 L.size);
00703
00704
00705 for (il = 0; il < L.size; il++) {
00706 (*this)[il][im][ia] = slicePSD1D[il];
00707 }
00708 }
00709 }
00710 }
00711
00712
00738 void PSD::Diffusion_pc_alpha(double dt, double Lpp,
00739 DiffusionCoefficient &Dpcpc, DiffusionCoefficient &DpcpcLpp,
00740 DiffusionCoefficient &Daa, DiffusionCoefficient &DaaLpp,
00741 DiffusionCoefficient &Dpca, DiffusionCoefficient &DpcaLpp,
00742 GridElement &L, GridElement &pc, GridElement &alpha, Matrix3D<double> Jacobian,
00743 Matrix2D<double> pc_lowerBoundaryCondition,
00744 Matrix2D<double> pc_upperBoundaryCondition,
00745 Matrix2D<double> alpha_lowerBoundaryCondition,
00746 Matrix2D<double> alpha_upperBoundaryCondition,
00747 string pc_lowerBoundaryCondition_calculationType,
00748 string pc_upperBoundaryCondition_calculationType,
00749 string alpha_lowerBoundaryCondition_calculationType,
00750 string alpha_upperBoundaryCondition_calculationType)
00751 {
00752
00753
00754
00755 static CalculationMatrix
00756 matr_A(alpha.size, pc.size, 1, 1),
00757 matr_B(alpha.size, pc.size, 1, 0),
00758 matr_C(alpha.size, pc.size, 1, 0);
00759
00760
00761 static Matrix3D<double> Zero_Matrix(1, pc.size, alpha.size);
00762 static Matrix2D<double> Zero_Matrix_2d(1, 1);
00763 static Matrix1D<double> Zero_Matrix_1d(1);
00764 Zero_Matrix = 0;
00765 Zero_Matrix_2d = 0;
00766 Zero_Matrix_1d = 0;
00767
00768
00769
00770 static Matrix1D<double> slicePSD1D(matr_A.total_size), RHS(matr_A.total_size);
00771
00772
00773 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);
00774 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);
00775
00776 static Matrix2D<double> pc_lowerBoundaryCondition_slice(1, alpha.size), pc_upperBoundaryCondition_slice(1, alpha.size);
00777 static Matrix2D<double> alpha_lowerBoundaryCondition_slice(1, pc.size), alpha_upperBoundaryCondition_slice(1, pc.size);
00778
00779
00780 int il, im, ia, in;
00781 DiagMatrix::iterator it;
00782 for (il = 0; il < L.size; il++) {
00783
00784 for (im = 0; im < pc.size; im++) {
00785 for (ia = 0; ia < alpha.size; ia++) {
00786
00787 in = matr_A.index1d(ia, im);
00788 slicePSD1D[in] = (*this)[il][im][ia];
00789
00790 Dpcpc_slice[0][im][ia] = Dpcpc[il][im][ia];
00791 DpcpcLpp_slice[0][im][ia] = DpcpcLpp[il][im][ia];
00792 Daa_slice[0][im][ia] = Daa[il][im][ia];
00793 DaaLpp_slice[0][im][ia] = DaaLpp[il][im][ia];
00794 Dpca_slice[0][im][ia] = Dpca[il][im][ia];
00795 DpcaLpp_slice[0][im][ia] = DpcaLpp[il][im][ia];
00796
00797 L_slice[0][im][ia] = L[il][im][ia];
00798 pc_slice[0][im][ia] = pc[il][im][ia];
00799 alpha_slice[0][im][ia] = alpha[il][im][ia];
00800
00801 Jacobian_slice[0][im][ia] = Jacobian[il][im][ia];
00802 }
00803 }
00804
00805 for (ia = 0; ia < alpha.size; ia++) {
00806 pc_upperBoundaryCondition_slice[0][ia] = pc_upperBoundaryCondition[il][ia];
00807 pc_lowerBoundaryCondition_slice[0][ia] = pc_lowerBoundaryCondition[il][ia];
00808 }
00809 for (im = 0; im < pc.size; im++) {
00810 alpha_upperBoundaryCondition_slice[0][im] = alpha_upperBoundaryCondition[il][im];
00811 alpha_lowerBoundaryCondition_slice[0][im] = alpha_lowerBoundaryCondition[il][im];
00812 }
00813
00814
00815 MakeModelMatrix_3D (
00816 matr_A, matr_B, matr_C,
00817 L_slice, pc_slice, alpha_slice,
00818 1, pc.size, alpha.size,
00819 Zero_Matrix_2d, Zero_Matrix_2d,
00820 pc_lowerBoundaryCondition_slice, pc_upperBoundaryCondition_slice,
00821 alpha_lowerBoundaryCondition_slice, alpha_upperBoundaryCondition_slice,
00822 "BCT_CONSTANT_VALUE", "BCT_CONSTANT_VALUE",
00823 pc_lowerBoundaryCondition_calculationType, pc_upperBoundaryCondition_calculationType,
00824 alpha_lowerBoundaryCondition_calculationType, alpha_upperBoundaryCondition_calculationType,
00825 Zero_Matrix,
00826 Dpcpc_slice, DpcpcLpp_slice,
00827 Daa_slice, DaaLpp_slice,
00828 Dpca_slice, DpcaLpp_slice,
00829 Jacobian_slice, dt, Lpp);
00830
00831
00832
00833
00834
00835
00836
00837
00838 DiagMatrix::iterator it;
00839 for (im = 0; im < pc.size; im++) {
00840 for (ia = 0; ia < alpha.size; ia++) {
00841 in = matr_B.index1d(ia, im);
00842
00843
00844 RHS[in] = matr_C[0][in];
00845 for (it = matr_B.begin(); it != matr_B.end(); it++) {
00846
00847 if (in + it->first >= 0 && in + it->first < matr_B.total_size) {
00848 RHS[in] += it->second[in] * slicePSD1D[in + it->first];
00849 }
00850 }
00851
00852 }
00853 }
00854
00855
00856 RHS.writeToFile("./Debug/pc_alpha_RHS.dat");
00857
00858
00859 int n = pc.size * alpha.size;
00860
00861 static Matrix2D<double> A;
00862 static Matrix1D<double> B;
00863 int modelMatrixLineNumber;
00864 int k, l, in;
00865 double sum, error, error_max;
00866 static bool recalculate_A = true;
00867 if (PSD_parameters.solutionMethod == "SM_Gauss") {
00868
00869 if (recalculate_A) {
00870 B = RHS;
00871 if (!A.initialized) {
00872 A.AllocateMemory(n, n);
00873 }
00874 A = 0;
00875
00876
00877
00878
00879 Output::echo(5, "Inverting A matrix (gauss)... ");
00880
00881 for (it = matr_A.begin(); it != matr_A.end(); it++) {
00882 for (modelMatrixLineNumber = 0; modelMatrixLineNumber < n; modelMatrixLineNumber++) {
00883
00884
00885 if (modelMatrixLineNumber + it->first >= 0 && modelMatrixLineNumber + it->first < n) {
00886
00887
00888
00889 A[modelMatrixLineNumber][modelMatrixLineNumber + it->first] = it->second[modelMatrixLineNumber];
00890 }
00891 }
00892 }
00893
00894
00895
00896
00897
00898 gauss_solve(&A[0][0], &B[0], n);
00899 recalculate_A = true;
00900 Output::echo(5, "done.\n");
00901 }
00902
00903 Output::echo(5, "Calculation of new PSD... ");
00904 slicePSD1D[n-1] = B[n-1] / A[n-1][n-1];
00905 for (k = n - 2; k >= 0; k--) {
00906 sum = 0;
00907 for (l = k + 1; l < n; l++)
00908 sum += A[k][l] * slicePSD1D[l];
00909 slicePSD1D[k] = B[k] - sum;
00910 }
00911
00912
00913
00914
00915 error_max = 0.;
00916 error = 0;
00917 for (in = 0; in < n; in++) {
00918 error = RHS[in];
00919 for (it = matr_A.begin(); it != matr_A.end(); it++) {
00920
00921 if (in + it->first >= 0 && in + it->first < n) {
00922 error -= it->second[in] * slicePSD1D[in + it->first];
00923 }
00924 }
00925 error_max = VF::max( error_max, error );
00926 }
00927 Output::echo(5, "done.\n");
00928 Output::echo(5, "Error_norm = %e \n", error_max);
00929
00930 } else if (PSD_parameters.solutionMethod == "SM_Relaxation") {
00931
00932
00933 Output::echo(5, "Inverting A matrix (relaxation)... ");
00934
00935 over_relaxation_diag(matr_A, RHS, slicePSD1D);
00936 Output::echo(5, "done.");
00937
00938 } else if (PSD_parameters.solutionMethod == "SM_Lapack") {
00939
00940
00941 Output::echo(5, "Inverting A matrix (Lapack)... ");
00942
00943 Lapack(matr_A, RHS, slicePSD1D);
00944 Output::echo(5, "done.");
00945
00946 } else if (PSD_parameters.solutionMethod == "SM_GMRES") {
00947
00948
00949 Output::echo(5, "Inverting A matrix (GMRES)... ");
00950
00951 gmres_wrapout(matr_A, RHS, slicePSD1D, PSD_parameters.SOL_maxiter, PSD_parameters.SOL_i_max, PSD_parameters.SOL_max_iter_err);
00952 Output::echo(5, "done.");
00953
00954 } else if (PSD_parameters.solutionMethod == "SM_Tridiag") {
00955
00956 throw error_msg("SOLUTIONERR", "Tridiagonal method can't be used");
00957 } else {
00958 throw error_msg("SOLUTIONERR", "Unknown solution method");
00959 }
00960
00961
00962 for (im = 0; im < pc.size; im++) {
00963 for (ia = 0; ia < alpha.size; ia++) {
00964 in = matr_A.index1d(ia, im);
00965 (*this)[il][im][ia] = slicePSD1D[in];
00966 }
00967 }
00968 }
00969
00970 }
00971
00972
00981 void checkInf(Matrix3D<double> arr, int il, int im, int ia, double maxNum = 1e99) {
00982
00983
00984 if (!(arr[il][im][ia] >= 0 || arr[il][im][ia] <= 0) || arr[il][im][ia] > maxNum || arr[il][im][ia] < -maxNum){
00985 throw error_msg("INTERPOLATION_ERROR", "Interpolation error at [%d][%d][%d]", il, im, ia);
00986 }
00987 }
00988
00997 void checkInf(Matrix1D<double> arr, int il, int im, int ia, double maxNum = 1e99) {
00998
00999 if (!(arr[im] >= 0 || arr[im] <= 0) || arr[im] > maxNum || arr[im] < -maxNum){
01000 throw error_msg("INTERPOLATION_ERROR", "Interpolation error at [%d][%d][%d]", il, im, ia);
01001 }
01002 }
01003
01015 void PSD::Interpolate(PSD &otherPSD, ParamStructure::Interpolation interpolationParamStructure, Grid &oldGrid, Grid &newGrid, Matrix2D<double> newGrid_pc_lowerBoundaryCondition, Matrix2D<double> newGrid_pc_upperBoundaryCondition) {
01016
01017
01018 int il, im, ia;
01019 double lowerBoundary1d=0, upperBoundary1d=0;
01020
01021
01022
01023 if (interpolationParamStructure.type == "IT_NONE") {
01024 return;
01025 } else if (oldGrid.type == newGrid.type) {
01026
01027
01028 for (il = 0; il < oldGrid.L.size; il++) {
01029 for (ia = 0; ia < oldGrid.alpha.size; ia++) {
01030 for (im = 0; im < oldGrid.pc.size; im++) {
01031
01032 (*this)[il][im][ia] = otherPSD[il][im][ia];
01033
01034 if ((*this)[il][im][ia] < VC::zero_f) {
01035 (*this)[il][im][ia] = VC::zero_f;
01036 }
01037 }
01038 }
01039 }
01040 return;
01041 }
01042
01043 Matrix1D<double> matrix_array_1d_old(oldGrid.pc.size);
01044 Matrix1D<double> matrix_array_1d_new(newGrid.pc.size);
01045 Matrix1D<double> grid_1d_old(oldGrid.pc.size);
01046 Matrix1D<double> grid_1d_new(newGrid.pc.size);
01047
01048 for (il = 0; il < oldGrid.L.size; il++) {
01049 for (ia = 0; ia < oldGrid.alpha.size; ia++) {
01050 for (im = 0; im < oldGrid.pc.size; im++) {
01051
01052 if (newGrid.L[il][im][ia] != oldGrid.L[il][im][ia] || newGrid.alpha[il][im][ia] != oldGrid.alpha[il][im][ia]) {
01053
01054 throw error_msg("INTERPOLATION_ERROR", "Interpolation error: grids must have the same L and Alpha");
01055 }
01056 #ifdef DEBUG_MODE
01057
01058 checkInf((*this), il, im, ia);
01059 #endif
01060
01061
01062 if (interpolationParamStructure.useLog == "Log_10") {
01063
01064 if (otherPSD[il][im][ia] >= VC::zero_f) {
01065 matrix_array_1d_old[im] = log10(otherPSD[il][im][ia]);
01066 } else {
01067
01068 matrix_array_1d_old[im] = log10(VC::zero_f);
01069 }
01070
01071
01072
01073 lowerBoundary1d = log10(VF::max(newGrid_pc_lowerBoundaryCondition[il][ia], VC::zero_f));
01074 upperBoundary1d = log10(VF::max(newGrid_pc_upperBoundaryCondition[il][ia], VC::zero_f));
01075 } else if (interpolationParamStructure.useLog == "Log_2") {
01076
01077 if (otherPSD[il][im][ia] >= VC::zero_f) {
01078 matrix_array_1d_old[im] = log(otherPSD[il][im][ia])/log(2.0);
01079 } else {
01080
01081 matrix_array_1d_old[im] = log(VC::zero_f)/log(2.0);
01082 }
01083
01084
01085
01086 lowerBoundary1d = log10(VF::max(newGrid_pc_lowerBoundaryCondition[il][ia], VC::zero_f));
01087 upperBoundary1d = log10(VF::max(newGrid_pc_upperBoundaryCondition[il][ia], VC::zero_f));
01088 } else if (interpolationParamStructure.useLog == "Log_E") {
01089
01090 if (otherPSD[il][im][ia] >= VC::zero_f) {
01091 matrix_array_1d_old[im] = log(otherPSD[il][im][ia]);
01092 } else {
01093
01094 matrix_array_1d_old[im] = log(VC::zero_f);
01095 }
01096
01097
01098
01099 lowerBoundary1d = log10(VF::max(newGrid_pc_lowerBoundaryCondition[il][ia], VC::zero_f));
01100 upperBoundary1d = log10(VF::max(newGrid_pc_upperBoundaryCondition[il][ia], VC::zero_f));
01101 } else if (interpolationParamStructure.useLog == "NoLog") {
01102 matrix_array_1d_old[im] = otherPSD[il][im][ia];
01103 lowerBoundary1d = VF::max(newGrid_pc_lowerBoundaryCondition[il][ia], VC::zero_f);
01104 upperBoundary1d = VF::max(newGrid_pc_upperBoundaryCondition[il][ia], VC::zero_f);
01105 } else {
01106 throw error_msg("INTERPOLATION_ERROR", "Usage of log(PSD) is not define properly: %s", interpolationParamStructure.useLog.c_str());
01107 break;
01108 }
01109
01110 #ifdef DEBUG_MODE
01111
01112 checkInf(matrix_array_1d_old, il, im, ia, 99);
01113 #endif
01114
01115 grid_1d_old[im] = oldGrid.pc[il][im][ia];
01116 grid_1d_new[im] = newGrid.pc[il][im][ia];
01117 }
01118
01119 if (interpolationParamStructure.type == "IT_LINEAR") {
01120
01121 matrix_array_1d_new.Interpolate(matrix_array_1d_old, grid_1d_old, grid_1d_new);
01122 } else if (interpolationParamStructure.type == "IT_POLYNOMIAL") {
01123
01124 matrix_array_1d_new.Polilinear(matrix_array_1d_old, grid_1d_old, grid_1d_new, lowerBoundary1d, upperBoundary1d);
01125 } else if (interpolationParamStructure.type == "IT_SPLINE") {
01126
01127 matrix_array_1d_new.Spline(matrix_array_1d_old, grid_1d_old, grid_1d_new, lowerBoundary1d, upperBoundary1d, interpolationParamStructure.linearSplineCoef, interpolationParamStructure.maxSecondDerivative);
01128 } else if (interpolationParamStructure.type == "IT_COPY") {
01129
01130 matrix_array_1d_new = matrix_array_1d_old;
01131 } else {
01132 throw error_msg("INTERPOLATION_TYPE_ERROR", "Unknown interpolation type");
01133 }
01134
01135 for (im = 0; im < newGrid.pc.size; im++) {
01136
01137 #ifdef DEBUG_MODE
01138 checkInf(matrix_array_1d_new, il, im, ia, 99);
01139 #endif
01140
01141
01142 if (interpolationParamStructure.useLog == "Log_10") (*this)[il][im][ia] = pow(10, matrix_array_1d_new[im]);
01143 else if (interpolationParamStructure.useLog == "Log_2") (*this)[il][im][ia] = pow(2, matrix_array_1d_new[im]);
01144 else if (interpolationParamStructure.useLog == "Log_E") (*this)[il][im][ia] = pow(VC::exp, matrix_array_1d_new[im]);
01145 else (*this)[il][im][ia] = matrix_array_1d_new[im];
01146
01147 #ifdef DEBUG_MODE
01148 checkInf((*this), il, im, ia);
01149 #endif
01150
01151 }
01152 }
01153 }
01154
01155
01156 }
01157
01165 void PSD::LoadInitialValue(ParamStructure::PSD parameters, Grid &grid, Matrix2D<double> L_UpperBoundaryCondition) {
01166
01167 int il, im, ia;
01168
01169
01170 if (parameters.initial_PSD_Type == "IPSDT_ORBIT_FLUX_2D") {
01171
01172 if (grid.L.size != 1) {
01173 throw error_msg("INITIAL_PSD_2D_METHOD_FOR_3D_GRID", "Initial PSD error: imposibble to use ORBIT_FLUX_2D if L.size > 1");
01174 }
01175 Load_initial_f_2d(grid.L, grid.pc, grid.alpha, parameters.initial_PSD_fileName.c_str());
01176 } else if (parameters.initial_PSD_Type == "IPSDT_STEADY_STATE") {
01177
01178 if (grid.L.size == 1) {
01179 throw error_msg("INITIAL_PSD_3D_METHOD_FOR_2D_GRID", "Initial PSD error: imposibble to use STEADY_STATE if L.size == 1");
01180 }
01181
01183 Load_initial_f(grid.L, grid.pc, grid.alpha, parameters.initial_PSD_tauSteadyState, parameters.initial_PSD_Kp0, parameters.initial_PSD_some_constant_value, parameters.initial_PSD_J_L7_function, parameters.initial_PSD_outer_psd, parameters.initial_PSD_inner_psd);
01184 } else if (parameters.initial_PSD_Type == "IPSDT_STEADY_STATE_FROM_BOUNDARY") {
01185
01186 if (grid.L.size == 1) {
01187 throw error_msg("INITIAL_PSD_3D_METHOD_FOR_2D_GRID", "Initial PSD error: imposibble to use STEADY_STATE if L.size == 1");
01188 }
01189 if (!L_UpperBoundaryCondition.initialized) {
01190 throw error_msg("INITIAL_PSD_ERROR", "Initial PSD error: upper boundary is not initialized, can not be used for steady state calculation.");
01191 }
01192
01194 Load_initial_f(grid.L, grid.pc, grid.alpha, parameters.initial_PSD_tauSteadyState, parameters.initial_PSD_Kp0, L_UpperBoundaryCondition, parameters.initial_PSD_some_constant_value, parameters.initial_PSD_outer_psd, parameters.initial_PSD_inner_psd);
01195 } else if (parameters.initial_PSD_Type == "IPSDT_INTERPOLATION") {
01196
01197
01198 throw error_msg("INITIAL_PSD_TYPE_ERROR", "Using interpolation initial PSD type without second PSD to interpolate from");
01199 } else if (parameters.initial_PSD_Type == "IPSDT_CONSTANT") {
01200
01201 for (il = 0; il < grid.L.size; il++)
01202 for (im = 0; im < grid.pc.size; im++)
01203 for (ia = 0; ia < grid.alpha.size; ia++)
01204 (*this)[il][im][ia] = parameters.initial_PSD_some_constant_value;
01205 } else if (parameters.initial_PSD_Type == "IPSDT_FILE") {
01206
01207 Load_initial_f_file(grid.L, grid.epc, grid.alpha, parameters.initial_PSD_fileName.c_str(), false);
01208 } else if (parameters.initial_PSD_Type == "IPSDT_FILE_GRID") {
01209
01210 Load_initial_f_file(grid.L, grid.epc, grid.alpha, parameters.initial_PSD_fileName.c_str(), true);
01211 } else {
01212 throw error_msg("INITIAL_PSD_TYPE_UNKNOWN", "Unknown initial PSD type");
01213 }
01214
01215
01216 if (parameters.initial_PSD_Type == "IPSDT_ORBIT_FLUX_2D" ||
01217 parameters.initial_PSD_Type == "IPSDT_STEADY_STATE" ||
01218 parameters.initial_PSD_Type == "IPSDT_STEADY_STATE_FROM_BOUNDARY") {
01219 for (il = 0; il < grid.L.size; il++)
01220 for (im = 0; im < grid.pc.size; im++)
01221 for (ia = 0; ia < grid.alpha.size; ia++)
01222 if (grid.alpha[il][im][ia] < VF::alc(grid.L[il][im][ia]) || grid.alpha[il][im][ia] > VC::pi - VF::alc(grid.L[il][im][ia]))
01223 (*this)[il][im][ia] = VC::zero_f;
01224 }
01225
01226
01227 }
01228
01229
01230
01231
01232
01238 void PSD::Output_without_grid(double time) {
01239 int il, im, ia;
01240 (*output_without_grid_file) << "ZONE T=\"" << internal << time << "\" I=" << (*this).size_z << ", J=" << (*this).size_y << ", K=" << (*this).size_x << endl;
01241 for (il = 0; il < (*this).size_x; il++) {
01242 for (im = 0; im < (*this).size_y; im++) {
01243 for (ia = 0; ia < (*this).size_z; ia++) {
01244 (*output_without_grid_file) << (*this)[il][im][ia] << endl;
01245 }
01246 }
01247 }
01248
01249 }
01250
01251
01260 void PSD::Load_initial_f_file(GridElement &L, GridElement &epc, GridElement &alpha, const char *filename, bool withGrid) {
01261 int il, im, ia;
01262 string inBuf;
01263 double Load_L=-1, Load_epc=-1, Load_alpha=-1;
01264 double err = 1.e-3;
01265
01266
01267 ifstream input_f(filename);
01268 if (input_f != NULL && !input_f.eof()) {
01269
01270 getline(input_f, inBuf);
01271 getline(input_f, inBuf);
01272 for (il = 0; il < L.size; il++) {
01273 for (im = 0; im < epc.size; im++) {
01274 for (ia = 0; ia < alpha.size; ia++) {
01275 if (withGrid) input_f >> Load_L >> Load_epc >> Load_alpha;
01276 input_f >> (*this)[il][im][ia];
01277 if (withGrid)
01278 if (fabs(Load_L - L[il][im][ia]) > err || fabs(Load_epc - epc[il][im][ia]) > err || fabs(Load_alpha - alpha[il][im][ia]) > err) {
01279 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[il][im][ia], epc[il][im][ia], alpha[il][im][ia]);
01280 }
01281 }
01282 }
01283 }
01284 } else {
01285 throw error_msg("INITIAL_PSD_ERROR", "Error reading PSD input file %s.\n", filename);
01286
01287
01288
01289 }
01290 input_f.close();
01291 }
01292
01307 void PSD::Load_initial_f(GridElement &L, GridElement &pc, GridElement &alpha, double tau, double Kp, double min_f, string J_L7_function, double fb_out, double fb_in) {
01308 int il, im, ia;
01309 double p1;
01310 Matrix1D<double> L_1d(L.size), fL(L.size);
01311 Matrix2D<double> fLm(L.size, pc.size);
01312
01313
01314 for (il = 0; il < L.size; il++) L_1d[il] = L[il][0][alpha.size-1];
01315 steady_state(fL, tau, Kp, L.size, L_1d, fb_out, fb_in);
01316 for (int i = 0; i < fL.size_x; i++) if (fL[i] < min_f) fL[i] = min_f;
01317
01318
01319
01320
01321
01322
01323
01324 #ifdef DEBUG_MODE
01325 fL.writeToFile("./Debug/fL.dat");
01326 #endif
01327
01328 for (im = 0; im < pc.size; im++) {
01329 for (il = 0; il < L.size; il++) {
01330 for (ia = 0; ia < alpha.size; ia++) {
01331
01332
01333
01334
01335
01336
01337
01338
01339 int ia_90=alpha.size-1;
01340 while (alpha[il][im][ia_90] > VC::pi/2 && ia_90 > 0) ia_90--;
01341 p1 = pc[il][im][ia_90] * sqrt( VF::B(7.)/VF::B(L[il][im][ia_90]));
01342 if (J_L7_function == "J_L7") {
01343 (*this)[il][im][ia] = sin(alpha[il][im][ia]) * fL[il] * VF::J_L7( VF::Kfunc(p1),0.2,2.e3,1,7 ) * pow(1.0/p1,2);
01344 } else if (J_L7_function == "J_L7_corrected") {
01345 (*this)[il][im][ia] = sin(alpha[il][im][ia]) * fL[il] * VF::J_L7_corrected( VF::Kfunc(p1) ) * pow(1.0/p1,2);
01346 } else {
01347 throw error_msg("INITIAL_PSD", "Unknown boundary function J_L7: %s. Should be 'J_L7' or 'J_L7_corrected'.", J_L7_function.c_str());
01348 }
01349
01350 if (alpha[il][im][ia] < VF::alc(L[il][im][ia]) || alpha[il][im][ia] > VC::pi - VF::alc(L[il][im][ia])) {
01351
01352
01353 (*this)[il][im][ia] = (*this)[il][im][ia] * pow(sinh(alpha[il][im][ia]/VF::alc(L[il][im][ia])),2)/pow(sinh(1.0),2);
01354 }
01355
01356
01357
01358 }
01359 }
01360 }
01361
01362 #ifdef DEBUG_MODE
01363 (*this).writeToFile("./Debug/f_0.plt", L, pc, alpha);
01364 #endif
01365
01366 }
01367
01368
01382 void PSD::Load_initial_f(GridElement &L, GridElement &pc, GridElement &alpha, double tau, double Kp, Matrix2D<double> L_UpperBoundaryCondition, double min_f, double fb_out, double fb_in) {
01383 int il, im, ia;
01384 double p1;
01385 Matrix1D<double> L_1d(L.size), fL(L.size);
01386 Matrix2D<double> fLm(L.size, pc.size);
01387
01388
01389
01390
01391 for (il = 0; il < L.size; il++) L_1d[il] = L[il][0][alpha.size-1];
01392 steady_state(fL, tau, Kp, L.size, L_1d, fb_out, fb_in);
01393 for (int i = 0; i < fL.size_x; i++) if (fL[i] < min_f) fL[i] = min_f;
01394
01395 #ifdef DEBUG_MODE
01396 fL.writeToFile("./Debug/fL.dat");
01397 #endif
01398
01399
01400 for (im = 0; im < pc.size; im++) {
01401 for (il = 0; il < L.size; il++) {
01402 for (ia = 0; ia < alpha.size; ia++) {
01403 p1 = sqrt( VF::B(7.)/VF::B(L[il][im][ia]));
01404 (*this)[il][im][ia] = sin(alpha[il][im][ia]) * fL[il] * L_UpperBoundaryCondition[im][ia] * pow(1.0/p1,2);
01405 }
01406 }
01407 }
01408
01409 #ifdef DEBUG_MODE
01410 (*this).writeToFile("./Debug/f_0.plt", L, pc, alpha);
01411 #endif
01412
01413 }
01414
01415
01424 void PSD::Load_initial_f_2d(GridElement &L, GridElement &pc, GridElement &alpha, const char *filename) {
01425 int i, j, k;
01426 int il, im, ia;
01427 int jmax, lmax;
01428
01429 lmax=8;
01430
01431 jmax=16;
01432
01433 Matrix2D<double> jperp(lmax+1, jmax+1);
01434 Matrix1D<double> energy(jmax+2), jcopy(jmax+2);
01435 Matrix1D<double> j_0m(pc.size), f_0m(pc.size);
01436 Matrix1D<double> dayno(lmax+1), eventtime(lmax+1), orbit(lmax+1);
01437
01438
01439 ifstream input_f(filename);
01440 if (input_f != NULL && !input_f.eof()) {
01441 for (k = 0; k <= lmax; k++) {
01442 input_f >> orbit[k] >> dayno[k] >> eventtime[k];
01443 for (j = 0; j <= jmax; j++) {
01444 input_f >> i >> energy[j] >> jperp[k][j];
01445 energy[j] = energy[j]/1000;
01446 }
01447 }
01448 } else {
01449 throw error_msg("INITIAL_PSD_ERROR", "Error reading initial conditions input file %s.\n", filename);
01450
01451
01452
01453 }
01454 input_f.close();
01455
01456
01457 for (i = 0; i <= jmax; i++) jcopy[i] = jperp[2][i];
01458
01459
01460 energy[jmax+1] = VF::Kfunc(pc.GridElement_parameters.max);
01461 jcopy[jmax+1] = VC::zero_f;
01462
01463 #ifdef DEBUG_MODE
01464 jcopy.writeToFile("./Debug/jcopy.plt", energy);
01465 #endif
01466
01467 if (VF::Kfunc(pc.GridElement_parameters.min) < energy[0]) printf("Warning: Emin less than minimum data in input file!\n");
01468
01469
01470 Matrix1D<double> epc1D(pc.size);
01471 for (i = 0; i < pc.size; i++)
01472 epc1D[i] = VF::Kfunc(pc[0][i][0]);
01473
01474
01475 for (i = 0; i < jmax+2; i++) {
01476 jcopy[i] = log10(jcopy[i]);
01477 }
01478 j_0m.Spline(jcopy, energy, epc1D, jcopy[0], log10(VC::zero_f), 0, 1);
01479 for (i = 0; i < pc.size; i++) {
01480 j_0m[i] = pow(10, j_0m[i]);
01481 }
01482
01483
01484 for (i = 0; i < pc.size; i++) {
01485 if (j_0m[i] < 0) j_0m[i] = VC::zero_f;
01486 f_0m[i] = j_0m[i]/pow(pc[0][i][0], 2);
01487 }
01488
01489 #ifdef DEBUG_MODE
01490 j_0m.writeToFile("./Debug/j_0m.plt", epc1D);
01491 f_0m.writeToFile("./Debug/f_0m.plt", epc1D);
01492 #endif
01493
01494 for (il = 0; il < L.size; il++) {
01495 for (im = 0; im < pc.size; im++) {
01496 for (ia = 0; ia < alpha.size; ia++) {
01497 (*this)[il][im][ia] = f_0m[im] * sin( alpha[il][im][ia] );
01498 if (alpha[il][im][ia] < VF::alc(L[il][im][ia]) || alpha[il][im][ia] > VC::pi - VF::alc(L[il][im][ia])) {
01499 (*this)[il][im][ia] = (*this)[il][im][ia] * pow(sinh(alpha[il][im][ia]/VF::alc(L[il][im][ia])),2)/pow(sinh(1.0),2);
01500 }
01501 }
01502 }
01503 }
01504 }
01505
01506
01507
01508
01509
01510
01511
01512
01513
01514
01515
01516
01517
01518
01519
01520
01521
01522
01534 void steady_state(Matrix1D<double> &f, double tau, double Kp, int nx, Matrix1D<double> &L, double f_bnd_out, double f_bnd_in) {
01535
01536
01537
01538 Matrix1D<double> C_1(nx);
01539 Matrix1D<double> C_2(nx);
01540 Matrix1D<double> C_3(nx);
01541 Matrix1D<double> alfa(nx);
01542 Matrix1D<double> A(nx);
01543 Matrix1D<double> B(nx);
01544 Matrix1D<double> C(nx);
01545 Matrix1D<double> D(nx);
01546 Matrix1D<double> Dxx(nx);
01547
01548 int i;
01549
01550 for (i = 0; i <= nx-1; i++) {
01551 Dxx[i] = VF::Df(L[i], Kp);
01552 }
01553
01554 for (i = 1; i < nx-1; i++) {
01555 alfa[i] = 1.0/((L[i+1] - L[i-1])/2)*pow(L[i],2);
01556
01557 C_1[i] = ((Dxx[i] + Dxx[i-1])/2)/(L[i] - L[i-1])/(pow(((L[i] + L[i-1])/2), 2));
01558 C_2[i] = ((Dxx[i+1] + Dxx[i])/2)/(L[i+1] - L[i])/(pow(((L[i+1] + L[i])/2), 2));
01559 C_3[i] = 1./tau;
01560
01561 }
01562
01563 for (i = 1; i < nx-1; i++) {
01564 A[i] = alfa[i] * C_1[i];
01565 B[i] = -alfa[i] * (C_1[i] + C_2[i]) - C_3[i];
01566 C[i] = alfa[i] * C_2[i];
01567 }
01568
01569 for (i = 0; i < nx; i++) D[i] = 0;;
01570
01571 i = 0;
01572 A[i] = 0.;
01573 B[i] = 1;
01574 C[i] = 0;
01575 D[i] = f_bnd_in;
01576
01577 i=nx-1;
01578 A[i] = 0;
01579 B[i] = 1;
01580 C[i] = 0;
01581 D[i] = f_bnd_out;
01582
01583 #ifdef DEBUG_MODE
01584 Dxx.writeToFile("./Debug/SS_DLL.plt");
01585 A.writeToFile("./Debug/SS_A.plt");
01586 B.writeToFile("./Debug/SS_B.plt");
01587 C.writeToFile("./Debug/SS_C.plt");
01588 D.writeToFile("./Debug/SS_D.plt");
01589 #endif
01590
01591 tridag(&A[0], &B[0], &C[0], &D[0], &f[0], nx);
01592
01593
01594
01595 }