VERB_code_2.2  2
 All Classes Namespaces Files Functions Variables Typedefs Macros Pages
Main.cpp
Go to the documentation of this file.
1 /**
2  * \mainpage
3  *
4  * The document <a href="../Algorithms.pdf">Algorithms.pdf</a> may be useful for understanding the structure of the code.
5  *
6  * This document generated by Doxygen system. It has a description of API, used my main program code which run main() function.
7  * If you whant to understand how the VERB code works, please, read the owerview and after read the code and coments itself.
8  *
9  * # Overview.
10  *
11  * The code is designed to solve the Fokker-Planck diffusion equation: \f[
12  \frac{\partial f}{\partial t} = \frac{1}{p^2} \frac{\partial}{\partial p} p^2 D_{pp} \frac{\partial f}{\partial p}\Big|_{\textrm{const }\alpha, L} +
13 + \frac{1}{\cos(\alpha) \sin(alpha) T(\alpha)} \frac{\partial}{\partial \alpha} \cos(\alpha) \sin(\alpha) T(\alpha) D_{\alpha \alpha} \frac{\partial f}{\partial \alpha}\Big|_{\textrm{const } p, L} +
14 + L^2 \frac{\partial}{\partial L} \frac{1}{L^2} D_{LL} \frac{\partial f}{\partial L}\Big|_{\textrm{const } J_1, J_2} + Q - S
15  * \f],
16  * describs the evolution of phase space density (PSD) in time and phase space.
17  * It has three diffusion terms: _Radial_, _Energy_, and _Pitch-angle_. _Radial_ diffusion needs different grid, than two other (local, _Energy_, and _Pitch-angle_) diffusions. Q - local sources. S - local losses.
18  * In the code we split the equation into \b three parts and solve them one after anouther with interpolation between radial and local grids:
19  * - \e Radial diffusion
20  * - Interpolation from radial grid to local grid.
21  * - \e Energy diffusion.
22  * - \e Pitch-angle diffusion.
23  * - Interpolation from local to radial grid.
24  *
25  * We have \b two grids in the code and phase space density arrays on each of them. %PSD values and grids are stored in the class PSD. The structure of the classes is following:
26  *
27  * - Matrix3D - basic class describe matrix type of variable to manipulate them. This class has to be considered to understand the logic of the code.
28  * - PSD - class, holds phase space density values and does diffusions.
29  * - Grid - grid, in which all that diffusions are (This class stored grid variables). Operate with variables of GridElement class.
30  * - GridElement. In the code you may find next variables:
31  * - L( radial ditance),
32  * - pc ( momentum * speed of light in units of MeV)
33  * - &alpha; (pitch angle in rad)
34  * - epc( energy in MeV )
35  * - three GridElement%s (since we have a 3D grid), and forth one, epc, is a function of pc and is just useful enough to keep it.
36  * - BoundaryCondition - upper, lower - upper and lower boundary conditions for diffusion on each GridElement%s
37  *
38  * PSD as general class, contains methods to compute diffusion. It keeps PSD values in the parent class Matrix3D.
39  * Grid has 3-dimensions, so it has 3 GridElement%s, and one additional element epc which is just useful for the output.
40  * Each GridElement%s is a 3D-Array of values of coordinates. We have to do it this way, because grids are irregular.
41  * BoundaryCondition class is used for the boundary conditions for diffusion calculation.
42  *
43  * PSD needs diffusion coefficients, grid, and boundary conditions as an input parameter for diffusion functions PSD::Diffusion_L(), PSD::Diffusion_pc(), PSD::Diffusion_alpha(), PSD::Diffusion_pc_alpha().
44  * Diffusion coefficients for each diffusion are stored inside DiffusionCoefficientsGroup class. The class structure:
45  *
46  * - DiffusionCoefficientsGroup Daa, Dpcpc etc. - summary of diffusion coefficients for one direction. For example, when we have different waves and different Daa coefficients respectively we have to combine them indo one diffusion coefficient for the calculation.
47  * - DiffusionCoefficient *_chorus, *_hiss, *_EMIC etc - diffusion coefficients by each type of wave separetly. The coefficients load from file or calculate and stored in array of DiffusionCoefficient elements which combine into DiffusionCoefficientsGroup.
48 
49  * DiffusionCoefficient%s, produced by each wave type, grouped by directions in the DiffusionCoefficientsGroup class. That class can turn the waves on and off and scale their effect,
50  * and store the summation as Daa or Dpp.
51  *
52  * Initial conditions defines by PSD::Create_Initial_PSD.
53  *
54  * ## Classes constructions.
55  * The code has empty constructor and functions, like `AllocateMemory()` and `Initialize()` which define an object. We also have a constructor, which call these functions during execution,
56  * just for shortening of the code. And it is like that for almost each class.
57  *
58  * ## Abbreviations
59  * - PSD - Phase space density.
60  * - RCM - Rice convection model (http://rocco.rice.edu/~toffo/research/rcm/rcm.html, electron seed population)
61  * - GMRES - Generalized minimal residual method. For matrix inversion.
62  * - LAPACK - Linear Algebra PACKage. For linear algebra methods.
63  *
64  * \todo
65  * - make a special page for mail.cpp to describe algoritms and sections
66  * - specify page on link for parameters
67  * - why we need load boundary in two different sections? LoadBoundaryCondition | MakeBoundaryCondition
68  * - is MixTerms in a AM_Block method? - Yes
69  * - why psdRadialDiffusion - dosn't check for === 0?
70  * - why local->radial interpolation after write the Output? -(for debug)
71  * - how to define LAPACK in Doxygen?
72  * - * - what to exclude from Doxigen about VERB code??
73  * - how to exclude "Constant groups" from namespace defenition
74  */
75 
76 /**
77 * \file Main.cpp
78 * \author Developed by Dmitry Subbotin under supervision of the PI Yuri Shprits.
79 * \brief This is the main program file.
80 *
81 * \details Code below correspond to the commented part. It works only in uncommented state.
82 *
83 * \code
84 * // \def DEBUG
85 * // \brief Reserved for debug purposes if defined
86 * // #define DEBUG
87 *
88 * // \def DEBUG_MODE
89 * // \brief define DEBUG_MODE value. Checks if memory was alocated for a matrix at each use.
90 * // \warning Additional checks and outputs for debbugging if defined. MUCH SLOWER, FOR DEBUG ONLY!
91 * // #define DEBUG_MODE
92 *
93 * \endcode
94 */
95 
96 /**
97 * \def _CRT_SECURE_NO_DEPRECATE
98 * \brief No "depreciated" warnings
99 *
100 * \def __VERB_VERSION_NUMBER__
101 * \brief Number of current VERB code version
102 */
103 
104 #define _CRT_SECURE_NO_DEPRECATE
105 #define __VERB_VERSION_NUMBER__ "2.03.000"
106 
107 
108 #include <valarray>
109 
110 // for usind specific matematical functions and constants
111 #include <math.h>
112 
113 // for using time in the code
114 #include <ctime>
115 #include <time.h>
116 
117 // class to make output to the log file and to the screen
118 #include "../Logging/Output.h"
119 
120 // various functions (converting, calculation of magnetic field, etc)
121 #include "../VariousFunctions/variousFunctions.h"
122 
123 // class opperatind matrix of phase space density. Also, do diffusion
124 #include "../Diffusion/PSD.h"
125 
126 // class initialized the calculation grid and allow to use it
127 #include "../Grid/Grid.h"
128 
129 // class load and give as acces to boundary conditions
130 #include "../Grid/BoundaryConditions.h"
131 
132 // class to operate by sources and losses processes
133 #include "../Grid/AdditionalSourcesAndLosses.h"
134 
135 // class load and give an access to parameters from parameter.ini file
136 #include "../Parameters/Parameters.h"
137 
138 
139 // Memory leaks ???
140 #if defined(_WIN32) || defined(_WIN64)
141 #define _CRTDBG_MAP_ALLOC
142 #include <stdlib.h>
143 #include <crtdbg.h>
144 
145 /** \code
146 * //* #ifdef _DEBUG
147 * //* #define DEBUG_NEW new(_NORMAL_BLOCK, __FILE__, __LINE__)
148 * //* #define new DEBUG_NEW
149 * //* #endif
150 * \endcode
151 */
152 
153 #endif
154 
155 /**
156 * \namespace std
157 * \brief General namespace
158 */
159 using namespace std;
160 
161 
162 
163 /// \brief Parameters structure, with all parameters from the parameters.ini file. The default parameters defined in the constructor
165 
166 /**
167 * \brief Main code.
168  * \param[in] argv[] The name of the parameter file
169  * \param argc Number of arguments, mandatory automatic parameter (does not need to be specified).
170  */
171 int main(int argc, char* argv[]) {
172 
173  //* \deprecated:
174  //* memory leaks
175  //* _CrtSetDbgFlag ( _CRTDBG_ALLOC_MEM_DF | _CRTDBG_LEAK_CHECK_DF );
176  //* _crtBreakAlloc = 24800; // use the number from memory leack debug output to find the leack, if any exists
177 
178  // \deprecated:
179  // try-catch block for exception handeling. Any exeption from the code are handeled at the end of the block
180  try {
181 
182  // extern - once again we define the variable parameters (see above)
183  // The default parameters will be defined and the constructor of this class
185 
186  // code uses 2 grids for calculations Radial and Local diffusion
187  Grid
188  radialDiffusionGrid, localDiffusionsGrid;
189 
190  // for each calculations code use dufferent PSD objects ofr Radial and Local
191  PSD
192  psdRadialDiffusion, psdLocalDiffusions;
193 
194  // This is a basic matrix used for check purposes calculation
195  // (to check the solution at each time step_
196  Matrix3D<double> psdPrevious;
197 
198  // Each boundary condition used special class object
200  L_LowerBoundaryCondition, L_UpperBoundaryCondition,
201  pc_LowerBoundaryCondition, pc_UpperBoundaryCondition,
202  alpha_LowerBoundaryCondition, alpha_UpperBoundaryCondition;
203 
204  // Diffusion coefficient for Radial diffusion
206 
207  // Diffusion coefficients for Local diffusion
209  Dpcpc, DpcpcLpp,
210  Daa, DaaLpp,
211  Dpca, DpcaLpp;
212 
213  // Extra losses which could be indluded in the simulation (t/tau)
215 
216  // Global counter for simulation ???
217  double simulation_time = 0;
218 
219  //----------------------------------------------------------------------------------------------------
220  // Setting up and loading parameters
221  //----------------------------------------------------------------------------------------------------
222 
223  //=== timer, to measure wall clock time ===
224  // This time required for check expiration time for using the code.
225  std::clock_t zero_time, start_time, current_time;
226  zero_time = clock();
227  time_t now = time(NULL);
228 
229  printf("%d\n", localtime(&now)->tm_year);
230  if (localtime(&now)->tm_year + 1900 >= 2014) {
231  cout << "The program version is outdated, please update." << endl;
232  return 0;
233  }
234  //=== end of timer section ===
235 
236  // Initialize log-file, save log-file name for future use
237  // open log file from [namespace Output defined in Output.h]
238  Output::open_log_file("logfile.log");
239 
240  time(&now); // current time into variable now
241  Output::echo(0, "Version: %s\n", __VERB_VERSION_NUMBER__);
242  Output::echo(0, "Compiled time: %s %s\n", __DATE__, __TIME__); // Predefined compiler macros used to save compilation time in the file - very useful to figure out later which version of the code was used
243  Output::echo(0, "Starting at %s\n", ctime(&now) ); // ctime - return char ofr time
244 
245  // if we have more than one argument in calling of main program (that means we have some arguments)
246  if (argc > 1) {
247  // in this case use first argument as the name of file for loading different variables and parameters
248  // Using procedure (function) "Load_parameters", that is described in file "Parameters.cpp" for loading parameters from file
249  parameters.Load_parameters(argv[1]);
250  } else {
251  // in other case try using Parameters.ini as the name of file with variables and constants to load
252  parameters.Load_parameters("Parameters.ini");
253  }
254 
255  // Output debug info or not.
256  // parameters.outputLvl - contain a level of debug information
257  Output::set_output_lvl(parameters.outputLvl);
258 
259  // output1D - file for writing 1d stuff (like Kp, Lpp etc). Open new stream for this filename
260  ofstream output1D((parameters.general_Output_parameters.folderName + parameters.general_Output_parameters.fileName1D).c_str());
261 
262  //----------------------------------------------------------------------------------------------------
263  // Creating grids. Parameters should be loaded from ini file
264  //----------------------------------------------------------------------------------------------------
265 
266  // Parameters.ini define the Radial grid
267  Output::echo(1, "Creating radial diffusion gird... ");
268  radialDiffusionGrid.Create_Grid(parameters.radialDiffusionGrid_L,
269  parameters.radialDiffusionGrid_pc,
270  parameters.radialDiffusionGrid_alpha,
271  parameters.radialDiffusionGrid_epc,
272  parameters.radialDiffusionGrid_filename,
273  parameters.radialDiffusionGrid_type);
274 
275  Output::echo(1, "done.\n");
276 
277  // Parameters.ini define the Local grid
278  Output::echo(1, "Creating local diffusions gird... ");
279  localDiffusionsGrid.Create_Grid(parameters.localDiffusionsGrid_L,
280  parameters.localDiffusionsGrid_pc,
281  parameters.localDiffusionsGrid_alpha,
282  parameters.localDiffusionsGrid_epc,
283  parameters.localDiffusionsGrid_filename,
284  parameters.localDiffusionsGrid_type, radialDiffusionGrid);
285  Output::echo(1, "done!\n");
286 
287 
288  // Output grids
289  // After grids defined they write down to file.
290  // See Grid::Output
291  localDiffusionsGrid.Output((parameters.general_Output_parameters.folderName + "perp_grid.plt").c_str());
292  radialDiffusionGrid.Output((parameters.general_Output_parameters.folderName + "rad_grid.plt").c_str());
293 
294 
295  //----------------------------------------------------------------------------------------------------
296  // (LoadBoundaryCondition) Setting up boundary conditions that should be loaded from file
297  //----------------------------------------------------------------------------------------------------
298 
299  // Loading the boundaries that should be load from a file
300  // Check the pareameters, if the type of boundary is defined we can load the boundary.
301  // Else the default value ofr boudary type is BCT_CONSTANT_PSD
302  // The boundary defined only on the slice of the grid See Matrix3D::xSlice,ySlice
303  if (parameters.L_UpperBoundaryCondition.type == "BCT_FILE" || parameters.L_UpperBoundaryCondition.type == "BCT_FILE_GRID" || parameters.L_UpperBoundaryCondition.type == "BCT_MULTIPLE_FILES") {
304  L_UpperBoundaryCondition.LoadBoundaryCondition(
305  parameters.L_UpperBoundaryCondition,
306  radialDiffusionGrid.epc.arr.xSlice(radialDiffusionGrid.L.size-1),
307  radialDiffusionGrid.alpha.arr.xSlice(radialDiffusionGrid.L.size-1));
308  }
309  if (parameters.pc_LowerBoundaryCondition.type == "BCT_FILE" || parameters.pc_LowerBoundaryCondition.type == "BCT_FILE_GRID" || parameters.pc_LowerBoundaryCondition.type == "BCT_MULTIPLE_FILES") {
310  pc_LowerBoundaryCondition.LoadBoundaryCondition(
311  parameters.pc_LowerBoundaryCondition,
312  localDiffusionsGrid.L.arr.ySlice(0),
313  localDiffusionsGrid.alpha.arr.ySlice(0));
314  }
315 
316  //----------------------------------------------------------------------------------------------------
317  // (Radial) Initializing Phase Space Densities AND loading initial values
318  //----------------------------------------------------------------------------------------------------
319 
320  // see PSD::Create_Initial_PSD for details
321  Output::echo(1, "Radial diffusion initialization... ");
322  psdRadialDiffusion.Create_Initial_PSD(parameters.psdRadialDiffusion, radialDiffusionGrid, L_UpperBoundaryCondition);
323  Output::echo(1, "done!\n");
324 
325  //----------------------------------------------------------------------------------------------------
326  // (MakeBoundaryCondition) Setting up boundary conditions that should be set up using PSD values
327  //----------------------------------------------------------------------------------------------------
328 
329 
330  // Making boundary conditions. Same as before
331  // Parameters are PSD and grid values at the boundary, mad by 'Slice' procedure
332  Output::echo(1, "Making boundary conditions... ");
333  L_LowerBoundaryCondition.MakeBoundaryCondition(
334  parameters.L_LowerBoundaryCondition,
335  psdRadialDiffusion.arr.xSlice(0),
336  radialDiffusionGrid.epc.arr.xSlice(0),
337  radialDiffusionGrid.alpha.arr.xSlice(0));
338  L_UpperBoundaryCondition.MakeBoundaryCondition(
339  parameters.L_UpperBoundaryCondition,
340  psdRadialDiffusion.arr.xSlice(radialDiffusionGrid.L.size-1),
341  radialDiffusionGrid.epc.arr.xSlice(radialDiffusionGrid.L.size-1),
342  radialDiffusionGrid.alpha.arr.xSlice(radialDiffusionGrid.L.size-1));
343 
344  // psdRadialDiffusion everywhere because we don't have local PSD yet
345  pc_LowerBoundaryCondition.MakeBoundaryCondition(
346  parameters.pc_LowerBoundaryCondition,
347  psdRadialDiffusion.arr.ySlice(0),
348  localDiffusionsGrid.L.arr.ySlice(0),
349  localDiffusionsGrid.alpha.arr.ySlice(0));
350  pc_UpperBoundaryCondition.MakeBoundaryCondition(
351  parameters.pc_UpperBoundaryCondition,
352  psdRadialDiffusion.arr.ySlice(radialDiffusionGrid.pc.size-1),
353  localDiffusionsGrid.L.arr.ySlice(radialDiffusionGrid.pc.size-1),
354  localDiffusionsGrid.alpha.arr.ySlice(radialDiffusionGrid.pc.size-1));
355 
356  alpha_LowerBoundaryCondition.MakeBoundaryCondition(
357  parameters.alpha_LowerBoundaryCondition,
358  psdRadialDiffusion.arr.zSlice(0),
359  localDiffusionsGrid.L.arr.zSlice(0),
360  localDiffusionsGrid.epc.arr.zSlice(0));
361  alpha_UpperBoundaryCondition.MakeBoundaryCondition(
362  parameters.alpha_UpperBoundaryCondition,
363  psdRadialDiffusion.arr.zSlice(radialDiffusionGrid.alpha.size-1),
364  localDiffusionsGrid.L.arr.zSlice(radialDiffusionGrid.alpha.size-1),
365  localDiffusionsGrid.epc.arr.zSlice(radialDiffusionGrid.alpha.size-1));
366  Output::echo(1, "done!\n");
367 
368  //----------------------------------------------------------------------------------------------------
369  // (Local) Initializing Phase Space Densities AND loading initial values
370  //----------------------------------------------------------------------------------------------------
371 
372  Output::echo(1, "Local diffusions initialization... ");
373  // Create local diffusion PSD
374  // see PSD::Create_Initial_PSD for details
375  psdLocalDiffusions.Create_Initial_PSD(parameters.psdLocalDiffusions, localDiffusionsGrid, L_UpperBoundaryCondition);
376  // Radial to Local grids initial PSD interpolation. It is nessesarry to have same PSD on different grids
377  psdLocalDiffusions.Interpolate(psdRadialDiffusion, parameters.interpolation, radialDiffusionGrid, localDiffusionsGrid, pc_LowerBoundaryCondition.arr, pc_UpperBoundaryCondition.arr);
378  Output::echo(1, "done!\n");
379 
380  //----------------------------------------------------------------------------------------------------
381  // Initializing PSD complite
382  //----------------------------------------------------------------------------------------------------
383 
384  // After grids defined we can work with Matrix3d::psdPrevious required for check calculation
385  psdPrevious.AllocateMemory(localDiffusionsGrid.L.size, localDiffusionsGrid.pc.size, localDiffusionsGrid.alpha.size);
386 
387  // First data writing for output files.
388  // output initial conditions into the file as "picture" at the 0-step
389  psdLocalDiffusions.Output_without_grid(0);
390  psdRadialDiffusion.Output_without_grid(0);
391 
392  // Create any additional sources and losses (load from file or something).
393  // Right now it's empty, but is useful for future development of the code.
394  SL.Create_SL(parameters.SL, localDiffusionsGrid, parameters.totalIterationsNumber, parameters.timeStep);
395 
396  //----------------------------------------------------------------------------------------------------
397  // Creating diffusion coefficients
398  //----------------------------------------------------------------------------------------------------
399 
400  // see DiffusionCoefficient.h
401  CreateAllDiffusionCoefficients(DLL, Daa, Dpcpc, Dpca, DaaLpp, DpcpcLpp, DpcaLpp, parameters, radialDiffusionGrid, localDiffusionsGrid);
402 
403  //----------------------------------------------------------------------------------------------------
404  // Writing header into the 1d output file
405  //----------------------------------------------------------------------------------------------------
406  // Writing headers for all diffusion coefficients scaling coefficients
407  Output1DHeaders(output1D, Daa, Dpcpc, Dpca, DaaLpp, DpcpcLpp, DpcaLpp);
408  // write information about 1d variables into the file with 1d variables
409  Output1DValues(output1D, Daa, Dpcpc, Dpca, DaaLpp, DpcpcLpp, DpcaLpp, simulation_time, parameters, 0);
410 
411  //----------------------------------------------------------------------------------------------------
412  //
413  // Main loop. Making loop from 1 to total number of iterations
414  //
415  //----------------------------------------------------------------------------------------------------
416 
417  // You don't believe, but there are counters for loops
418  // il - for L, im - fro pc, ia - for alpha
419  int il, im, ia;
420 
421  Output::echo(0, "Starting...\n");
422  // totalIterationsNumber - calculated while parameters loaded
423  for (int iteration = 1; iteration < parameters.totalIterationsNumber; iteration++) {
424  // calculating time (in days) by multiplication time step (dt) to iteration number (it)
425  simulation_time = parameters.timeStep * iteration;
426 
427  // Output information about current calculations on screen
428  Output::echo(1, "Step %d of %d. Time = %d days, %.0f hours.\n", iteration, parameters.totalIterationsNumber, int(simulation_time), (simulation_time-int(simulation_time))*24);
429 
430  //----------------------------------------------------------------------------------------------------
431  // Radial diffusion
432  //----------------------------------------------------------------------------------------------------
433  if (radialDiffusionGrid.L.size > 1) {
434 
435  // Update boundary conditions, if needed see BoundaryCondition::Update
436  Output::echo(3, "Updating L-boundaries...\n");
437  L_LowerBoundaryCondition.Update(iteration, psdRadialDiffusion.arr.xSlice(0));
438  L_UpperBoundaryCondition.Update(iteration, psdRadialDiffusion.arr.xSlice(radialDiffusionGrid.L.size-1), simulation_time, parameters.timeStep);
439  Output::echo(3, "Done.\n");
440 
441  // Update DLL. DLL calculates on each step by defined formula or other method (see DiffusionCoefficients::MakeDLL)
442  if (parameters.useRadialDiffusion) {
443  Output::echo(3, "Updating radial diffusion coefficient...\n");
444  DLL.MakeDLL(radialDiffusionGrid.L, radialDiffusionGrid.pc, radialDiffusionGrid.alpha, parameters.Kp[iteration], parameters.DLLType);
445  Output::echo(3, "Done (Max(DLL) = %e).\n", DLL.arr.max());
446  }
447 
448  // If Debug folder exists - put calculated DLL into file
449  // Why it moved away of the if(parameters.useRadialDiffusion) ???
450  DLL.arr.writeToFile("./Debug/DLL.dat", radialDiffusionGrid.L.arr, radialDiffusionGrid.epc.arr, radialDiffusionGrid.alpha.arr);
451 
452  if (parameters.useRadialDiffusion) Output::echo(3, "Calculating radial diffusion... \n");
453 
454  //---------------------------------
455  // Radial diffusion PSD calculation
456  //---------------------------------
457  psdRadialDiffusion.Diffusion_L(parameters.timeStep, parameters.Lpp[iteration],
458  DLL,
459  radialDiffusionGrid.L, radialDiffusionGrid.pc, radialDiffusionGrid.alpha, radialDiffusionGrid.Jacobian,
460  L_LowerBoundaryCondition.arr,
461  L_UpperBoundaryCondition.arr * parameters.Bf[iteration],
462  L_LowerBoundaryCondition.calculationType,
463  L_UpperBoundaryCondition.calculationType,
464  parameters.tau[iteration],
465  parameters.tauLpp[iteration]);
466  /*
467  How to calculate BF - it is a boundary condition file
468  μ = μ(L-boundary)
469  Bf = Flux-data(L-boundary, μ) / pc(L-boundary, μ) /[ J_L7 / pc(L=7, μ) * steady_state(L-boundary) / steady_state(L=7)]
470  */
471  if (parameters.useRadialDiffusion) Output::echo(3, "Done. \n");
472 
473  // Interpolation from radial diffusion grin onto local diffusion grid
474  // radial->local
475  Output::echo(3, "Interpolation radial->local... \n");
476 
477  //* For time-dependent boundary conditions (used for RCM coupling) use next
478  //* *.arr.xSlice(iteration) create a 2D array of boundary values at the current time
479  //* psdLocalDiffusions.Interpolate(psdRadialDiffusion, parameters.interpolation, radialDiffusionGrid, localDiffusionsGrid, pc_LowerBoundaryCondition.arr.xSlice(iteration), pc_UpperBoundaryCondition.arr.xSlice(iteration));
480 
481  //* Constant in time boundary conditions (a typical simulation)
482  //* psdLocalDiffusions.Interpolate(psdRadialDiffusion, parameters.interpolation, radialDiffusionGrid, localDiffusionsGrid, pc_LowerBoundaryCondition.arr.xSlice(0), pc_UpperBoundaryCondition.arr.xSlice(0));
483 
484  // Constant in time low-energy boundary conditions, but flexible high-Energy boundary.
485  // Created to fix interpolation problems.
486  psdLocalDiffusions.Interpolate(psdRadialDiffusion, parameters.interpolation, radialDiffusionGrid, localDiffusionsGrid, pc_LowerBoundaryCondition.arr, psdRadialDiffusion.arr.ySlice(radialDiffusionGrid.pc.size-1));
487  Output::echo(3, "Done. \n");
488  } // end if for L-grid size checking > 1;
489 
490  //----------------------------------------------------------------------------------------------------
491  // Local Diffusion
492  //----------------------------------------------------------------------------------------------------
493 
494  // Update everything (analogy of radal diffusion)
495  // Update boundary conditions (see BoundaryCondition::Update)
496  if (localDiffusionsGrid.pc.size > 1) {
497  Output::echo(3, "Updating pc-boundaries... \n");
498  pc_LowerBoundaryCondition.Update(iteration, psdLocalDiffusions.arr.ySlice(0));
499  pc_UpperBoundaryCondition.Update(iteration, psdLocalDiffusions.arr.ySlice(radialDiffusionGrid.pc.size-1));
500  Output::echo(3, "Done. \n");
501  }
502  if (localDiffusionsGrid.alpha.size > 1) {
503  Output::echo(3, "Updating alpha-boundaries... \n");
504  alpha_LowerBoundaryCondition.Update(iteration, psdLocalDiffusions.arr.zSlice(0));
505  alpha_UpperBoundaryCondition.Update(iteration, psdLocalDiffusions.arr.zSlice(radialDiffusionGrid.alpha.size-1));
506  Output::echo(3, "Done. \n");
507  }
508 
509  // Update (scale) diffusion coefficients (see DiffusionCoefficientsGroup::ActivateAndScale)
510  // ??? Dpcpc.CurrentDxx.arr.max(), DpcpcLpp.CurrentDxx.arr.max() -- %%
511  if (parameters.useEnergyDiffusion) {
512  Output::echo(3, "Updating momentum diffusion coefficients...\n");
513  Dpcpc.ActivateAndScale(simulation_time, parameters.Kp[iteration]);
514  DpcpcLpp.ActivateAndScale(simulation_time, parameters.Kp[iteration]);
515  Output::echo(3, "Done (Max(Dpp) = %e / Max(DppLpp) = %e).\n", Dpcpc.CurrentDxx.arr.max(), DpcpcLpp.CurrentDxx.arr.max());
516  }
517  if (parameters.useAlphaDiffusion) {
518  Output::echo(3, "Updating pitch-angle diffusion coefficients...\n");
519  Daa.ActivateAndScale(simulation_time, parameters.Kp[iteration]);
520  DaaLpp.ActivateAndScale(simulation_time, parameters.Kp[iteration]);
521  Output::echo(3, "Done (Max(Daa) = %e / Max(DaaLpp) = %e).\n", Daa.CurrentDxx.arr.max(), DaaLpp.CurrentDxx.arr.max());
522  }
523  if (parameters.useAlphaDiffusion && parameters.useEnergyDiffusion) {
524  Output::echo(3, "Updating mixed diffusion coefficients... \n");
525  Dpca.ActivateAndScale(simulation_time, parameters.Kp[iteration]);
526  DpcaLpp.ActivateAndScale(simulation_time, parameters.Kp[iteration]);
527  Output::echo(3, "Done (Max(Dpa) = %e / Max(DpaLpp) = %e).\n", Dpca.CurrentDxx.arr.max(), DpcaLpp.CurrentDxx.arr.max());
528  }
529 
530 
531  //---------------------------------
532  // Local diffusion PSD calculation
533  //---------------------------------
534  // Local diffusions according to derivatives approximation method (block or split simply)
535  if (parameters.psdLocalDiffusions.approximationMethod == "AM_Block_C" || parameters.psdLocalDiffusions.approximationMethod == "AM_Block_LR") {
536  // "_C" and "_LR" are "central" and "left-right" approximation methods
537 
538  if (localDiffusionsGrid.pc.size > 1 && localDiffusionsGrid.alpha.size > 1) {
539  Output::echo(3, "Calculating local diffusions using BLOCK method... \n");
540  // local diffusions Momentum + Pitch angle
541  // ??? mixterms dosn't work here at all ???
542  psdLocalDiffusions.Diffusion_pc_alpha(parameters.timeStep,
543  parameters.Lpp[iteration],
544  Dpcpc.CurrentDxx, DpcpcLpp.CurrentDxx,
545  Daa.CurrentDxx, DaaLpp.CurrentDxx,
546  Dpca.CurrentDxx, DpcaLpp.CurrentDxx,
547  localDiffusionsGrid.L, localDiffusionsGrid.pc, localDiffusionsGrid.alpha, localDiffusionsGrid.Jacobian,
548  pc_LowerBoundaryCondition.arr,
549  pc_UpperBoundaryCondition.arr,
550  alpha_LowerBoundaryCondition.arr,
551  alpha_UpperBoundaryCondition.arr,
552  pc_LowerBoundaryCondition.calculationType,
553  pc_UpperBoundaryCondition.calculationType,
554  alpha_LowerBoundaryCondition.calculationType,
555  alpha_UpperBoundaryCondition.calculationType);
556  Output::echo(3, "Done. \n");
557 
558  if (iteration == 1 && parameters.outputModelMatrix) { // Id Debug folder exists and outputModelMatrix - write down the Matrix from the Local diffusion computation
559  psdLocalDiffusions.matr_A.writeToFile("./Debug/Matr_A.dat");
560  psdLocalDiffusions.matr_B.writeToFile("./Debug/Matr_B.dat");
561  psdLocalDiffusions.matr_C.writeToFile("./Debug/Matr_C.dat");
562  }
563  } // grid check
564 
565  } // approximationMethod
566  else if (parameters.psdLocalDiffusions.approximationMethod == "AM_Split_C"
567  || parameters.psdLocalDiffusions.approximationMethod == "AM_Split_LR"
568  || parameters.psdLocalDiffusions.approximationMethod == "AM_Split_C_oldway"
569  || parameters.psdLocalDiffusions.approximationMethod == "AM_Split_LR_oldway") {
570 
571  // local diffusion
572  // Momentum diffusion
573  if (parameters.useEnergyDiffusion && localDiffusionsGrid.pc.size > 1) {
574  Output::echo(3, "Calculating energy diffusions using SPLIT method... \n");
575  psdLocalDiffusions.Diffusion_pc(parameters.timeStep,
576  parameters.Lpp[iteration],
577  Dpcpc.CurrentDxx, DpcpcLpp.CurrentDxx,
578  localDiffusionsGrid.L, localDiffusionsGrid.pc, localDiffusionsGrid.alpha, localDiffusionsGrid.Jacobian,
579  pc_LowerBoundaryCondition.arr,
580  pc_UpperBoundaryCondition.arr,
581  pc_LowerBoundaryCondition.calculationType,
582  pc_UpperBoundaryCondition.calculationType);
583  Output::echo(3, "Done. \n");
584 
585  if (iteration == 1 && parameters.outputModelMatrix) { /* TBD: Model matrix output here */ }
586  } // grid check
587 
588  // local diffusion
589  // Pitch angle diffusion
590  if (parameters.useAlphaDiffusion && localDiffusionsGrid.alpha.size > 1) {
591  Output::echo(3, "Calculating pitch angle diffusions using SPLIT method... \n");
592  psdLocalDiffusions.Diffusion_alpha(parameters.timeStep,
593  parameters.Lpp[iteration],
594  Daa.CurrentDxx, DaaLpp.CurrentDxx,
595  localDiffusionsGrid.L, localDiffusionsGrid.pc, localDiffusionsGrid.alpha, localDiffusionsGrid.Jacobian,
596  alpha_LowerBoundaryCondition.arr,
597  alpha_UpperBoundaryCondition.arr,
598  alpha_LowerBoundaryCondition.calculationType,
599  alpha_UpperBoundaryCondition.calculationType);
600  Output::echo(3, "Done. \n");
601 
602  if (iteration == 1 && parameters.outputModelMatrix) { /* TBD: Model matrix output here */ }
603  } // grid check
604 
605  // Explicit mixed terms
606  if (parameters.useEnergyAlphaMixedTerms && localDiffusionsGrid.pc.size > 1 && localDiffusionsGrid.alpha.size > 1) {
607 
608  Output::echo(3, "Calculating mixed terms EXPLICITLY ... \n");
609  psdLocalDiffusions.DiffusionMixTermExplicit(parameters.timeStep,
610  parameters.Lpp[iteration],
611  Dpca.CurrentDxx, DpcaLpp.CurrentDxx,
612  localDiffusionsGrid.L, localDiffusionsGrid.pc, localDiffusionsGrid.alpha, localDiffusionsGrid.Jacobian,
613  pc_LowerBoundaryCondition.arr,
614  pc_UpperBoundaryCondition.arr,
615  alpha_LowerBoundaryCondition.arr,
616  alpha_UpperBoundaryCondition.arr,
617  pc_LowerBoundaryCondition.calculationType,
618  pc_UpperBoundaryCondition.calculationType,
619  alpha_LowerBoundaryCondition.calculationType,
620  alpha_UpperBoundaryCondition.calculationType);
621  Output::echo(3, "Done. \n");
622  } // grid check
623 
624  } // approximationMethod
625  else if (parameters.psdLocalDiffusions.approximationMethod == "AM_Kyung-Chan") { // Undocumented Kyung-Chan method
626 
627  if (localDiffusionsGrid.pc.size > 1 && localDiffusionsGrid.alpha.size > 1) {
628  Output::echo(3, "Calculating local diffusions using Kyung-Chan's approximation... \n");
629  // do local diffusions
630  psdLocalDiffusions.Diffusion_pc_alpha_KC(parameters.timeStep,
631  parameters.Lpp[iteration],
632  Dpcpc.CurrentDxx, DpcpcLpp.CurrentDxx,
633  Daa.CurrentDxx, DaaLpp.CurrentDxx,
634  Dpca.CurrentDxx, DpcaLpp.CurrentDxx,
635  localDiffusionsGrid.L, localDiffusionsGrid.pc, localDiffusionsGrid.alpha, localDiffusionsGrid.Jacobian,
636  pc_LowerBoundaryCondition.arr,
637  pc_UpperBoundaryCondition.arr,
638  alpha_LowerBoundaryCondition.arr,
639  alpha_UpperBoundaryCondition.arr,
640  pc_LowerBoundaryCondition.calculationType,
641  pc_UpperBoundaryCondition.calculationType,
642  alpha_LowerBoundaryCondition.calculationType,
643  alpha_UpperBoundaryCondition.calculationType);
644  Output::echo(3, "Done. \n");
645 
646  if (iteration == 1 && parameters.outputModelMatrix) { // Id Debug folder exists and outputModelMatrix - write down the Matrix from the Local diffusion computation
647  psdLocalDiffusions.matr_A.writeToFile("./Debug/Matr_A.dat");
648  psdLocalDiffusions.matr_B.writeToFile("./Debug/Matr_B.dat");
649  psdLocalDiffusions.matr_C.writeToFile("./Debug/Matr_C.dat");
650  }
651  } // grid check
652 
653  } else {
654  throw error_msg("APPROX_ERROR", "Approximation unknown");
655  } // end of if(approximationMethod)
656 
657 
658  //-----------------
659  // Sources & Losses
660  //-----------------
661  // Loss-Cone losses, electron lifetime losses (if any), sources from coupling with RCM, etc
662  psdLocalDiffusions.SourcesAndLosses(
663  parameters,
664  localDiffusionsGrid.L, localDiffusionsGrid.pc, localDiffusionsGrid.alpha,
665  SL,
666  parameters.timeStep,
667  parameters.Lpp[iteration],
668  parameters.tau[iteration],
669  parameters.tauLpp[iteration], parameters.Kp[iteration]);
670 
671 
672  //----------------------------------------------------------------------------------------------------
673  // Make PSD positive
674  //----------------------------------------------------------------------------------------------------
675  if (parameters.NoNegative) // This check needs for correcd strange calculation when PSD === 0
676  for (il = 0; il < localDiffusionsGrid.L.size; il++)
677  for (im = 0; im < localDiffusionsGrid.pc.size; im++)
678  for (ia = 0; ia < localDiffusionsGrid.alpha.size; ia++)
679  if (psdLocalDiffusions.arr[il][im][ia] < VC::zero_f)
680  psdLocalDiffusions.arr[il][im][ia] = VC::zero_f; // keep negative value to see where the negative values appears...
681 
682  //----------------------------------------------------------------------------------------------------
683  // Checking the answer (probably pretty meaningless, since we check the inversion of the model matrix)
684  //----------------------------------------------------------------------------------------------------
685  // ??? \deprecated ?
686  if (false) {
687  // check the local solution
688  // error = f(t+1) - f(t) - dt* d/dp Dpp d/dp f - dt*d/da Daa d/da f - dt* d/dp Dpa d/da f - dt*d/da Dap d/dp f - dt*f/tau
689  static Matrix3D<double> PSD_Der_a_L(localDiffusionsGrid.L.size, localDiffusionsGrid.pc.size, localDiffusionsGrid.alpha.size);
690  static Matrix3D<double> PSD_Der_a_R(localDiffusionsGrid.L.size, localDiffusionsGrid.pc.size, localDiffusionsGrid.alpha.size);
691  static Matrix3D<double> PSD_Der_m_L(localDiffusionsGrid.L.size, localDiffusionsGrid.pc.size, localDiffusionsGrid.alpha.size);
692  static Matrix3D<double> PSD_Der_m_R(localDiffusionsGrid.L.size, localDiffusionsGrid.pc.size, localDiffusionsGrid.alpha.size);
693 
694  static Matrix3D<double> error(localDiffusionsGrid.L.size, localDiffusionsGrid.pc.size, localDiffusionsGrid.alpha.size);
695 
696  // Some variables
697  double LaaLR, LaaRL, LmmLR, LmmRL;
698  double LmaLL, LmaLR, LmaRL, LmaRR;
699  double LamLL, LamLR, LamRL, LamRR;
700  double taulc;
701 
702  error = 0;
703  double maxerror = 0, aveerror = 0;
704 
705  for (il = 0; il < localDiffusionsGrid.L.size; il++) {
706  for (im = 0; im < localDiffusionsGrid.pc.size; im++) {
707  ia = 0;
708  PSD_Der_a_R[il][im][ia] = localDiffusionsGrid.Jacobian[il][im][ia] * Dpca.CurrentDxx.arr[il][im][ia] * (psdLocalDiffusions.arr[il][im][ia+1] - psdLocalDiffusions.arr[il][im][ia])/(localDiffusionsGrid.alpha.arr[il][im][ia+1] - localDiffusionsGrid.alpha.arr[il][im][ia]);
709  for (ia = 1; ia < localDiffusionsGrid.alpha.size-1; ia++) {
710  PSD_Der_a_L[il][im][ia] = localDiffusionsGrid.Jacobian[il][im][ia] * Dpca.CurrentDxx.arr[il][im][ia] * (psdLocalDiffusions.arr[il][im][ia] - psdLocalDiffusions.arr[il][im][ia-1])/(localDiffusionsGrid.alpha.arr[il][im][ia] - localDiffusionsGrid.alpha.arr[il][im][ia-1]);
711  PSD_Der_a_R[il][im][ia] = localDiffusionsGrid.Jacobian[il][im][ia] * Dpca.CurrentDxx.arr[il][im][ia] * (psdLocalDiffusions.arr[il][im][ia+1] - psdLocalDiffusions.arr[il][im][ia])/(localDiffusionsGrid.alpha.arr[il][im][ia+1] - localDiffusionsGrid.alpha.arr[il][im][ia]);
712  }
713  ia = localDiffusionsGrid.alpha.size-1;
714  PSD_Der_a_L[il][im][ia] = localDiffusionsGrid.Jacobian[il][im][ia] * Dpca.CurrentDxx.arr[il][im][ia] * (psdLocalDiffusions.arr[il][im][ia] - psdLocalDiffusions.arr[il][im][ia-1])/(localDiffusionsGrid.alpha.arr[il][im][ia] - localDiffusionsGrid.alpha.arr[il][im][ia-1]);
715  }
716  for (ia = 0; ia < localDiffusionsGrid.alpha.size; ia++) {
717  im = 0;
718  PSD_Der_m_R[il][im][ia] = localDiffusionsGrid.Jacobian[il][im][ia] * Dpca.CurrentDxx.arr[il][im][ia]*(psdLocalDiffusions.arr[il][im+1][ia] - psdLocalDiffusions.arr[il][im][ia])/(localDiffusionsGrid.pc.arr[il][im+1][ia] - localDiffusionsGrid.pc.arr[il][im][ia]);
719  for (im = 1; im < localDiffusionsGrid.pc.size-1; im++) {
720  PSD_Der_m_L[il][im][ia] = localDiffusionsGrid.Jacobian[il][im][ia] * Dpca.CurrentDxx.arr[il][im][ia]*(psdLocalDiffusions.arr[il][im][ia] - psdLocalDiffusions.arr[il][im-1][ia])/(localDiffusionsGrid.pc.arr[il][im][ia] - localDiffusionsGrid.pc.arr[il][im-1][ia]);
721  PSD_Der_m_R[il][im][ia] = localDiffusionsGrid.Jacobian[il][im][ia] * Dpca.CurrentDxx.arr[il][im][ia]*(psdLocalDiffusions.arr[il][im+1][ia] - psdLocalDiffusions.arr[il][im][ia])/(localDiffusionsGrid.pc.arr[il][im+1][ia] - localDiffusionsGrid.pc.arr[il][im][ia]);
722  }
723  im = localDiffusionsGrid.pc.size-1;
724  PSD_Der_m_L[il][im][ia] = localDiffusionsGrid.Jacobian[il][im][ia] * Dpca.CurrentDxx.arr[il][im][ia]*(psdLocalDiffusions.arr[il][im][ia] - psdLocalDiffusions.arr[il][im-1][ia])/(localDiffusionsGrid.pc.arr[il][im][ia] - localDiffusionsGrid.pc.arr[il][im-1][ia]);
725  }
726  }
727 
728  for (il = 0; il < localDiffusionsGrid.L.size; il++) {
729  for (im = 1; im < localDiffusionsGrid.pc.size-1; im++) {
730  for (ia = 1; ia < localDiffusionsGrid.alpha.size-1; ia++) {
731 
732  // dda D d/da f
733  LaaLR = 1.0/(localDiffusionsGrid.Jacobian[il][im][ia]) *
734  (PSD_Der_a_R[il][im][ia] - PSD_Der_a_R[il][im][ia-1])/
735  (localDiffusionsGrid.alpha.arr[il][im][ia] - localDiffusionsGrid.alpha.arr[il][im][ia-1]);
736  LaaRL = 1.0/(localDiffusionsGrid.Jacobian[il][im][ia]) *
737  (PSD_Der_a_L[il][im][ia+1] - PSD_Der_a_L[il][im][ia])/
738  (localDiffusionsGrid.alpha.arr[il][im][ia+1] - localDiffusionsGrid.alpha.arr[il][im][ia]);
739 
740  // dda D d/da f
741  LmmLR = 1.0/(localDiffusionsGrid.Jacobian[il][im][ia]) *
742  (PSD_Der_m_R[il][im][ia] - PSD_Der_m_R[il][im-1][ia])/
743  (localDiffusionsGrid.pc.arr[il][im][ia] - localDiffusionsGrid.pc.arr[il][im-1][ia]);
744  LmmRL = 1.0/(localDiffusionsGrid.Jacobian[il][im][ia]) *
745  (PSD_Der_m_L[il][im+1][ia] - PSD_Der_m_L[il][im][ia])/
746  (localDiffusionsGrid.pc.arr[il][im+1][ia] - localDiffusionsGrid.pc.arr[il][im][ia]);
747 
748  // d/dp D d/da f
749  LmaRR = 1.0/(localDiffusionsGrid.Jacobian[il][im][ia]) *
750  (PSD_Der_a_R[il][im+1][ia] - PSD_Der_a_R[il][im][ia])/
751  (localDiffusionsGrid.pc.arr[il][im+1][ia] - localDiffusionsGrid.pc.arr[il][im][ia]);
752  LmaRL = 1.0/(localDiffusionsGrid.Jacobian[il][im][ia]) *
753  (PSD_Der_a_L[il][im+1][ia] - PSD_Der_a_L[il][im][ia])/
754  (localDiffusionsGrid.pc.arr[il][im+1][ia] - localDiffusionsGrid.pc.arr[il][im][ia]);
755  LmaLR = 1.0/(localDiffusionsGrid.Jacobian[il][im][ia]) *
756  (PSD_Der_a_R[il][im][ia] - PSD_Der_a_R[il][im-1][ia])/
757  (localDiffusionsGrid.pc.arr[il][im][ia] - localDiffusionsGrid.pc.arr[il][im-1][ia]);
758  LmaLL = 1.0/(localDiffusionsGrid.Jacobian[il][im][ia]) *
759  (PSD_Der_a_L[il][im][ia] - PSD_Der_a_L[il][im-1][ia])/
760  (localDiffusionsGrid.pc.arr[il][im][ia] - localDiffusionsGrid.pc.arr[il][im-1][ia]);
761 
762 
763  // d/da D d/dp f
764  LamRR = 1.0/(localDiffusionsGrid.Jacobian[il][im][ia]) *
765  (PSD_Der_m_R[il][im][ia+1] - PSD_Der_m_R[il][im][ia])/
766  (localDiffusionsGrid.alpha.arr[il][im][ia+1] - localDiffusionsGrid.alpha.arr[il][im][ia]);
767  LamRL = 1.0/(localDiffusionsGrid.Jacobian[il][im][ia]) *
768  (PSD_Der_m_L[il][im][ia+1] - PSD_Der_m_L[il][im][ia])/
769  (localDiffusionsGrid.alpha.arr[il][im][ia+1] - localDiffusionsGrid.alpha.arr[il][im][ia]);
770  LamLR = 1.0/(localDiffusionsGrid.Jacobian[il][im][ia]) *
771  (PSD_Der_m_R[il][im][ia] - PSD_Der_m_R[il][im][ia-1])/
772  (localDiffusionsGrid.alpha.arr[il][im][ia] - localDiffusionsGrid.alpha.arr[il][im][ia-1]);
773  LamLL = 1.0/(localDiffusionsGrid.Jacobian[il][im][ia]) *
774  (PSD_Der_m_L[il][im][ia] - PSD_Der_m_L[il][im][ia-1])/
775  (localDiffusionsGrid.alpha.arr[il][im][ia] - localDiffusionsGrid.alpha.arr[il][im][ia-1]);
776 
777  error[il][im][ia] = psdLocalDiffusions.arr[il][im][ia] - psdPrevious[il][im][ia]
778  - (LaaLR + LaaRL) / 2 * parameters.timeStep // d/da D d/da f
779  - (LmmLR + LmmRL) / 2 * parameters.timeStep // d/da D d/da f
780  - (LmaRR + LmaRL + LmaLR + LmaLL) / 4 * parameters.timeStep // d/dm D d/da f
781  - (LamRR + LamRL + LamLR + LamLL) / 4 * parameters.timeStep; // d/dm D d/da f
782 
783  if (parameters.useLossCone && ((localDiffusionsGrid.alpha.arr[il][im][ia] <= VF::alc(localDiffusionsGrid.L.arr[il][im][0])) || (localDiffusionsGrid.alpha.arr[il][im][ia] >= (VC::pi - VF::alc(localDiffusionsGrid.L.arr[il][im][0])))) ) {
784  taulc = 0.25 * VF::bounce_time_new(localDiffusionsGrid.L.arr[il][im][ia], localDiffusionsGrid.pc.arr[il][im][ia], localDiffusionsGrid.alpha.arr[il][im][ia]);
785  error[il][im][ia] += psdLocalDiffusions.arr[il][im][ia]/taulc * parameters.timeStep;
786  }
787 
788  aveerror += fabs(error[il][im][ia]);
789  maxerror = (fabs(maxerror) > fabs(error[il][im][ia]))?maxerror:error[il][im][ia];
790  }
791  }
792  }
793  aveerror = aveerror / (localDiffusionsGrid.L.size) / (localDiffusionsGrid.epc.size-2) / (localDiffusionsGrid.alpha.size-2);
794  Output::echo(20, "Max calculation error = %e; ave calculation error = %e\n", maxerror, aveerror);
795 
796  stringstream filename;
797  filename << "./Debug/Calc_error_" << iteration << ".plt";
798 
799  error.writeToFile(filename.str().c_str(), localDiffusionsGrid.L.arr, localDiffusionsGrid.epc.arr, localDiffusionsGrid.alpha.arr);
800 
801  psdPrevious = psdLocalDiffusions.arr;
802  }
803 
804 
805  //----------------------------------------------------------------------------------------------------
806  // Writing output
807  //----------------------------------------------------------------------------------------------------
808  // Write PSD into a file
809  // namespace VF defined in variousFunctions.cpp
810  // If this calculation step correspond to write time - do write to output file
811  if (VF::check_time(iteration, parameters.general_Output_parameters.iterStep)) {
812  Output::echo(0, "Step %d of %d. Time = %d days, %.0f hours. ", iteration, parameters.totalIterationsNumber, int(simulation_time), (simulation_time-int(simulation_time))*24);
813  Output::echo(0, "Writing data files.\n");
814 
815  // Write information about 1d variables into the file with 1d variables
816  Output1DValues(output1D, Daa, Dpcpc, Dpca, DaaLpp, DpcpcLpp, DpcaLpp, simulation_time, parameters, iteration);
817 
818  //
819  // Most important function with write down the calculated PSD
820  //
821  psdLocalDiffusions.Output_without_grid(simulation_time);
822  if (radialDiffusionGrid.L.size > 1) psdRadialDiffusion.Output_without_grid(simulation_time);
823 
824  // 'Flushing' of log files (saving from memory to disk). Otherwise it saves it only at the end, which is not convenient.
826 
827  // Save the last PSD so it can be used as an initial condition to continue the calculation
828  psdRadialDiffusion.arr.writeToFile(parameters.general_Output_parameters.folderName + "/PSD0.dat", radialDiffusionGrid.L.arr, radialDiffusionGrid.epc.arr, radialDiffusionGrid.alpha.arr);
829  }
830 
831  if (radialDiffusionGrid.L.size > 1) {
832  // Interpolation from Local diffusions grid to Radial diffusion grid
833  // local->radial
834  Output::echo(3, "Interpolation local->radial... ");
835 
836  //* For time-dependent boundary conditions (used for RCM coupling) use next
837  //* *.arr.xSlice(iteration) create a 2D array of boundary values at the current time
838  //* psdRadialDiffusion.Interpolate(psdLocalDiffusions, parameters.interpolation, localDiffusionsGrid, radialDiffusionGrid, pc_LowerBoundaryCondition.arr.xSlice(iteration), pc_UpperBoundaryCondition.arr.xSlice(iteration));
839 
840  //* Constant in time boundary conditions (a typical simulation)
841  //* psdRadialDiffusion.Interpolate(psdLocalDiffusions, parameters.interpolation, localDiffusionsGrid, radialDiffusionGrid, pc_LowerBoundaryCondition.arr.xSlice(0), pc_UpperBoundaryCondition.arr.xSlice(0));
842 
843  // Constant in time low-energy boundary conditions, but flexible high-Energy boundary.
844  // Created to fix interpolation problems.
845  psdRadialDiffusion.Interpolate(psdLocalDiffusions, parameters.interpolation, localDiffusionsGrid, radialDiffusionGrid, pc_LowerBoundaryCondition.arr, psdLocalDiffusions.arr.ySlice(radialDiffusionGrid.pc.size-1));
846  Output::echo(3, "done!\n");
847  }
848 
849  } // end if for Main loop.
850 
851  //----------------------------------
852  // End of simulation
853  //----------------------------------
854 
855  // Last write the log
856  std::time(&now);
857  current_time = std::clock();
858  Output::echo(0, "End of simulation %s", ctime(&now));
859  Output::echo(0, "Total time %0.2f sec\n", ( current_time - zero_time ) / (double)CLOCKS_PER_SEC );
860 
861  // Close the log
863  return 0; // end program if there wasn't any exeption
864 
865 
866  //----------------------------------
867  // Exeptions
868  //----------------------------------
869  // \deprecated?
870 
871  // try block ends.
872  } catch (error_msg err) {
873  Output::echo(0, "\nCRITICAL ERROR:\n%s\n", err.what().c_str());
875  return 1;
876  } catch (exception exp) {
877  Output::echo(0, "\nCRITICAL ERROR:\n%s\n", exp.what());
879  return 1;
880  }
881 
882  //* #ifndef DEBUG_MODE
883  //* catch (...) {
884  //* Output::echo(0, "Unknown exception!");
885  //* }
886  //* #endif
887 
888  Output::echo(0, "UNEXPECTED EXIT\n");
890 
891  //* Memory leaks
892  //* _CrtDumpMemoryLeaks();
893  return 0;
894 }