VERB_code_2.3
PSD.cpp
Go to the documentation of this file.
1 
16 //#define DEBUG_MODE
17 
18 // Header file
19 #include "PSD.h"
20 
21 // Mathematical function
22 #include <math.h>
23 
24 // Various function requred for Radiation Belts simulation
25 #include "../VariousFunctions/variousFunctions.h"
26 
27 
28 
29 //
30 #include "../Exceptions/error.h"
31 
32 // General c++ headers
33 #include <iostream>
34 #include <string>
35 #include <ctime>
36 
37 #include <stdio.h>
38 #include <stdlib.h>
39 
40 using namespace std;
41 
42 
48 {
49 
50  if (!arr.initialized) {
51  arr.AllocateMemory(grid.L.size, grid.pc.size, grid.alpha.size);
52  }
53 
54  // Copying from parameters to the class variables
55  PSD_parameters = parameters;
56 
57  // Loading initial values, from time-dependent boundary conditions takes only first time step values
58  //LoadInitialValue(parameters, grid, L_UpperBoundaryCondition.xSlice(0));
59  // iterators
60  int il, im, ia;
61 
62  // loading (calculating) initial values
63  if (parameters.initial_PSD_Type == "IPSDT_ORBIT_FLUX_2D") {
64  // loading from orbit flux (for 2d)
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");
67  }
68  Load_initial_f_2d(grid.L, grid.pc, grid.alpha, parameters.initial_PSD_fileName.c_str());
70  }
71  else if (parameters.initial_PSD_Type == "IPSDT_FLUX_2D_MAXWELLIAN") {
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");
74  }
75  cout<<"kim................................"<<endl;
76  Load_initial_f_maxwell(grid.L, grid.pc, grid.alpha);
77 
79  }
80  else if (parameters.initial_PSD_Type == "IPSDT_STEADY_STATE") {
81  // calculating as a steady state (for 3d)
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");
84  }
85  // if tauSteadyState == 0 - using 6.0/Kp
87  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, false); // false indicating old method
89  }
90  // new version of steady state solution, that requires additional boundary conditions using L7 - ADrozdov
91  else if (parameters.initial_PSD_Type == "IPSDT_STEADY_STATE_L7") {
92  // calculating as a steady state (for 3d)
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");
95  }
96  // if tauSteadyState == 0 - using 6.0/Kp
98  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, true); // true indicating new method
100  }
101  else if (parameters.initial_PSD_Type == "IPSDT_STEADY_STATE_TWO_ZONE") {
102  // calculating as a steady state (for 3d)
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");
105  }
106  Load_initial_f_two_zone(grid.L, grid.pc, grid.alpha, parameters.initial_PSD_tauSteadyState, parameters.initial_PSD_tauLppSteadyState, 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);
108  }
109  else if (parameters.initial_PSD_Type == "IPSDT_STEADY_STATE_FROM_BOUNDARY") {
110  // calculating as a steady state (for 3d)
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");
113  }
114  if (!L_UpperBoundaryCondition.arr.initialized) {
115  throw error_msg("INITIAL_PSD_ERROR", "Initial PSD error: upper boundary is not initialized, can not be used for steady state calculation.");
116  }
117  // if tauSteadyState == 0 - using 6.0/Kp
119  Load_initial_f_from_outer_L(grid.L, grid.pc, grid.alpha,
120  parameters.initial_PSD_tauSteadyState, parameters.initial_PSD_Kp0,
121  L_UpperBoundaryCondition.arr, parameters.initial_PSD_some_constant_value,
122  parameters.initial_PSD_outer_psd, parameters.initial_PSD_inner_psd);
123  }
124  //not used
125  //else if (parameters.initial_PSD_Type == "IPSDT_INTERPOLATION") {
126  // Initial values for local grid as an interpolation from radial grid, or opposite
127  // throw error_msg("INITIAL_PSD_TYPE_ERROR", "Using interpolation initial PSD type without second PSD to interpolate from");
128  //}
129  else if (parameters.initial_PSD_Type == "IPSDT_CONSTANT") {
130  // just constant values
131  for (il = 0; il < grid.L.size; il++)
132  for (im = 0; im < grid.pc.size; im++)
133  for (ia = 0; ia < grid.alpha.size; ia++)
134  arr[il][im][ia] = parameters.initial_PSD_some_constant_value;
135  }
136  else if (parameters.initial_PSD_Type == "IPSDT_FILE") {
137  // loading from orbit flux (for 2d)
138  Load_initial_f_file(grid.L, grid.epc, grid.alpha, parameters.initial_PSD_fileName.c_str(), false);
139  }
140  else if (parameters.initial_PSD_Type == "IPSDT_FILE_GRID") {
141  // loading from orbit flux (for 2d)
142  Load_initial_f_file(grid.L, grid.epc, grid.alpha, parameters.initial_PSD_fileName.c_str(), true);
143  }
144  else {
145  throw error_msg("INITIAL_PSD_TYPE_UNKNOWN", "Unknown initial PSD type");
146  }
147 
148  // Emptying loss cone
149  if (parameters.initial_PSD_Type == "IPSDT_ORBIT_FLUX_2D" ||
150  parameters.initial_PSD_Type == "IPSDT_STEADY_STATE" ||
151  parameters.initial_PSD_Type == "IPSDT_STEADY_STATE_L7" ||
152  parameters.initial_PSD_Type == "IPSDT_STEADY_STATE_FROM_BOUNDARY") {
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++) {
156  // lc = loss cone calculation result in rad
157  double lc = VF::alc(grid.L.arr[il][im][ia]);
158  double a = grid.alpha.arr[il][im][ia];
159  if (a < lc || a > VC::pi - lc)
160  // arr[il][im][ia] = VC::zero_f;
161  // alpha dependence inside the loss cone
162  //if (grid.alpha.arr[il][im][ia] < lc || grid.alpha.arr[il][im][ia] > VC::pi - VF::alc(grid.L.arr[il][im][ia])) {
163  //fLmp[*,i,*]=fLmp[*,i,*]*sinh(i/8.)^2/sinh(1)^2
164  // ???????? arr[il][im][ia] = arr[il][im][ia] * pow(sinh(alpha[il][im][ia]/8.),2)/pow(sinh(1.0),2);
165  arr[il][im][ia] = arr[il][im][ia] * pow(sinh(a/lc),2)/pow(sinh(1.0),2);
166  //}
167  }
168 
169  }
170 
171 
172 
173  // Output - create stream
174  output_without_grid_file.open((parameters.output_PSD_folderName + parameters.output_PSD_fileName4D).c_str());
175  output_without_grid_file << "VARIABLES = \"Phase_Space_Density\" " << endl;
176  output_without_grid_file.setf(ios_base::scientific, ios_base::floatfield);
177 
178  // initialized = true;
179 }
180 
181 
185 /*
186 void PSD::ScaleBoundaryFlux(Parameters_structure::PSD parameters, GridElement &L, GridElement &pc, GridElement &alpha)
187 {
188 
189  Matrix1D<double> L_1d(61),fL(61),tau_0(61),Ke1(61);
190  Matrix2D<double> scale_fm(61,pc.size);
191  int il,im,ia;
192  double L_upper, scale_f, alpha1,Lpp,tau_c,tau_p;
193  double Kp = parameters.initial_PSD_Kp0;
194  double tau = parameters.initial_PSD_tauSteadyState;
195  double tauLpp = parameters.initial_PSD_tauLppSteadyState;
196  double fb_in,fb_out,Lmin,Lmax,Lsize,dL;
197 
198 
199  if (parameters.initial_PSD_Type == "IPSDT_STEADY_STATE") { // || parameters.initial_PSD_Type == "IPSDT_STEADY_STATE_TWO_ZONE") {
200  // solve 1d steady-state eq. from L=1 to L=7
201 
202  fb_in = 0;
203  fb_out = 1;
204  Lmin = 1.0;
205  Lmax = 7.0;
206  dL = 0.1;
207  Lsize = (Lmax-Lmin)/dL+1;
208 
209  for (il=0; il<Lsize; il++) L_1d[il] = 1+il*dL; // make 1d grid
210  steady_state(fL, tau, Kp, Lsize, L_1d, fb_out, fb_in); // solve 1d steady state eq. for L=1 to 7
211 
212  L_upper = L[L.size-1][0][0]; // L-upper boundary we want to change, which is given in parameters.ini file
213 
214  if(L_upper != Lmax) { // if L is not equal to 7
215  unsigned i=0; // interpolate PSD using Lagragian
216  while (L_1d[++i] <= L_upper);
217  scale_f = fL[i-1] * ( (L_upper-L_1d[i])/(L_1d[i-1]-L_1d[i]) )
218  + fL[i] * ( (L_upper-L_1d[i-1])/(L_1d[i]-L_1d[i-1]) );
219  } else { // if L is equal to 7
220  scale_f = fb_out; // scale_f = 1;
221  }
222 
223 
224  for (il = 0; il < L.size; il++) {
225  for (im = 0; im < pc.size; im++) {
226  for (ia = 0; ia < alpha.size; ia++) {
227  arr[il][im][ia] = arr[il][im][ia] * scale_f;
228 
229  }
230  }
231  }
232 
233 
234 
235  }
236 
237 
239 
240  if (parameters.initial_PSD_Type == "IPSDT_STEADY_STATE_TWO_ZONE") {
241 
242  fb_in = 0;
243  fb_out = 1;
244  Lmin = 1.0;
245  Lmax = 7.0;
246  dL = 0.1;
247  Lsize = (Lmax-Lmin)/dL+1;
248  //Lpp = (5.6 - 0.46*Kp);
249 
250  for (il=0; il<Lsize; il++) L_1d[il] = 1+il*dL; // make 1d grid
251 
252  for (im = 0; im < pc.size; im++){
253  for (il = 0; il < L.size; il++) {
254 
255  Ke1[il] = VF::Kfunc(pc[il][im][alpha.size-1]); // Energy in MeV
256  alpha1 = alpha[il][im][alpha.size-1]; // ~90 deg alpha
257  // cout<<Ke1[il]<<" "<<il<<Lsize<<endl;
258  // make a array of tau
259  //if(L_1d[il] >= Lpp) { // calculate tau
260  // tau_0[il] = tau;
261  //} else { // calculate tauLpp
262  if(tauLpp == 1e99) {
263  tau_0[il] = VF::Coulomb(L_1d[il],Ke1[il]);
264  } else if (tauLpp == 2e99) {
265  tau_0[il] = VF::Precipitation(Kp,L_1d[il],Ke1[il]);
266  } else if (tauLpp == 3e99) {
267  tau_c = VF::Coulomb(L_1d[il],Ke1[il]);
268  tau_p = VF::Precipitation(Kp,L_1d[il],Ke1[il]);
269  tau_0[il] = 1.0 / (1.0/tau_p+1.0/tau_c);
270  } else {
271  tau_0[il] = tauLpp;
272  }
273  //cout<<tau_0[il]<<" "<<il<<endl;
274  //}
275  }
276  steady_state_two_zone(fL, tau_0, Kp, alpha1, Ke1, L.size, L_1d, fb_out, fb_in);
277 
278  L_upper = L[L.size-1][0][0]; // L-upper boundary we want to change, which is given in parameters.ini file
279 
280  if(L_upper != Lmax) { // if L is not equal to 7
281  unsigned i=0; // interpolate PSD using Lagragian
282  while (L_1d[++i] <= L_upper);
283  scale_fm[il][im] = fL[i-1] * ( (L_upper-L_1d[i])/(L_1d[i-1]-L_1d[i]) )
284  + fL[i] * ( (L_upper-L_1d[i-1])/(L_1d[i]-L_1d[i-1]) );
285  } else { // if L is equal to 7
286  scale_fm[il][im] = fb_out; // scale_f = 1;
287  }
288  }
289  // scaling PSD ///////////////////////
290  for (il = 0; il < L.size; il++) {
291  for (im = 0; im < pc.size; im++) {
292  for (ia = 0; ia < alpha.size; ia++) {
293  arr[il][im][ia] = arr[il][im][ia] * scale_fm[il][im];
294  }
295  }
296  }
298  }
299 
300 
301 }
302 */
303 
304 
321  GridElement &L, GridElement &pc, GridElement &alpha,
323  double dt, double Lpp,
324  double tau, double tauLpp, double Kp)
325 {
326  // iterators
327  int il, im, ia;
328  // lifetime in loss cone
329  double taulc;
331  double tau_p,tau_c,tauLpp0;
333 
334  for (il = 0; il < L.size; il++) {
335  for (im = 0; im < pc.size; im++) {
336  for (ia = 0; ia < alpha.size; ia++) {
337 
338  // Loss cone. - does not work here, needs too small time step to be split from pitch-angle diffusion
339  // df/dt = -f/taulc
340 
341 
342 /*
343  if (parameters.useLossCone && ((alpha.arr[il][im][ia] <= VF::alc(L.arr[il][im][0])) || (alpha.arr[il][im][ia] >= (VC::pi - VF::alc(L.arr[il][im][0])))) ) {
344  // taulc = VF::bounce_time_res(VF::Kfunc(pc.arr[il][im][ia]), L.arr[il][im][ia])/60./60./24.;
345  // Quarter bounce time
346  taulc = 0.25 * VF::bounce_time_new(L.arr[il][im][ia], pc.arr[il][im][ia], alpha.arr[il][im][ia]);
347  arr[il][im][ia] = arr[il][im][ia] / dt / (1/dt + 1/taulc); // implicit
348  // arr[il][im][ia] = arr[il][im][ia] * exp(-dt/taulc); // some other
349  // arr[il][im][ia] = exp(log(arr[il][im][ia]) - dt/taulc); // to not have PSD=0
350  }
351 */
352  // Additional losses, if any
353  // df/dt = -f/tau
354  if (L.arr[il][im][ia] >= Lpp) {
355  // arr[il][im][ia] = arr[il][im][ia] / dt / (1/dt + 1/tau); // implicit
356  arr[il][im][ia] = arr[il][im][ia] - arr[il][im][ia] * (dt / tau); // explicit
357 
358  } else {
359  // (*this)[il][im][ia] = (*this)[il][im][ia] / dt / (1/dt + 1/tauLpp); // implicit
360 
361  // kinetic energy computed
362  double k = VF::Kfunc(pc.arr[il][im][ia]);
363  if(parameters.usetauLpp == "coulomb") {
364  // tau in days - initial
365  tauLpp0 = VF::Coulomb(L.arr[il][im][ia], k);
366  } else if (parameters.usetauLpp == "precipitaion") {
367  tauLpp0 = VF::Precipitation(Kp, L.arr[il][im][ia], k);
368  } else if (parameters.usetauLpp == "combined") {
369  tau_p = VF::Precipitation(Kp, L.arr[il][im][ia], k);
370  tau_c = VF::Coulomb(L.arr[il][im][ia], k);
371  tauLpp0 = 1.0/(1.0/tau_p+1.0/tau_c);
372  } else {
373  tauLpp0 = tauLpp;
374  }
375  arr[il][im][ia] = arr[il][im][ia] - arr[il][im][ia] * (dt / tauLpp0); // explicit, probably should be implicit! I.e.:
376  //arr[il][im][ia] = arr[il][im][ia] / dt / (1/dt + 1/tauLpp0); // <-- implicit
377  }
378 
379 
380 
381  // Add losses according to Magnetopause shadowing
382 
383  if (parameters.SL.useLmp == true)
384  arr[il][im][ia] = arr[il][im][ia] * SL.arr[il][im][ia];
385 
386 
387  }
388  }
389  } // end of loop over all grids
390 
391 }
392 
393 
414 void PSD::DiffusionMixTermExplicit(double dt, double Lpp,
416  GridElement &L, GridElement &pc, GridElement &alpha, Matrix3D<double> Jacobian,
417  Matrix2D<double> pc_lowerBoundaryCondition,
418  Matrix2D<double> pc_upperBoundaryCondition,
419  Matrix2D<double> alpha_lowerBoundaryCondition,
420  Matrix2D<double> alpha_upperBoundaryCondition,
421  string pc_lowerBoundaryCondition_calculationType,
422  string pc_upperBoundaryCondition_calculationType,
423  string alpha_lowerBoundaryCondition_calculationType,
424  string alpha_upperBoundaryCondition_calculationType)
425 {
426 
427  int il, im, ia;
428 
429  // Declaration of arrays for derivatives
430  static Matrix3D<double> PSD_Der_a_L(L.size, pc.size, alpha.size);
431  static Matrix3D<double> PSD_Der_a_R(L.size, pc.size, alpha.size);
432  static Matrix3D<double> PSD_Der_m_L(L.size, pc.size, alpha.size);
433  static Matrix3D<double> PSD_Der_m_R(L.size, pc.size, alpha.size);
434 
435  // Some variables
436  // center, alpha left, alpha right, pc left, pc right
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;
439  // alpha pc left, pc alpha left, alpha pc right, pc alpha right
440  double Laml, Lmal, Lamr, Lmar;
441  double PSD_m_r, PSD_m_l, PSD_a_r, PSD_a_l, Lam, Lma;
442 
443  // Different way of derivatives approximation
452  if (PSD_parameters.approximationMethod == "AM_Split_LR") {
453 
454  // Derivatives calculation
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++) {
458  // LHS derivative for alpha
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]);
460  // RHS derivative for alpha
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]);
462  }
463  }
464  for (im = 1; im < pc.size-1; im++) {
465  for (ia = 0; ia < alpha.size; ia++) {
466  // LHS derivative for pc
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]);
468  // RHS derivative for pc
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]);
470  }
471  }
472  }
473 
474  // Approximation and calculation
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++) {
478  // for debugging...
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];
489 
490 
491  // d/dp D d/da f
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]);
498 
499  // d/da D d/dp f
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]);
506 
507  // approximation calculation using L's
509  arr[il][im][ia] = arr[il][im][ia] + dt * 0.5 * (Lmar + Lamr + Lmal + Laml);
510 
511  /* L12p = 1/J*(Jacobian[il][im+1][ia]*Dpca[il][im+1][ia]*PSD_Der_a_L[il][im+1][ia] -
512  Jacobian[il][im ][ia]*Dpca[il][im ][ia]*PSD_Der_a_L[il][im ][ia])/(pc[il][im+1][ia] - pc[il][im ][ia]);
513  L12m = 1/J*(Jacobian[il][im ][ia]*Dpca[il][im ][ia]*PSD_Der_a_R[il][im ][ia] -
514  Jacobian[il][im-1][ia]*Dpca[il][im-1][ia]*PSD_Der_a_R[il][im-1][ia])/(pc[il][im ][ia] - pc[il][im-1][ia]);
515 
516  L21p = 1/J*(Jacobian[il][im][ia+1]*Dpca[il][im][ia+1]*PSD_Der_m_L[il][im][ia+1] -
517  Jacobian[il][im][ia ]*Dpca[il][im][ia ]*PSD_Der_m_L[il][im][ia ])/(alpha[il][im][ia+1] - alpha[il][im][ia]);
518  L21m = 1/J*(Jacobian[il][im][ia ]*Dpca[il][im][ia ]*PSD_Der_m_R[il][im][ia ] -
519  Jacobian[il][im][ia-1]*Dpca[il][im][ia-1]*PSD_Der_m_R[il][im][ia-1])/(alpha[il][im][ia] - alpha[il][im][ia-1]);*/
520  }
521  }
522  }
523  }
537  else if (PSD_parameters.approximationMethod == "AM_Split_C") { // was not fully checked, may not work
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++) {
541  // for debugging...
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];
552 
553  // d/dp D d/da f
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]);
557 
558  // d/da D d/dm f
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]);
562 
563  // cout << ((Lma + Lam) - (Lmar + Lamr + Lmal + Laml))/((Lma + Lam) + (Lmar + Lamr + Lmal + Laml)) << "%" << endl;
564  // approximation calculation using L's
565  arr[il][im][ia] = arr[il][im][ia] + dt * (Lma + Lam);
566  }
567  }
568  }
569 
570  } else {
571  throw error_msg("APPROXIMATION", "Unknown approximation method for mixed terms");
572  }
573 
574 }
575 
576 
578 
579 
599 void PSD::Diffusion_alpha(double dt, double Lpp,
601  GridElement &L, GridElement &pc, GridElement &alpha, Matrix3D<double> Jacobian,
602  Matrix2D<double> alpha_lowerBoundaryCondition,
603  Matrix2D<double> alpha_upperBoundaryCondition,
604  string alpha_lowerBoundaryCondition_calculationType,
605  string alpha_upperBoundaryCondition_calculationType)
606 {
607 
608  // seems to be no longer used
609  //char out_name_matr_a[256], out_name_matr_b[256], out_name_matr_c[256];
610 
611  // Calculation matrices A, B, and C for 1d diffusion problem
612  // A * PSD_1D(t+1) = B * PSD_1D(t) + C, for each L and Energy
613 
614  // Creation is moved to PSD.h, so it can be accessed from Main.cpp
615  //static CalculationMatrix
616  // matr_A(alpha.size, 1, 1, 1),
617  // matr_B(alpha.size, 1, 1, 0),
618  // matr_C(alpha.size, 1, 1, 0);
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);
622 
623 
624  // Create zero matrices, temporary matrices to be used later
625  static Matrix3D<double> Zero_Matrix(1, 1, alpha.size);
626  static Matrix2D<double> Zero_Matrix_2d(1, 1);
627  static Matrix1D<double> Zero_Matrix_1d(alpha.size);
628  Zero_Matrix = 0;
629  Zero_Matrix_2d = 0;
630  Zero_Matrix_1d = 0;
631 
632  // Create 1d arrays
633  // Creating 1d arrays for PSD and right hand side (RHS = B * PSD + C)
634  static Matrix1D<double> slicePSD1D(alpha.size), RHS(alpha.size);
635 
636  // Creating 1d arrays which will be sent to the MakeMatrix function
637  // Creating 1d matrixes, but with type of Matrix3D, cause we need that type for MakeMatrix function
638  // Create 3d and 2d matrices instead of 1d and single values so that it can meet parameters in MakeModelMatrix_3D()
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);
640  // Creating temporary arrays size 1x1 for boundary conditions, which are just values
641  static Matrix2D<double> alpha_lowerBoundaryCondition_1d(1, 1), alpha_upperBoundaryCondition_1d(1, 1);
642 
643  // iterators
644  int il, im, ia, in; // matrix element numbers on 3D grid
645  DiagMatrix::iterator it;
646  for (il = 0; il < L.size; il++) {
647  for (im = 0; im < pc.size; im++) {
648 
649  // Copying 3D arrays into 1d arrays
650  for (ia = 0; ia < alpha.size; ia++) {
651  slicePSD1D[ia] = arr[il][im][ia];
652 
653  Daa_1d[0][0][ia] = Daa.arr[il][im][ia];
654  DaaLpp_1d[0][0][ia] = DaaLpp.arr[il][im][ia];
655 
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];
659 
660  Jacobian_1d[0][0][ia] = Jacobian[il][im][ia];
661  }
662  // boundary conditions, these are arrays size 1x1, i.e. values
663  alpha_upperBoundaryCondition_1d[0][0] = alpha_upperBoundaryCondition[il][im];
664  alpha_lowerBoundaryCondition_1d[0][0] = alpha_lowerBoundaryCondition[il][im];
665 
666  // Make A, B, C matrices
668  matr_A, matr_B, matr_C, // matrices A, B, C
669  L_1d, pc_1d, alpha_1d, // L, pc, alpha grids. We use L crossection here
670  1, 1, alpha.size, // L-size, pc-size, alpha-size
671  Zero_Matrix_2d, Zero_Matrix_2d, // L-boundaries
672  Zero_Matrix_2d, Zero_Matrix_2d, // pc-boundaries
673  alpha_lowerBoundaryCondition_1d, alpha_upperBoundaryCondition_1d, // alpha-boundaries
674  "BCT_CONSTANT_VALUE", "BCT_CONSTANT_VALUE", // Boundary conditions at L-boundaries
675  "BCT_CONSTANT_VALUE", "BCT_CONSTANT_VALUE", // Boundary conditions at pc-boundaries
676  alpha_lowerBoundaryCondition_calculationType, alpha_upperBoundaryCondition_calculationType, // Boundary conditions at alpha-boundaries
677  Zero_Matrix, // DLL - use 0 matrix
678  Zero_Matrix, Zero_Matrix, // Dpp, DppLpp - use 0 matrix
679  Daa_1d, DaaLpp_1d, // Daa, DaaLpp
680  Zero_Matrix, Zero_Matrix, // Dpca, Dpca_Lpp - use 0 matrix
681  Jacobian_1d, dt, Lpp);
682 
683  /*Output::echo(0, "Writing calculation matrices. ");
684  sprintf(out_name_matr_a, "./Debug/alpha_matr_A_%d_%d_0.dat", il, im);
685  sprintf(out_name_matr_b, "./Debug/alpha_matr_B_%d_%d_0.dat", il, im);
686  sprintf(out_name_matr_c, "./Debug/alpha_matr_C_%d_%d_0.dat", il, im);
687  matr_A.writeToFile(out_name_matr_a);
688  matr_B.writeToFile(out_name_matr_b);
689  matr_C.writeToFile(out_name_matr_c);
690  Output::echo(0, "done.\n");*/
691 
692  // make RHS = B*f + C
693  for (in = 0; in < alpha.size; in++) {
694  // If assume matr_B is 1d and has only +1, 0, -1 diagonals, can use this:
695  // RHS[in] = matr_B[-1][in] * slicePSD1D[in-1] + matr_B[0][in] * slicePSD1D[in] + matr_B[+1][in] * slicePSD1D[in+1] + matr_C[0][in];
696  RHS[in] = matr_C[0][in];
697  for (it = matr_B.begin(); it != matr_B.end(); it++) {
698  // multiplication B * f
699  if (in + it->first >= 0 && in + it->first < matr_B.total_size) {
700  RHS[in] += it->second[in] * slicePSD1D[in + it->first];
701  }
702  }
703  }
704 
705  // right hand side to files
706  // RHS.writeToFile("alpha_RHS.dat");
707 
708  // Solving model matrix, obtaining PSD
709  // Assume matr_A has only +1, 0, -1 diagonals
710  // TODO: Check if matr_A has only +1, 0, -1 diagonals, i.e. is tridiagonal
711  tridag(
712  &matr_A[-1][0], // just a way to get a pointer to first element of an array
713  &matr_A[0][0],
714  &matr_A[+1][0],
715  &RHS[0],
716  &slicePSD1D[0],
717  alpha.size);
718 
719  // copy PSD back, from 1d array to 3d array
720  for (ia = 0; ia < alpha.size; ia++) {
721  arr[il][im][ia] = slicePSD1D[ia];
722  }
723  }
724  }
725 }
726 
727 
747 void PSD::Diffusion_pc(double dt, double Lpp,
748  DiffusionCoefficient &Dpcpc, DiffusionCoefficient &DpcpcLpp,
749  GridElement &L, GridElement &pc, GridElement &alpha, Matrix3D<double> Jacobian,
750  Matrix2D<double> pc_lowerBoundaryCondition,
751  Matrix2D<double> pc_upperBoundaryCondition,
752  string pc_lowerBoundaryCondition_calculationType,
753  string pc_upperBoundaryCondition_calculationType)
754 {
755  // no longer used
756  // char out_name_matr_a[256], out_name_matr_b[256], out_name_matr_c[256];
757 
758 
759  // Calculation matrices A, B, and C for 1d diffusion problem
760  // A * PSD_1D(t+1) = B * PSD_1D(t) + C, for each L and Energy
761 
762  // Creation is moved to PSD.h, so it can be accessed from Main.cpp
763  //static CalculationMatrix
764  // matr_A(1, pc.size, 1, 1),
765  // matr_B(1, pc.size, 1, 0),
766  // matr_C(1, pc.size, 1, 0);
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);
770 
771  // Create zero matrices, temporary matrices to be used later
772  static Matrix3D<double> Zero_Matrix(1, pc.size, 1);
773  static Matrix2D<double> Zero_Matrix_2d(1, 1);
774  static Matrix1D<double> Zero_Matrix_1d(alpha.size);
775  Zero_Matrix = 0;
776  Zero_Matrix_2d = 0;
777  Zero_Matrix_1d = 0;
778 
779  // Create 1d arrays
780  // Creating 1d arrays for PSD and right hand side (RHS = B * PSD + C)
781  static Matrix1D<double> slicePSD1D(pc.size), RHS(pc.size);
782 
783  // Creating 1d arrays which will be sent to the MakeMatrix function
784  // Creating 1d matrixes, but with type of Matrix3D, cause we need that type for MakeMatrix function
785  // Create 3d and 2d matrices instead of 1d and single values so that it can meet parameters in MakeModelMatrix_3D()
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);
787  // Creating temporary arrays size 1x1 for boundary conditions, which are just values
788  static Matrix2D<double> pc_lowerBoundaryCondition_1d(1, 1), pc_upperBoundaryCondition_1d(1, 1);
789 
790  // iterators
791  int il, im, ia, in; // matrix element numbers on 3D grid
792  DiagMatrix::iterator it;
793  for (il = 0; il < L.size; il++) {
794  for (ia = 0; ia < alpha.size; ia++) {
795 
796  // Copying 3D arrays into 1d arrays
797  for (im = 0; im < pc.size; im++) {
798  // Copying 3D PSD to 1d PSD
799  slicePSD1D[im] = arr[il][im][ia];
800 
801  Dpcpc_1d[0][im][0] = Dpcpc.arr[il][im][ia];
802  DpcpcLpp_1d[0][im][0] = DpcpcLpp.arr[il][im][ia];
803 
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];
807 
808  Jacobian_1d[0][im][0] = Jacobian[il][im][ia];
809  }
810  // boundary conditions, these are arrays size 1x1, i.e. values
811  pc_upperBoundaryCondition_1d[0][0] = pc_upperBoundaryCondition[il][ia];
812  pc_lowerBoundaryCondition_1d[0][0] = pc_lowerBoundaryCondition[il][ia];
813 
814  // Make A, B, C matrices
816  matr_A, matr_B, matr_C, // matrices A, B, C
817  L_1d, pc_1d, alpha_1d, // L, pc, alpha grids.
818  1, pc.size, 1, // L-size, pc-size, alpha-size
819  Zero_Matrix_2d, Zero_Matrix_2d, // L-boundaries
820  pc_lowerBoundaryCondition_1d, pc_upperBoundaryCondition_1d, // pc-boundaries
821  Zero_Matrix_2d, Zero_Matrix_2d, // alpha-boundaries
822  "BCT_CONSTANT_VALUE", "BCT_CONSTANT_VALUE", // Boundary conditions at L-boundaries
823  pc_lowerBoundaryCondition_calculationType, pc_upperBoundaryCondition_calculationType, // Boundary conditions at pc-boundaries
824  "BCT_CONSTANT_VALUE", "BCT_CONSTANT_VALUE", // Boundary conditions at alpha-boundaries
825  Zero_Matrix, // DLL - use 0 matrix
826  Dpcpc_1d, DpcpcLpp_1d, // Dpp, DppLpp
827  Zero_Matrix, Zero_Matrix, // Daa, DaaLpp - use 0 matrix
828  Zero_Matrix, Zero_Matrix, // Dpca, Dpca_Lpp - use 0 matrix
829  Jacobian_1d, dt, Lpp);
830 
831  /*Output::echo(0, "Writing calculation matrices. ");
832  sprintf(out_name_matr_a, "./Debug/pc_matr_A_%d_0_%d.dat", il, ia);
833  sprintf(out_name_matr_b, "./Debug/pc_matr_B_%d_0_%d.dat", il, ia);
834  sprintf(out_name_matr_c, "./Debug/pc_matr_C_%d_0_%d.dat", il, ia);
835  matr_A.writeToFile(out_name_matr_a);
836  matr_B.writeToFile(out_name_matr_b);
837  matr_C.writeToFile(out_name_matr_c);
838  Output::echo(0, "done.\n");*/
839 
840  // make RHS = B*f + C
841  DiagMatrix::iterator it;
842  for (in = 0; in < pc.size; in++) {
843  // I guess matr_B has only main diagonal in implicit method, so it can be used like this:
844  // RHS[in] = matr_B[0][in] * slicePSD1D[in] + matr_C[0][in];
845  RHS[in] = matr_C[0][in];
846  for (it = matr_B.begin(); it != matr_B.end(); it++) {
847  // multiplication B * f
848  if (in + it->first >= 0 && in + it->first < matr_B.total_size) {
849  RHS[in] += it->second[in] * slicePSD1D[in + it->first];
850  }
851  }
852  }
853 
854  // right hand side to files
855  // RHS.writeToFile("pc_RHS.dat");
856 
857  // Solving model matrix, obtaining PSD
858  // Assume matr_A has only +1, 0, -1 diagonals
859  // TODO: Check if matr_A has only +1, 0, -1 diagonals, i.e. is tridiagonal
860  tridag(
861  &matr_A[-1][0], // just a way to get a pointer to first element of an array
862  &matr_A[0][0],
863  &matr_A[+1][0],
864  &RHS[0],
865  &slicePSD1D[0],
866  pc.size);
867 
868  // copying PSD back, from 1d array to 3d array
869  for (im = 0; im < pc.size; im++) {
870  arr[il][im][ia] = slicePSD1D[im];
871  }
872  }
873  }
874 }
875 
876 
895 void PSD::Diffusion_L( double dt, double Lpp,
897  GridElement &L, GridElement &pc, GridElement &alpha, Matrix3D<double> Jacobian,
898  Matrix2D<double> L_lowerBoundaryCondition,
899  Matrix2D<double> L_upperBoundaryCondition,
900  string L_lowerBoundaryCondition_calculationType,
901  string L_upperBoundaryCondition_calculationType,
902  double tau, double tauLpp)
903 {
904  // no longer used
905  // char out_name_matr_a[256], out_name_matr_b[256], out_name_matr_c[256];
906 
907 
908  // Calculation matrices A, B, and C for 1d diffusion problem
909  // A * PSD_1D(t+1) = B * PSD_1D(t) + C, for each L and Energy
910 
911  // Creation is moved to PSD.h, so it can be accessed from Main.cpp
912  //static CalculationMatrix
913  // matr_A(1, 1, L.size, 1),
914  // matr_B(1, 1, L.size, 0),
915  // matr_C(1, 1, L.size, 0);
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);
919 
920  // Create zero matrices, temporary matrices to be used later
921  static Matrix3D<double> Zero_Matrix(L.size, 1, 1);
922  static Matrix2D<double> Zero_Matrix_2d(1, 1);
923  static Matrix1D<double> Zero_Matrix_1d(1);
924  Zero_Matrix = 0;
925  Zero_Matrix_2d = 0;
926  Zero_Matrix_1d = 0;
927 
928  // Create 1d arrays
929  // Creating 1d arrays for PSD and right hand side (RHS = B * PSD + C)
930  static Matrix1D<double> slicePSD1D(L.size), RHS(L.size);
931 
932  // Creating 1d arrays which will be send to the MakeMatrix function
933  // Creating 1d matrixes, but with type of Matrix3D, cause we need that type for MakeMatrix function
934  // Create 3d and 2d matrices instead of 1d and single values so that it can meet parameters in MakeModelMatrix_3D()
935  static Matrix3D<double> DLL_1d(L.size, 1, 1), L_1d(L.size, 1, 1), Jacobian_1d(L.size, 1, 1); // pc_1d(L.size, 1, 1), alpha_1d(L.size, 1, 1)
936  // Creating temporary arrays size 1x1 for boundary conditions, which are just values
937  static Matrix2D<double> L_lowerBoundaryCondition_1d(1, 1), L_upperBoundaryCondition_1d(1, 1);
938 
939  // iterators
940  int il, im, ia, in; // matrix element numbers on 3D grid
941  DiagMatrix::iterator it;
942  for (ia = 0; ia < alpha.size; ia++) {
943  for (im = 0; im < pc.size; im++) {
944 
945  // Copying 3D arrays into 1d arrays
946  for (il = 0; il < L.size; il++) {
947  // Copying 3D PSD to 1d PSD
948  slicePSD1D[il] = arr[il][im][ia];
949 
950  DLL_1d[il][0][0] = DLL.arr[il][im][ia];
951 
952  L_1d[il][0][0] = L.arr[il][im][ia];
953  // pc_1d[il][0][0] = pc[il][im][ia];
954  // alpha_1d[il][0][0] = alpha[il][im][ia];
955 
956  Jacobian_1d[il][0][0] = Jacobian[il][im][ia];
957  }
958  // boundary conditions, these are arrays size 1x1, i.e. values
959  L_upperBoundaryCondition_1d[0][0] = L_upperBoundaryCondition[im][ia];
960  L_lowerBoundaryCondition_1d[0][0] = L_lowerBoundaryCondition[im][ia];
961 
962  // Make A, B, C matrices
964  matr_A, matr_B, matr_C, // matrices A, B, C
965  L_1d, Zero_Matrix, Zero_Matrix, // L, pc, alpha grids. We use L grid only here
966  L.size, 1, 1, // L-size, pc-size, alpha-size
967  L_lowerBoundaryCondition_1d, L_upperBoundaryCondition_1d, // L-boundaries
968  Zero_Matrix_2d, Zero_Matrix_2d, // pc-boundaries
969  Zero_Matrix_2d, Zero_Matrix_2d, // alpha-boundaries
970  L_lowerBoundaryCondition_calculationType, L_upperBoundaryCondition_calculationType, // Boundary conditions at L-boundaries
971  "BCT_CONSTANT_VALUE", "BCT_CONSTANT_VALUE", // Boundary conditions at pc-boundaries
972  "BCT_CONSTANT_VALUE", "BCT_CONSTANT_VALUE", // Boundary conditions at alpha-boundaries
973  DLL_1d, // DLL
974  Zero_Matrix, Zero_Matrix, // Dpp, DppLpp - use 0
975  Zero_Matrix, Zero_Matrix, // Daa, DaaLpp - use 0
976  Zero_Matrix, Zero_Matrix, // Dpca, Dpca_Lpp - use 0
977  Jacobian_1d, dt, Lpp,
978  tau, tauLpp);
979 
980  /*
981  Output::echo(0, "Writing calculation matrices. ");
982  sprintf(out_name_matr_a, "./Debug/L_matr_A_0_%d_%d.dat", im, ia);
983  sprintf(out_name_matr_b, "./Debug/L_matr_B_0_%d_%d.dat", im, ia);
984  sprintf(out_name_matr_c, "./Debug/L_matr_C_0_%d_%d.dat", im, ia);
985  matr_A.writeToFile(out_name_matr_a);
986  matr_B.writeToFile(out_name_matr_b);
987  matr_C.writeToFile(out_name_matr_c);
988  Output::echo(0, "done.\n");
989  */
990 
991  // make RHS = B*f + C
992  DiagMatrix::iterator it;
993  for (in = 0; in < L.size; in++) {
994  // I guess matr_B has only main diagonal in implicit method, so it can be used like this:
995  // RHS[in] = matr_B[0][in] * slicePSD1D[in] + matr_C[0][in];
996  RHS[in] = matr_C[0][in];
997  for (it = matr_B.begin(); it != matr_B.end(); it++) {
998  // multiplication B * f
999  if (in + it->first >= 0 && in + it->first < matr_B.total_size) {
1000  RHS[in] += it->second[in] * slicePSD1D[in + it->first];
1001  }
1002  }
1003  }
1004 
1005  // right hand side to files
1006  // RHS.writeToFile("pc_RHS.dat");
1007 
1008  // Solving model matrix, obtaining PSD
1009  // Assume matr_A has only +1, 0, -1 diagonals
1010  // TODO: Check if matr_A has only +1, 0, -1 diagonals, i.e. is tridiagonal
1011  tridag(
1012  &matr_A[-1][0], // just a way to get a pointer to first element of an array
1013  &matr_A[0][0],
1014  &matr_A[+1][0],
1015  &RHS[0],
1016  &slicePSD1D[0],
1017  L.size);
1018 
1019  // copying PSD back, from 1d array to 3d array
1020  for (il = 0; il < L.size; il++) {
1021  arr[il][im][ia] = slicePSD1D[il];
1022  }
1023  }
1024  }
1025 }
1026 
1027 
1053 void PSD::Diffusion_pc_alpha(double dt, double Lpp,
1054  DiffusionCoefficient &Dpcpc, DiffusionCoefficient &DpcpcLpp,
1057  GridElement &L, GridElement &pc, GridElement &alpha, Matrix3D<double> Jacobian,
1058  Matrix2D<double> pc_lowerBoundaryCondition,
1059  Matrix2D<double> pc_upperBoundaryCondition,
1060  Matrix2D<double> alpha_lowerBoundaryCondition,
1061  Matrix2D<double> alpha_upperBoundaryCondition,
1062  string pc_lowerBoundaryCondition_calculationType,
1063  string pc_upperBoundaryCondition_calculationType,
1064  string alpha_lowerBoundaryCondition_calculationType,
1065  string alpha_upperBoundaryCondition_calculationType)
1066 {
1067 
1068  // Calculation matrices A, B, and C
1069  // A * PSD_1D(t+1) = B * PSD_1D(t) + C, for each L and Energy
1070 
1071  // Creation is moved to PSD.h, so it can be accessed from Main.cpp
1072  //static CalculationMatrix
1073  // matr_A(alpha.size, pc.size, 1, 1),
1074  // matr_B(alpha.size, pc.size, 1, 0),
1075  // matr_C(alpha.size, pc.size, 1, 0);
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);
1079 
1080  // Create zero matrices, temporary matrices to be used later
1081  static Matrix3D<double> Zero_Matrix(1, pc.size, alpha.size);
1082  static Matrix2D<double> Zero_Matrix_2d(1, 1);
1083  static Matrix1D<double> Zero_Matrix_1d(1);
1084  Zero_Matrix = 0;
1085  Zero_Matrix_2d = 0;
1086  Zero_Matrix_1d = 0;
1087 
1088  // Create 1d arrays
1089  // Creating 1d arrays for PSD and right hand side (RHS = B * PSD + C)
1090  static Matrix1D<double> slicePSD1D(matr_A.total_size), RHS(matr_A.total_size);
1091  // Creating 1d arrays which will be send to the MakeMatrix function
1092  // Creating 1d matrixes, but with type of Matrix3D, cause we need that type for MakeMatrix function
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);
1095  // Creating temporary arrays for boundary conditions
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);
1098 
1099  // iterators
1100  int il, im, ia, in; // matrix element numbers on 3D grid
1101  DiagMatrix::iterator it;
1102  for (il = 0; il < L.size; il++) {
1103  // Copying 3D arrays into 1d arrays
1104  for (im = 0; im < pc.size; im++) {
1105  for (ia = 0; ia < alpha.size; ia++) {
1106  // Copying 3D PSD to 1d PSD
1107  in = matr_A.index1d(ia, im);
1108  slicePSD1D[in] = arr[il][im][ia];
1109 
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];
1116 
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];
1120 
1121  Jacobian_slice[0][im][ia] = Jacobian[il][im][ia];
1122  }
1123  }
1124  // boundary conditions
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];
1128  }
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];
1132  }
1133 
1134  // Make
1136  matr_A, matr_B, matr_C, // matrices A, B, C
1137  L_slice, pc_slice, alpha_slice, // L, pc, alpha grids. We use pc and alpha crossections here
1138  1, pc.size, alpha.size, // L-size, pc-size, alpha-size
1139  Zero_Matrix_2d, Zero_Matrix_2d, // L-boundaries
1140  pc_lowerBoundaryCondition_slice, pc_upperBoundaryCondition_slice, // pc-boundaries
1141  alpha_lowerBoundaryCondition_slice, alpha_upperBoundaryCondition_slice, // alpha-boundaries
1142  "BCT_CONSTANT_VALUE", "BCT_CONSTANT_VALUE", // Boundary conditions at L-boundaries
1143  pc_lowerBoundaryCondition_calculationType, pc_upperBoundaryCondition_calculationType, // Boundary conditions at pc-boundaries
1144  alpha_lowerBoundaryCondition_calculationType, alpha_upperBoundaryCondition_calculationType, // Boundary conditions at alpha-boundaries
1145  Zero_Matrix, // DLL
1146  Dpcpc_slice, DpcpcLpp_slice, // Dpp, DppLpp - use 0
1147  Daa_slice, DaaLpp_slice, // Daa, DaaLpp - use 0
1148  Dpca_slice, DpcaLpp_slice, // Dpca, Dpca_Lpp - use 0
1149  Jacobian_slice, dt, Lpp);
1150 
1151  /* Output::echo(0, "Writing calculation matrices for L[%d] = %f... ", il, L[il][0][0]);
1152  matr_A.writeToFile("./Debug/pc_alpha_matr_A.dat");
1153  matr_B.writeToFile("./Debug/pc_alpha_matr_B.dat");
1154  matr_C.writeToFile("./Debug/pc_alpha_matr_C.dat");
1155  Output::echo(0, "done.\n"); */
1156 
1157  // make RHS = B*f + C
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);
1162  // I guess matr_B has only main diagonal in implicit method, so it can be used like this:
1163  // RHS[in] = matr_B[0][in] * slicePSD1D[in] + matr_C[0][in];
1164  RHS[in] = matr_C[0][in]; // RHS = C
1165  for (it = matr_B.begin(); it != matr_B.end(); it++) {
1166  // multiplication B * f
1167  if (in + it->first >= 0 && in + it->first < matr_B.total_size) {
1168  RHS[in] += it->second[in] * slicePSD1D[in + it->first]; // RHS = C + Bf
1169  }
1170  }
1171 
1172  }
1173  }
1174 
1175  // right hand side to files
1176  // RHS.writeToFile("./Debug/pc_alpha_RHS.dat");
1177 
1178  // Solving model matrix
1179  int n = pc.size * alpha.size;
1180 
1181  if (PSD_parameters.solutionMethod == "SM_Gauss") {
1182  // external gauss function
1183  gauss(matr_A, RHS, slicePSD1D, n);
1184 
1185  Output::echo(5, "done.\n");
1186 
1187  /* moved this block of code to helper function gauss() located in MatrixSolver */
1188 
1189  } else if (PSD_parameters.solutionMethod == "SM_Relaxation") {
1190  // Relaxation method
1191 
1192  Output::echo(5, "Inverting A matrix (relaxation)... ");
1193  // runs the over_relaxation_diag solution method
1194  over_relaxation_diag(matr_A, RHS, slicePSD1D);
1195  Output::echo(5, "done.");
1196 
1197  } else if (PSD_parameters.solutionMethod == "SM_Lapack") {
1198  // Relaxation method
1199 
1200  Output::echo(5, "Inverting A matrix (Lapack)... ");
1201  // runs the over_relaxation_diag solution method
1202  Lapack(matr_A, RHS, slicePSD1D);
1203  Output::echo(5, "done.");
1204 
1205  } else if (PSD_parameters.solutionMethod == "SM_GMRES") {
1206  // Relaxation method
1207 
1208  Output::echo(5, "Inverting A matrix (GMRES)... ");
1209  // runs the over_relaxation_diag solution method
1210  //gmres_wrapout(matr_A, RHS, slicePSD1D, PSD_parameters.SOL_maxiter, PSD_parameters.SOL_i_max, PSD_parameters.SOL_max_iter_err);
1211  gmres_wrapout(matr_A, RHS, slicePSD1D, PSD_parameters.GMRES_parameters);
1212  Output::echo(5, "done.");
1213 
1214  } else if (PSD_parameters.solutionMethod == "SM_GMRES2") {
1215  // Relaxation method
1216 
1217  Output::echo(5, "Inverting A matrix (GMRES2)... ");
1218  // runs the over_relaxation_diag solution method
1219  //gmres_wrapout(matr_A, RHS, slicePSD1D, PSD_parameters.SOL_maxiter, PSD_parameters.SOL_i_max, PSD_parameters.SOL_max_iter_err);
1220  gmres2_wrapout(matr_A, RHS, slicePSD1D, PSD_parameters.GMRES_parameters);
1221  Output::echo(5, "done.");
1222 
1223 
1224  } else if (PSD_parameters.solutionMethod == "SM_Tridiag") {
1225  // Tridiagonal method
1226  throw error_msg("SOLUTIONERR", "Tridiagonal method can't be used");
1227  } else { // Unknown method
1228  throw error_msg("SOLUTIONERR", "Unknown solution method");
1229  }
1230 
1231  // copying PSD back, from 1d array to 3d array
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];
1236  }
1237  }
1238  }
1239 
1240 }
1241 
1242 
1269 void PSD::Diffusion_pc_alpha_KC(double dt, double Lpp,
1270  DiffusionCoefficient &Dpcpc, DiffusionCoefficient &DpcpcLpp,
1273  GridElement &L, GridElement &pc, GridElement &alpha, Matrix3D<double> Jacobian,
1274  Matrix2D<double> pc_lowerBoundaryCondition,
1275  Matrix2D<double> pc_upperBoundaryCondition,
1276  Matrix2D<double> alpha_lowerBoundaryCondition,
1277  Matrix2D<double> alpha_upperBoundaryCondition,
1278  string pc_lowerBoundaryCondition_calculationType,
1279  string pc_upperBoundaryCondition_calculationType,
1280  string alpha_lowerBoundaryCondition_calculationType,
1281  string alpha_upperBoundaryCondition_calculationType)
1282 {
1283 
1284  // Calculation matrices A, B, and C
1285  // A * PSD_1D(t+1) = B * PSD_1D(t) + C, for each L and Energy
1286 
1287  // Creation is moved to PSD.h, so it can be accessed from Main.cpp
1288  //static CalculationMatrix
1289  // matr_A(alpha.size, pc.size, 1, 1),
1290  // matr_B(alpha.size, pc.size, 1, 0),
1291  // matr_C(alpha.size, pc.size, 1, 0);
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);
1295 
1296 
1297 
1298 
1299  // Create zero matrices, temporary matrices to be used later
1300  static Matrix3D<double> Zero_Matrix(1, pc.size, alpha.size);
1301  static Matrix2D<double> Zero_Matrix_2d(1, 1);
1302  static Matrix1D<double> Zero_Matrix_1d(1);
1303  Zero_Matrix = 0;
1304  Zero_Matrix_2d = 0;
1305  Zero_Matrix_1d = 0;
1306 
1307  // Create 1d arrays
1308  // Creating 1d arrays for PSD and right hand side (RHS = B * PSD + C)
1309  static Matrix1D<double> slicePSD1D(matr_A.total_size), RHS(matr_A.total_size);
1310  // Creating 1d arrays which will be send to the MakeMatrix function
1311  // Creating 1d matrixes, but with type of Matrix3D, cause we need that type for MakeMatrix function
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);
1314  // Creating temporary arrays for boundary conditions
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);
1317 
1318  // iterators
1319  int il, im, ia, in; // matrix element numbers on 3D grid
1320  DiagMatrix::iterator it;
1321  for (il = 0; il < L.size; il++) {
1322  // Copying 3D arrays into 1d arrays
1323  for (im = 0; im < pc.size; im++) {
1324  for (ia = 0; ia < alpha.size; ia++) {
1325  // Copying 3D PSD to 1d PSD
1326  in = matr_A.index1d(ia, im);
1327  slicePSD1D[in] = arr[il][im][ia];
1328 
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];
1335 
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];
1339 
1340  Jacobian_slice[0][im][ia] = Jacobian[il][im][ia];
1341  }
1342  }
1343  // boundary conditions
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];
1347  }
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];
1351  }
1352 
1353  // Make
1355  matr_A, matr_B, matr_C, // matrices A, B, C
1356  L_slice, pc_slice, alpha_slice, // L, pc, alpha grids. We use pc and alpha crossections here
1357  1, pc.size, alpha.size, // L-size, pc-size, alpha-size
1358  Zero_Matrix_2d, Zero_Matrix_2d, // L-boundaries
1359  pc_lowerBoundaryCondition_slice, pc_upperBoundaryCondition_slice, // pc-boundaries
1360  alpha_lowerBoundaryCondition_slice, alpha_upperBoundaryCondition_slice, // alpha-boundaries
1361  "BCT_CONSTANT_VALUE", "BCT_CONSTANT_VALUE", // Boundary conditions at L-boundaries
1362  pc_lowerBoundaryCondition_calculationType, pc_upperBoundaryCondition_calculationType, // Boundary conditions at pc-boundaries
1363  alpha_lowerBoundaryCondition_calculationType, alpha_upperBoundaryCondition_calculationType, // Boundary conditions at alpha-boundaries
1364  Zero_Matrix, // DLL
1365  Dpcpc_slice, DpcpcLpp_slice, // Dpp, DppLpp - use 0
1366  Daa_slice, DaaLpp_slice, // Daa, DaaLpp - use 0
1367  Dpca_slice, DpcaLpp_slice, // Dpca, Dpca_Lpp - use 0
1368  Jacobian_slice, dt, Lpp);
1369 
1370  /* Output::echo(0, "Writing calculation matrices for L[%d] = %f... ", il, L[il][0][0]);
1371  matr_A.writeToFile("./Debug/pc_alpha_matr_A.dat");
1372  matr_B.writeToFile("./Debug/pc_alpha_matr_B.dat");
1373  matr_C.writeToFile("./Debug/pc_alpha_matr_C.dat");
1374  Output::echo(0, "done.\n"); */
1375 
1376  // make RHS = B*f + C
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);
1381  // I guess matr_B has only main diagonal in implicit method, so it gan be used like this:
1382  // RHS[in] = matr_B[0][in] * slicePSD1D[in] + matr_C[0][in];
1383  RHS[in] = matr_C[0][in];
1384  for (it = matr_B.begin(); it != matr_B.end(); it++) {
1385  // multiplication B * f
1386  if (in + it->first >= 0 && in + it->first < matr_B.total_size) {
1387  RHS[in] += it->second[in] * slicePSD1D[in + it->first];
1388  }
1389  }
1390 
1391  }
1392  }
1393 
1394  // right hand side to files
1395  // RHS.writeToFile("./Debug/pc_alpha_RHS.dat");
1396 
1397  // Solving model matrix
1398  int n = pc.size * alpha.size;
1399  // need A and B for gauss method
1400  static Matrix2D<double> A;
1401  static Matrix1D<double> B;
1402  int modelMatrixLineNumber;
1403  int k, l, in;
1404  double sum, error, error_max;
1405  static bool recalculate_A = true;
1406  if (PSD_parameters.solutionMethod == "SM_Gauss") {
1407  // the routine does not actually invert matrix, but just solve AX=B eq
1408  if (recalculate_A) {
1409  B = RHS;
1410  if (!A.initialized) {
1411  A.AllocateMemory(n, n);
1412  }
1413  A = 0;
1414  // Output::echo("Writing A0 matrix... ");
1415  // A.writeToFile("A0.dat");
1416  // Output::echo("done.\n");
1417 
1418  Output::echo(5, "Inverting A matrix (gauss)... ");
1419  // Gauss method
1420  for (it = matr_A.begin(); it != matr_A.end(); it++) {
1421  for (modelMatrixLineNumber = 0; modelMatrixLineNumber < n; modelMatrixLineNumber++) {
1422  // check, if the element at line (modelMatrixLineNumber) and diagonal (it->first) is inside the matrix.
1423  // cause it can be outside, if line number is 1 and diagonal number is -5, for example.
1424  if (modelMatrixLineNumber + it->first >= 0 && modelMatrixLineNumber + it->first < n) {
1425  // converting matrix, stored as diagonals, into normal 2D matrix
1426  // cause we need 2D matrix to solve by Gauss method
1427  // cout << it->first << " - " << it->second[modelMatrixLineNumber] << endl;
1428  A[modelMatrixLineNumber][modelMatrixLineNumber + it->first] = it->second[modelMatrixLineNumber];
1429  }
1430  }
1431  }
1432 
1433  // inversion of A
1434  // sending address of the first element of the plane 2d array
1435  // it's not actually matrix inversion
1436  // it just makes A-triangle, so it's easy to solve A*X = B
1437  gauss_solve(&A[0][0], &B[0], n);
1438  recalculate_A = true;
1439  Output::echo(5, "done.\n");
1440  }
1441 
1442  Output::echo(5, "Calculation of new PSD... ");
1443  slicePSD1D[n-1] = B[n-1] / A[n-1][n-1];
1444  for (k = n - 2; k >= 0; k--) {
1445  sum = 0;
1446  for (l = k + 1; l < n; l++)
1447  sum += A[k][l] * slicePSD1D[l];
1448  slicePSD1D[k] = B[k] - sum;
1449  }
1450  /* end calc */
1451 
1452  /* control */
1453  //mult_vect2(&A[0], &slicePSD1D[0], &RHS[0], n);
1454  error_max = 0.;
1455  error = 0;
1456  for (in = 0; in < n; in++) {
1457  error = RHS[in];
1458  for (it = matr_A.begin(); it != matr_A.end(); it++) {
1459  // multiplying all diagonals to PSD to check error = RHS - A*PSD = 0
1460  if (in + it->first >= 0 && in + it->first < n) {
1461  error -= it->second[in] * slicePSD1D[in + it->first];
1462  }
1463  }
1464  error_max = VF::max( error_max, error );
1465  }
1466  Output::echo(5, "done.\n");
1467  Output::echo(5, "Error_norm = %e \n", error_max);
1468 
1469  } else if (PSD_parameters.solutionMethod == "SM_Relaxation") {
1470  // Relaxation method
1471 
1472  Output::echo(5, "Inverting A matrix (relaxation)... ");
1473  // runs the over_relaxation_diag solution method
1474  over_relaxation_diag(matr_A, RHS, slicePSD1D);
1475  Output::echo(5, "done.");
1476 
1477  } else if (PSD_parameters.solutionMethod == "SM_Lapack") {
1478  // Relaxation method
1479 
1480  Output::echo(5, "Inverting A matrix (Lapack)... ");
1481  // runs the over_relaxation_diag solution method
1482  Lapack(matr_A, RHS, slicePSD1D);
1483  Output::echo(5, "done.");
1484 
1485  } else if (PSD_parameters.solutionMethod == "SM_GMRES") {
1486  // Relaxation method
1487 
1488  Output::echo(5, "Inverting A matrix (GMRES)... ");
1489  // runs the over_relaxation_diag solution method
1490  //gmres_wrapout(matr_A, RHS, slicePSD1D, PSD_parameters.SOL_maxiter, PSD_parameters.SOL_i_max, PSD_parameters.SOL_max_iter_err);
1491  gmres_wrapout(matr_A, RHS, slicePSD1D, PSD_parameters.GMRES_parameters);
1492  Output::echo(5, "done.");
1493 
1494  } else if (PSD_parameters.solutionMethod == "SM_GMRES2") {
1495  // Relaxation method
1496 
1497  Output::echo(5, "Inverting A matrix (GMRES2)... ");
1498  // runs the over_relaxation_diag solution method
1499  //gmres_wrapout(matr_A, RHS, slicePSD1D, PSD_parameters.SOL_maxiter, PSD_parameters.SOL_i_max, PSD_parameters.SOL_max_iter_err);
1500  gmres2_wrapout(matr_A, RHS, slicePSD1D, PSD_parameters.GMRES_parameters);
1501  Output::echo(5, "done.");
1502 
1503 
1504  } else if (PSD_parameters.solutionMethod == "SM_Tridiag") {
1505  // Tridiagonal method
1506  throw error_msg("SOLUTIONERR", "Tridiagonal method can't be used");
1507  } else { // Unknown method
1508  throw error_msg("SOLUTIONERR", "Unknown solution method");
1509  }
1510 
1511  // copying PSD back, from 1d array to 3d array
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];
1516  }
1517  }
1518  }
1519 
1520 }
1521 
1522 
1523 
1532 void checkInf_3D(Matrix3D<double> arr, int il, int im, int ia, double maxNum = 1e99) {
1533  // Some weird condition... It tries to check for Not-A-Number which gives false for all comparisons and same time for infinite values... But I guess there is a bug here.
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);
1537  }
1538 }
1539 
1548 void checkInf_1D(Matrix1D<double> arr, int il, int im, int ia, double maxNum = 1e99) {
1549  // Some weird condition... It tries to check for Not-A-Number which gives false for all comparisons and same time for infinite values... But I guess there is a bug here.
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);
1553  }
1554 }
1555 
1557 void PSD::Interpolate(PSD &otherPSD, Parameters_structure::Interpolation interpolationParameters_structure, Grid &oldGrid, Grid &newGrid, Matrix2D<double> newGrid_pc_lowerBoundaryCondition, Matrix2D<double> newGrid_pc_upperBoundaryCondition) {
1558  // check...
1559 
1560  int il, im, ia;
1561  double lowerBoundary1d=0, upperBoundary1d=0;
1562 
1563  // if interpolation method is IT_NONE - exit
1564  // if grids have the same type, we just copy the values
1565  if (interpolationParameters_structure.type == "IT_NONE") {
1566  return;
1567  } else if (oldGrid.type == newGrid.type) {
1568  // arr.Matrix3D<double>::operator =(otherPSD);
1569 
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++) {
1573  // copying of the values from one grid to anouther
1574  arr[il][im][ia] = otherPSD.arr[il][im][ia];
1575  // some weird bug. Model matrix does not work correctly without this. Don't know why.
1576  if (arr[il][im][ia] < VC::zero_f) {
1577  arr[il][im][ia] = VC::zero_f;
1578  }
1579  }
1580  }
1581  }
1582  return;
1583  }
1584 
1585  // matrices to be interpolated
1586  Matrix1D<double> matrix_array_1d_old(oldGrid.pc.size);
1587  Matrix1D<double> matrix_array_1d_new(newGrid.pc.size);
1588  Matrix1D<double> grid_1d_old(oldGrid.pc.size);
1589  Matrix1D<double> grid_1d_new(newGrid.pc.size);
1590 
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++) {
1594 
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]) {
1596  // if L-size or alpha-size are different, then we can not interpolate by that interpolation function
1597  throw error_msg("INTERPOLATION_ERROR", "Interpolation error: grids must have the same L and Alpha");
1598  }
1599 #ifdef DEBUG_MODE
1600  // check for infinity and NAN
1601  checkInf_3D(arr, il, im, ia);
1602 #endif
1603 
1604  // some work around boundary conditions
1605  if (interpolationParameters_structure.useLog == "Log_10") {
1606  // if we are interpolating log of values, then store log of values to 1d array
1607  if (otherPSD.arr[il][im][ia] >= VC::zero_f) {
1608  matrix_array_1d_old[im] = log10(otherPSD.arr[il][im][ia]);
1609  } else {
1610  // if PSD lower then lowest possible value (it should not be zero, because we are interpolating log of values for example)
1611  matrix_array_1d_old[im] = log10(VC::zero_f);
1612  }
1613  // theoretically we don't need to make the boundary conditions lower then 10^21
1614  // because we check that in its constructor
1615  // But just in case...
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));
1618  }
1619  else if (interpolationParameters_structure.useLog == "Log_2") {
1620  // if we are interpolating log of values, then store log of values to 1d array
1621  if (otherPSD.arr[il][im][ia] >= VC::zero_f) {
1622  matrix_array_1d_old[im] = log(otherPSD.arr[il][im][ia])/log(2.0);
1623  } else {
1624  // if PSD lower then lowest possible value (it should not be zero, because we are interpolating log of values for example)
1625  matrix_array_1d_old[im] = log(VC::zero_f)/log(2.0);
1626  }
1627  // theoretically we don't need to make the boundary conditions lower then 10^21
1628  // because we check that in its constructor
1629  // But just in case...
1630 
1631  // adrozdov: Bug report. log10 replaced by log2
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);
1634  }
1635  else if (interpolationParameters_structure.useLog == "Log_E") {
1636  // if we are interpolating log of values, then store log of values to 1d array
1637  if (otherPSD.arr[il][im][ia] >= VC::zero_f) {
1638  matrix_array_1d_old[im] = log(otherPSD.arr[il][im][ia]);
1639  } else {
1640  // if PSD lower then lowest possible value (it should not be zero, because we are interpolating log of values for example)
1641  matrix_array_1d_old[im] = log(VC::zero_f);
1642  }
1643  // theoretically we don't need to make the boundary conditions lower then 10^21
1644  // because we check that in its constructor
1645  // But just in case...
1646 
1647  // adrozdov: Bug report. log10 replaced by log
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));
1650  }
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);
1655  }
1656  else {
1657  throw error_msg("INTERPOLATION_ERROR", "Usage of log(PSD) is not define properly: %s", interpolationParameters_structure.useLog.c_str());
1658  break;
1659  }
1660 
1661 #ifdef DEBUG_MODE
1662  // check for infinity and NAN
1663  checkInf_1D(matrix_array_1d_old, il, im, ia, 99);
1664 #endif
1665 
1666  grid_1d_old[im] = oldGrid.pc.arr[il][im][ia];
1667  grid_1d_new[im] = newGrid.pc.arr[il][im][ia];
1668  //cout<<grid_1d_old[im]<<"kim......."<<grid_1d_new[im]<<"..........."<<matrix_array_1d_old[im]<<endl;
1669  }
1670 
1671  // interpolation based on type
1672  if (interpolationParameters_structure.type == "IT_LINEAR") {
1673  // do linear interpolation
1674  matrix_array_1d_new.Interpolate(matrix_array_1d_old, grid_1d_old, grid_1d_new);
1675  } else if (interpolationParameters_structure.type == "IT_POLYNOMIAL") {
1676  // do polynomial - linear interpolation (linear close to the boundaries)
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") {
1679  // do spline interpolation
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") {
1682  // just copy
1683  matrix_array_1d_new = matrix_array_1d_old;
1684  } else { // unknown
1685  throw error_msg("INTERPOLATION_TYPE_ERROR", "Unknown interpolation type");
1686  }
1687 
1688  for (im = 0; im < newGrid.pc.size; im++) {
1689 
1690 #ifdef DEBUG_MODE
1691  checkInf_1D(matrix_array_1d_new, il, im, ia, 99);
1692 #endif
1693 
1694  // copying values back to 3d arrays
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];
1699 
1700 #ifdef DEBUG_MODE
1701  checkInf_3D(arr, il, im, ia);
1702 #endif
1703 
1704  }
1705  }
1706  }
1707 
1708  // check result...
1709 }
1710 
1711 /*-------------------------------------------------------------------------------------------------------------------------------
1712 Outputs
1713 -------------------------------------------------------------------------------------------------------------------------------*/
1714 
1720 void PSD::Output_without_grid(double time) {
1721  int il, im, ia;
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;
1727  }
1728  }
1729  }
1730 
1731 }
1732 
1733 
1742 void PSD::Load_initial_f_file(GridElement &L, GridElement &epc, GridElement &alpha, const char *filename, bool withGrid) {
1743  int il, im, ia;
1744  string inBuf;
1745  double Load_L=-1, Load_epc=-1, Load_alpha=-1;// temporary variables for the grid point
1746  double err = 1.e-3;
1747 
1748  // Load j from input file
1749  ifstream input_f(filename);
1750  if (input_f != NULL && !input_f.eof()) {
1751  // Skipping first two lines. Should search for 'ZONE' and read from following line better.
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];
1759  if (withGrid)
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]);
1762  }
1763  }
1764  }
1765  }
1766  } else {
1767  throw error_msg("INITIAL_PSD_ERROR", "Error reading PSD input file %s.\n", filename);
1768  //cout << "Error reading initial conditions input file " << filename << endl;
1769  //getchar();
1770  //exit(-1);
1771  }
1772  input_f.close();
1773 }
1774 
1775 
1790 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, bool using_ia_90) {
1791  int il, im, ia;
1792  double p1, Lmax, Lmin, scale_factor, Lgrid;
1793  int Lsize_7, addGrid; // These two numbers should be of integer type.
1794  Matrix1D<double> L_1d, fL;
1795  //Matrix2D<double> fLm(L.size, pc.size);
1796 
1798  // Steady state solution for any outer-L boundary, D. Subbotin implimentation
1799  // Create 1D grid at min energy and 90 deg alpha
1800  // First, find a closest index of 90 deg pitch-angle
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--;
1803 
1804  // Make sure we have a point at L=7
1805  // Our boundary values J_L7() are at L=7, so we need to know the steady state solution at L=7.
1806  if (L.arr[L.size-1][0][ia_90] < 7) {
1807  // If maximum L value is less than 7, we add another points from Lmax to 7
1808 
1809  // estimate grid points to add
1810  if (using_ia_90) {
1811  // L values set based on ia_90
1812  Lmax = L.arr[L.size-1][0][ia_90];
1813  Lmin = L.arr[0][0][ia_90];
1814  } else {
1815  Lmax = L.arr[L.size-1][0][0];
1816  Lmin = L.arr[0][0][0];
1817  }
1818 
1819  Lgrid = (Lmax-Lmin)/(L.size-1);
1820  Lsize_7 = (7.0-Lmin)/Lgrid + 1;
1821  addGrid = Lsize_7 - L.size;
1822 
1823  // construct new L and psd grid
1824  L_1d.AllocateMemory(L.size+addGrid);
1825  fL.AllocateMemory(L.size+addGrid);
1826 
1827  for (il = 0; il < L.size; il++) {
1828  L_1d[il] = L.arr[il][0][ia_90];
1829  }
1830  for (il = L.size; il < Lsize_7; il++) { // add another points from Lmax to 7
1831  L_1d[il] = Lmin+Lgrid*il;
1832  }
1833 
1834  } else {
1835  Lsize_7 = L.size;
1836  L_1d.AllocateMemory(L.size);
1837  fL.AllocateMemory(L.size);
1838  for (il = 0; il < L.size; il++) L_1d[il] = L.arr[il][0][ia_90];
1839 
1840  }
1841 
1842  steady_state(fL, tau, Kp, Lsize_7, L_1d, fb_out, fb_in);
1843  if (using_ia_90) {
1844  scale_factor = 1; // Bugfix. In terms of usage bf.txt as Flux/MeanFlux PSD should be 1 at L==7
1845  } else {
1846  scale_factor = fb_out/fL[L.size-1]; // scale factor to make psd to be 1 at the outer boundary (=Lmax(<7))
1847  }
1848 
1849  // If max(L)>7, normalize, so the steady state fL(L=7) = 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));
1854  fL = fL / fL(il_7);
1855  }
1856  //
1858 
1859  // if PSD is lower then min_f - make it min_f
1860  for (int i = 0; i < fL.size_x; i++) if (fL[i] < min_f) fL[i] = min_f;
1861 
1862 
1863  // inner belt simulation
1864  // int il_2=L.size-1;
1865  // while (L[il_2][0][0] > 2 && il_2 > 0) il_2--;
1866  // for (int i = il_2; i >= 0; i--) fL[i] = min_f/pc[i][0][0]/pc[i][0][0]*pc[il_2][0][0]*pc[il_2][0][0];
1867 
1868 
1869  //#ifdef DEBUG_MODE
1870  fL.writeToFile("./Debug/fL.dat");
1871  //#endif
1872 
1873  for (im = 0; im < pc.size; im++) {
1874  for (il = 0; il < L.size; il++) {
1875  for (ia = 0; ia < alpha.size; ia++) {
1876  // filling for all pc
1877  // Why do we do it for alpha=90?:
1878  // Each point is calculated at the equator and multiplied by sin(pitch angle).
1879  // BUT! The equatorial point should have the SAME ENERGY.
1880  // So, we shouldn't use pc[...][ia_90], but should just use pc[...][ia]
1881  // From the other hand, LATER we use this boundary for the ortogonal grid WITHOUD CHANGES,
1882  // so, calculation with pc[...][ia_90] is surprisingly more correct at this point.
1883  // p1 = pc[il][im][alpha.size-1] * sqrt( VF::B(7.)/VF::B(L[il][im][alpha.size-1]));
1884  if (!using_ia_90) {
1885  ia_90=alpha.size-1;
1886  while (alpha.arr[il][im][ia_90] > VC::pi/2 && ia_90 > 0) ia_90--;
1887  }
1888 
1889  p1 = pc.arr[il][im][ia_90] * sqrt( VF::B(7.)/VF::B(L.arr[il][im][ia_90]));
1890 
1891  arr[il][im][ia] = sin(alpha.arr[il][im][ia]) * scale_factor * fL[il] * VF::Outer_Boundary_Flux_at_L7(VF::Kfunc(p1),J_L7_function) * pow(1.0/p1,2);
1892 
1893  // alpha dependence inside the loss cone
1894  if (alpha.arr[il][im][ia] < VF::alc(L.arr[il][im][ia]) || alpha.arr[il][im][ia] > VC::pi - VF::alc(L.arr[il][im][ia])) {
1895  //fLmp[*,i,*]=fLmp[*,i,*]*sinh(i/8.)^2/sinh(1)^2
1896  // ???????? arr[il][im][ia] = arr[il][im][ia] * pow(sinh(alpha[il][im][ia]/8.),2)/pow(sinh(1.0),2);
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);
1898  }
1899  // fLm[*,im]=fL[*]*f1L7(epc[im])*(1/pc[im])^2/fL[index7]
1900  // fLm[iL,im]=fL[iL]
1901  //cout << fLm[il][im] << endl;
1902  }
1903  }
1904  }
1905 
1906 }
1907 
1908 
1909 
1911 
1926 void PSD::Load_initial_f_two_zone(GridElement &L, GridElement &pc, GridElement &alpha, double tau, double tauLpp, double Kp, double min_f, string J_L7_function, double fb_out, double fb_in) {
1927  int il, im, ia;
1928  double Lmax, Lmin, Lgrid, Lpp;
1929  int Lsize_7, addGrid; // These two numbers should be of integer type.
1930  double tau_p, tau_c, alpha1, p1;
1931  Matrix1D<double> L_1d_i, fL_i, tau_0, Ke1, scale_factor(pc.size);
1932  Matrix2D<double> pc_i, fLm;
1933 
1934  addGrid = 0;
1935  Lsize_7 = L.size;
1936 
1938  // Steady state solution for any outer-L boundary
1939  // Create 1D grid at min energy and 90 deg alpha
1940  // First, find a closest index of 90 deg pitch-angle
1941 
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--;
1944 
1945  // If maximum L value is less than 7, we add another points from Lmax(<7) to 7
1946  // and construct pc grid at each added point by moving adiabatically pc values at Lmax(<7)
1947  if (L.arr[L.size - 1][0][ia_90] < 7.0) {
1948 
1949  // construct new L grid
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; // grid points to add
1955  L_1d_i.AllocateMemory(L.size + addGrid); // construct new L grid
1956  fL_i.AllocateMemory(L.size + addGrid); // and psd grid
1957 
1958  // construct new pc grid at each added L point
1959  pc_i.AllocateMemory(L.size + addGrid, pc.size);
1960 
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];
1965  }
1966  }
1967 
1968  // add pc and L values at each added point
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;
1972  pc_i[il][im] = pc.arr[L.size - 1][im][ia_90] * sqrt(VF::B(L_1d_i[il]) / VF::B(L.arr[L.size - 1][im][ia_90]));
1973  }
1974  }
1975  } else {
1976 
1977  // construct new L grid
1978  addGrid = 0; // grid points to add
1979  L_1d_i.AllocateMemory(L.size); // construct new L grid
1980  fL_i.AllocateMemory(L.size); // and psd grid
1981 
1982  // construct new pc grid at each added L point
1983  pc_i.AllocateMemory(L.size, pc.size);
1984 
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];
1989  }
1990  }
1991 
1992  }
1993 
1994  // solve steady state equation based on the newly constructed grid
1995  Lpp = (5.6 - 0.46 * Kp);
1996 
1997  tau_0.AllocateMemory(L.size + addGrid);
1998  Ke1.AllocateMemory(L.size + addGrid);
1999  fLm.AllocateMemory(L.size + addGrid, pc.size);
2000 
2001  for (im = 0; im < pc.size; im++) {
2002  for (il = 0; il < Lsize_7; il++) {
2003 
2004  Ke1[il] = VF::Kfunc(pc_i[il][im]); // Energy in MeV
2005  alpha1 = alpha.arr[L.size - 1][0][alpha.size - 1]; // ~90 deg alpha
2006 
2007  if(PSD_parameters.usetauLpp_SteadyState == "coulomb") { // tau by coulomb scattering
2008  tau_0[il] = VF::Coulomb(L_1d_i[il],Ke1[il]);
2009  } else if(PSD_parameters.usetauLpp_SteadyState == "precipitation") { // tau by precipitaion
2010  tau_0[il] = VF::Precipitation(Kp,L_1d_i[il],Ke1[il]);
2011  } else if (PSD_parameters.usetauLpp_SteadyState == "combined") { // tau by combined scattering
2012  tau_c = VF::Coulomb(L_1d_i[il],Ke1[il]);
2013  tau_p = VF::Precipitation(Kp,L_1d_i[il],Ke1[il]);
2014  tau_0[il] = 1.0 / (1.0/tau_p+1.0/tau_c);
2015  } else { // user specificed tau
2016  tau_0[il] = tauLpp;
2017  }
2018 
2019  }
2020 
2021  L_1d_i.writeToFile("L_1d_i.dat");
2022  Ke1.writeToFile("Ke1.dat");
2023  tau_0.writeToFile("tau_0.dat");
2024 
2025  steady_state_two_zone(fL_i, tau_0, Kp, alpha1, Ke1, Lsize_7, L_1d_i, fb_out, fb_in);
2026 
2027  fL_i.writeToFile("fL_i.dat");
2028 
2029  for (il = 0; il < Lsize_7; il++) {
2030  if (fL_i[il] < min_f) {
2031  fLm[il][im] = min_f;
2032  } else {
2033  fLm[il][im] = fL_i[il];
2034  }
2035  }
2036  }
2037 
2038  // estimate scale factor that would be used to scale down psd at the outer boundary (=Lmax(<7))
2039  for (im = 0; im < pc.size; im++) {
2040  scale_factor[im] = fb_out / fLm[L.size - 1][im];
2041  }
2042 
2043  for (im = 0; im < pc.size; im++) {
2044  for (il = 0; il < L.size; il++) {
2045  for (ia = 0; ia < alpha.size; ia++) {
2046  // filling for all pc
2047  // Why do we do it for alpha=90?:
2048  // Each point is calculated at the equator and multiplied by sin(pitch angle).
2049  // BUT! The equatorial point should have the SAME ENERGY.
2050  // So, we shouldn't use pc[...][ia_90], but should just use pc[...][ia]
2051  // From the other hand, LATER we use this boundary for the ortogonal grid WITHOUD CHANGES,
2052  // so, calculation with pc[...][ia_90] is surprisingly more correct at this point.
2053  // p1 = pc[il][im][alpha.size-1] * sqrt( VF::B(7.)/VF::B(L[il][im][alpha.size-1]));
2054  int ia_90 = alpha.size - 1;
2055  while (alpha.arr[il][im][ia_90] > VC::pi / 2 && ia_90 > 0)
2056  ia_90--;
2057 
2058  p1 = pc.arr[il][im][ia_90] * sqrt(VF::B(7.) / VF::B(
2059  L.arr[il][im][ia_90]));
2060  arr[il][im][ia] = sin(alpha.arr[il][im][ia]) * scale_factor[im]
2061  * fLm[il][im] * VF::Outer_Boundary_Flux_at_L7(
2062  VF::Kfunc(p1), J_L7_function) * pow(1.0 / p1, 2);
2063 
2064  // alpha dependence inside the loss cone
2065  if (alpha.arr[il][im][ia] < VF::alc(L.arr[il][im][ia])
2066  || alpha.arr[il][im][ia] > VC::pi - VF::alc(
2067  L.arr[il][im][ia])) {
2068  arr[il][im][ia] = arr[il][im][ia] * pow(
2069  sinh(alpha.arr[il][im][ia] / VF::alc(
2070  L.arr[il][im][ia])), 2) / pow(sinh(1.0), 2);
2071  }
2072  }
2073  }
2074  }
2075 
2076 #ifdef DEBUG_MODE
2077  //arr.writeToFile("./Debug/f_0.plt", L, pc, alpha);
2078 #endif
2079 
2080 }
2081 
2082 
2083 
2084 
2085 
2099 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) {//, int nit
2100 
2101  //Matrix1D<double> f(nx);
2102  //Matrix1D<double> L(nx);
2103  Matrix1D<double> C_1(nx);
2104  Matrix1D<double> C_2(nx);
2105  Matrix1D<double> C_3(nx);
2106  Matrix1D<double> alfa(nx);
2107  Matrix1D<double> A(nx);
2108  Matrix1D<double> B(nx);
2109  Matrix1D<double> C(nx);
2110  Matrix1D<double> D(nx);
2111  Matrix1D<double> Dxx(nx);
2112 
2113  int i;
2114 
2115  // calculate Dxx
2116  for (i = 0; i <= nx-1; i++) {
2117  //Dxx[i] = VF::Df(L[i], Kp) + VF::Dfe(alpha,Ke[i],L[i],Kp);
2118  Dxx[i] = pow(10, 0.506*Kp-9.325)*pow(L[i],10) + VF::Dfe(alpha,Ke[i],L[i],Kp);
2119  }
2120 
2121  Dxx.writeToFile("Dxx.dat");
2122 
2123  for (i = 1; i < nx-1; i++) {
2124  alfa[i] = 1.0/((L[i+1] - L[i-1])/2)*pow(L[i],2);
2125  //C_1[i] = (Dxx[i] + Dxx[i-1])/(pow(((L[i] + L[i-1])/2), 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));
2128  C_3[i] = 1./tau[i];
2129  //alfa[i] = pow(L[i],2)/(2.*pow((L[i+1] - L[i-1])/2, 2));
2130  }
2131 
2132  // compute values for the diagonals
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];
2137  }
2138 
2139  for (i = 0; i < nx; i++) D[i] = 0;;
2140 
2141  // setting upper and lower boundaries and values
2142  i = 0;
2143  A[i] = 0.;
2144  B[i] = 1;
2145  C[i] = 0;
2146  D[i] = f_bnd_in;
2147 
2148  i=nx-1;
2149  A[i] = 0;
2150  B[i] = 1;
2151  C[i] = 0;
2152  D[i] = f_bnd_out;
2153 
2154 #ifdef DEBUG_MODE
2155  Dxx.writeToFile("./Debug/SS_DLL.plt");
2156  A.writeToFile("./Debug/SS_A.plt");
2157  B.writeToFile("./Debug/SS_B.plt");
2158  C.writeToFile("./Debug/SS_C.plt");
2159  D.writeToFile("./Debug/SS_D.plt");
2160 #endif
2161 
2162  tridag(&A[0], &B[0], &C[0], &D[0], &f[0], nx);
2163 
2164  // for (i = 0; i < nx; i++) f[i] = f[i]/f_bnd_out;
2165 
2166 }
2167 
2168 
2183 void PSD::Load_initial_f_from_outer_L(GridElement &L, GridElement &pc, GridElement &alpha, double tau, double Kp, Matrix2D<double> L_UpperBoundaryCondition, double min_f, double fb_out, double fb_in) {
2184  int il, im, ia;
2185  double p1;
2186  Matrix1D<double> L_1d(L.size), fL(L.size);
2187  Matrix2D<double> fLm(L.size, pc.size);
2188 
2189  // L_UpperBoundaryCondition.writeToFile("./Debug/L_UpperBoundaryCondition.dat");
2190 
2191  // steady state
2192  for (il = 0; il < L.size; il++) L_1d[il] = L.arr[il][0][alpha.size-1]; // create 1D grid at min energy and 90 deg alpha
2193  steady_state(fL, tau, Kp, L.size, L_1d, fb_out, fb_in);
2194  for (int i = 0; i < fL.size_x; i++) if (fL[i] < min_f) fL[i] = min_f;
2195 
2196 #ifdef DEBUG_MODE
2197  fL.writeToFile("./Debug/fL.dat");
2198 #endif
2199 
2200  // Steady state from L-upper boundary
2201  for (im = 0; im < pc.size; im++) {
2202  for (il = 0; il < L.size; il++) {
2203  for (ia = 0; ia < alpha.size; ia++) {
2204  p1 = sqrt( VF::B(7.)/VF::B(L.arr[il][im][ia])); // why 7?
2205  arr[il][im][ia] = sin(alpha.arr[il][im][ia]) * fL[il] * L_UpperBoundaryCondition[im][ia] * pow(1.0/p1,2);
2206  }
2207  }
2208  }
2209 
2210 #ifdef DEBUG_MODE
2211  //arr.writeToFile("./Debug/f_0.plt", L, pc, alpha);
2212 #endif
2213 
2214 }
2215 
2216 
2217 
2227 void PSD::Load_initial_f_2d(GridElement &L, GridElement &pc, GridElement &alpha, const char *filename) {
2228  int i, j, k;
2229  int il, im, ia;
2230  int jmax, lmax;
2231  // number of orbits
2232  lmax=8;
2233  // number of j on orbit
2234  jmax=16;
2235 
2236  Matrix2D<double> jperp(lmax+1, jmax+1);
2237  Matrix1D<double> energy(jmax+2), jcopy(jmax+2);
2238  Matrix1D<double> j_0m(pc.size), f_0m(pc.size);
2239  Matrix1D<double> dayno(lmax+1), eventtime(lmax+1), orbit(lmax+1);
2240 
2241  // Load j from input file
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;
2249  }
2250  }
2251  } else {
2252  throw error_msg("INITIAL_PSD_ERROR", "Error reading initial conditions input file %s.\n", filename);
2253  //cout << "Error reading initial conditions input file " << filename << endl;
2254  //getchar();
2255  //exit(-1);
2256  }
2257  input_f.close();
2258 
2259  // Copying j into 1D matrix_array for interpolation
2260  for (i = 0; i <= jmax; i++) jcopy[i] = jperp[2][i];
2261 
2262  // Adding upper border point and border condition
2263  energy[jmax+1] = VF::Kfunc(pc.GridElement_parameters.max);
2264  jcopy[jmax+1] = VC::zero_f;
2265 
2266 #ifdef DEBUG_MODE
2267  jcopy.writeToFile("./Debug/jcopy.plt", energy);
2268 #endif
2269 
2270  if (VF::Kfunc(pc.GridElement_parameters.min) < energy[0]) printf("Warning: Emin less than minimum data in input file!\n");
2271 
2272  // log all values in jcopy
2273 // for (i = 0; i < jmax+2; i++) {
2274 // jcopy[i] = log10(jcopy[i]);
2275 // }
2276 
2277  // Interpolating(old grid, new grid);
2278  Matrix1D<double> epc1D(pc.size); // temporary 1D energy-grid
2279  j_0m.Spline(jcopy, energy, epc1D, jcopy[0], log10(VC::zero_f), 0, 1); // spline into epc1D
2280  for (i = 0; i < pc.size; i++) {
2281  epc1D[i] = VF::Kfunc(pc.arr[0][i][0]);
2282  // j_0m.Interpolate(jcopy, energy, epc1D);
2283 // j_0m[i] = pow(10, j_0m[i]);
2284 //
2285 // // Calculating from j_0m to f_0m.
2286 // if (j_0m[i] < 0)
2287 // j_0m[i] = VC::zero_f;
2288 // f_0m[i] = j_0m[i]/pow(pc.arr[0][i][0], 2);
2289  }
2290  for (i = 0; i < jmax+2; i++) {
2291  jcopy[i] = log10(jcopy[i]);
2292  }
2293  j_0m.Spline(jcopy, energy, epc1D, jcopy[0], log10(VC::zero_f), 0, 1); // spline into epc1D
2294 
2295 
2296  for (i = 0; i < pc.size; i++) {
2297  j_0m[i] = pow(10, j_0m[i]);
2298  }
2299 
2300  // Calculating from j_0m to f_0m.
2301  for (i = 0; i < pc.size; i++) {
2302  if (j_0m[i] < 0) j_0m[i] = VC::zero_f;
2303  f_0m[i] = j_0m[i]/pow(pc.arr[0][i][0], 2);
2304  }
2305 
2306 #ifdef DEBUG_MODE
2307  j_0m.writeToFile("./Debug/j_0m.plt", epc1D);
2308  f_0m.writeToFile("./Debug/f_0m.plt", epc1D);
2309 #endif
2310 
2312  for (il = 0; il < L.size; il++) {
2313  for (im = 0; im < pc.size; im++) {
2314  for (ia = 0; ia < alpha.size; ia++) {
2315  // copying over interpolated values
2316  arr[il][im][ia] = f_0m[im] * sin( alpha.arr[il][im][ia] );
2317  if (alpha.arr[il][im][ia] < VF::alc(L.arr[il][im][ia]) || alpha.arr[il][im][ia] > VC::pi - VF::alc(L.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);
2319  }
2320  }
2321  }
2322  }
2323 }
2324 
2325 
2326 
2336  int il, im, ia;
2337  double Ke, alf, flux;
2338 
2339  for (il = 0; il < L.size; il++) {
2340  for (im = 0; im < pc.size; im++) {
2341  for (ia = 0; ia < alpha.size; ia++) {
2342 
2343  Ke = VF::Kfunc(pc.arr[0][im][0]); // Energy in MeV
2344  alf = alpha.arr[il][im][ia]; // pitch angle in radian
2345  flux = exp(-(Ke-0.2)/0.1); //* ( sin(alf)-sin(VF::alc(L[il][im][ia])) );
2346  //flux = exp(-(Ke-0.2)/0.1)* ( sin(alf)-sin(VF::alc(L[il][im][ia])) ); */
2347  arr[il][im][ia] = flux / pow(pc.arr[0][im][0], 2);
2348 
2349  }
2350  }
2351  }
2352 }
2353 
2354 
2355 
2367 void steady_state(Matrix1D<double> &f, double tau, double Kp, int nx, Matrix1D<double> &L, double f_bnd_out, double f_bnd_in) {//, int nit
2368 
2369  //Matrix1D<double> f(nx);
2370  //Matrix1D<double> L(nx);
2371  Matrix1D<double> C_1(nx);
2372  Matrix1D<double> C_2(nx);
2373  Matrix1D<double> C_3(nx);
2374  Matrix1D<double> alfa(nx);
2375  Matrix1D<double> A(nx);
2376  Matrix1D<double> B(nx);
2377  Matrix1D<double> C(nx);
2378  Matrix1D<double> D(nx);
2379  Matrix1D<double> Dxx(nx);
2380 
2381  int i;
2382 
2383  for (i = 0; i <= nx-1; i++) {
2384  // Radial diffusion coefficient approximation by Brautigam and Albert [2000]
2385  Dxx[i] = pow(10, 0.506*Kp-9.325)*pow(L[i],10);
2386  }
2387 
2388  for (i = 1; i < nx-1; i++) {
2389  alfa[i] = 1.0/((L[i+1] - L[i-1])/2)*pow(L[i],2);
2390  //C_1[i] = (Dxx[i] + Dxx[i-1])/(pow(((L[i] + L[i-1])/2), 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));
2393  C_3[i] = 1./tau;
2394  //alfa[i] = pow(L[i],2)/(2.*pow((L[i+1] - L[i-1])/2, 2));
2395  }
2396 
2397  // compute values for the diagonals
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];
2402  }
2403 
2404  for (i = 0; i < nx; i++) D[i] = 0;;
2405 
2406  // setting upper and lower boundary values
2407  i = 0;
2408  A[i] = 0.;
2409  B[i] = 1;
2410  C[i] = 0;
2411  D[i] = f_bnd_in;
2412 
2413  i=nx-1;
2414  A[i] = 0;
2415  B[i] = 1;
2416  C[i] = 0;
2417  D[i] = f_bnd_out;
2418 
2419 #ifdef DEBUG_MODE
2420  Dxx.writeToFile("./Debug/SS_DLL.plt");
2421  A.writeToFile("./Debug/SS_A.plt");
2422  B.writeToFile("./Debug/SS_B.plt");
2423  C.writeToFile("./Debug/SS_C.plt");
2424  D.writeToFile("./Debug/SS_D.plt");
2425 #endif
2426 
2427  tridag(&A[0], &B[0], &C[0], &D[0], &f[0], nx);
2428 
2429  // for (i = 0; i < nx; i++) f[i] = f[i]/f_bnd_out;
2430 
2431 }
bool initialized
Flag, equal true if initialized.
Definition: Matrix.h:138
GridElement alpha
grid elements to be used
Definition: Grid.h:59
bool useLmp
magnetopause position flag
Definition: Parameters.h:246
Matrix3D< double > arr
array of grid points
Definition: Grid.h:30
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)
Definition: Matrix.cpp:1455
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)
Definition: PSD.cpp:895
void Interpolate(PSD &otherPSD, Parameters_structure::Interpolation interpolationParameters, Grid &oldGrid, Grid &newGrid, Matrix2D< double > newGrid_pc_lowerBoundaryCondition, Matrix2D< double > newGrid_pc_upperBoundaryCondition)
Interpolation function.
Definition: PSD.cpp:1557
Interpolation parameters structure
Definition: Parameters.h:253
void Polilinear(Matrix1D< T > &old_function, Matrix1D< T > &old_grid, Matrix1D< T > &new_grid, double lb, double ub)
Definition: Matrix.cpp:1673
void SourcesAndLosses(Parameters_structure parameters, GridElement &L, GridElement &pc, GridElement &alpha, AdditionalSourcesAndLosses &SL, double dt, double Lpp, double tau, double tauLpp, double Kp)
Definition: PSD.cpp:319
double initial_PSD_inner_psd
Inner PSD boundary value for steady state, default 0.
Definition: Parameters.h:207
Array of values of coordinate axes.
Definition: Grid.h:28
void Load_initial_f_file(GridElement &L, GridElement &pc, GridElement &alpha, const char *filename, bool with_grid)
Definition: PSD.cpp:1742
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.
Definition: Parameters.h:212
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)
Definition: PSD.cpp:2367
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)
Definition: PSD.cpp:2099
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)
Definition: PSD.cpp:599
Holds diffusion coefficient matrix and routines to load and calculate it.
string usetauLpp
key word to define if and how to use tau/tauLpp
Definition: Parameters.h:120
double initial_PSD_Kp0
Kp for steady state?
Definition: Parameters.h:203
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)
Definition: PSD.cpp:1053
void AllocateMemory(int size_x, int size_y)
Definition: Matrix.cpp:524
void Create_Initial_PSD(Parameters_structure::PSD parameters, Grid &grid, BoundaryCondition L_UpperBoundaryCondition)
Create inital PSD (Steady State)
Definition: PSD.cpp:47
void Output_without_grid(double time)
Definition: PSD.cpp:1720
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.
Definition: Parameters.h:195
string type
grid type
Definition: Grid.h:65
Parameters_structure::GridElement GridElement_parameters
parameters for grid element
Definition: Grid.h:36
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...
Definition: Parameters.h:193
void echo(int msglvl, const char *format,...)
Basic function! Write the message to the file.
Definition: Output.cpp:44
void AllocateMemory(int x_size)
Definition: Matrix.cpp:173
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)
Definition: PSD.cpp:1926
string type
Interpolation type: spline, linear etc.
Definition: Parameters.h:255
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)
Definition: PSD.cpp:2335
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...
Definition: PSD.cpp:2183
static const double pi
π
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)
Definition: PSD.cpp:414
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)
Definition: PSD.cpp:747
Main parameters structure that holds smaller structures for individual parameters.
Definition: Parameters.h:94
double Coulomb(double L, double energy)
Coulomb scattering electron lifetime.
int size
size of grid element
Definition: Grid.h:39
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)
Definition: PSD.cpp:1790
double initial_PSD_some_constant_value
Some psd value. Used as initial PSD value everywhere if initial_PSD_Type = IPSDT_CONSTANT or as minim...
Definition: Parameters.h:204
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.
Definition: Parameters.h:199
double B(double Lparam)
Dipole magnetic field.
void writeToFile(string filename)
Definition: Matrix.cpp:366
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)
Definition: PSD.cpp:1269
Matrix3D< double > arr
array of PSD values
Definition: PSD.h:53
void Interpolate(Matrix1D< T > &old_function, Matrix1D< T > &old_grid, Matrix1D< T > &new_grid)
Definition: Matrix.cpp:1414
Phase Space Density class.
Definition: PSD.h:48
GridElement L
grid elements to be used
Definition: Grid.h:59
double alc(double L)
Loss cone calculations.
void checkInf_3D(Matrix3D< double > arr, int il, int im, int ia, double maxNum=1e99)
Definition: PSD.cpp:1532
Holds upper and lower boundary conditions.
string initial_PSD_fileName
File name, if we need to load initial values from file.
Definition: Parameters.h:194
string output_PSD_folderName
Folder name for psd output.
Definition: Parameters.h:211
Parameters_structure parameters
Parameters structure, with all parameters from the parameters.ini file. The default parameters define...
Definition: Main.cpp:187
string useLog
Flag, if we should interpolate logarithm of the values.
Definition: Parameters.h:257
GridElement epc
grid elements to be used
Definition: Grid.h:59
Error message - stack of single_errors.
Definition: error.h:54
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'.
Definition: Parameters.h:205
double initial_PSD_outer_psd
Outer PSD boundary value for steady state, default 1.
Definition: Parameters.h:206
Computational grid composed of 3 different GridElement.
Definition: Grid.h:53
PSD parameters structure.
Definition: Parameters.h:190
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)
Definition: PSD.cpp:1548
static const double zero_f
Minimum value for PSD.
int size_x
size x
Definition: Matrix.h:46
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)
Definition: PSD.cpp:2227
GridElement pc
grid elements to be used
Definition: Grid.h:59
double G(double alpha)
G-function.