VERB_code_2.2  2
 All Classes Namespaces Files Functions Variables Typedefs Macros Pages
PSD.cpp
Go to the documentation of this file.
1 /**
2 * \file PSD.cpp
3 *
4 * \brief Makes operations with %PSD (Phase Space Density (%PSD).) (like, diffusion).
5 *
6 * General view of Fokker-Planck diffusion equation in terms of linearization:
7 * \f[ A \times x = B \times x + C \f]
8 * A, B, C refered as MatrixA, MatrixB and MatrixC in the code
9 *
10 * \todo
11 * - A lot of corrections should be done in PSD class to make it more logical and less spread.
12 *
13 * \author Developed under supervision of the PI Yuri Shprits
14 */
15 
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 #include "../Exceptions/error.h"
29 
30 // General c++ headers
31 #include <iostream>
32 #include <string>
33 #include <ctime>
34 
35 #include <stdio.h>
36 #include <stdlib.h>
37 
38 using namespace std;
39 
40 
41 /**
42 * Initializing: storing parameters, loading initial values, making boundary conditions, initializing output parameters
43 * Simply, it is a creation of the object.
44 */
46 {
47 
48  if (!arr.initialized) {
49  arr.AllocateMemory(grid.L.size, grid.pc.size, grid.alpha.size);
50  }
51 
52  // Copying from parameters to the class variables
53  PSD_parameters = parameters;
54 
55  // Loading initial values, from time-dependent boundary conditions takes only first time step values
56  //LoadInitialValue(parameters, grid, L_UpperBoundaryCondition.xSlice(0));
57  // iterators
58  int il, im, ia;
59 
60  // loading (calculating) initial values
61  if (parameters.initial_PSD_Type == "IPSDT_ORBIT_FLUX_2D") {
62  // loading from orbit flux (for 2d)
63  if (grid.L.size != 1) {
64  throw error_msg("INITIAL_PSD_2D_METHOD_FOR_3D_GRID", "Initial PSD error: imposibble to use ORBIT_FLUX_2D if L.size > 1");
65  }
66  Load_initial_f_2d(grid.L, grid.pc, grid.alpha, parameters.initial_PSD_fileName.c_str());
67  ////////////////////////// kckim (mawellian) /////////////////////////
68  }
69  else if (parameters.initial_PSD_Type == "IPSDT_FLUX_2D_MAXWELLIAN") {
70  if (grid.L.size != 1) {
71  throw error_msg("INITIAL_PSD_2D_METHOD_FOR_3D_GRID", "Initial PSD error: imposibble to use ORBIT_FLUX_2D if L.size > 1");
72  }
73  cout<<"kim................................"<<endl;
74  Load_initial_f_maxwell(grid.L, grid.pc, grid.alpha);
75 
76  //////////////////////////////////////////////////////////////////////
77  }
78  else if (parameters.initial_PSD_Type == "IPSDT_STEADY_STATE") {
79  // calculating as a steady state (for 3d)
80  if (grid.L.size == 1) {
81  throw error_msg("INITIAL_PSD_3D_METHOD_FOR_2D_GRID", "Initial PSD error: imposibble to use STEADY_STATE if L.size == 1");
82  }
83  // if tauSteadyState == 0 - using 6.0/Kp
84  /// moved to parameters.cpp if (parameters.initial_PSD_tauSteadyState <= 1.e-99) parameters.initial_PSD_tauSteadyState = 4.0/parameters.initial_PSD_Kp0;
85  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);
86  /////////////////////// kckim /////////////////////////////////////////////////
87  }
88  else if (parameters.initial_PSD_Type == "IPSDT_STEADY_STATE_TWO_ZONE") {
89  // calculating as a steady state (for 3d)
90  if (grid.L.size == 1) {
91  throw error_msg("INITIAL_PSD_3D_METHOD_FOR_2D_GRID", "Initial PSD error: imposibble to use STEADY_STATE if L.size == 1");
92  }
93  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);
94  //////////////////////////////////////////////////////////////////////////////////
95  }
96  else if (parameters.initial_PSD_Type == "IPSDT_STEADY_STATE_FROM_BOUNDARY") {
97  // calculating as a steady state (for 3d)
98  if (grid.L.size == 1) {
99  throw error_msg("INITIAL_PSD_3D_METHOD_FOR_2D_GRID", "Initial PSD error: imposibble to use STEADY_STATE if L.size == 1");
100  }
101  if (!L_UpperBoundaryCondition.arr.initialized) {
102  throw error_msg("INITIAL_PSD_ERROR", "Initial PSD error: upper boundary is not initialized, can not be used for steady state calculation.");
103  }
104  // if tauSteadyState == 0 - using 6.0/Kp
105  /// moved to parameters.cpp if (parameters.initial_PSD_tauSteadyState <= 1.e-99) parameters.initial_PSD_tauSteadyState = 4.0/parameters.initial_PSD_Kp0;
106  Load_initial_f_from_outer_L(grid.L, grid.pc, grid.alpha,
107  parameters.initial_PSD_tauSteadyState, parameters.initial_PSD_Kp0,
108  L_UpperBoundaryCondition.arr, parameters.initial_PSD_some_constant_value,
109  parameters.initial_PSD_outer_psd, parameters.initial_PSD_inner_psd);
110  }
111  else if (parameters.initial_PSD_Type == "IPSDT_INTERPOLATION") {
112  // Initial values for local grid as an interpolation from radial grid, or opposite
113  // throw error_msg("INITIAL_PSD_TYPE_ERROR", "Using interpolation initial PSD type without second PSD to interpolate from");
114  }
115  else if (parameters.initial_PSD_Type == "IPSDT_CONSTANT") {
116  // just constant values
117  for (il = 0; il < grid.L.size; il++)
118  for (im = 0; im < grid.pc.size; im++)
119  for (ia = 0; ia < grid.alpha.size; ia++)
120  arr[il][im][ia] = parameters.initial_PSD_some_constant_value;
121  }
122  else if (parameters.initial_PSD_Type == "IPSDT_FILE") {
123  // loading from orbit flux (for 2d)
124  Load_initial_f_file(grid.L, grid.epc, grid.alpha, parameters.initial_PSD_fileName.c_str(), false);
125  }
126  else if (parameters.initial_PSD_Type == "IPSDT_FILE_GRID") {
127  // loading from orbit flux (for 2d)
128  Load_initial_f_file(grid.L, grid.epc, grid.alpha, parameters.initial_PSD_fileName.c_str(), true);
129  }
130  else {
131  throw error_msg("INITIAL_PSD_TYPE_UNKNOWN", "Unknown initial PSD type");
132  }
133 
134  // Empting loss cone
135  if (parameters.initial_PSD_Type == "IPSDT_ORBIT_FLUX_2D" ||
136  parameters.initial_PSD_Type == "IPSDT_STEADY_STATE" ||
137  parameters.initial_PSD_Type == "IPSDT_STEADY_STATE_FROM_BOUNDARY") {
138  for (il = 0; il < grid.L.size; il++)
139  for (im = 0; im < grid.pc.size; im++)
140  for (ia = 0; ia < grid.alpha.size; ia++) {
141  if (grid.alpha.arr[il][im][ia] < VF::alc(grid.L.arr[il][im][ia]) || grid.alpha.arr[il][im][ia] > VC::pi - VF::alc(grid.L.arr[il][im][ia]))
142  // arr[il][im][ia] = VC::zero_f;
143  // alpha dependence inside the loss cone
144  if (grid.alpha.arr[il][im][ia] < VF::alc(grid.L.arr[il][im][ia]) || grid.alpha.arr[il][im][ia] > VC::pi - VF::alc(grid.L.arr[il][im][ia])) {
145  //fLmp[*,i,*]=fLmp[*,i,*]*sinh(i/8.)^2/sinh(1)^2
146  // ???????? arr[il][im][ia] = arr[il][im][ia] * pow(sinh(alpha[il][im][ia]/8.),2)/pow(sinh(1.0),2);
147  arr[il][im][ia] = arr[il][im][ia] * pow(sinh(grid.alpha.arr[il][im][ia]/VF::alc(grid.L.arr[il][im][ia])),2)/pow(sinh(1.0),2);
148  }
149  }
150 
151  }
152 
153 
154 
155  // Output - create stream
156  output_without_grid_file.open((parameters.output_PSD_folderName + parameters.output_PSD_fileName4D).c_str());
157  output_without_grid_file << "VARIABLES = \"Phase_Space_Density\" " << endl;
158  output_without_grid_file.setf(ios_base::scientific, ios_base::floatfield);
159 
160  // initialized = true;
161 }
162 
163 
164 /**
165 * Scaling initial PSD if changed upper L-boundary (not 7Re)
166 */
167 /*
168 void PSD::ScaleBoundaryFlux(Parameters_structure::PSD parameters, GridElement &L, GridElement &pc, GridElement &alpha)
169 {
170 
171  Matrix1D<double> L_1d(61),fL(61),tau_0(61),Ke1(61);
172  Matrix2D<double> scale_fm(61,pc.size);
173  int il,im,ia;
174  double L_upper, scale_f, alpha1,Lpp,tau_c,tau_p;
175  double Kp = parameters.initial_PSD_Kp0;
176  double tau = parameters.initial_PSD_tauSteadyState;
177  double tauLpp = parameters.initial_PSD_tauLppSteadyState;
178  double fb_in,fb_out,Lmin,Lmax,Lsize,dL;
179 
180 
181  if (parameters.initial_PSD_Type == "IPSDT_STEADY_STATE") { // || parameters.initial_PSD_Type == "IPSDT_STEADY_STATE_TWO_ZONE") {
182  // solve 1d steady-state eq. from L=1 to L=7
183 
184  fb_in = 0;
185  fb_out = 1;
186  Lmin = 1.0;
187  Lmax = 7.0;
188  dL = 0.1;
189  Lsize = (Lmax-Lmin)/dL+1;
190 
191  for (il=0; il<Lsize; il++) L_1d[il] = 1+il*dL; // make 1d grid
192  steady_state(fL, tau, Kp, Lsize, L_1d, fb_out, fb_in); // solve 1d steady state eq. for L=1 to 7
193 
194  L_upper = L[L.size-1][0][0]; // L-upper boundary we want to change, which is given in parameters.ini file
195 
196  if(L_upper != Lmax) { // if L is not equal to 7
197  unsigned i=0; // interpolate PSD using Lagragian
198  while (L_1d[++i] <= L_upper);
199  scale_f = fL[i-1] * ( (L_upper-L_1d[i])/(L_1d[i-1]-L_1d[i]) )
200  + fL[i] * ( (L_upper-L_1d[i-1])/(L_1d[i]-L_1d[i-1]) );
201  } else { // if L is equal to 7
202  scale_f = fb_out; // scale_f = 1;
203  }
204 
205 
206  for (il = 0; il < L.size; il++) {
207  for (im = 0; im < pc.size; im++) {
208  for (ia = 0; ia < alpha.size; ia++) {
209  arr[il][im][ia] = arr[il][im][ia] * scale_f;
210 
211  }
212  }
213  }
214 
215 
216 
217  }
218 
219 
220 ///
221 
222  if (parameters.initial_PSD_Type == "IPSDT_STEADY_STATE_TWO_ZONE") {
223 
224  fb_in = 0;
225  fb_out = 1;
226  Lmin = 1.0;
227  Lmax = 7.0;
228  dL = 0.1;
229  Lsize = (Lmax-Lmin)/dL+1;
230  //Lpp = (5.6 - 0.46*Kp);
231 
232  for (il=0; il<Lsize; il++) L_1d[il] = 1+il*dL; // make 1d grid
233 
234  for (im = 0; im < pc.size; im++){
235  for (il = 0; il < L.size; il++) {
236 
237  Ke1[il] = VF::Kfunc(pc[il][im][alpha.size-1]); // Energy in MeV
238  alpha1 = alpha[il][im][alpha.size-1]; // ~90 deg alpha
239  // cout<<Ke1[il]<<" "<<il<<Lsize<<endl;
240  // make a array of tau
241  //if(L_1d[il] >= Lpp) { // calculate tau
242  // tau_0[il] = tau;
243  //} else { // calculate tauLpp
244  if(tauLpp == 1e99) {
245  tau_0[il] = VF::Coulomb(L_1d[il],Ke1[il]);
246  } else if (tauLpp == 2e99) {
247  tau_0[il] = VF::Precipitation(Kp,L_1d[il],Ke1[il]);
248  } else if (tauLpp == 3e99) {
249  tau_c = VF::Coulomb(L_1d[il],Ke1[il]);
250  tau_p = VF::Precipitation(Kp,L_1d[il],Ke1[il]);
251  tau_0[il] = 1.0 / (1.0/tau_p+1.0/tau_c);
252  } else {
253  tau_0[il] = tauLpp;
254  }
255  //cout<<tau_0[il]<<" "<<il<<endl;
256  //}
257  }
258  steady_state_two_zone(fL, tau_0, Kp, alpha1, Ke1, L.size, L_1d, fb_out, fb_in);
259 
260  L_upper = L[L.size-1][0][0]; // L-upper boundary we want to change, which is given in parameters.ini file
261 
262  if(L_upper != Lmax) { // if L is not equal to 7
263  unsigned i=0; // interpolate PSD using Lagragian
264  while (L_1d[++i] <= L_upper);
265  scale_fm[il][im] = fL[i-1] * ( (L_upper-L_1d[i])/(L_1d[i-1]-L_1d[i]) )
266  + fL[i] * ( (L_upper-L_1d[i-1])/(L_1d[i]-L_1d[i-1]) );
267  } else { // if L is equal to 7
268  scale_fm[il][im] = fb_out; // scale_f = 1;
269  }
270  }
271  // scaling PSD ///////////////////////
272  for (il = 0; il < L.size; il++) {
273  for (im = 0; im < pc.size; im++) {
274  for (ia = 0; ia < alpha.size; ia++) {
275  arr[il][im][ia] = arr[il][im][ia] * scale_fm[il][im];
276  }
277  }
278  }
279  //////////////////////////////////////
280  }
281 
282 
283 }
284 */
285 
286 
287 /**
288 * Sources and losses term from the FP eq.
289 * \param &L - L-grid array
290 * \param &pc - pc-grid array
291 * \param &alpha - alpha-grid array
292 * \param &SL - sources/losses 3D matrix
293 * \param dt - time sted
294 * \param Lpp - plasmapause location
295 * \param tau - lifetime outside of the plasmapause
296 * \param tauLpp - lifetime inside of the plasmapause
297 *
298 */
301  GridElement &L, GridElement &pc, GridElement &alpha,
302  Matrix3D<double> &SL,
303  double dt, double Lpp,
304  double tau, double tauLpp, double Kp)
305 {
306  int il, im, ia;
307  double taulc;
308  ///////////////// kckim /////////////////
309  double tau_p,tau_c,tauLpp0;
310  /////////////////////////////////////////
311 
312  for (il = 0; il < L.size; il++) {
313  for (im = 0; im < pc.size; im++) {
314  for (ia = 0; ia < alpha.size; ia++) {
315 
316  // Loss cone. - does not work here, needs too small time step to be splitted from pitch-angle diffusion
317  // df/dt = -f/taulc
318 
319 
320 /*
321  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])))) ) {
322  // taulc = VF::bounce_time_res(VF::Kfunc(pc.arr[il][im][ia]), L.arr[il][im][ia])/60./60./24.;
323  // Quarter bounce time
324  taulc = 0.25 * VF::bounce_time_new(L.arr[il][im][ia], pc.arr[il][im][ia], alpha.arr[il][im][ia]);
325  arr[il][im][ia] = arr[il][im][ia] / dt / (1/dt + 1/taulc); // implicit
326  // arr[il][im][ia] = arr[il][im][ia] * exp(-dt/taulc); // some other
327  // arr[il][im][ia] = exp(log(arr[il][im][ia]) - dt/taulc); // to not have PSD=0
328  }
329 */
330  // Additional losses, if any
331  // df/dt = -f/tau
332  if (L.arr[il][im][ia] >= Lpp) {
333  // arr[il][im][ia] = arr[il][im][ia] / dt / (1/dt + 1/tau); // implicit
334  arr[il][im][ia] = arr[il][im][ia] - arr[il][im][ia] * (dt / tau); // explicit
335 
336  } else {
337  // (*this)[il][im][ia] = (*this)[il][im][ia] / dt / (1/dt + 1/tauLpp); // implicit
338  if(parameters.usetauLpp == "coulomb") {
339  tauLpp0 = VF::Coulomb(L.arr[il][im][ia], VF::Kfunc(pc.arr[il][im][ia]));
340  } else if (parameters.usetauLpp == "precipitaion") {
341  tauLpp0 = VF::Precipitation(Kp, L.arr[il][im][ia], VF::Kfunc(pc.arr[il][im][ia]));
342  } else if (parameters.usetauLpp == "combined") {
343  tau_p = VF::Precipitation(Kp, L.arr[il][im][ia], VF::Kfunc(pc.arr[il][im][ia]));
344  tau_c = VF::Coulomb(L.arr[il][im][ia], VF::Kfunc(pc.arr[il][im][ia]));
345  tauLpp0 = 1.0/(1.0/tau_p+1.0/tau_c);
346  } else {
347  tauLpp0 = tauLpp;
348  }
349  arr[il][im][ia] = arr[il][im][ia] - arr[il][im][ia] * (dt / tauLpp0); // explicit, probably should be implicit! I.e.:
350  //arr[il][im][ia] = arr[il][im][ia] / dt / (1/dt + 1/tauLpp0); // <-- implicit
351  }
352  }
353  }
354  }
355 
356 }
357 
358 
359 /**
360 * Mixed terms calculation by explicit method.
361 *
362 * \param dt - time step
363 * \param Lpp - plasma pause location
364 * \param &Dpca - pc-alpha diffusion coefficient
365 * \param &DpcaLpp - pc-alpha diffusion coefficient under plasma pause location
366 * \param &L - grid element L
367 * \param &pc - grid element pc
368 * \param &alpha - grid element alpha
369 * \param &Jacobian - jacobian
370 * \param &pc_lowerBoundaryCondition - lower boundary condition on pc
371 * \param &pc_upperBoundaryCondition - upper boundary condition on pc
372 * \param &alpha_lowerBoundaryCondition - lower boundary condition on alpha
373 * \param &alpha_upperBoundaryCondition - upper boundary condition on alpha
374 * \param pc_lowerBoundaryCondition_calculationType - pc lower boundary condition type,
375 * \param pc_upperBoundaryCondition_calculationType - pc upper boundary condition type,
376 * \param alpha_lowerBoundaryCondition_calculationType - alpha lower boundary condition type,
377 * \param alpha_upperBoundaryCondition_calculationType - alpha upper boundary condition type.
378 */
379 void PSD::DiffusionMixTermExplicit(double dt, double Lpp,
381  GridElement &L, GridElement &pc, GridElement &alpha, Matrix3D<double> Jacobian,
382  Matrix2D<double> pc_lowerBoundaryCondition,
383  Matrix2D<double> pc_upperBoundaryCondition,
384  Matrix2D<double> alpha_lowerBoundaryCondition,
385  Matrix2D<double> alpha_upperBoundaryCondition,
386  string pc_lowerBoundaryCondition_calculationType,
387  string pc_upperBoundaryCondition_calculationType,
388  string alpha_lowerBoundaryCondition_calculationType,
389  string alpha_upperBoundaryCondition_calculationType)
390 {
391 
392  int il, im, ia;
393 
394  // Declaration of arrays for derivatives
395  static Matrix3D<double> PSD_Der_a_L(L.size, pc.size, alpha.size);
396  static Matrix3D<double> PSD_Der_a_R(L.size, pc.size, alpha.size);
397  static Matrix3D<double> PSD_Der_m_L(L.size, pc.size, alpha.size);
398  static Matrix3D<double> PSD_Der_m_R(L.size, pc.size, alpha.size);
399 
400  // Some variables
401  double Dpca_, Dpca_a_L, Dpca_a_R, Dpca_m_L, Dpca_m_R;
402  double J_, J_a_L, J_a_R, J_m_L, J_m_R;
403  double Laml, Lmal, Lamr, Lmar;
404  double PSD_m_r, PSD_m_l, PSD_a_r, PSD_a_l, Lam, Lma;
405 
406  // Different way of derivatives approximation
407  if (PSD_parameters.approximationMethod == "AM_Split_LR") {
408 
409  // Derivatives calculation
410  for (il = 0; il < L.size; il++) {
411  for (im = 0; im < pc.size; im++) {
412  for (ia = 1; ia < alpha.size-1; ia++) {
413  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]);
414  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]);
415  }
416  }
417  for (im = 1; im < pc.size-1; im++) {
418  for (ia = 0; ia < alpha.size; ia++) {
419  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]);
420  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]);
421  }
422  }
423  }
424 
425  // Approximation and calculation
426  for (il = 0; il < L.size; il++) {
427  for (im = 1; im < pc.size-1; im++) {
428  for (ia = 1; ia < alpha.size-1; ia++) {
429  // for debugging...
430  Dpca_ = Dpca.arr[il][im][ia];
431  Dpca_a_L = Dpca.arr[il][im][ia-1];
432  Dpca_a_R = Dpca.arr[il][im][ia+1];
433  Dpca_m_L = Dpca.arr[il][im-1][ia];
434  Dpca_m_R = Dpca.arr[il][im+1][ia];
435  J_ = Jacobian[il][im][ia];
436  J_a_L = Jacobian[il][im][ia-1];
437  J_a_R = Jacobian[il][im][ia+1];
438  J_m_L = Jacobian[il][im-1][ia];
439  J_m_R = Jacobian[il][im+1][ia];
440 
441 
442  // d/dp D d/da f
443  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]);
444  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]);
445 
446  // d/da D d/dp f
447  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]);
448  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]);
449 
450  arr[il][im][ia] = arr[il][im][ia] + dt * 0.5 * (Lmar + Lamr + Lmal + Laml);
451 
452  /* L12p = 1/J*(Jacobian[il][im+1][ia]*Dpca[il][im+1][ia]*PSD_Der_a_L[il][im+1][ia] -
453  Jacobian[il][im ][ia]*Dpca[il][im ][ia]*PSD_Der_a_L[il][im ][ia])/(pc[il][im+1][ia] - pc[il][im ][ia]);
454  L12m = 1/J*(Jacobian[il][im ][ia]*Dpca[il][im ][ia]*PSD_Der_a_R[il][im ][ia] -
455  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]);
456 
457  L21p = 1/J*(Jacobian[il][im][ia+1]*Dpca[il][im][ia+1]*PSD_Der_m_L[il][im][ia+1] -
458  Jacobian[il][im][ia ]*Dpca[il][im][ia ]*PSD_Der_m_L[il][im][ia ])/(alpha[il][im][ia+1] - alpha[il][im][ia]);
459  L21m = 1/J*(Jacobian[il][im][ia ]*Dpca[il][im][ia ]*PSD_Der_m_R[il][im][ia ] -
460  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]);*/
461  }
462  }
463  }
464  } else if (PSD_parameters.approximationMethod == "AM_Split_C") { // was not fully checked, may not work
465  for (il = 0; il < L.size; il++) {
466  for (im = 1; im < pc.size-1; im++) {
467  for (ia = 1; ia < alpha.size-1; ia++) {
468  // for debugging...
469  Dpca_ = Dpca.arr[il][im][ia];
470  Dpca_a_L = Dpca.arr[il][im][ia-1];
471  Dpca_a_R = Dpca.arr[il][im][ia+1];
472  Dpca_m_L = Dpca.arr[il][im-1][ia];
473  Dpca_m_R = Dpca.arr[il][im+1][ia];
474  J_ = Jacobian[il][im][ia];
475  J_a_L = Jacobian[il][im][ia-1];
476  J_a_R = Jacobian[il][im][ia+1];
477  J_m_L = Jacobian[il][im-1][ia];
478  J_m_R = Jacobian[il][im+1][ia];
479 
480  // d/dp D d/da f
481  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]);
482  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]);
483  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]);
484 
485  // d/da D d/dm f
486  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]);
487  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]);
488  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]);
489 
490  // cout << ((Lma + Lam) - (Lmar + Lamr + Lmal + Laml))/((Lma + Lam) + (Lmar + Lamr + Lmal + Laml)) << "%" << endl;
491  arr[il][im][ia] = arr[il][im][ia] + dt * (Lma + Lam);
492  }
493  }
494  }
495 
496  } else {
497  throw error_msg("APPROXIMATION", "Unknown approximation method for mixed terms");
498  }
499 
500 }
501 
502 
503 /** Pitch angle diffusion calculation function.
504 * Takes diffusion coefficients, boundary conditions and grid elements as input.
505 * It calculates model matrix by MakeModelMatrix_3D for 3D problem: A_3D * PSD_3D(t+1) = B_3D * PSD_3D(t) + C_3D
506 * then splits it for 1D problems: A * PSD_1D(t+1) = B * PSD_1D(t) + C
507 * and solve by tridiag method.
508 *
509 * \param dt - time step
510 * \param Lpp - plasma pause location
511 * \param &Daa - L diffusion coefficient
512 * \param &DaaLpp - L diffusion coefficient
513 * \param &L - grid element L
514 * \param &pc - grid elelent pc
515 * \param &alpha - grid element alpha
516 * \param &Jacobian - jacobian
517 * \param &alpha_lowerBoundaryCondition - lower boundary condition
518 * \param &alpha_upperBoundaryCondition - upper boundary condition
519 * \param &alpha_lowerBoundaryCondition_calculationType - lower boundary condition type (on value/on derivative)
520 * \param &alpha_upperBoundaryCondition_calculationType - upper boundary condition type (on value/on derivative)
521 */
522 void PSD::Diffusion_alpha(double dt, double Lpp,
524  GridElement &L, GridElement &pc, GridElement &alpha, Matrix3D<double> Jacobian,
525  Matrix2D<double> alpha_lowerBoundaryCondition,
526  Matrix2D<double> alpha_upperBoundaryCondition,
527  string alpha_lowerBoundaryCondition_calculationType,
528  string alpha_upperBoundaryCondition_calculationType)
529 {
530 
531  char out_name_matr_a[256], out_name_matr_b[256], out_name_matr_c[256];
532 
533  // Calculation matrices A, B, and C for 1d diffusion problem
534  // A * PSD_1D(t+1) = B * PSD_1D(t) + C, for each L and Energy
535 
536  // Creation is moved to PSD.h, so it can be accessed from Main.cpp
537  //static CalculationMatrix
538  // matr_A(alpha.size, 1, 1, 1),
539  // matr_B(alpha.size, 1, 1, 0),
540  // matr_C(alpha.size, 1, 1, 0);
541  matr_A.Initialize(alpha.size, 1, 1, 1);
542  matr_B.Initialize(alpha.size, 1, 1, 0),
543  matr_C.Initialize(alpha.size, 1, 1, 0);
544 
545 
546  // Create zero matrices, temporary matrices to be used later
547  static Matrix3D<double> Zero_Matrix(1, 1, alpha.size);
548  static Matrix2D<double> Zero_Matrix_2d(1, 1);
549  static Matrix1D<double> Zero_Matrix_1d(alpha.size);
550  Zero_Matrix = 0;
551  Zero_Matrix_2d = 0;
552  Zero_Matrix_1d = 0;
553 
554  // Create 1d arrays
555  // Creating 1d arrays for PSD and right hand side (RHS = B * PSD + C)
556  static Matrix1D<double> slicePSD1D(alpha.size), RHS(alpha.size);
557 
558  // Creating 1d arrays which will be send to the MakeMatrix function
559  // Creating 1d matrixes, but with type of Matrix3D, cause we need that type for MakeMatrix function
560  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);
561  // Creating temporary arrays size 1x1 for boundary conditions, which are just values
562  static Matrix2D<double> alpha_lowerBoundaryCondition_1d(1, 1), alpha_upperBoundaryCondition_1d(1, 1);
563 
564  // iterators
565  int il, im, ia, in; // matrix element numbers on 3D grid
566  DiagMatrix::iterator it;
567  for (il = 0; il < L.size; il++) {
568  for (im = 0; im < pc.size; im++) {
569 
570  // Copying 3D arrays into 1d arrays
571  for (ia = 0; ia < alpha.size; ia++) {
572  slicePSD1D[ia] = arr[il][im][ia];
573 
574  Daa_1d[0][0][ia] = Daa.arr[il][im][ia];
575  DaaLpp_1d[0][0][ia] = DaaLpp.arr[il][im][ia];
576 
577  L_1d[0][0][ia] = L.arr[il][im][ia];
578  pc_1d[0][0][ia] = pc.arr[il][im][ia];
579  alpha_1d[0][0][ia] = alpha.arr[il][im][ia];
580 
581  Jacobian_1d[0][0][ia] = Jacobian[il][im][ia];
582  }
583  // boundary conditions, these are arrays size 1x1, i.e. values
584  alpha_upperBoundaryCondition_1d[0][0] = alpha_upperBoundaryCondition[il][im];
585  alpha_lowerBoundaryCondition_1d[0][0] = alpha_lowerBoundaryCondition[il][im];
586 
587  // Make A, B, C matrices
589  matr_A, matr_B, matr_C, // matrices A, B, C
590  L_1d, pc_1d, alpha_1d, // L, pc, alpha grids. We use L crossection here
591  1, 1, alpha.size, // L-size, pc-size, alpha-size
592  Zero_Matrix_2d, Zero_Matrix_2d, // L-boundaries
593  Zero_Matrix_2d, Zero_Matrix_2d, // pc-boundaries
594  alpha_lowerBoundaryCondition_1d, alpha_upperBoundaryCondition_1d, // alpha-boundaries
595  "BCT_CONSTANT_VALUE", "BCT_CONSTANT_VALUE", // Boundary conditions at L-boundaries
596  "BCT_CONSTANT_VALUE", "BCT_CONSTANT_VALUE", // Boundary conditions at pc-boundaries
597  alpha_lowerBoundaryCondition_calculationType, alpha_upperBoundaryCondition_calculationType, // Boundary conditions at alpha-boundaries
598  Zero_Matrix, // DLL - use 0 matrix
599  Zero_Matrix, Zero_Matrix, // Dpp, DppLpp - use 0 matrix
600  Daa_1d, DaaLpp_1d, // Daa, DaaLpp
601  Zero_Matrix, Zero_Matrix, // Dpca, Dpca_Lpp - use 0 matrix
602  Jacobian_1d, dt, Lpp);
603 
604  /*Output::echo(0, "Writing calculation matrices. ");
605  sprintf(out_name_matr_a, "./Debug/alpha_matr_A_%d_%d_0.dat", il, im);
606  sprintf(out_name_matr_b, "./Debug/alpha_matr_B_%d_%d_0.dat", il, im);
607  sprintf(out_name_matr_c, "./Debug/alpha_matr_C_%d_%d_0.dat", il, im);
608  matr_A.writeToFile(out_name_matr_a);
609  matr_B.writeToFile(out_name_matr_b);
610  matr_C.writeToFile(out_name_matr_c);
611  Output::echo(0, "done.\n");*/
612 
613  // make RHS = B*f + C
614  for (in = 0; in < alpha.size; in++) {
615  // If assume matr_B is 1d and has only +1, 0, -1 diagonals, can use this:
616  // 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];
617  RHS[in] = matr_C[0][in];
618  for (it = matr_B.begin(); it != matr_B.end(); it++) {
619  // multiplication B * f
620  if (in + it->first >= 0 && in + it->first < matr_B.total_size) {
621  RHS[in] += it->second[in] * slicePSD1D[in + it->first];
622  }
623  }
624  }
625 
626  // right hand side to files
627  // RHS.writeToFile("alpha_RHS.dat");
628 
629  // Solving model matrix, obtaining PSD
630  // Assume matr_A has only +1, 0, -1 diagonals
631  // TODO: Check if matr_A has only +1, 0, -1 diagonals, i.e. is tridiagonal
632  tridag(
633  &matr_A[-1][0], // just a way to get a pointer to first element of an array
634  &matr_A[0][0],
635  &matr_A[+1][0],
636  &RHS[0],
637  &slicePSD1D[0],
638  alpha.size);
639 
640  // copy PSD back, from 1d array to 3d array
641  for (ia = 0; ia < alpha.size; ia++) {
642  arr[il][im][ia] = slicePSD1D[ia];
643  }
644  }
645  }
646 }
647 
648 
649 /** Energy diffusion calculation function.
650 * Takes diffusion coefficients, boundary conditions and grid elements as input.
651 * It calculates model matrix by MakeModelMatrix_3D for 3D problem: A_3D * PSD_3D(t+1) = B_3D * PSD_3D(t) + C_3D
652 * then splits it for 1D problems: A * PSD_1D(t+1) = B * PSD_1D(t) + C
653 * and solve by tridiag method.
654 *
655 * \param dt - time step
656 * \param Lpp - plasma pause location
657 * \param &Dpcpc - L diffusion coefficient
658 * \param &DpcpcLpp - L diffusion coefficient
659 * \param &L - grid element L
660 * \param &pc - grid elelent pc
661 * \param &alpha - grid element alpha
662 * \param &Jacobian - jacobian
663 * \param &pc_lowerBoundaryCondition - lower boundary condition
664 * \param &pc_upperBoundaryCondition - upper boundary condition
665 * \param &pc_lowerBoundaryCondition_calculationType - lower boundary condition type (on value/on derivative)
666 * \param &pc_upperBoundaryCondition_calculationType - upper boundary condition type (on value/on derivative)
667 */
668 void PSD::Diffusion_pc(double dt, double Lpp,
669  DiffusionCoefficient &Dpcpc, DiffusionCoefficient &DpcpcLpp,
670  GridElement &L, GridElement &pc, GridElement &alpha, Matrix3D<double> Jacobian,
671  Matrix2D<double> pc_lowerBoundaryCondition,
672  Matrix2D<double> pc_upperBoundaryCondition,
673  string pc_lowerBoundaryCondition_calculationType,
674  string pc_upperBoundaryCondition_calculationType)
675 {
676  char out_name_matr_a[256], out_name_matr_b[256], out_name_matr_c[256];
677 
678 
679  // Calculation matrices A, B, and C for 1d diffusion problem
680  // A * PSD_1D(t+1) = B * PSD_1D(t) + C, for each L and Energy
681 
682  // Creation is moved to PSD.h, so it can be accessed from Main.cpp
683  //static CalculationMatrix
684  // matr_A(1, pc.size, 1, 1),
685  // matr_B(1, pc.size, 1, 0),
686  // matr_C(1, pc.size, 1, 0);
687  matr_A.Initialize(1, pc.size, 1, 1);
688  matr_B.Initialize(1, pc.size, 1, 0),
689  matr_C.Initialize(1, pc.size, 1, 0);
690 
691  // Create zero matrices, temporary matrices to be used later
692  static Matrix3D<double> Zero_Matrix(1, pc.size, 1);
693  static Matrix2D<double> Zero_Matrix_2d(1, 1);
694  static Matrix1D<double> Zero_Matrix_1d(alpha.size);
695  Zero_Matrix = 0;
696  Zero_Matrix_2d = 0;
697  Zero_Matrix_1d = 0;
698 
699  // Create 1d arrays
700  // Creating 1d arrays for PSD and right hand side (RHS = B * PSD + C)
701  static Matrix1D<double> slicePSD1D(pc.size), RHS(pc.size);
702 
703  // Creating 1d arrays which will be send to the MakeMatrix function
704  // Creating 1d matrixes, but with type of Matrix3D, cause we need that type for MakeMatrix function
705  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);
706  // Creating temporary arrays size 1x1 for boundary conditions, which are just values
707  static Matrix2D<double> pc_lowerBoundaryCondition_1d(1, 1), pc_upperBoundaryCondition_1d(1, 1);
708 
709  // iterators
710  int il, im, ia, in; // matrix element numbers on 3D grid
711  DiagMatrix::iterator it;
712  for (il = 0; il < L.size; il++) {
713  for (ia = 0; ia < alpha.size; ia++) {
714 
715  // Copying 3D arrays into 1d arrays
716  for (im = 0; im < pc.size; im++) {
717  // Copying 3D PSD to 1d PSD
718  slicePSD1D[im] = arr[il][im][ia];
719 
720  Dpcpc_1d[0][im][0] = Dpcpc.arr[il][im][ia];
721  DpcpcLpp_1d[0][im][0] = DpcpcLpp.arr[il][im][ia];
722 
723  L_1d[0][im][0] = L.arr[il][im][ia];
724  pc_1d[0][im][0] = pc.arr[il][im][ia];
725  alpha_1d[0][im][0] = alpha.arr[il][im][ia];
726 
727  Jacobian_1d[0][im][0] = Jacobian[il][im][ia];
728  }
729  // boundary conditions, these are arrays size 1x1, i.e. values
730  pc_upperBoundaryCondition_1d[0][0] = pc_upperBoundaryCondition[il][ia];
731  pc_lowerBoundaryCondition_1d[0][0] = pc_lowerBoundaryCondition[il][ia];
732 
733  // Make A, B, C matrices
735  matr_A, matr_B, matr_C, // matrices A, B, C
736  L_1d, pc_1d, alpha_1d, // L, pc, alpha grids.
737  1, pc.size, 1, // L-size, pc-size, alpha-size
738  Zero_Matrix_2d, Zero_Matrix_2d, // L-boundaries
739  pc_lowerBoundaryCondition_1d, pc_upperBoundaryCondition_1d, // pc-boundaries
740  Zero_Matrix_2d, Zero_Matrix_2d, // alpha-boundaries
741  "BCT_CONSTANT_VALUE", "BCT_CONSTANT_VALUE", // Boundary conditions at L-boundaries
742  pc_lowerBoundaryCondition_calculationType, pc_upperBoundaryCondition_calculationType, // Boundary conditions at pc-boundaries
743  "BCT_CONSTANT_VALUE", "BCT_CONSTANT_VALUE", // Boundary conditions at alpha-boundaries
744  Zero_Matrix, // DLL - use 0 matrix
745  Dpcpc_1d, DpcpcLpp_1d, // Dpp, DppLpp
746  Zero_Matrix, Zero_Matrix, // Daa, DaaLpp - use 0 matrix
747  Zero_Matrix, Zero_Matrix, // Dpca, Dpca_Lpp - use 0 matrix
748  Jacobian_1d, dt, Lpp);
749 
750  /*Output::echo(0, "Writing calculation matrices. ");
751  sprintf(out_name_matr_a, "./Debug/pc_matr_A_%d_0_%d.dat", il, ia);
752  sprintf(out_name_matr_b, "./Debug/pc_matr_B_%d_0_%d.dat", il, ia);
753  sprintf(out_name_matr_c, "./Debug/pc_matr_C_%d_0_%d.dat", il, ia);
754  matr_A.writeToFile(out_name_matr_a);
755  matr_B.writeToFile(out_name_matr_b);
756  matr_C.writeToFile(out_name_matr_c);
757  Output::echo(0, "done.\n");*/
758 
759  // make RHS = B*f + C
760  DiagMatrix::iterator it;
761  for (in = 0; in < pc.size; in++) {
762  // I guess matr_B has only main diagonal in implicit method, so it gan be used like this:
763  // RHS[in] = matr_B[0][in] * slicePSD1D[in] + matr_C[0][in];
764  RHS[in] = matr_C[0][in];
765  for (it = matr_B.begin(); it != matr_B.end(); it++) {
766  // multiplication B * f
767  if (in + it->first >= 0 && in + it->first < matr_B.total_size) {
768  RHS[in] += it->second[in] * slicePSD1D[in + it->first];
769  }
770  }
771  }
772 
773  // right hand side to files
774  // RHS.writeToFile("pc_RHS.dat");
775 
776  // Solving model matrix, obtaining PSD
777  // Assume matr_A has only +1, 0, -1 diagonals
778  // TODO: Check if matr_A has only +1, 0, -1 diagonals, i.e. is tridiagonal
779  tridag(
780  &matr_A[-1][0], // just a way to get a pointer to first element of an array
781  &matr_A[0][0],
782  &matr_A[+1][0],
783  &RHS[0],
784  &slicePSD1D[0],
785  pc.size);
786 
787  // copying PSD back, from 1d array to 3d array
788  for (im = 0; im < pc.size; im++) {
789  arr[il][im][ia] = slicePSD1D[im];
790  }
791  }
792  }
793 }
794 
795 
796 /** Radial diffusion calculation function.
797 * Takes diffusion coefficients, boundary conditions and grid elements as parameters.
798 * It caltulates model matrix by makeMatrix and then solve it by SolveMatrix method.
799 *
800 * \param dt - time step
801 * \param Lpp - plasma pause location
802 * \param &DLL - L diffusion coefficient
803 * \param &L - grid element L
804 * \param &pc - grid elelent pc
805 * \param &alpha - grid element alpha
806 * \param &Jacobian - jacobian
807 * \param &L_lowerBoundaryCondition - lower boundary condition on L
808 * \param &L_upperBoundaryCondition - upper boundary condition on L
809 * \param &L_lowerBoundaryCondition_calculationType - lower boundary condition on L calculation type
810 * \param &L_upperBoundaryCondition_calculationType - upper boundary condition on L calculation type
811 * \param tau - life time upper location of the plasma pause
812 * \param tauLpp - life time lower location of the plasma pause
813 */
814 void PSD::Diffusion_L( double dt, double Lpp,
816  GridElement &L, GridElement &pc, GridElement &alpha, Matrix3D<double> Jacobian,
817  Matrix2D<double> L_lowerBoundaryCondition,
818  Matrix2D<double> L_upperBoundaryCondition,
819  string L_lowerBoundaryCondition_calculationType,
820  string L_upperBoundaryCondition_calculationType,
821  double tau, double tauLpp)
822 {
823  char out_name_matr_a[256], out_name_matr_b[256], out_name_matr_c[256];
824 
825 
826  // Calculation matrices A, B, and C for 1d diffusion problem
827  // A * PSD_1D(t+1) = B * PSD_1D(t) + C, for each L and Energy
828 
829  // Creation is moved to PSD.h, so it can be accessed from Main.cpp
830  //static CalculationMatrix
831  // matr_A(1, 1, L.size, 1),
832  // matr_B(1, 1, L.size, 0),
833  // matr_C(1, 1, L.size, 0);
834  matr_A.Initialize(1, 1, L.size, 1);
835  matr_B.Initialize(1, 1, L.size, 0),
836  matr_C.Initialize(1, 1, L.size, 0);
837 
838  // Create zero matrices, temporary matrices to be used later
839  static Matrix3D<double> Zero_Matrix(L.size, 1, 1);
840  static Matrix2D<double> Zero_Matrix_2d(1, 1);
841  static Matrix1D<double> Zero_Matrix_1d(1);
842  Zero_Matrix = 0;
843  Zero_Matrix_2d = 0;
844  Zero_Matrix_1d = 0;
845 
846  // Create 1d arrays
847  // Creating 1d arrays for PSD and right hand side (RHS = B * PSD + C)
848  static Matrix1D<double> slicePSD1D(L.size), RHS(L.size);
849 
850  // Creating 1d arrays which will be send to the MakeMatrix function
851  // Creating 1d matrixes, but with type of Matrix3D, cause we need that type for MakeMatrix function
852  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)
853  // Creating temporary arrays size 1x1 for boundary conditions, which are just values
854  static Matrix2D<double> L_lowerBoundaryCondition_1d(1, 1), L_upperBoundaryCondition_1d(1, 1);
855 
856  // iterators
857  int il, im, ia, in; // matrix element numbers on 3D grid
858  DiagMatrix::iterator it;
859  for (ia = 0; ia < alpha.size; ia++) {
860  for (im = 0; im < pc.size; im++) {
861 
862  // Copying 3D arrays into 1d arrays
863  for (il = 0; il < L.size; il++) {
864  // Copying 3D PSD to 1d PSD
865  slicePSD1D[il] = arr[il][im][ia];
866 
867  DLL_1d[il][0][0] = DLL.arr[il][im][ia];
868 
869  L_1d[il][0][0] = L.arr[il][im][ia];
870  // pc_1d[il][0][0] = pc[il][im][ia];
871  // alpha_1d[il][0][0] = alpha[il][im][ia];
872 
873  Jacobian_1d[il][0][0] = Jacobian[il][im][ia];
874  }
875  // boundary conditions, these are arrays size 1x1, i.e. values
876  L_upperBoundaryCondition_1d[0][0] = L_upperBoundaryCondition[im][ia];
877  L_lowerBoundaryCondition_1d[0][0] = L_lowerBoundaryCondition[im][ia];
878 
879  // Make A, B, C matrices
881  matr_A, matr_B, matr_C, // matrices A, B, C
882  L_1d, Zero_Matrix, Zero_Matrix, // L, pc, alpha grids. We use L grid only here
883  L.size, 1, 1, // L-size, pc-size, alpha-size
884  L_lowerBoundaryCondition_1d, L_upperBoundaryCondition_1d, // L-boundaries
885  Zero_Matrix_2d, Zero_Matrix_2d, // pc-boundaries
886  Zero_Matrix_2d, Zero_Matrix_2d, // alpha-boundaries
887  L_lowerBoundaryCondition_calculationType, L_upperBoundaryCondition_calculationType, // Boundary conditions at L-boundaries
888  "BCT_CONSTANT_VALUE", "BCT_CONSTANT_VALUE", // Boundary conditions at pc-boundaries
889  "BCT_CONSTANT_VALUE", "BCT_CONSTANT_VALUE", // Boundary conditions at alpha-boundaries
890  DLL_1d, // DLL
891  Zero_Matrix, Zero_Matrix, // Dpp, DppLpp - use 0
892  Zero_Matrix, Zero_Matrix, // Daa, DaaLpp - use 0
893  Zero_Matrix, Zero_Matrix, // Dpca, Dpca_Lpp - use 0
894  Jacobian_1d, dt, Lpp,
895  tau, tauLpp);
896 
897  /*
898  Output::echo(0, "Writing calculation matrices. ");
899  sprintf(out_name_matr_a, "./Debug/L_matr_A_0_%d_%d.dat", im, ia);
900  sprintf(out_name_matr_b, "./Debug/L_matr_B_0_%d_%d.dat", im, ia);
901  sprintf(out_name_matr_c, "./Debug/L_matr_C_0_%d_%d.dat", im, ia);
902  matr_A.writeToFile(out_name_matr_a);
903  matr_B.writeToFile(out_name_matr_b);
904  matr_C.writeToFile(out_name_matr_c);
905  Output::echo(0, "done.\n");
906  */
907 
908  // make RHS = B*f + C
909  DiagMatrix::iterator it;
910  for (in = 0; in < L.size; in++) {
911  // I guess matr_B has only main diagonal in implicit method, so it can be used like this:
912  // RHS[in] = matr_B[0][in] * slicePSD1D[in] + matr_C[0][in];
913  RHS[in] = matr_C[0][in];
914  for (it = matr_B.begin(); it != matr_B.end(); it++) {
915  // multiplication B * f
916  if (in + it->first >= 0 && in + it->first < matr_B.total_size) {
917  RHS[in] += it->second[in] * slicePSD1D[in + it->first];
918  }
919  }
920  }
921 
922  // right hand side to files
923  // RHS.writeToFile("pc_RHS.dat");
924 
925  // Solving model matrix, obtaining PSD
926  // Assume matr_A has only +1, 0, -1 diagonals
927  // TODO: Check if matr_A has only +1, 0, -1 diagonals, i.e. is tridiagonal
928  tridag(
929  &matr_A[-1][0], // just a way to get a pointer to first element of an array
930  &matr_A[0][0],
931  &matr_A[+1][0],
932  &RHS[0],
933  &slicePSD1D[0],
934  L.size);
935 
936  // copying PSD back, from 1d array to 3d array
937  for (il = 0; il < L.size; il++) {
938  arr[il][im][ia] = slicePSD1D[il];
939  }
940  }
941  }
942 }
943 
944 
945 /** Radial diffusion calculation function.
946 * Takes diffusion coefficients, boundary conditions and grid elements as parameters.
947 * It caltulates model matrix by makeMatrix and then solve it by SolveMatrix method.
948 *
949 * \param dt - time step
950 * \param Lpp - plasma pause location
951 * \param &Dpcpc - diffusion coefficient
952 * \param &DpcpcLpp - diffusion coefficient
953 * \param &Daa - diffusion coefficient
954 * \param &DaaLpp - diffusion coefficient
955 * \param &Dpca - diffusion coefficient
956 * \param &DpcaLpp - diffusion coefficient
957 * \param &L - grid element L
958 * \param &pc - grid elelent pc
959 * \param &alpha - grid element alpha
960 * \param &Jacobian - jacobian
961 * \param &pc_lowerBoundaryCondition - lower boundary condition
962 * \param &pc_upperBoundaryCondition - upper boundary condition
963 * \param &alpha_lowerBoundaryCondition - lower boundary condition
964 * \param &alpha_upperBoundaryCondition - upper boundary condition
965 * \param &pc_lowerBoundaryCondition_calculationType - lower boundary condition calculation type
966 * \param &pc_upperBoundaryCondition_calculationType - upper boundary condition calculation type
967 * \param &alpha_lowerBoundaryCondition_calculationType - lower boundary condition calculation type
968 * \param &alpha_upperBoundaryCondition_calculationType - upper boundary condition calculation type
969 */
970 void PSD::Diffusion_pc_alpha(double dt, double Lpp,
971  DiffusionCoefficient &Dpcpc, DiffusionCoefficient &DpcpcLpp,
974  GridElement &L, GridElement &pc, GridElement &alpha, Matrix3D<double> Jacobian,
975  Matrix2D<double> pc_lowerBoundaryCondition,
976  Matrix2D<double> pc_upperBoundaryCondition,
977  Matrix2D<double> alpha_lowerBoundaryCondition,
978  Matrix2D<double> alpha_upperBoundaryCondition,
979  string pc_lowerBoundaryCondition_calculationType,
980  string pc_upperBoundaryCondition_calculationType,
981  string alpha_lowerBoundaryCondition_calculationType,
982  string alpha_upperBoundaryCondition_calculationType)
983 {
984 
985  // Calculation matrices A, B, and C
986  // A * PSD_1D(t+1) = B * PSD_1D(t) + C, for each L and Energy
987 
988  // Creation is moved to PSD.h, so it can be accessed from Main.cpp
989  //static CalculationMatrix
990  // matr_A(alpha.size, pc.size, 1, 1),
991  // matr_B(alpha.size, pc.size, 1, 0),
992  // matr_C(alpha.size, pc.size, 1, 0);
993  matr_A.Initialize(alpha.size, pc.size, 1, 1);
994  matr_B.Initialize(alpha.size, pc.size, 1, 0),
995  matr_C.Initialize(alpha.size, pc.size, 1, 0);
996 
997  // Create zero matrices, temporary matrices to be used later
998  static Matrix3D<double> Zero_Matrix(1, pc.size, alpha.size);
999  static Matrix2D<double> Zero_Matrix_2d(1, 1);
1000  static Matrix1D<double> Zero_Matrix_1d(1);
1001  Zero_Matrix = 0;
1002  Zero_Matrix_2d = 0;
1003  Zero_Matrix_1d = 0;
1004 
1005  // Create 1d arrays
1006  // Creating 1d arrays for PSD and right hand side (RHS = B * PSD + C)
1007  static Matrix1D<double> slicePSD1D(matr_A.total_size), RHS(matr_A.total_size);
1008  // Creating 1d arrays which will be send to the MakeMatrix function
1009  // Creating 1d matrixes, but with type of Matrix3D, cause we need that type for MakeMatrix function
1010  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);
1011  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);
1012  // Creating temporary arrays for boundary conditions
1013  static Matrix2D<double> pc_lowerBoundaryCondition_slice(1, alpha.size), pc_upperBoundaryCondition_slice(1, alpha.size);
1014  static Matrix2D<double> alpha_lowerBoundaryCondition_slice(1, pc.size), alpha_upperBoundaryCondition_slice(1, pc.size);
1015 
1016  // iterators
1017  int il, im, ia, in; // matrix element numbers on 3D grid
1018  DiagMatrix::iterator it;
1019  for (il = 0; il < L.size; il++) {
1020  // Copying 3D arrays into 1d arrays
1021  for (im = 0; im < pc.size; im++) {
1022  for (ia = 0; ia < alpha.size; ia++) {
1023  // Copying 3D PSD to 1d PSD
1024  in = matr_A.index1d(ia, im);
1025  slicePSD1D[in] = arr[il][im][ia];
1026 
1027  Dpcpc_slice[0][im][ia] = Dpcpc.arr[il][im][ia];
1028  DpcpcLpp_slice[0][im][ia] = DpcpcLpp.arr[il][im][ia];
1029  Daa_slice[0][im][ia] = Daa.arr[il][im][ia];
1030  DaaLpp_slice[0][im][ia] = DaaLpp.arr[il][im][ia];
1031  Dpca_slice[0][im][ia] = Dpca.arr[il][im][ia];
1032  DpcaLpp_slice[0][im][ia] = DpcaLpp.arr[il][im][ia];
1033 
1034  L_slice[0][im][ia] = L.arr[il][im][ia];
1035  pc_slice[0][im][ia] = pc.arr[il][im][ia];
1036  alpha_slice[0][im][ia] = alpha.arr[il][im][ia];
1037 
1038  Jacobian_slice[0][im][ia] = Jacobian[il][im][ia];
1039  }
1040  }
1041  // boundary conditions
1042  for (ia = 0; ia < alpha.size; ia++) {
1043  pc_upperBoundaryCondition_slice[0][ia] = pc_upperBoundaryCondition[il][ia];
1044  pc_lowerBoundaryCondition_slice[0][ia] = pc_lowerBoundaryCondition[il][ia];
1045  }
1046  for (im = 0; im < pc.size; im++) {
1047  alpha_upperBoundaryCondition_slice[0][im] = alpha_upperBoundaryCondition[il][im];
1048  alpha_lowerBoundaryCondition_slice[0][im] = alpha_lowerBoundaryCondition[il][im];
1049  }
1050 
1051  // Make
1053  matr_A, matr_B, matr_C, // matrices A, B, C
1054  L_slice, pc_slice, alpha_slice, // L, pc, alpha grids. We use pc and alpha crossections here
1055  1, pc.size, alpha.size, // L-size, pc-size, alpha-size
1056  Zero_Matrix_2d, Zero_Matrix_2d, // L-boundaries
1057  pc_lowerBoundaryCondition_slice, pc_upperBoundaryCondition_slice, // pc-boundaries
1058  alpha_lowerBoundaryCondition_slice, alpha_upperBoundaryCondition_slice, // alpha-boundaries
1059  "BCT_CONSTANT_VALUE", "BCT_CONSTANT_VALUE", // Boundary conditions at L-boundaries
1060  pc_lowerBoundaryCondition_calculationType, pc_upperBoundaryCondition_calculationType, // Boundary conditions at pc-boundaries
1061  alpha_lowerBoundaryCondition_calculationType, alpha_upperBoundaryCondition_calculationType, // Boundary conditions at alpha-boundaries
1062  Zero_Matrix, // DLL
1063  Dpcpc_slice, DpcpcLpp_slice, // Dpp, DppLpp - use 0
1064  Daa_slice, DaaLpp_slice, // Daa, DaaLpp - use 0
1065  Dpca_slice, DpcaLpp_slice, // Dpca, Dpca_Lpp - use 0
1066  Jacobian_slice, dt, Lpp);
1067 
1068  /* Output::echo(0, "Writing calculation matrices for L[%d] = %f... ", il, L[il][0][0]);
1069  matr_A.writeToFile("./Debug/pc_alpha_matr_A.dat");
1070  matr_B.writeToFile("./Debug/pc_alpha_matr_B.dat");
1071  matr_C.writeToFile("./Debug/pc_alpha_matr_C.dat");
1072  Output::echo(0, "done.\n"); */
1073 
1074  // make RHS = B*f + C
1075  DiagMatrix::iterator it;
1076  for (im = 0; im < pc.size; im++) {
1077  for (ia = 0; ia < alpha.size; ia++) {
1078  in = matr_B.index1d(ia, im);
1079  // I guess matr_B has only main diagonal in implicit method, so it can be used like this:
1080  // RHS[in] = matr_B[0][in] * slicePSD1D[in] + matr_C[0][in];
1081  RHS[in] = matr_C[0][in];
1082  for (it = matr_B.begin(); it != matr_B.end(); it++) {
1083  // multiplication B * f
1084  if (in + it->first >= 0 && in + it->first < matr_B.total_size) {
1085  RHS[in] += it->second[in] * slicePSD1D[in + it->first];
1086  }
1087  }
1088 
1089  }
1090  }
1091 
1092  // right hand side to files
1093  // RHS.writeToFile("./Debug/pc_alpha_RHS.dat");
1094 
1095  // Solving model matrix
1096  int n = pc.size * alpha.size;
1097  // need A and B for gauss method
1098  static Matrix2D<double> A;
1099  static Matrix1D<double> B;
1100  int modelMatrixLineNumber;
1101  int k, l, in;
1102  double sum, error, error_max;
1103  static bool recalculate_A = true;
1104  if (PSD_parameters.solutionMethod == "SM_Gauss") {
1105  // the routine does not actually invert matrix, but just solve AX=B eq
1106  if (recalculate_A) {
1107  B = RHS;
1108  if (!A.initialized) {
1109  A.AllocateMemory(n, n);
1110  }
1111  A = 0;
1112  // Output::echo("Writing A0 matrix... ");
1113  // A.writeToFile("A0.dat");
1114  // Output::echo("done.\n");
1115 
1116  Output::echo(5, "Inverting A matrix (gauss)... ");
1117  // Gauss method
1118  for (it = matr_A.begin(); it != matr_A.end(); it++) {
1119  for (modelMatrixLineNumber = 0; modelMatrixLineNumber < n; modelMatrixLineNumber++) {
1120  // check, if the element at line (modelMatrixLineNumber) and diagonal (it->first) is inside the matrix.
1121  // cause it can be outside, if line number is 1 and diagonal number is -5, for example.
1122  if (modelMatrixLineNumber + it->first >= 0 && modelMatrixLineNumber + it->first < n) {
1123  // converting matrix, stored as diagonals, into normal 2D matrix
1124  // cause we need 2D matrix to solve by Gauss method
1125  // cout << it->first << " - " << it->second[modelMatrixLineNumber] << endl;
1126  A[modelMatrixLineNumber][modelMatrixLineNumber + it->first] = it->second[modelMatrixLineNumber];
1127  }
1128  }
1129  }
1130 
1131  // inversion of A
1132  // sending address of the first element of the plane 2d array
1133  // it's not actually matrix inversion
1134  // it just makes A-triangle, so it's easy to solve A*X = B
1135  gauss_solve(&A[0][0], &B[0], n);
1136  recalculate_A = true;
1137  Output::echo(5, "done.\n");
1138  }
1139 
1140  Output::echo(5, "Calculation of new PSD... ");
1141  slicePSD1D[n-1] = B[n-1] / A[n-1][n-1];
1142  for (k = n - 2; k >= 0; k--) {
1143  sum = 0;
1144  for (l = k + 1; l < n; l++)
1145  sum += A[k][l] * slicePSD1D[l];
1146  slicePSD1D[k] = B[k] - sum;
1147  }
1148  /* end calc */
1149 
1150  /* control */
1151  //mult_vect2(&A[0], &slicePSD1D[0], &RHS[0], n);
1152  error_max = 0.;
1153  error = 0;
1154  for (in = 0; in < n; in++) {
1155  error = RHS[in];
1156  for (it = matr_A.begin(); it != matr_A.end(); it++) {
1157  // multiplying all diagonals to PSD to check error = RHS - A*PSD = 0
1158  if (in + it->first >= 0 && in + it->first < n) {
1159  error -= it->second[in] * slicePSD1D[in + it->first];
1160  }
1161  }
1162  error_max = VF::max( error_max, error );
1163  }
1164  Output::echo(5, "done.\n");
1165  Output::echo(5, "Error_norm = %e \n", error_max);
1166 
1167  } else if (PSD_parameters.solutionMethod == "SM_Relaxation") {
1168  // Relaxation method
1169 
1170  Output::echo(5, "Inverting A matrix (relaxation)... ");
1171  // runs the over_relaxation_diag solution method
1172  over_relaxation_diag(matr_A, RHS, slicePSD1D);
1173  Output::echo(5, "done.");
1174 
1175  } else if (PSD_parameters.solutionMethod == "SM_Lapack") {
1176  // Relaxation method
1177 
1178  Output::echo(5, "Inverting A matrix (Lapack)... ");
1179  // runs the over_relaxation_diag solution method
1180  Lapack(matr_A, RHS, slicePSD1D);
1181  Output::echo(5, "done.");
1182 
1183  } else if (PSD_parameters.solutionMethod == "SM_GMRES") {
1184  // Relaxation method
1185 
1186  Output::echo(5, "Inverting A matrix (GMRES)... ");
1187  // runs the over_relaxation_diag solution method
1188  //gmres_wrapout(matr_A, RHS, slicePSD1D, PSD_parameters.SOL_maxiter, PSD_parameters.SOL_i_max, PSD_parameters.SOL_max_iter_err);
1189  gmres_wrapout(matr_A, RHS, slicePSD1D, PSD_parameters.GMRES_parameters);
1190  Output::echo(5, "done.");
1191 
1192  } else if (PSD_parameters.solutionMethod == "SM_GMRES2") {
1193  // Relaxation method
1194 
1195  Output::echo(5, "Inverting A matrix (GMRES2)... ");
1196  // runs the over_relaxation_diag solution method
1197  //gmres_wrapout(matr_A, RHS, slicePSD1D, PSD_parameters.SOL_maxiter, PSD_parameters.SOL_i_max, PSD_parameters.SOL_max_iter_err);
1198  gmres2_wrapout(matr_A, RHS, slicePSD1D, PSD_parameters.GMRES_parameters);
1199  Output::echo(5, "done.");
1200 
1201 
1202  } else if (PSD_parameters.solutionMethod == "SM_Tridiag") {
1203  // Tridiagonal method
1204  throw error_msg("SOLUTIONERR", "Tridiagonal method can't be used");
1205  } else { // Unknown method
1206  throw error_msg("SOLUTIONERR", "Unknown solution method");
1207  }
1208 
1209  // copying PSD back, from 1d array to 3d array
1210  for (im = 0; im < pc.size; im++) {
1211  for (ia = 0; ia < alpha.size; ia++) {
1212  in = matr_A.index1d(ia, im);
1213  arr[il][im][ia] = slicePSD1D[in];
1214  }
1215  }
1216  }
1217 
1218 }
1219 
1220 
1221 void PSD::Diffusion_pc_alpha_KC(double dt, double Lpp,
1222  DiffusionCoefficient &Dpcpc, DiffusionCoefficient &DpcpcLpp,
1225  GridElement &L, GridElement &pc, GridElement &alpha, Matrix3D<double> Jacobian,
1226  Matrix2D<double> pc_lowerBoundaryCondition,
1227  Matrix2D<double> pc_upperBoundaryCondition,
1228  Matrix2D<double> alpha_lowerBoundaryCondition,
1229  Matrix2D<double> alpha_upperBoundaryCondition,
1230  string pc_lowerBoundaryCondition_calculationType,
1231  string pc_upperBoundaryCondition_calculationType,
1232  string alpha_lowerBoundaryCondition_calculationType,
1233  string alpha_upperBoundaryCondition_calculationType)
1234 {
1235 
1236  // Calculation matrices A, B, and C
1237  // A * PSD_1D(t+1) = B * PSD_1D(t) + C, for each L and Energy
1238 
1239  // Creation is moved to PSD.h, so it can be accessed from Main.cpp
1240  //static CalculationMatrix
1241  // matr_A(alpha.size, pc.size, 1, 1),
1242  // matr_B(alpha.size, pc.size, 1, 0),
1243  // matr_C(alpha.size, pc.size, 1, 0);
1244  matr_A.Initialize(alpha.size, pc.size, 1, 1);
1245  matr_B.Initialize(alpha.size, pc.size, 1, 0),
1246  matr_C.Initialize(alpha.size, pc.size, 1, 0);
1247 
1248 
1249 
1250 
1251  // Create zero matrices, temporary matrices to be used later
1252  static Matrix3D<double> Zero_Matrix(1, pc.size, alpha.size);
1253  static Matrix2D<double> Zero_Matrix_2d(1, 1);
1254  static Matrix1D<double> Zero_Matrix_1d(1);
1255  Zero_Matrix = 0;
1256  Zero_Matrix_2d = 0;
1257  Zero_Matrix_1d = 0;
1258 
1259  // Create 1d arrays
1260  // Creating 1d arrays for PSD and right hand side (RHS = B * PSD + C)
1261  static Matrix1D<double> slicePSD1D(matr_A.total_size), RHS(matr_A.total_size);
1262  // Creating 1d arrays which will be send to the MakeMatrix function
1263  // Creating 1d matrixes, but with type of Matrix3D, cause we need that type for MakeMatrix function
1264  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);
1265  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);
1266  // Creating temporary arrays for boundary conditions
1267  static Matrix2D<double> pc_lowerBoundaryCondition_slice(1, alpha.size), pc_upperBoundaryCondition_slice(1, alpha.size);
1268  static Matrix2D<double> alpha_lowerBoundaryCondition_slice(1, pc.size), alpha_upperBoundaryCondition_slice(1, pc.size);
1269 
1270  // iterators
1271  int il, im, ia, in; // matrix element numbers on 3D grid
1272  DiagMatrix::iterator it;
1273  for (il = 0; il < L.size; il++) {
1274  // Copying 3D arrays into 1d arrays
1275  for (im = 0; im < pc.size; im++) {
1276  for (ia = 0; ia < alpha.size; ia++) {
1277  // Copying 3D PSD to 1d PSD
1278  in = matr_A.index1d(ia, im);
1279  slicePSD1D[in] = arr[il][im][ia];
1280 
1281  Dpcpc_slice[0][im][ia] = Dpcpc.arr[il][im][ia];
1282  DpcpcLpp_slice[0][im][ia] = DpcpcLpp.arr[il][im][ia];
1283  Daa_slice[0][im][ia] = Daa.arr[il][im][ia];
1284  DaaLpp_slice[0][im][ia] = DaaLpp.arr[il][im][ia];
1285  Dpca_slice[0][im][ia] = Dpca.arr[il][im][ia];
1286  DpcaLpp_slice[0][im][ia] = DpcaLpp.arr[il][im][ia];
1287 
1288  L_slice[0][im][ia] = L.arr[il][im][ia];
1289  pc_slice[0][im][ia] = pc.arr[il][im][ia];
1290  alpha_slice[0][im][ia] = alpha.arr[il][im][ia];
1291 
1292  Jacobian_slice[0][im][ia] = Jacobian[il][im][ia];
1293  }
1294  }
1295  // boundary conditions
1296  for (ia = 0; ia < alpha.size; ia++) {
1297  pc_upperBoundaryCondition_slice[0][ia] = pc_upperBoundaryCondition[il][ia];
1298  pc_lowerBoundaryCondition_slice[0][ia] = pc_lowerBoundaryCondition[il][ia];
1299  }
1300  for (im = 0; im < pc.size; im++) {
1301  alpha_upperBoundaryCondition_slice[0][im] = alpha_upperBoundaryCondition[il][im];
1302  alpha_lowerBoundaryCondition_slice[0][im] = alpha_lowerBoundaryCondition[il][im];
1303  }
1304 
1305  // Make
1307  matr_A, matr_B, matr_C, // matrices A, B, C
1308  L_slice, pc_slice, alpha_slice, // L, pc, alpha grids. We use pc and alpha crossections here
1309  1, pc.size, alpha.size, // L-size, pc-size, alpha-size
1310  Zero_Matrix_2d, Zero_Matrix_2d, // L-boundaries
1311  pc_lowerBoundaryCondition_slice, pc_upperBoundaryCondition_slice, // pc-boundaries
1312  alpha_lowerBoundaryCondition_slice, alpha_upperBoundaryCondition_slice, // alpha-boundaries
1313  "BCT_CONSTANT_VALUE", "BCT_CONSTANT_VALUE", // Boundary conditions at L-boundaries
1314  pc_lowerBoundaryCondition_calculationType, pc_upperBoundaryCondition_calculationType, // Boundary conditions at pc-boundaries
1315  alpha_lowerBoundaryCondition_calculationType, alpha_upperBoundaryCondition_calculationType, // Boundary conditions at alpha-boundaries
1316  Zero_Matrix, // DLL
1317  Dpcpc_slice, DpcpcLpp_slice, // Dpp, DppLpp - use 0
1318  Daa_slice, DaaLpp_slice, // Daa, DaaLpp - use 0
1319  Dpca_slice, DpcaLpp_slice, // Dpca, Dpca_Lpp - use 0
1320  Jacobian_slice, dt, Lpp);
1321 
1322  /* Output::echo(0, "Writing calculation matrices for L[%d] = %f... ", il, L[il][0][0]);
1323  matr_A.writeToFile("./Debug/pc_alpha_matr_A.dat");
1324  matr_B.writeToFile("./Debug/pc_alpha_matr_B.dat");
1325  matr_C.writeToFile("./Debug/pc_alpha_matr_C.dat");
1326  Output::echo(0, "done.\n"); */
1327 
1328  // make RHS = B*f + C
1329  DiagMatrix::iterator it;
1330  for (im = 0; im < pc.size; im++) {
1331  for (ia = 0; ia < alpha.size; ia++) {
1332  in = matr_B.index1d(ia, im);
1333  // I guess matr_B has only main diagonal in implicit method, so it gan be used like this:
1334  // RHS[in] = matr_B[0][in] * slicePSD1D[in] + matr_C[0][in];
1335  RHS[in] = matr_C[0][in];
1336  for (it = matr_B.begin(); it != matr_B.end(); it++) {
1337  // multiplication B * f
1338  if (in + it->first >= 0 && in + it->first < matr_B.total_size) {
1339  RHS[in] += it->second[in] * slicePSD1D[in + it->first];
1340  }
1341  }
1342 
1343  }
1344  }
1345 
1346  // right hand side to files
1347  // RHS.writeToFile("./Debug/pc_alpha_RHS.dat");
1348 
1349  // Solving model matrix
1350  int n = pc.size * alpha.size;
1351  // need A and B for gauss method
1352  static Matrix2D<double> A;
1353  static Matrix1D<double> B;
1354  int modelMatrixLineNumber;
1355  int k, l, in;
1356  double sum, error, error_max;
1357  static bool recalculate_A = true;
1358  if (PSD_parameters.solutionMethod == "SM_Gauss") {
1359  // the routine does not actually invert matrix, but just solve AX=B eq
1360  if (recalculate_A) {
1361  B = RHS;
1362  if (!A.initialized) {
1363  A.AllocateMemory(n, n);
1364  }
1365  A = 0;
1366  // Output::echo("Writing A0 matrix... ");
1367  // A.writeToFile("A0.dat");
1368  // Output::echo("done.\n");
1369 
1370  Output::echo(5, "Inverting A matrix (gauss)... ");
1371  // Gauss method
1372  for (it = matr_A.begin(); it != matr_A.end(); it++) {
1373  for (modelMatrixLineNumber = 0; modelMatrixLineNumber < n; modelMatrixLineNumber++) {
1374  // check, if the element at line (modelMatrixLineNumber) and diagonal (it->first) is inside the matrix.
1375  // cause it can be outside, if line number is 1 and diagonal number is -5, for example.
1376  if (modelMatrixLineNumber + it->first >= 0 && modelMatrixLineNumber + it->first < n) {
1377  // converting matrix, stored as diagonals, into normal 2D matrix
1378  // cause we need 2D matrix to solve by Gauss method
1379  // cout << it->first << " - " << it->second[modelMatrixLineNumber] << endl;
1380  A[modelMatrixLineNumber][modelMatrixLineNumber + it->first] = it->second[modelMatrixLineNumber];
1381  }
1382  }
1383  }
1384 
1385  // inversion of A
1386  // sending address of the first element of the plane 2d array
1387  // it's not actually matrix inversion
1388  // it just makes A-triangle, so it's easy to solve A*X = B
1389  gauss_solve(&A[0][0], &B[0], n);
1390  recalculate_A = true;
1391  Output::echo(5, "done.\n");
1392  }
1393 
1394  Output::echo(5, "Calculation of new PSD... ");
1395  slicePSD1D[n-1] = B[n-1] / A[n-1][n-1];
1396  for (k = n - 2; k >= 0; k--) {
1397  sum = 0;
1398  for (l = k + 1; l < n; l++)
1399  sum += A[k][l] * slicePSD1D[l];
1400  slicePSD1D[k] = B[k] - sum;
1401  }
1402  /* end calc */
1403 
1404  /* control */
1405  //mult_vect2(&A[0], &slicePSD1D[0], &RHS[0], n);
1406  error_max = 0.;
1407  error = 0;
1408  for (in = 0; in < n; in++) {
1409  error = RHS[in];
1410  for (it = matr_A.begin(); it != matr_A.end(); it++) {
1411  // multiplying all diagonals to PSD to check error = RHS - A*PSD = 0
1412  if (in + it->first >= 0 && in + it->first < n) {
1413  error -= it->second[in] * slicePSD1D[in + it->first];
1414  }
1415  }
1416  error_max = VF::max( error_max, error );
1417  }
1418  Output::echo(5, "done.\n");
1419  Output::echo(5, "Error_norm = %e \n", error_max);
1420 
1421  } else if (PSD_parameters.solutionMethod == "SM_Relaxation") {
1422  // Relaxation method
1423 
1424  Output::echo(5, "Inverting A matrix (relaxation)... ");
1425  // runs the over_relaxation_diag solution method
1426  over_relaxation_diag(matr_A, RHS, slicePSD1D);
1427  Output::echo(5, "done.");
1428 
1429  } else if (PSD_parameters.solutionMethod == "SM_Lapack") {
1430  // Relaxation method
1431 
1432  Output::echo(5, "Inverting A matrix (Lapack)... ");
1433  // runs the over_relaxation_diag solution method
1434  Lapack(matr_A, RHS, slicePSD1D);
1435  Output::echo(5, "done.");
1436 
1437  } else if (PSD_parameters.solutionMethod == "SM_GMRES") {
1438  // Relaxation method
1439 
1440  Output::echo(5, "Inverting A matrix (GMRES)... ");
1441  // runs the over_relaxation_diag solution method
1442  //gmres_wrapout(matr_A, RHS, slicePSD1D, PSD_parameters.SOL_maxiter, PSD_parameters.SOL_i_max, PSD_parameters.SOL_max_iter_err);
1443  gmres_wrapout(matr_A, RHS, slicePSD1D, PSD_parameters.GMRES_parameters);
1444  Output::echo(5, "done.");
1445 
1446  } else if (PSD_parameters.solutionMethod == "SM_GMRES2") {
1447  // Relaxation method
1448 
1449  Output::echo(5, "Inverting A matrix (GMRES2)... ");
1450  // runs the over_relaxation_diag solution method
1451  //gmres_wrapout(matr_A, RHS, slicePSD1D, PSD_parameters.SOL_maxiter, PSD_parameters.SOL_i_max, PSD_parameters.SOL_max_iter_err);
1452  gmres2_wrapout(matr_A, RHS, slicePSD1D, PSD_parameters.GMRES_parameters);
1453  Output::echo(5, "done.");
1454 
1455 
1456  } else if (PSD_parameters.solutionMethod == "SM_Tridiag") {
1457  // Tridiagonal method
1458  throw error_msg("SOLUTIONERR", "Tridiagonal method can't be used");
1459  } else { // Unknown method
1460  throw error_msg("SOLUTIONERR", "Unknown solution method");
1461  }
1462 
1463  // copying PSD back, from 1d array to 3d array
1464  for (im = 0; im < pc.size; im++) {
1465  for (ia = 0; ia < alpha.size; ia++) {
1466  in = matr_A.index1d(ia, im);
1467  arr[il][im][ia] = slicePSD1D[in];
1468  }
1469  }
1470  }
1471 
1472 }
1473 
1474 
1475 
1476 /** Checking for infinity value appears during interpolation for 3D arrays.
1477 *
1478 * \param arr - array of values
1479 * \param il - index
1480 * \param im - index
1481 * \param ia - index
1482 * \param maxNum = 1e99 - max number (to compare with)
1483 */
1484 void checkInf_3D(Matrix3D<double> arr, int il, int im, int ia, double maxNum = 1e99) {
1485  // Some wired 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.
1486  // \todo Check checkInf.
1487  if (!(arr[il][im][ia] >= 0 || arr[il][im][ia] <= 0) || arr[il][im][ia] > maxNum || arr[il][im][ia] < -maxNum){
1488  throw error_msg("INTERPOLATION_ERROR", "Interpolation error at [%d][%d][%d]", il, im, ia);
1489  }
1490 }
1491 
1492 /** Checking for infinity value appears during interpolation for 1D arrays.
1493 *
1494 * \param arr - array of values
1495 * \param il - index
1496 * \param im - index
1497 * \param ia - index
1498 * \param maxNum = 1e99 - max number (to compare with)
1499 */
1500 void checkInf_1D(Matrix1D<double> arr, int il, int im, int ia, double maxNum = 1e99) {
1501  // Some wired 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.
1502  if (!(arr[im] >= 0 || arr[im] <= 0) || arr[im] > maxNum || arr[im] < -maxNum){
1503  throw error_msg("INTERPOLATION_ERROR", "Interpolation error at [%d][%d][%d]", il, im, ia);
1504  }
1505 }
1506 
1507 /// It does interpolation... for log(function) or for just function, depends on parameters.
1508 void PSD::Interpolate(PSD &otherPSD, Parameters_structure::Interpolation interpolationParameters_structure, Grid &oldGrid, Grid &newGrid, Matrix2D<double> newGrid_pc_lowerBoundaryCondition, Matrix2D<double> newGrid_pc_upperBoundaryCondition) {
1509  // check...
1510 
1511  int il, im, ia;
1512  double lowerBoundary1d=0, upperBoundary1d=0;
1513 
1514  // if interoplation method is IT_NONE - exit
1515  // if grids have the same type, we just copy the values
1516  if (interpolationParameters_structure.type == "IT_NONE") {
1517  return;
1518  } else if (oldGrid.type == newGrid.type) { // \todo If grid types are the same, it's not necessary that grids are the same
1519  // arr.Matrix3D<double>::operator =(otherPSD);
1520 
1521  for (il = 0; il < oldGrid.L.size; il++) {
1522  for (ia = 0; ia < oldGrid.alpha.size; ia++) {
1523  for (im = 0; im < oldGrid.pc.size; im++) {
1524  // copying of the values from one grid to anouther
1525  arr[il][im][ia] = otherPSD.arr[il][im][ia];
1526  // some wired bug. Model matrix does not work correctly without this. Don't know why.
1527  if (arr[il][im][ia] < VC::zero_f) {
1528  arr[il][im][ia] = VC::zero_f;
1529  }
1530  }
1531  }
1532  }
1533  return;
1534  }
1535 
1536  Matrix1D<double> matrix_array_1d_old(oldGrid.pc.size);
1537  Matrix1D<double> matrix_array_1d_new(newGrid.pc.size);
1538  Matrix1D<double> grid_1d_old(oldGrid.pc.size);
1539  Matrix1D<double> grid_1d_new(newGrid.pc.size);
1540 
1541  for (il = 0; il < oldGrid.L.size; il++) {
1542  for (ia = 0; ia < oldGrid.alpha.size; ia++) {
1543  for (im = 0; im < oldGrid.pc.size; im++) {
1544 
1545  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]) {
1546  // if L-size or alpha-size are different, then we can not interpolate by that interpolation function
1547  throw error_msg("INTERPOLATION_ERROR", "Interpolation error: grids must have the same L and Alpha");
1548  }
1549 #ifdef DEBUG_MODE
1550  // check for infinity and NAN
1551  checkInf_3D(arr, il, im, ia);
1552 #endif
1553 
1554  // some work around boundary conditions
1555  if (interpolationParameters_structure.useLog == "Log_10") {
1556  // if we are interpolating log of values, then store log of values to 1d array
1557  if (otherPSD.arr[il][im][ia] >= VC::zero_f) {
1558  matrix_array_1d_old[im] = log10(otherPSD.arr[il][im][ia]);
1559  } else {
1560  // if PSD lower then lowest possible value (it should not be zero, because we are interpolating log of values for example)
1561  matrix_array_1d_old[im] = log10(VC::zero_f);
1562  }
1563  // theoretically we don't need to make the boundary conditions not lover then 10^21
1564  // because we check that in its constructor
1565  // But just in case...
1566  lowerBoundary1d = log10(VF::max(newGrid_pc_lowerBoundaryCondition[il][ia], VC::zero_f));
1567  upperBoundary1d = log10(VF::max(newGrid_pc_upperBoundaryCondition[il][ia], VC::zero_f));
1568  }
1569  else if (interpolationParameters_structure.useLog == "Log_2") {
1570  // if we are interpolating log of values, then store log of values to 1d array
1571  if (otherPSD.arr[il][im][ia] >= VC::zero_f) {
1572  matrix_array_1d_old[im] = log(otherPSD.arr[il][im][ia])/log(2.0);
1573  } else {
1574  // if PSD lower then lowest possible value (it should not be zero, because we are interpolating log of values for example)
1575  matrix_array_1d_old[im] = log(VC::zero_f)/log(2.0);
1576  }
1577  // theoretically we don't need to make the boundary conditions not lover then 10^21
1578  // because we check that in its constructor
1579  // But just in case...
1580 
1581  // adrozdov: Bug report. log10 replaced by log2
1582  lowerBoundary1d = log(VF::max(newGrid_pc_lowerBoundaryCondition[il][ia], VC::zero_f)) / log(2.0);
1583  upperBoundary1d = log(VF::max(newGrid_pc_upperBoundaryCondition[il][ia], VC::zero_f)) / log(2.0);
1584  }
1585  else if (interpolationParameters_structure.useLog == "Log_E") {
1586  // if we are interpolating log of values, then store log of values to 1d array
1587  if (otherPSD.arr[il][im][ia] >= VC::zero_f) {
1588  matrix_array_1d_old[im] = log(otherPSD.arr[il][im][ia]);
1589  } else {
1590  // if PSD lower then lowest possible value (it should not be zero, because we are interpolating log of values for example)
1591  matrix_array_1d_old[im] = log(VC::zero_f);
1592  }
1593  // theoretically we don't need to make the boundary conditions not lover then 10^21
1594  // because we check that in its constructor
1595  // But just in case...
1596 
1597  // adrozdov: Bug report. log10 replaced by log
1598  lowerBoundary1d = log(VF::max(newGrid_pc_lowerBoundaryCondition[il][ia], VC::zero_f));
1599  upperBoundary1d = log(VF::max(newGrid_pc_upperBoundaryCondition[il][ia], VC::zero_f));
1600  }
1601  else if (interpolationParameters_structure.useLog == "NoLog") {
1602  matrix_array_1d_old[im] = otherPSD.arr[il][im][ia];
1603  lowerBoundary1d = VF::max(newGrid_pc_lowerBoundaryCondition[il][ia], VC::zero_f);
1604  upperBoundary1d = VF::max(newGrid_pc_upperBoundaryCondition[il][ia], VC::zero_f);
1605  }
1606  else {
1607  throw error_msg("INTERPOLATION_ERROR", "Usage of log(PSD) is not define properly: %s", interpolationParameters_structure.useLog.c_str());
1608  break;
1609  }
1610 
1611 #ifdef DEBUG_MODE
1612  // check for infinity and NAN
1613  checkInf_1D(matrix_array_1d_old, il, im, ia, 99);
1614 #endif
1615 
1616  grid_1d_old[im] = oldGrid.pc.arr[il][im][ia];
1617  grid_1d_new[im] = newGrid.pc.arr[il][im][ia];
1618  //cout<<grid_1d_old[im]<<"kim......."<<grid_1d_new[im]<<"..........."<<matrix_array_1d_old[im]<<endl;
1619  }
1620 
1621  if (interpolationParameters_structure.type == "IT_LINEAR") {
1622  // do linear interpolation
1623  matrix_array_1d_new.Interpolate(matrix_array_1d_old, grid_1d_old, grid_1d_new);
1624  } else if (interpolationParameters_structure.type == "IT_POLYNOMIAL") {
1625  // do polinomial - linear interpolation (linear close to the boundaryes)
1626  matrix_array_1d_new.Polilinear(matrix_array_1d_old, grid_1d_old, grid_1d_new, lowerBoundary1d, upperBoundary1d);
1627  } else if (interpolationParameters_structure.type == "IT_SPLINE") {
1628  // do spline interpolation
1629  matrix_array_1d_new.Spline(matrix_array_1d_old, grid_1d_old, grid_1d_new, lowerBoundary1d, upperBoundary1d, interpolationParameters_structure.linearSplineCoef, interpolationParameters_structure.maxSecondDerivative);
1630  } else if (interpolationParameters_structure.type == "IT_COPY") {
1631  // just copy
1632  matrix_array_1d_new = matrix_array_1d_old;
1633  } else { // unknown
1634  throw error_msg("INTERPOLATION_TYPE_ERROR", "Unknown interpolation type");
1635  }
1636 
1637  for (im = 0; im < newGrid.pc.size; im++) {
1638 
1639 #ifdef DEBUG_MODE
1640  checkInf_1D(matrix_array_1d_new, il, im, ia, 99);
1641 #endif
1642 
1643  // copying values back to 3d arrays
1644  if (interpolationParameters_structure.useLog == "Log_10") arr[il][im][ia] = pow(10, matrix_array_1d_new[im]);
1645  else if (interpolationParameters_structure.useLog == "Log_2") arr[il][im][ia] = pow(2, matrix_array_1d_new[im]);
1646  else if (interpolationParameters_structure.useLog == "Log_E") arr[il][im][ia] = pow(VC::exp, matrix_array_1d_new[im]);
1647  else arr[il][im][ia] = matrix_array_1d_new[im];
1648 
1649 #ifdef DEBUG_MODE
1650  checkInf_3D(arr, il, im, ia);
1651 #endif
1652 
1653  }
1654  }
1655  }
1656 
1657  // check result...
1658 }
1659 
1660 /*-------------------------------------------------------------------------------------------------------------------------------
1661 Outputs
1662 -------------------------------------------------------------------------------------------------------------------------------*/
1663 
1664 /**
1665 * %PSD output.
1666 *
1667 * \param time - time
1668 */
1669 void PSD::Output_without_grid(double time) {
1670  int il, im, ia;
1671  output_without_grid_file << "ZONE T=\"" << internal << time << "\" I=" << arr.size_z << ", J=" << arr.size_y << ", K=" << arr.size_x << endl;
1672  for (il = 0; il < arr.size_x; il++) {
1673  for (im = 0; im < arr.size_y; im++) {
1674  for (ia = 0; ia < arr.size_z; ia++) {
1675  output_without_grid_file << arr[il][im][ia] << endl;
1676  }
1677  }
1678  }
1679 
1680 }
1681 
1682 
1683 /**
1684 * Load initial %PSD from file
1685 * \param &L - grid element L
1686 * \param &epc - grid element pc
1687 * \param &alpha - grid element alpha
1688 * \param filename - filename
1689 * \param withGrid - indicate if the file contains grid values
1690 */
1691 void PSD::Load_initial_f_file(GridElement &L, GridElement &epc, GridElement &alpha, const char *filename, bool withGrid) {
1692  int il, im, ia;
1693  string inBuf;
1694  double Load_L=-1, Load_epc=-1, Load_alpha=-1;// temporaraly variables for the grid point
1695  double err = 1.e-3;
1696 
1697  // Load j from input file
1698  ifstream input_f(filename);
1699  if (input_f != NULL && !input_f.eof()) {
1700  // Skipping first two lines. Should serch for 'ZONE' and read from following line better.
1701  getline(input_f, inBuf);
1702  getline(input_f, inBuf);
1703  for (il = 0; il < L.size; il++) {
1704  for (im = 0; im < epc.size; im++) {
1705  for (ia = 0; ia < alpha.size; ia++) {
1706  if (withGrid) input_f >> Load_L >> Load_epc >> Load_alpha;
1707  input_f >> arr[il][im][ia];
1708  if (withGrid)
1709  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) {
1710  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]);
1711  }
1712  }
1713  }
1714  }
1715  } else {
1716  throw error_msg("INITIAL_PSD_ERROR", "Error reading PSD input file %s.\n", filename);
1717  //cout << "Error reading initial conditions input file " << filename << endl;
1718  //getchar();
1719  //exit(-1);
1720  }
1721  input_f.close();
1722 }
1723 
1724 ///
1725 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) {
1726  int il, im, ia;
1727  double p1, Lmax, Lmin, scale_factor, Lgrid;
1728  int Lsize_7, addGrid; // These two numbers should be of integer type.
1729  Matrix1D<double> L_1d, fL;
1730  //Matrix2D<double> fLm(L.size, pc.size);
1731 
1732  //////////
1733  // Steady state solution for any outer-L boundary, D. Subbotin implimentation
1734  // Create 1D grid at min energy and 90 deg alpha
1735  // First, find a closest index of 90 deg pitch-angle
1736  int ia_90=alpha.size-1;
1737  while (alpha.arr[L.size-1][0][ia_90] > VC::pi/2 && ia_90 > 0) ia_90--;
1738 
1739  // Make sure we have a point at L=7
1740  // Our boundary values J_L7() are at L=7, so we need to know the steady state solution at L=7.
1741  if (L.arr[L.size-1][0][ia_90] < 7) {
1742  // If maximum L value is less than 7, we add another points from Lmax to 7
1743 
1744  // estimate grid points to add
1745  Lmax = L.arr[L.size-1][0][ia_90];
1746  Lmin = L.arr[0][0][ia_90];
1747  Lgrid = (Lmax-Lmin)/(L.size-1);
1748  Lsize_7 = (7.0-Lmin)/Lgrid + 1;
1749  addGrid = Lsize_7 - L.size;
1750 
1751  // construct new L and psd grid
1752  L_1d.AllocateMemory(L.size+addGrid);
1753  fL.AllocateMemory(L.size+addGrid);
1754 
1755  for (il = 0; il < L.size; il++) {
1756  L_1d[il] = L.arr[il][0][ia_90];
1757  }
1758  for (il = L.size; il < Lsize_7; il++) { // add another points from Lmax to 7
1759  L_1d[il] = Lmin+Lgrid*il;
1760  }
1761 
1762  } else {
1763  Lsize_7 = L.size;
1764  L_1d.AllocateMemory(L.size);
1765  fL.AllocateMemory(L.size);
1766  for (il = 0; il < L.size; il++) L_1d[il] = L.arr[il][0][ia_90];
1767 
1768  }
1769 
1770  steady_state(fL, tau, Kp, Lsize_7, L_1d, fb_out, fb_in);
1771  //scale_factor = fb_out/fL[L.size-1]; // scale factor to make psd to be 1 at the outer boundary (=Lmax(<7))
1772  scale_factor = 1; // Bugfix. In terms of usage bf.txt as Flux/MeanFlux PSD should be 1 at L==7
1773 
1774  // If max(L)>7, normalize, so the steady state fL(L=7) = 1;
1775  if (L_1d[L_1d.size_x-1] > 7) {
1776  int il_7 = L.size-1;
1777  while (L_1d[il_7] > 7 && il_7 > 0) il_7--;
1778  Output::echo(5, "Steady state normalized by fL(L[%d] = %f) = %f.\n", il_7, L_1d[il_7], fL(il_7));
1779  fL = fL / fL(il_7);
1780  }
1781  //
1782  ///////////////
1783 
1784  // if PSD is lower then min_f - make it min_f
1785  for (int i = 0; i < fL.size_x; i++) if (fL[i] < min_f) fL[i] = min_f;
1786 
1787 
1788  // inner belt simulation
1789 // int il_2=L.size-1;
1790 // while (L[il_2][0][0] > 2 && il_2 > 0) il_2--;
1791 // 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];
1792 
1793 
1794 //#ifdef DEBUG_MODE
1795  fL.writeToFile("./Debug/fL.dat");
1796 //#endif
1797 
1798  for (im = 0; im < pc.size; im++) {
1799  for (il = 0; il < L.size; il++) {
1800  for (ia = 0; ia < alpha.size; ia++) {
1801  // filling for all pc
1802  // Why do we do it for alpha=90?:
1803  // Each point is calculated at the equator and multiplied by sin(pitch angle).
1804  // BUT! The equatorial point should have the SAME ENERGY.
1805  // So, we shouldn't use pc[...][ia_90], but should just use pc[...][ia]
1806  // From the other hand, LATER we use this boundary for the ortogonal grid WITHOUD CHANGES,
1807  // so, calculation with pc[...][ia_90] is surprisingly more correct at this point.
1808  // p1 = pc[il][im][alpha.size-1] * sqrt( VF::B(7.)/VF::B(L[il][im][alpha.size-1]));
1809 
1810  // ia_90=alpha.size-1;
1811  // while (alpha.arr[il][im][ia_90] > VC::pi/2 && ia_90 > 0) ia_90--;
1812 
1813  p1 = pc.arr[il][im][ia_90] * sqrt( VF::B(7.)/VF::B(L.arr[il][im][ia_90]));
1814 
1815  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);
1816 
1817  // alpha dependence inside the loss cone
1818  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])) {
1819  //fLmp[*,i,*]=fLmp[*,i,*]*sinh(i/8.)^2/sinh(1)^2
1820  // ???????? arr[il][im][ia] = arr[il][im][ia] * pow(sinh(alpha[il][im][ia]/8.),2)/pow(sinh(1.0),2);
1821  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);
1822  }
1823  // fLm[*,im]=fL[*]*f1L7(epc[im])*(1/pc[im])^2/fL[index7]
1824  // fLm[iL,im]=fL[iL]
1825  //cout << fLm[il][im] << endl;
1826  }
1827  }
1828  }
1829 
1830 }
1831 
1832 
1833 
1834 /// two zone initial profile (kckim)
1835 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) {
1836  int il, im, ia;
1837  double Lmax, Lmin, Lgrid, Lpp;
1838  int Lsize_7, addGrid; // These two numbers should be of integer type.
1839  double tau_p, tau_c, alpha1, p1;
1840  Matrix1D<double> L_1d_i, fL_i, tau_0, Ke1, scale_factor(pc.size);
1841  Matrix2D<double> pc_i, fLm;
1842 
1843  addGrid = 0;
1844  Lsize_7 = L.size;
1845 
1846  //////////
1847  // Steady state solution for any outer-L boundary
1848  // Create 1D grid at min energy and 90 deg alpha
1849  // First, find a closest index of 90 deg pitch-angle
1850 
1851  int ia_90 = alpha.size - 1;
1852  while (alpha.arr[L.size - 1][0][ia_90] > VC::pi / 2 && ia_90 > 0) ia_90--;
1853 
1854  // If maximum L value is less than 7, we add another points from Lmax(<7) to 7
1855  // and construct pc grid at each added point by moving adiabatically pc values at Lmax(<7)
1856  if (L.arr[L.size - 1][0][ia_90] < 7.0) {
1857 
1858  // construct new L grid
1859  Lmax = L.arr[L.size - 1][0][ia_90];
1860  Lmin = L.arr[0][0][ia_90];
1861  Lgrid = (Lmax - Lmin) / (L.size - 1);
1862  Lsize_7 = (7.0 - Lmin) / Lgrid + 1;
1863  addGrid = Lsize_7 - L.size; // grid points to add
1864  L_1d_i.AllocateMemory(L.size + addGrid); // construct new L grid
1865  fL_i.AllocateMemory(L.size + addGrid); // and psd grid
1866 
1867  // construct new pc grid at each added L point
1868  pc_i.AllocateMemory(L.size + addGrid, pc.size);
1869 
1870  for (il = 0; il < L.size; il++) {
1871  for (im = 0; im < pc.size; im++) {
1872  L_1d_i[il] = L.arr[il][0][ia_90];
1873  pc_i[il][im] = pc.arr[il][im][ia_90];
1874  }
1875  }
1876 
1877  // add pc and L values at each added point
1878  for (il = L.size; il < Lsize_7; il++) {
1879  for (im = 0; im < pc.size; im++) {
1880  L_1d_i[il] = Lmin + Lgrid * il;
1881  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]));
1882  }
1883  }
1884  } else {
1885 
1886  // construct new L grid
1887  addGrid = 0; // grid points to add
1888  L_1d_i.AllocateMemory(L.size); // construct new L grid
1889  fL_i.AllocateMemory(L.size); // and psd grid
1890 
1891  // construct new pc grid at each added L point
1892  pc_i.AllocateMemory(L.size, pc.size);
1893 
1894  for (il = 0; il < L.size; il++) {
1895  for (im = 0; im < pc.size; im++) {
1896  L_1d_i[il] = L.arr[il][0][ia_90];
1897  pc_i[il][im] = pc.arr[il][im][ia_90];
1898  }
1899  }
1900 
1901  }
1902 
1903  // solve steady state equation based on the newly constructed grid
1904  Lpp = (5.6 - 0.46 * Kp);
1905 
1906  tau_0.AllocateMemory(L.size + addGrid);
1907  Ke1.AllocateMemory(L.size + addGrid);
1908  fLm.AllocateMemory(L.size + addGrid, pc.size);
1909 
1910  for (im = 0; im < pc.size; im++) {
1911  for (il = 0; il < Lsize_7; il++) {
1912 
1913  Ke1[il] = VF::Kfunc(pc_i[il][im]); // Energy in MeV
1914  alpha1 = alpha.arr[L.size - 1][0][alpha.size - 1]; // ~90 deg alpha
1915 
1916  if(PSD_parameters.usetauLpp_SteadyState == "coulomb") { // tau by coulomb scattering
1917  tau_0[il] = VF::Coulomb(L_1d_i[il],Ke1[il]);
1918  } else if(PSD_parameters.usetauLpp_SteadyState == "precipitation") { // tau by precipitaion
1919  tau_0[il] = VF::Precipitation(Kp,L_1d_i[il],Ke1[il]);
1920  } else if (PSD_parameters.usetauLpp_SteadyState == "combined") { // tau by combined scattering
1921  tau_c = VF::Coulomb(L_1d_i[il],Ke1[il]);
1922  tau_p = VF::Precipitation(Kp,L_1d_i[il],Ke1[il]);
1923  tau_0[il] = 1.0 / (1.0/tau_p+1.0/tau_c);
1924  } else { // user specificed tau
1925  tau_0[il] = tauLpp;
1926  }
1927 
1928  }
1929 
1930  L_1d_i.writeToFile("L_1d_i.dat");
1931  Ke1.writeToFile("Ke1.dat");
1932  tau_0.writeToFile("tau_0.dat");
1933 
1934  steady_state_two_zone(fL_i, tau_0, Kp, alpha1, Ke1, Lsize_7, L_1d_i, fb_out, fb_in);
1935 
1936  fL_i.writeToFile("fL_i.dat");
1937 
1938  for (il = 0; il < Lsize_7; il++) {
1939  if (fL_i[il] < min_f) {
1940  fLm[il][im] = min_f;
1941  } else {
1942  fLm[il][im] = fL_i[il];
1943  }
1944  }
1945  }
1946 
1947  // estimate scale factor that would be used to scale down psd at the outer boundary (=Lmax(<7))
1948  for (im = 0; im < pc.size; im++) {
1949  scale_factor[im] = fb_out / fLm[L.size - 1][im];
1950  }
1951 
1952  for (im = 0; im < pc.size; im++) {
1953  for (il = 0; il < L.size; il++) {
1954  for (ia = 0; ia < alpha.size; ia++) {
1955  // filling for all pc
1956  // Why do we do it for alpha=90?:
1957  // Each point is calculated at the equator and multiplied by sin(pitch angle).
1958  // BUT! The equatorial point should have the SAME ENERGY.
1959  // So, we shouldn't use pc[...][ia_90], but should just use pc[...][ia]
1960  // From the other hand, LATER we use this boundary for the ortogonal grid WITHOUD CHANGES,
1961  // so, calculation with pc[...][ia_90] is surprisingly more correct at this point.
1962  // p1 = pc[il][im][alpha.size-1] * sqrt( VF::B(7.)/VF::B(L[il][im][alpha.size-1]));
1963  int ia_90 = alpha.size - 1;
1964  while (alpha.arr[il][im][ia_90] > VC::pi / 2 && ia_90 > 0)
1965  ia_90--;
1966 
1967  p1 = pc.arr[il][im][ia_90] * sqrt(VF::B(7.) / VF::B(
1968  L.arr[il][im][ia_90]));
1969  arr[il][im][ia] = sin(alpha.arr[il][im][ia]) * scale_factor[im]
1970  * fLm[il][im] * VF::Outer_Boundary_Flux_at_L7(
1971  VF::Kfunc(p1), J_L7_function) * pow(1.0 / p1, 2);
1972 
1973  // alpha dependence inside the loss cone
1974  if (alpha.arr[il][im][ia] < VF::alc(L.arr[il][im][ia])
1975  || alpha.arr[il][im][ia] > VC::pi - VF::alc(
1976  L.arr[il][im][ia])) {
1977  arr[il][im][ia] = arr[il][im][ia] * pow(
1978  sinh(alpha.arr[il][im][ia] / VF::alc(
1979  L.arr[il][im][ia])), 2) / pow(sinh(1.0), 2);
1980  }
1981  }
1982  }
1983  }
1984 
1985 #ifdef DEBUG_MODE
1986  arr.writeToFile("./Debug/f_0.plt", L, pc, alpha);
1987 #endif
1988 
1989 }
1990 
1991 
1992 
1993 
1994 
1995 /**
1996 * Calculate steady state.
1997 *
1998 * \param &f - function
1999 * \param tau - life time
2000 * \param Kp - Kp index value
2001 * \param nx - number of points
2002 * \param &CL - grid
2003 */
2004 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
2005 
2006  //Matrix1D<double> f(nx);
2007  //Matrix1D<double> L(nx);
2008  Matrix1D<double> C_1(nx);
2009  Matrix1D<double> C_2(nx);
2010  Matrix1D<double> C_3(nx);
2011  Matrix1D<double> alfa(nx);
2012  Matrix1D<double> A(nx);
2013  Matrix1D<double> B(nx);
2014  Matrix1D<double> C(nx);
2015  Matrix1D<double> D(nx);
2016  Matrix1D<double> Dxx(nx);
2017 
2018  int i;
2019 
2020  for (i = 0; i <= nx-1; i++) {
2021  //Dxx[i] = VF::Df(L[i], Kp) + VF::Dfe(alpha,Ke[i],L[i],Kp);
2022  Dxx[i] = pow(10, 0.506*Kp-9.325)*pow(L[i],10) + VF::Dfe(alpha,Ke[i],L[i],Kp);
2023  }
2024 
2025  Dxx.writeToFile("Dxx.dat");
2026 
2027  for (i = 1; i < nx-1; i++) {
2028  alfa[i] = 1.0/((L[i+1] - L[i-1])/2)*pow(L[i],2);
2029  //C_1[i] = (Dxx[i] + Dxx[i-1])/(pow(((L[i] + L[i-1])/2), 2));
2030  C_1[i] = ((Dxx[i] + Dxx[i-1])/2)/(L[i] - L[i-1])/(pow(((L[i] + L[i-1])/2), 2));
2031  C_2[i] = ((Dxx[i+1] + Dxx[i])/2)/(L[i+1] - L[i])/(pow(((L[i+1] + L[i])/2), 2));
2032  C_3[i] = 1./tau[i];
2033  //alfa[i] = pow(L[i],2)/(2.*pow((L[i+1] - L[i-1])/2, 2));
2034  }
2035 
2036  for (i = 1; i < nx-1; i++) {
2037  A[i] = alfa[i] * C_1[i];
2038  B[i] = -alfa[i] * (C_1[i] + C_2[i]) - C_3[i];
2039  C[i] = alfa[i] * C_2[i];
2040  }
2041 
2042  for (i = 0; i < nx; i++) D[i] = 0;;
2043 
2044  i = 0;
2045  A[i] = 0.;
2046  B[i] = 1;
2047  C[i] = 0;
2048  D[i] = f_bnd_in;
2049 
2050  i=nx-1;
2051  A[i] = 0;
2052  B[i] = 1;
2053  C[i] = 0;
2054  D[i] = f_bnd_out;
2055 
2056 #ifdef DEBUG_MODE
2057  Dxx.writeToFile("./Debug/SS_DLL.plt");
2058  A.writeToFile("./Debug/SS_A.plt");
2059  B.writeToFile("./Debug/SS_B.plt");
2060  C.writeToFile("./Debug/SS_C.plt");
2061  D.writeToFile("./Debug/SS_D.plt");
2062 #endif
2063 
2064  tridag(&A[0], &B[0], &C[0], &D[0], &f[0], nx);
2065 
2066  // for (i = 0; i < nx; i++) f[i] = f[i]/f_bnd_out;
2067 
2068 }
2069 
2070 
2071 
2072 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) {
2073  int il, im, ia;
2074  double p1;
2075  Matrix1D<double> L_1d(L.size), fL(L.size);
2076  Matrix2D<double> fLm(L.size, pc.size);
2077 
2078  // L_UpperBoundaryCondition.writeToFile("./Debug/L_UpperBoundaryCondition.dat");
2079 
2080  // steady state
2081  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
2082  steady_state(fL, tau, Kp, L.size, L_1d, fb_out, fb_in);
2083  for (int i = 0; i < fL.size_x; i++) if (fL[i] < min_f) fL[i] = min_f;
2084 
2085 #ifdef DEBUG_MODE
2086  fL.writeToFile("./Debug/fL.dat");
2087 #endif
2088 
2089  // Steady state from L-upper boundary
2090  for (im = 0; im < pc.size; im++) {
2091  for (il = 0; il < L.size; il++) {
2092  for (ia = 0; ia < alpha.size; ia++) {
2093  p1 = sqrt( VF::B(7.)/VF::B(L.arr[il][im][ia]));
2094  arr[il][im][ia] = sin(alpha.arr[il][im][ia]) * fL[il] * L_UpperBoundaryCondition[im][ia] * pow(1.0/p1,2);
2095  }
2096  }
2097  }
2098 
2099 #ifdef DEBUG_MODE
2100  arr.writeToFile("./Debug/f_0.plt", L, pc, alpha);
2101 #endif
2102 
2103 }
2104 
2105 
2106 
2107 /**
2108 * Load initial %PSD from 2d-file (for 2d calculations).
2109 *
2110 * \param &L - grid element L
2111 * \param &pc - grid element pc
2112 * \param &alpha - grid element alpha
2113 * \param *filename - file name
2114 */
2115 
2116 void PSD::Load_initial_f_2d(GridElement &L, GridElement &pc, GridElement &alpha, const char *filename) {
2117  int i, j, k;
2118  int il, im, ia;
2119  int jmax, lmax;
2120  // number of orbits
2121  lmax=8;
2122  // number of j on orbit
2123  jmax=16;
2124 
2125  Matrix2D<double> jperp(lmax+1, jmax+1);
2126  Matrix1D<double> energy(jmax+2), jcopy(jmax+2);
2127  Matrix1D<double> j_0m(pc.size), f_0m(pc.size);
2128  Matrix1D<double> dayno(lmax+1), eventtime(lmax+1), orbit(lmax+1);
2129 
2130  // Load j from input file
2131  ifstream input_f(filename);
2132  if (input_f != NULL && !input_f.eof()) {
2133  for (k = 0; k <= lmax; k++) {
2134  input_f >> orbit[k] >> dayno[k] >> eventtime[k];
2135  for (j = 0; j <= jmax; j++) {
2136  input_f >> i >> energy[j] >> jperp[k][j];
2137  energy[j] = energy[j]/1000;
2138  }
2139  }
2140  } else {
2141  throw error_msg("INITIAL_PSD_ERROR", "Error reading initial conditions input file %s.\n", filename);
2142  //cout << "Error reading initial conditions input file " << filename << endl;
2143  //getchar();
2144  //exit(-1);
2145  }
2146  input_f.close();
2147 
2148  // Copying j into 1D matrix_array for interpolation
2149  for (i = 0; i <= jmax; i++) jcopy[i] = jperp[2][i];
2150 
2151  // Adding upper border point and border condition
2152  energy[jmax+1] = VF::Kfunc(pc.GridElement_parameters.max);
2153  jcopy[jmax+1] = VC::zero_f;
2154 
2155 #ifdef DEBUG_MODE
2156  jcopy.writeToFile("./Debug/jcopy.plt", energy);
2157 #endif
2158 
2159  if (VF::Kfunc(pc.GridElement_parameters.min) < energy[0]) printf("Warning: Emin less than minimum data in input file!\n");
2160 
2161  // Interpolating(old grid, new grid);
2162  Matrix1D<double> epc1D(pc.size); // temporary 1D energy-grid
2163  for (i = 0; i < pc.size; i++)
2164  epc1D[i] = VF::Kfunc(pc.arr[0][i][0]);
2165  // j_0m.Interpolate(jcopy, energy, epc1D);
2166 
2167  for (i = 0; i < jmax+2; i++) {
2168  jcopy[i] = log10(jcopy[i]);
2169  }
2170  j_0m.Spline(jcopy, energy, epc1D, jcopy[0], log10(VC::zero_f), 0, 1);
2171  for (i = 0; i < pc.size; i++) {
2172  j_0m[i] = pow(10, j_0m[i]);
2173  }
2174 
2175  // Calculating from j_0m to f_0m.
2176  for (i = 0; i < pc.size; i++) {
2177  if (j_0m[i] < 0) j_0m[i] = VC::zero_f;
2178  f_0m[i] = j_0m[i]/pow(pc.arr[0][i][0], 2);
2179  }
2180 
2181 #ifdef DEBUG_MODE
2182  j_0m.writeToFile("./Debug/j_0m.plt", epc1D);
2183  f_0m.writeToFile("./Debug/f_0m.plt", epc1D);
2184 #endif
2185 
2186  for (il = 0; il < L.size; il++) {
2187  for (im = 0; im < pc.size; im++) {
2188  for (ia = 0; ia < alpha.size; ia++) {
2189  arr[il][im][ia] = f_0m[im] * sin( alpha.arr[il][im][ia] );
2190  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])) {
2191  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);
2192  }
2193  }
2194  }
2195  }
2196 }
2197 
2198 
2199 
2200 /**
2201 * Load initial %PSD from 2d-file (for 2d calculations).
2202 *
2203 * \param &L - grid element L
2204 * \param &pc - grid element pc
2205 * \param &alpha - grid element alpha
2206 * \param *filename - file name
2207 */
2208 
2210  int il, im, ia;
2211  double Ke, alf, flux;
2212 
2213  for (il = 0; il < L.size; il++) {
2214  for (im = 0; im < pc.size; im++) {
2215  for (ia = 0; ia < alpha.size; ia++) {
2216 
2217  Ke = VF::Kfunc(pc.arr[0][im][0]); // Energy in MeV
2218  alf = alpha.arr[il][im][ia]; // pitch angle in radian
2219  flux = exp(-(Ke-0.2)/0.1); //* ( sin(alf)-sin(VF::alc(L[il][im][ia])) );
2220  //flux = exp(-(Ke-0.2)/0.1)* ( sin(alf)-sin(VF::alc(L[il][im][ia])) ); */
2221  arr[il][im][ia] = flux / pow(pc.arr[0][im][0], 2);
2222 
2223  }
2224  }
2225  }
2226 }
2227 
2228 
2229 
2230 /**
2231 * Calculate steady state.
2232 *
2233 * \param &f - function
2234 * \param tau - life time
2235 * \param Kp - Kp index value
2236 * \param nx - number of points
2237 * \param &L - grid
2238 * \param f_bnd_out - the value on the outer boundary
2239 * \param f_bnd_in - the value on the inner boundary
2240 */
2241 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
2242 
2243  //Matrix1D<double> f(nx);
2244  //Matrix1D<double> L(nx);
2245  Matrix1D<double> C_1(nx);
2246  Matrix1D<double> C_2(nx);
2247  Matrix1D<double> C_3(nx);
2248  Matrix1D<double> alfa(nx);
2249  Matrix1D<double> A(nx);
2250  Matrix1D<double> B(nx);
2251  Matrix1D<double> C(nx);
2252  Matrix1D<double> D(nx);
2253  Matrix1D<double> Dxx(nx);
2254 
2255  int i;
2256 
2257  for (i = 0; i <= nx-1; i++) {
2258  // Radial diffusion coefficient approximation by Brautigam and Albert [2000]
2259  Dxx[i] = pow(10, 0.506*Kp-9.325)*pow(L[i],10);
2260  }
2261 
2262  for (i = 1; i < nx-1; i++) {
2263  alfa[i] = 1.0/((L[i+1] - L[i-1])/2)*pow(L[i],2);
2264  //C_1[i] = (Dxx[i] + Dxx[i-1])/(pow(((L[i] + L[i-1])/2), 2));
2265  C_1[i] = ((Dxx[i] + Dxx[i-1])/2)/(L[i] - L[i-1])/(pow(((L[i] + L[i-1])/2), 2));
2266  C_2[i] = ((Dxx[i+1] + Dxx[i])/2)/(L[i+1] - L[i])/(pow(((L[i+1] + L[i])/2), 2));
2267  C_3[i] = 1./tau;
2268  //alfa[i] = pow(L[i],2)/(2.*pow((L[i+1] - L[i-1])/2, 2));
2269  }
2270 
2271  for (i = 1; i < nx-1; i++) {
2272  A[i] = alfa[i] * C_1[i];
2273  B[i] = -alfa[i] * (C_1[i] + C_2[i]) - C_3[i];
2274  C[i] = alfa[i] * C_2[i];
2275  }
2276 
2277  for (i = 0; i < nx; i++) D[i] = 0;;
2278 
2279  i = 0;
2280  A[i] = 0.;
2281  B[i] = 1;
2282  C[i] = 0;
2283  D[i] = f_bnd_in;
2284 
2285  i=nx-1;
2286  A[i] = 0;
2287  B[i] = 1;
2288  C[i] = 0;
2289  D[i] = f_bnd_out;
2290 
2291 #ifdef DEBUG_MODE
2292  Dxx.writeToFile("./Debug/SS_DLL.plt");
2293  A.writeToFile("./Debug/SS_A.plt");
2294  B.writeToFile("./Debug/SS_B.plt");
2295  C.writeToFile("./Debug/SS_C.plt");
2296  D.writeToFile("./Debug/SS_D.plt");
2297 #endif
2298 
2299  tridag(&A[0], &B[0], &C[0], &D[0], &f[0], nx);
2300 
2301  // for (i = 0; i < nx; i++) f[i] = f[i]/f_bnd_out;
2302 
2303 }