VERB_code_2.3
Main.cpp
Go to the documentation of this file.
1 
123 #define _CRT_SECURE_NO_DEPRECATE
124 #define __VERB_VERSION_NUMBER__ "2.3.0.0" // Daniel Choi version
125 // presious __VERB_VERSION_NUMBER__ 2.03.000 is deprecated. New version starts with 2.2.3
126 // 2.2.0 - Ozeke Dll coefficient [version title never used]
127 // 2.2.1 - new BF [version title never used]
128 // 2.2.2 - magnetopause [version title never used]
129 // 2.2.3 - parametrization of Hiss
130 // 2.2.3.3 - New expiration date for the code
131 
132 
133 
134 #include <valarray>
135 
136 // for using specific mathematical functions and constants
137 #include <math.h>
138 
139 // for using time in the code
140 #include <ctime>
141 #include <time.h>
142 
143 // class to make output to the log file and to the screen
144 #include "../Logging/Output.h"
145 
146 // various functions (converting, calculation of magnetic field, etc)
147 #include "../VariousFunctions/variousFunctions.h"
148 
149 // class operating matrix of phase space density. Also, do diffusion
150 #include "../Diffusion/PSD.h"
151 
152 // class initializes the calculation grid and allows to use it
153 #include "../Grid/Grid.h"
154 
155 // class load and give as access to boundary conditions
156 #include "../Grid/BoundaryConditions.h"
157 
158 // class to operate by sources and losses processes
159 #include "../Grid/AdditionalSourcesAndLosses.h"
160 
161 // class load and give an access to parameters from parameter.ini file
162 #include "../Parameters/Parameters.h"
163 
164 
165 // Memory leaks ???
166 #if defined(_WIN32) || defined(_WIN64)
167 #define _CRTDBG_MAP_ALLOC
168 #include <stdlib.h>
169 #include <crtdbg.h>
170 
179 #endif
180 
181 
182 using namespace std;
183 
184 
185 
188 
195 int main(int argc, char* argv[]) {
196 
197  //* \deprecated:
198  //* memory leaks
199  //* _CrtSetDbgFlag ( _CRTDBG_ALLOC_MEM_DF | _CRTDBG_LEAK_CHECK_DF );
200  //* _crtBreakAlloc = 24800; // use the number from memory leak debug output to find the leak, if any exists
201 
202  // \deprecated:
203  // try-catch block for exception handling. Any exception from the code are handled at the end of the block
204  try {
205 
206  // extern - once again we define the variable parameters (see above)
207  // The default parameters will be defined and the constructor of this class
209 
210  // code uses 2 grids for calculations Radial and Local diffusion
211  Grid
212  radialDiffusionGrid, localDiffusionsGrid;
213 
214  // for each calculation code uses different PSD objects for Radial and Local
215  PSD
216  psdRadialDiffusion, psdLocalDiffusions;
217 
218  // This is a basic matrix used for check purposes calculation
219  // (to check the solution at each time step)
220  Matrix3D<double> psdPrevious;
221 
222  // Each boundary condition used special class object
224  L_LowerBoundaryCondition, L_UpperBoundaryCondition,
225  pc_LowerBoundaryCondition, pc_UpperBoundaryCondition,
226  alpha_LowerBoundaryCondition, alpha_UpperBoundaryCondition;
227 
228  // Diffusion coefficient for Radial diffusion
230 
231  // Diffusion coefficients for Local diffusion
233  Dpcpc, DpcpcLpp,
234  Daa, DaaLpp,
235  Dpca, DpcaLpp;
236 
237  // Extra losses which could be included in the simulation (t/tau)
239 
240  // Global counter for simulation ???
241  double simulation_time = 0;
242 
243  //----------------------------------------------------------------------------------------------------
244  // Setting up and loading parameters
245  //----------------------------------------------------------------------------------------------------
246 
247  //=== timer, to measure wall clock time ===
248  // This time required for checking expiration time for using the code.
249  std::clock_t zero_time, start_time, current_time;
250  zero_time = clock();
251  time_t now = time(NULL);
252 
253  printf("%d\n", localtime(&now)->tm_year);
254  // if current year is 2016 or later, needs to be updated
255  if (localtime(&now)->tm_year + 1900 >= 2016) {
256  cout << "The program version is outdated, please update." << endl;
257  return 0;
258  }
259  //=== end of timer section ===
260 
261  // Initialize log-file, save log-file name for future use
262  // open log file from [namespace Output defined in Output.h]
263  Output::open_log_file("logfile.log");
264 
265  time(&now); // current time into variable now
266  Output::echo(0, "Version: %s\n", __VERB_VERSION_NUMBER__);
267  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
268  Output::echo(0, "Starting at %s\n", ctime(&now) ); // ctime - return char for time
269 
270  // if we have more than one argument in calling of main program (that means we have some arguments)
271  if (argc > 1) {
272  // in this case use first argument as the name of file for loading different variables and parameters
273  // Using procedure (function) "Load_parameters", that is described in file "Parameters.cpp" for loading parameters from file
274  parameters.Load_parameters(argv[1]);
275  } else {
276  // in other case try using Parameters.ini as the name of file with variables and constants to load
277  parameters.Load_parameters("Parameters.ini");
278  }
279 
280  // Output debug info or not.
281  // parameters.outputLvl - contain a level of debug information
282  Output::set_output_lvl(parameters.outputLvl);
283 
284  // output1D - file for writing 1d stuff (like Kp, Lpp etc). Open new stream for this filename
285  ofstream output1D((parameters.general_Output_parameters.folderName + parameters.general_Output_parameters.fileName1D).c_str());
286 
287  //----------------------------------------------------------------------------------------------------
288  // Creating grids. Parameters should be loaded from ini file
289  //----------------------------------------------------------------------------------------------------
290 
301  // Parameters.ini define the Radial grid
302  Output::echo(1, "Creating radial diffusion grid... ");
303  radialDiffusionGrid.Create_Grid(parameters.radialDiffusionGrid_L,
304  parameters.radialDiffusionGrid_pc,
305  parameters.radialDiffusionGrid_alpha,
306  parameters.radialDiffusionGrid_epc,
307  parameters.radialDiffusionGrid_filename,
308  parameters.radialDiffusionGrid_type);
309 
310  Output::echo(1, "done.\n");
311 
312  // Parameters.ini define the Local grid
313  Output::echo(1, "Creating local diffusions grid... ");
314  localDiffusionsGrid.Create_Grid(parameters.localDiffusionsGrid_L,
315  parameters.localDiffusionsGrid_pc,
316  parameters.localDiffusionsGrid_alpha,
317  parameters.localDiffusionsGrid_epc,
318  parameters.localDiffusionsGrid_filename,
319  parameters.localDiffusionsGrid_type, radialDiffusionGrid);
320  Output::echo(1, "done!\n");
321 
322 
323  // Output grids
324  // After grids defined they write down to file.
325  // See Grid::Output
334  localDiffusionsGrid.Output((parameters.general_Output_parameters.folderName + "perp_grid.plt").c_str());
335  radialDiffusionGrid.Output((parameters.general_Output_parameters.folderName + "rad_grid.plt").c_str());
336 
337 
338  //----------------------------------------------------------------------------------------------------
339  // (LoadBoundaryCondition) Setting up boundary conditions that should be loaded from file
340  //----------------------------------------------------------------------------------------------------
341 
367  if (parameters.L_UpperBoundaryCondition.type == "BCT_FILE" || parameters.L_UpperBoundaryCondition.type == "BCT_FILE_GRID" || parameters.L_UpperBoundaryCondition.type == "BCT_MULTIPLE_FILES") {
368  L_UpperBoundaryCondition.LoadBoundaryCondition(
369  parameters.L_UpperBoundaryCondition,
370  radialDiffusionGrid.epc.arr.xSlice(radialDiffusionGrid.L.size-1),
371  radialDiffusionGrid.alpha.arr.xSlice(radialDiffusionGrid.L.size-1));
372  }
373  if (parameters.pc_LowerBoundaryCondition.type == "BCT_FILE" || parameters.pc_LowerBoundaryCondition.type == "BCT_FILE_GRID" || parameters.pc_LowerBoundaryCondition.type == "BCT_MULTIPLE_FILES") {
374  pc_LowerBoundaryCondition.LoadBoundaryCondition(
375  parameters.pc_LowerBoundaryCondition,
376  localDiffusionsGrid.L.arr.ySlice(0),
377  localDiffusionsGrid.alpha.arr.ySlice(0));
378  }
379 
380  //----------------------------------------------------------------------------------------------------
381  // (Radial) Initializing Phase Space Densities AND loading initial values
382  //----------------------------------------------------------------------------------------------------
383 
384  // see PSD::Create_Initial_PSD for details
385  Output::echo(1, "Radial diffusion initialization... ");
393  psdRadialDiffusion.Create_Initial_PSD(parameters.psdRadialDiffusion, radialDiffusionGrid, L_UpperBoundaryCondition);
394  Output::echo(1, "done!\n");
395 
396  //----------------------------------------------------------------------------------------------------
397  // (MakeBoundaryCondition) Setting up boundary conditions that should be set up using PSD values
398  //----------------------------------------------------------------------------------------------------
399 
400 
401  // Making boundary conditions. Same as before
402  // Parameters are PSD and grid values at the boundary, made by 'Slice' procedure
412  Output::echo(1, "Making boundary conditions... ");
413  L_LowerBoundaryCondition.MakeBoundaryCondition(
414  parameters.L_LowerBoundaryCondition,
415  psdRadialDiffusion.arr.xSlice(0),
416  radialDiffusionGrid.epc.arr.xSlice(0),
417  radialDiffusionGrid.alpha.arr.xSlice(0));
418  L_UpperBoundaryCondition.MakeBoundaryCondition(
419  parameters.L_UpperBoundaryCondition,
420  psdRadialDiffusion.arr.xSlice(radialDiffusionGrid.L.size-1),
421  radialDiffusionGrid.epc.arr.xSlice(radialDiffusionGrid.L.size-1),
422  radialDiffusionGrid.alpha.arr.xSlice(radialDiffusionGrid.L.size-1));
423 
424  // psdRadialDiffusion everywhere because we don't have local PSD yet
425  pc_LowerBoundaryCondition.MakeBoundaryCondition(
426  parameters.pc_LowerBoundaryCondition,
427  psdRadialDiffusion.arr.ySlice(0),
428  localDiffusionsGrid.L.arr.ySlice(0),
429  localDiffusionsGrid.alpha.arr.ySlice(0));
430  pc_UpperBoundaryCondition.MakeBoundaryCondition(
431  parameters.pc_UpperBoundaryCondition,
432  psdRadialDiffusion.arr.ySlice(radialDiffusionGrid.pc.size-1),
433  localDiffusionsGrid.L.arr.ySlice(radialDiffusionGrid.pc.size-1),
434  localDiffusionsGrid.alpha.arr.ySlice(radialDiffusionGrid.pc.size-1));
435 
436  alpha_LowerBoundaryCondition.MakeBoundaryCondition(
437  parameters.alpha_LowerBoundaryCondition,
438  psdRadialDiffusion.arr.zSlice(0),
439  localDiffusionsGrid.L.arr.zSlice(0),
440  localDiffusionsGrid.epc.arr.zSlice(0));
441  alpha_UpperBoundaryCondition.MakeBoundaryCondition(
442  parameters.alpha_UpperBoundaryCondition,
443  psdRadialDiffusion.arr.zSlice(radialDiffusionGrid.alpha.size-1),
444  localDiffusionsGrid.L.arr.zSlice(radialDiffusionGrid.alpha.size-1),
445  localDiffusionsGrid.epc.arr.zSlice(radialDiffusionGrid.alpha.size-1));
446  Output::echo(1, "done!\n");
447 
448  //----------------------------------------------------------------------------------------------------
449  // (Local) Initializing Phase Space Densities AND loading initial values
450  //----------------------------------------------------------------------------------------------------
451 
460  Output::echo(1, "Local diffusions initialization... ");
461  // Create local diffusion PSD
462  // see PSD::Create_Initial_PSD for details
463  psdLocalDiffusions.Create_Initial_PSD(parameters.psdLocalDiffusions, localDiffusionsGrid, L_UpperBoundaryCondition);
464  // Radial to Local grids initial PSD interpolation. It is necessary to have same PSD on different grids
465  psdLocalDiffusions.Interpolate(psdRadialDiffusion, parameters.interpolation, radialDiffusionGrid, localDiffusionsGrid, pc_LowerBoundaryCondition.arr, pc_UpperBoundaryCondition.arr);
466  Output::echo(1, "done!\n");
467 
468  //----------------------------------------------------------------------------------------------------
469  // Initializing PSD complete
470  //----------------------------------------------------------------------------------------------------
471 
472  // After grids defined we can work with Matrix3d::psdPrevious required for check calculation
473  psdPrevious.AllocateMemory(localDiffusionsGrid.L.size, localDiffusionsGrid.pc.size, localDiffusionsGrid.alpha.size);
474 
475  // First data writing for output files.
476  // output initial conditions into the file as "picture" at the 0-step
477  psdLocalDiffusions.Output_without_grid(0);
478  psdRadialDiffusion.Output_without_grid(0);
479 
480  // Create any additional sources and losses (load from file or something).
481  // Right now it's empty, but is useful for future development of the code.
482  SL.Create_SL(parameters.SL, localDiffusionsGrid, parameters.totalIterationsNumber, parameters.timeStep);
483 
484  //----------------------------------------------------------------------------------------------------
485  // Creating diffusion coefficients
486  //----------------------------------------------------------------------------------------------------
487 
495  // see DiffusionCoefficient.h
496  CreateAllDiffusionCoefficients(DLL, Daa, Dpcpc, Dpca, DaaLpp, DpcpcLpp, DpcaLpp, parameters, radialDiffusionGrid, localDiffusionsGrid);
497 
498  //----------------------------------------------------------------------------------------------------
499  // Writing header into the 1d output file
500  //----------------------------------------------------------------------------------------------------
501  // Writing headers for all diffusion coefficients scaling coefficients
502  Output1DHeaders(output1D, Daa, Dpcpc, Dpca, DaaLpp, DpcpcLpp, DpcaLpp);
503  // write information about 1d variables into the file with 1d variables
504  Output1DValues(output1D, Daa, Dpcpc, Dpca, DaaLpp, DpcpcLpp, DpcaLpp, simulation_time, parameters, 0);
505 
506  //----------------------------------------------------------------------------------------------------
507  //
508  // Main loop. Making loop from 1 to total number of iterations
509  //
510  //----------------------------------------------------------------------------------------------------
511 
512 
513  // You don't believe, but there are counters for loops
514  // il - for L, im - for pc, ia - for alpha
515  int il, im, ia;
516 
517  Output::echo(0, "Starting...\n");
518  // totalIterationsNumber - calculated while parameters loaded
519  for (int iteration = 1; iteration < parameters.totalIterationsNumber; iteration++) {
520  // calculating time (in days) by multiplication time step (dt) to iteration number (it)
521  simulation_time = parameters.timeStep * iteration;
522 
523  // Output information about current calculations on screen
524  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);
525 
526  //----------------------------------------------------------------------------------------------------
527  // Radial diffusion
528  //----------------------------------------------------------------------------------------------------
529  if (radialDiffusionGrid.L.size > 1) {
530 
531  // Update boundary conditions, if needed see BoundaryCondition::Update
540  Output::echo(3, "Updating L-boundaries...\n");
541  L_LowerBoundaryCondition.Update(iteration, psdRadialDiffusion.arr.xSlice(0));
542  L_UpperBoundaryCondition.Update(iteration, psdRadialDiffusion.arr.xSlice(radialDiffusionGrid.L.size-1), simulation_time, parameters.timeStep);
543  Output::echo(3, "Done.\n");
544 
554  // Update DLL. DLL calculates on each step by defined formula or other method (see DiffusionCoefficients::MakeDLL)
555  if (parameters.useRadialDiffusion) {
556  Output::echo(3, "Updating radial diffusion coefficient...\n");
557  DLL.MakeDLL(radialDiffusionGrid.L, radialDiffusionGrid.pc, radialDiffusionGrid.alpha, parameters.Kp[iteration], parameters.DLLType);
558  Output::echo(3, "Done (Max(DLL) = %e).\n", DLL.arr.max());
559  }
560 
561  // If Debug folder exists - put calculated DLL into file
562  // Why it moved away of the if(parameters.useRadialDiffusion) ???
563  DLL.arr.writeToFile("./Debug/DLL.dat", radialDiffusionGrid.L.arr, radialDiffusionGrid.epc.arr, radialDiffusionGrid.alpha.arr);
564 
565  if (parameters.useRadialDiffusion) Output::echo(3, "Calculating radial diffusion... \n");
566 
567  //---------------------------------
568  // Radial diffusion PSD calculation
569  //---------------------------------
570 
581  psdRadialDiffusion.Diffusion_L(parameters.timeStep, parameters.Lpp[iteration],
582  DLL,
583  radialDiffusionGrid.L, radialDiffusionGrid.pc, radialDiffusionGrid.alpha, radialDiffusionGrid.Jacobian,
584  L_LowerBoundaryCondition.arr,
585  L_UpperBoundaryCondition.arr * parameters.Bf[iteration],
586  L_LowerBoundaryCondition.calculationType,
587  L_UpperBoundaryCondition.calculationType,
588  parameters.tau[iteration],
589  parameters.tauLpp[iteration]);
598  if (parameters.useRadialDiffusion) Output::echo(3, "Done. \n");
599 
600  // Interpolation from radial diffusion grid onto local diffusion grid
601  // radial->local
602  Output::echo(3, "Interpolation radial->local... \n");
603 
604  //* For time-dependent boundary conditions (used for RCM coupling) use next
605  //* *.arr.xSlice(iteration) create a 2D array of boundary values at the current time
606  //* psdLocalDiffusions.Interpolate(psdRadialDiffusion, parameters.interpolation, radialDiffusionGrid, localDiffusionsGrid, pc_LowerBoundaryCondition.arr.xSlice(iteration), pc_UpperBoundaryCondition.arr.xSlice(iteration));
607 
608  //* Constant in time boundary conditions (a typical simulation)
609  //* psdLocalDiffusions.Interpolate(psdRadialDiffusion, parameters.interpolation, radialDiffusionGrid, localDiffusionsGrid, pc_LowerBoundaryCondition.arr.xSlice(0), pc_UpperBoundaryCondition.arr.xSlice(0));
610 
611  // Constant in time low-energy boundary conditions, but flexible high-Energy boundary.
612  // Created to fix interpolation problems.
620  psdLocalDiffusions.Interpolate(psdRadialDiffusion, parameters.interpolation, radialDiffusionGrid, localDiffusionsGrid, pc_LowerBoundaryCondition.arr, psdRadialDiffusion.arr.ySlice(radialDiffusionGrid.pc.size-1));
621  Output::echo(3, "Done. \n");
622  } // end if for L-grid size checking > 1;
623 
624  //----------------------------------------------------------------------------------------------------
625  // Local Diffusion
626  //----------------------------------------------------------------------------------------------------
627 
628  // Update everything (analogy of radial diffusion)
629  // Update boundary conditions (see BoundaryCondition::Update)
656  if (localDiffusionsGrid.pc.size > 1) {
657  Output::echo(3, "Updating pc-boundaries... \n");
658  pc_LowerBoundaryCondition.Update(iteration, psdLocalDiffusions.arr.ySlice(0));
659  pc_UpperBoundaryCondition.Update(iteration, psdLocalDiffusions.arr.ySlice(radialDiffusionGrid.pc.size-1));
660  Output::echo(3, "Done. \n");
661  }
662  if (localDiffusionsGrid.alpha.size > 1) {
663  Output::echo(3, "Updating alpha-boundaries... \n");
664  alpha_LowerBoundaryCondition.Update(iteration, psdLocalDiffusions.arr.zSlice(0));
665  alpha_UpperBoundaryCondition.Update(iteration, psdLocalDiffusions.arr.zSlice(radialDiffusionGrid.alpha.size-1));
666  Output::echo(3, "Done. \n");
667  }
668 
669  // Update (scale) diffusion coefficients (see DiffusionCoefficientsGroup::ActivateAndScale)
670  // ??? Dpcpc.CurrentDxx.arr.max(), DpcpcLpp.CurrentDxx.arr.max() -- %%
671  if (parameters.useEnergyDiffusion) {
672  Output::echo(3, "Updating momentum diffusion coefficients...\n");
673  Dpcpc.ActivateAndScale(simulation_time, parameters.Kp[iteration]);
674  DpcpcLpp.ActivateAndScale(simulation_time, parameters.Kp[iteration]);
675  Output::echo(3, "Done (Max(Dpp) = %e / Max(DppLpp) = %e).\n", Dpcpc.CurrentDxx.arr.max(), DpcpcLpp.CurrentDxx.arr.max());
676  }
677  if (parameters.useAlphaDiffusion) {
678  Output::echo(3, "Updating pitch-angle diffusion coefficients...\n");
679  Daa.ActivateAndScale(simulation_time, parameters.Kp[iteration]);
680  DaaLpp.ActivateAndScale(simulation_time, parameters.Kp[iteration]);
681  Output::echo(3, "Done (Max(Daa) = %e / Max(DaaLpp) = %e).\n", Daa.CurrentDxx.arr.max(), DaaLpp.CurrentDxx.arr.max());
682  }
683  if (parameters.useAlphaDiffusion && parameters.useEnergyDiffusion) {
684  Output::echo(3, "Updating mixed diffusion coefficients... \n");
685  Dpca.ActivateAndScale(simulation_time, parameters.Kp[iteration]);
686  DpcaLpp.ActivateAndScale(simulation_time, parameters.Kp[iteration]);
687  Output::echo(3, "Done (Max(Dpa) = %e / Max(DpaLpp) = %e).\n", Dpca.CurrentDxx.arr.max(), DpcaLpp.CurrentDxx.arr.max());
688  }
689 
690 
691  //---------------------------------
692  // Local diffusion PSD calculation
693  //---------------------------------
694  // Local diffusions according to derivatives approximation method (block or split simply)
710  if (parameters.psdLocalDiffusions.approximationMethod == "AM_Block_C" || parameters.psdLocalDiffusions.approximationMethod == "AM_Block_LR") {
711  // "_C" and "_LR" are "central" and "left-right" approximation methods
712 
713  if (localDiffusionsGrid.pc.size > 1 && localDiffusionsGrid.alpha.size > 1) {
714  Output::echo(3, "Calculating local diffusions using BLOCK method... \n");
715  // local diffusions Momentum + Pitch angle
716  // ??? mixterms dosn't work here at all ???
717  psdLocalDiffusions.Diffusion_pc_alpha(parameters.timeStep,
718  parameters.Lpp[iteration],
719  Dpcpc.CurrentDxx, DpcpcLpp.CurrentDxx,
720  Daa.CurrentDxx, DaaLpp.CurrentDxx,
721  Dpca.CurrentDxx, DpcaLpp.CurrentDxx,
722  localDiffusionsGrid.L, localDiffusionsGrid.pc, localDiffusionsGrid.alpha, localDiffusionsGrid.Jacobian,
723  pc_LowerBoundaryCondition.arr,
724  pc_UpperBoundaryCondition.arr,
725  alpha_LowerBoundaryCondition.arr,
726  alpha_UpperBoundaryCondition.arr,
727  pc_LowerBoundaryCondition.calculationType,
728  pc_UpperBoundaryCondition.calculationType,
729  alpha_LowerBoundaryCondition.calculationType,
730  alpha_UpperBoundaryCondition.calculationType);
731  Output::echo(3, "Done. \n");
732 
733  if (iteration == 1 && parameters.outputModelMatrix) { // Id Debug folder exists and outputModelMatrix - write down the Matrix from the Local diffusion computation
734  psdLocalDiffusions.matr_A.writeToFile("./Debug/Matr_A.dat");
735  psdLocalDiffusions.matr_B.writeToFile("./Debug/Matr_B.dat");
736  psdLocalDiffusions.matr_C.writeToFile("./Debug/Matr_C.dat");
737  }
738  } // grid check
739 
740  } // approximationMethod
760  else if (parameters.psdLocalDiffusions.approximationMethod == "AM_Split_C"
761  || parameters.psdLocalDiffusions.approximationMethod == "AM_Split_LR"
762  || parameters.psdLocalDiffusions.approximationMethod == "AM_Split_C_oldway"
763  || parameters.psdLocalDiffusions.approximationMethod == "AM_Split_LR_oldway") {
764 
765  // local diffusion
766  // Momentum diffusion
767  if (parameters.useEnergyDiffusion && localDiffusionsGrid.pc.size > 1) {
768  Output::echo(3, "Calculating energy diffusions using SPLIT method... \n");
769  psdLocalDiffusions.Diffusion_pc(parameters.timeStep,
770  parameters.Lpp[iteration],
771  Dpcpc.CurrentDxx, DpcpcLpp.CurrentDxx,
772  localDiffusionsGrid.L, localDiffusionsGrid.pc, localDiffusionsGrid.alpha, localDiffusionsGrid.Jacobian,
773  pc_LowerBoundaryCondition.arr,
774  pc_UpperBoundaryCondition.arr,
775  pc_LowerBoundaryCondition.calculationType,
776  pc_UpperBoundaryCondition.calculationType);
777  Output::echo(3, "Done. \n");
778 
779  if (iteration == 1 && parameters.outputModelMatrix) { /* TBD: Model matrix output here */ }
780  } // grid check
781 
782  // local diffusion
783  // Pitch angle diffusion
784  if (parameters.useAlphaDiffusion && localDiffusionsGrid.alpha.size > 1) {
785  Output::echo(3, "Calculating pitch angle diffusions using SPLIT method... \n");
786  psdLocalDiffusions.Diffusion_alpha(parameters.timeStep,
787  parameters.Lpp[iteration],
788  Daa.CurrentDxx, DaaLpp.CurrentDxx,
789  localDiffusionsGrid.L, localDiffusionsGrid.pc, localDiffusionsGrid.alpha, localDiffusionsGrid.Jacobian,
790  alpha_LowerBoundaryCondition.arr,
791  alpha_UpperBoundaryCondition.arr,
792  alpha_LowerBoundaryCondition.calculationType,
793  alpha_UpperBoundaryCondition.calculationType);
794  Output::echo(3, "Done. \n");
795 
796  if (iteration == 1 && parameters.outputModelMatrix) { /* TBD: Model matrix output here */ }
797  } // grid check
798 
799  // Explicit mixed terms
800  if (parameters.useEnergyAlphaMixedTerms && localDiffusionsGrid.pc.size > 1 && localDiffusionsGrid.alpha.size > 1) {
801 
802  Output::echo(3, "Calculating mixed terms EXPLICITLY ... \n");
803  psdLocalDiffusions.DiffusionMixTermExplicit(parameters.timeStep,
804  parameters.Lpp[iteration],
805  Dpca.CurrentDxx, DpcaLpp.CurrentDxx,
806  localDiffusionsGrid.L, localDiffusionsGrid.pc, localDiffusionsGrid.alpha, localDiffusionsGrid.Jacobian,
807  pc_LowerBoundaryCondition.arr,
808  pc_UpperBoundaryCondition.arr,
809  alpha_LowerBoundaryCondition.arr,
810  alpha_UpperBoundaryCondition.arr,
811  pc_LowerBoundaryCondition.calculationType,
812  pc_UpperBoundaryCondition.calculationType,
813  alpha_LowerBoundaryCondition.calculationType,
814  alpha_UpperBoundaryCondition.calculationType);
815  Output::echo(3, "Done. \n");
816  } // grid check
817 
818  } // approximationMethod
819  else if (parameters.psdLocalDiffusions.approximationMethod == "AM_Kyung-Chan") { // Undocumented Kyung-Chan method
820 
821  if (localDiffusionsGrid.pc.size > 1 && localDiffusionsGrid.alpha.size > 1) {
822  Output::echo(3, "Calculating local diffusions using Kyung-Chan's approximation... \n");
823  // do local diffusions
824  psdLocalDiffusions.Diffusion_pc_alpha_KC(parameters.timeStep,
825  parameters.Lpp[iteration],
826  Dpcpc.CurrentDxx, DpcpcLpp.CurrentDxx,
827  Daa.CurrentDxx, DaaLpp.CurrentDxx,
828  Dpca.CurrentDxx, DpcaLpp.CurrentDxx,
829  localDiffusionsGrid.L, localDiffusionsGrid.pc, localDiffusionsGrid.alpha, localDiffusionsGrid.Jacobian,
830  pc_LowerBoundaryCondition.arr,
831  pc_UpperBoundaryCondition.arr,
832  alpha_LowerBoundaryCondition.arr,
833  alpha_UpperBoundaryCondition.arr,
834  pc_LowerBoundaryCondition.calculationType,
835  pc_UpperBoundaryCondition.calculationType,
836  alpha_LowerBoundaryCondition.calculationType,
837  alpha_UpperBoundaryCondition.calculationType);
838  Output::echo(3, "Done. \n");
839 
840  if (iteration == 1 && parameters.outputModelMatrix) { // Id Debug folder exists and outputModelMatrix - write down the Matrix from the Local diffusion computation
841  psdLocalDiffusions.matr_A.writeToFile("./Debug/Matr_A.dat");
842  psdLocalDiffusions.matr_B.writeToFile("./Debug/Matr_B.dat");
843  psdLocalDiffusions.matr_C.writeToFile("./Debug/Matr_C.dat");
844  }
845  } // grid check
846 
847  } else {
848  throw error_msg("APPROX_ERROR", "Approximation unknown");
849  } // end of if(approximationMethod)
850 
851 
852  //-----------------
853  // Sources & Losses
854  //-----------------
855  // Loss-Cone losses, electron lifetime losses (if any), sources from coupling with RCM, etc
856 
865  SL.Update(iteration, parameters.SL, localDiffusionsGrid); // update array of SL
866 
867  psdLocalDiffusions.SourcesAndLosses(
868  parameters,
869  localDiffusionsGrid.L, localDiffusionsGrid.pc, localDiffusionsGrid.alpha,
870  SL,
871  parameters.timeStep,
872  parameters.Lpp[iteration],
873  parameters.tau[iteration],
874  parameters.tauLpp[iteration], parameters.Kp[iteration]);
875 
876 
877  //----------------------------------------------------------------------------------------------------
878  // Make PSD positive
879  //----------------------------------------------------------------------------------------------------
889  if (parameters.NoNegative) // This check needs for correct strange calculation when PSD === 0
890  for (il = 0; il < localDiffusionsGrid.L.size; il++)
891  for (im = 0; im < localDiffusionsGrid.pc.size; im++)
892  for (ia = 0; ia < localDiffusionsGrid.alpha.size; ia++)
893  if (psdLocalDiffusions.arr[il][im][ia] < VC::zero_f)
894  psdLocalDiffusions.arr[il][im][ia] = VC::zero_f; // keep negative value to see where the negative values appears...
895 
896  //----------------------------------------------------------------------------------------------------
897  // Checking the answer (probably pretty meaningless, since we check the inversion of the model matrix)
898  //----------------------------------------------------------------------------------------------------
899  // ??? \deprecated ?
900  if (false) {
901  // check the local solution
902  // 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
903  static Matrix3D<double> PSD_Der_a_L(localDiffusionsGrid.L.size, localDiffusionsGrid.pc.size, localDiffusionsGrid.alpha.size);
904  static Matrix3D<double> PSD_Der_a_R(localDiffusionsGrid.L.size, localDiffusionsGrid.pc.size, localDiffusionsGrid.alpha.size);
905  static Matrix3D<double> PSD_Der_m_L(localDiffusionsGrid.L.size, localDiffusionsGrid.pc.size, localDiffusionsGrid.alpha.size);
906  static Matrix3D<double> PSD_Der_m_R(localDiffusionsGrid.L.size, localDiffusionsGrid.pc.size, localDiffusionsGrid.alpha.size);
907 
908  static Matrix3D<double> error(localDiffusionsGrid.L.size, localDiffusionsGrid.pc.size, localDiffusionsGrid.alpha.size);
909 
910  // Some variables
911  double LaaLR, LaaRL, LmmLR, LmmRL;
912  double LmaLL, LmaLR, LmaRL, LmaRR;
913  double LamLL, LamLR, LamRL, LamRR;
914  double taulc;
915 
916  error = 0;
917  double maxerror = 0, aveerror = 0;
918 
919  for (il = 0; il < localDiffusionsGrid.L.size; il++) {
920  for (im = 0; im < localDiffusionsGrid.pc.size; im++) {
921  ia = 0;
922  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]);
923  for (ia = 1; ia < localDiffusionsGrid.alpha.size-1; ia++) {
924  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]);
925  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]);
926  }
927  ia = localDiffusionsGrid.alpha.size-1;
928  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]);
929  }
930  for (ia = 0; ia < localDiffusionsGrid.alpha.size; ia++) {
931  im = 0;
932  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]);
933  for (im = 1; im < localDiffusionsGrid.pc.size-1; im++) {
934  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]);
935  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]);
936  }
937  im = localDiffusionsGrid.pc.size-1;
938  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]);
939  }
940  }
941 
942  for (il = 0; il < localDiffusionsGrid.L.size; il++) {
943  for (im = 1; im < localDiffusionsGrid.pc.size-1; im++) {
944  for (ia = 1; ia < localDiffusionsGrid.alpha.size-1; ia++) {
945 
946  // dda D d/da f
947  LaaLR = 1.0/(localDiffusionsGrid.Jacobian[il][im][ia]) *
948  (PSD_Der_a_R[il][im][ia] - PSD_Der_a_R[il][im][ia-1])/
949  (localDiffusionsGrid.alpha.arr[il][im][ia] - localDiffusionsGrid.alpha.arr[il][im][ia-1]);
950  LaaRL = 1.0/(localDiffusionsGrid.Jacobian[il][im][ia]) *
951  (PSD_Der_a_L[il][im][ia+1] - PSD_Der_a_L[il][im][ia])/
952  (localDiffusionsGrid.alpha.arr[il][im][ia+1] - localDiffusionsGrid.alpha.arr[il][im][ia]);
953 
954  // dda D d/da f
955  LmmLR = 1.0/(localDiffusionsGrid.Jacobian[il][im][ia]) *
956  (PSD_Der_m_R[il][im][ia] - PSD_Der_m_R[il][im-1][ia])/
957  (localDiffusionsGrid.pc.arr[il][im][ia] - localDiffusionsGrid.pc.arr[il][im-1][ia]);
958  LmmRL = 1.0/(localDiffusionsGrid.Jacobian[il][im][ia]) *
959  (PSD_Der_m_L[il][im+1][ia] - PSD_Der_m_L[il][im][ia])/
960  (localDiffusionsGrid.pc.arr[il][im+1][ia] - localDiffusionsGrid.pc.arr[il][im][ia]);
961 
962  // d/dp D d/da f
963  LmaRR = 1.0/(localDiffusionsGrid.Jacobian[il][im][ia]) *
964  (PSD_Der_a_R[il][im+1][ia] - PSD_Der_a_R[il][im][ia])/
965  (localDiffusionsGrid.pc.arr[il][im+1][ia] - localDiffusionsGrid.pc.arr[il][im][ia]);
966  LmaRL = 1.0/(localDiffusionsGrid.Jacobian[il][im][ia]) *
967  (PSD_Der_a_L[il][im+1][ia] - PSD_Der_a_L[il][im][ia])/
968  (localDiffusionsGrid.pc.arr[il][im+1][ia] - localDiffusionsGrid.pc.arr[il][im][ia]);
969  LmaLR = 1.0/(localDiffusionsGrid.Jacobian[il][im][ia]) *
970  (PSD_Der_a_R[il][im][ia] - PSD_Der_a_R[il][im-1][ia])/
971  (localDiffusionsGrid.pc.arr[il][im][ia] - localDiffusionsGrid.pc.arr[il][im-1][ia]);
972  LmaLL = 1.0/(localDiffusionsGrid.Jacobian[il][im][ia]) *
973  (PSD_Der_a_L[il][im][ia] - PSD_Der_a_L[il][im-1][ia])/
974  (localDiffusionsGrid.pc.arr[il][im][ia] - localDiffusionsGrid.pc.arr[il][im-1][ia]);
975 
976 
977  // d/da D d/dp f
978  LamRR = 1.0/(localDiffusionsGrid.Jacobian[il][im][ia]) *
979  (PSD_Der_m_R[il][im][ia+1] - PSD_Der_m_R[il][im][ia])/
980  (localDiffusionsGrid.alpha.arr[il][im][ia+1] - localDiffusionsGrid.alpha.arr[il][im][ia]);
981  LamRL = 1.0/(localDiffusionsGrid.Jacobian[il][im][ia]) *
982  (PSD_Der_m_L[il][im][ia+1] - PSD_Der_m_L[il][im][ia])/
983  (localDiffusionsGrid.alpha.arr[il][im][ia+1] - localDiffusionsGrid.alpha.arr[il][im][ia]);
984  LamLR = 1.0/(localDiffusionsGrid.Jacobian[il][im][ia]) *
985  (PSD_Der_m_R[il][im][ia] - PSD_Der_m_R[il][im][ia-1])/
986  (localDiffusionsGrid.alpha.arr[il][im][ia] - localDiffusionsGrid.alpha.arr[il][im][ia-1]);
987  LamLL = 1.0/(localDiffusionsGrid.Jacobian[il][im][ia]) *
988  (PSD_Der_m_L[il][im][ia] - PSD_Der_m_L[il][im][ia-1])/
989  (localDiffusionsGrid.alpha.arr[il][im][ia] - localDiffusionsGrid.alpha.arr[il][im][ia-1]);
990 
991  error[il][im][ia] = psdLocalDiffusions.arr[il][im][ia] - psdPrevious[il][im][ia]
992  - (LaaLR + LaaRL) / 2 * parameters.timeStep // d/da D d/da f
993  - (LmmLR + LmmRL) / 2 * parameters.timeStep // d/da D d/da f
994  - (LmaRR + LmaRL + LmaLR + LmaLL) / 4 * parameters.timeStep // d/dm D d/da f
995  - (LamRR + LamRL + LamLR + LamLL) / 4 * parameters.timeStep; // d/dm D d/da f
996 
997  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])))) ) {
998  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]);
999  error[il][im][ia] += psdLocalDiffusions.arr[il][im][ia]/taulc * parameters.timeStep;
1000  }
1001 
1002  aveerror += fabs(error[il][im][ia]);
1003  maxerror = (fabs(maxerror) > fabs(error[il][im][ia]))?maxerror:error[il][im][ia];
1004  }
1005  }
1006  }
1007  aveerror = aveerror / (localDiffusionsGrid.L.size) / (localDiffusionsGrid.epc.size-2) / (localDiffusionsGrid.alpha.size-2);
1008  Output::echo(20, "Max calculation error = %e; ave calculation error = %e\n", maxerror, aveerror);
1009 
1010  stringstream filename;
1011  filename << "./Debug/Calc_error_" << iteration << ".plt";
1012 
1013  error.writeToFile(filename.str().c_str(), localDiffusionsGrid.L.arr, localDiffusionsGrid.epc.arr, localDiffusionsGrid.alpha.arr);
1014 
1015  psdPrevious = psdLocalDiffusions.arr;
1016  }
1017 
1018 
1019  //----------------------------------------------------------------------------------------------------
1020  // Writing output
1021  //----------------------------------------------------------------------------------------------------
1022  // Write PSD into a file
1023  // namespace VF defined in variousFunctions.cpp
1024  // If this calculation step correspond to write time - do write to output file
1025  if (VF::check_time(iteration, parameters.general_Output_parameters.iterStep)) {
1026  Output::echo(0, "Step %d of %d. Time = %d days, %.0f hours. ", iteration, parameters.totalIterationsNumber, int(simulation_time), (simulation_time-int(simulation_time))*24);
1027  Output::echo(0, "Writing data files.\n");
1028 
1029  // Write information about 1d variables into the file with 1d variables
1030  Output1DValues(output1D, Daa, Dpcpc, Dpca, DaaLpp, DpcpcLpp, DpcaLpp, simulation_time, parameters, iteration);
1031 
1032  //
1033  // Most important function with write down the calculated PSD
1034  //
1035  psdLocalDiffusions.Output_without_grid(simulation_time);
1036  if (radialDiffusionGrid.L.size > 1) psdRadialDiffusion.Output_without_grid(simulation_time);
1037 
1038  // 'Flushing' of log files (saving from memory to disk). Otherwise it saves it only at the end, which is not convenient.
1040 
1041  // Save the last PSD so it can be used as an initial condition to continue the calculation
1042  psdRadialDiffusion.arr.writeToFile(parameters.general_Output_parameters.folderName + "/PSD0.dat", radialDiffusionGrid.L.arr, radialDiffusionGrid.epc.arr, radialDiffusionGrid.alpha.arr);
1043  }
1044 
1052  if (radialDiffusionGrid.L.size > 1) {
1053  // Interpolation from Local diffusions grid to Radial diffusion grid
1054  // local->radial
1055  Output::echo(3, "Interpolation local->radial... ");
1056 
1057  //* For time-dependent boundary conditions (used for RCM coupling) use next
1058  //* *.arr.xSlice(iteration) create a 2D array of boundary values at the current time
1059  //* psdRadialDiffusion.Interpolate(psdLocalDiffusions, parameters.interpolation, localDiffusionsGrid, radialDiffusionGrid, pc_LowerBoundaryCondition.arr.xSlice(iteration), pc_UpperBoundaryCondition.arr.xSlice(iteration));
1060 
1061  //* Constant in time boundary conditions (a typical simulation)
1062  //* psdRadialDiffusion.Interpolate(psdLocalDiffusions, parameters.interpolation, localDiffusionsGrid, radialDiffusionGrid, pc_LowerBoundaryCondition.arr.xSlice(0), pc_UpperBoundaryCondition.arr.xSlice(0));
1063 
1064  // Constant in time low-energy boundary conditions, but flexible high-Energy boundary.
1065  // Created to fix interpolation problems.
1066  psdRadialDiffusion.Interpolate(psdLocalDiffusions, parameters.interpolation, localDiffusionsGrid, radialDiffusionGrid, pc_LowerBoundaryCondition.arr, psdLocalDiffusions.arr.ySlice(radialDiffusionGrid.pc.size-1));
1067  Output::echo(3, "done!\n");
1068  }
1069 
1070  } // end if for Main loop.
1071 
1072  //----------------------------------
1073  // End of simulation
1074  //----------------------------------
1075 
1076  // Last write the log
1077  std::time(&now);
1078  current_time = std::clock();
1079  Output::echo(0, "End of simulation %s", ctime(&now));
1080  Output::echo(0, "Total time %0.2f sec\n", ( current_time - zero_time ) / (double)CLOCKS_PER_SEC );
1081 
1082  // Close the log
1084  return 0; // end program if there wasn't any exeption
1085 
1086 
1087  //----------------------------------
1088  // Exceptions
1089  //----------------------------------
1090  // \deprecated?
1091 
1092  // try block ends.
1093  } catch (error_msg err) {
1094  Output::echo(0, "\nCRITICAL ERROR:\n%s\n", err.what().c_str());
1096  return 1;
1097  } catch (exception exp) {
1098  Output::echo(0, "\nCRITICAL ERROR:\n%s\n", exp.what());
1100  return 1;
1101  }
1102 
1103  //* #ifndef DEBUG_MODE
1104  //* catch (...) {
1105  //* Output::echo(0, "Unknown exception!");
1106  //* }
1107  //* #endif
1108 
1109  Output::echo(0, "UNEXPECTED EXIT\n");
1111 
1112  //* Memory leaks
1113  //* _CrtDumpMemoryLeaks();
1114  return 0;
1115 }
Holds list of instances of DiffusionCoefficient class of same type (like Daa, Dpp, etc), but produced by different waves (Daa_chorus, Daa_EMIC, etc).
void Update(int iteration, Parameters_structure::SL_structure &SL, Grid &grid)
GridElement alpha
grid elements to be used
Definition: Grid.h:59
Matrix3D< double > arr
array of grid points
Definition: Grid.h:30
Matrix1D< double > Kp
Kp array.
Definition: Parameters.h:106
void Diffusion_L(double dt, double Lpp, DiffusionCoefficient &DLL, GridElement &L, GridElement &pc, GridElement &alpha, Matrix3D< double > Jacobian, Matrix2D< double > lowerBoundaryCondition, Matrix2D< double > upperBoundaryCondition, string lowerBoundaryCondition_calculationType, string upperBoundaryCondition_calculationType, double tau, double tauLpp)
Definition: PSD.cpp:895
void Interpolate(PSD &otherPSD, Parameters_structure::Interpolation interpolationParameters, Grid &oldGrid, Grid &newGrid, Matrix2D< double > newGrid_pc_lowerBoundaryCondition, Matrix2D< double > newGrid_pc_upperBoundaryCondition)
Interpolation function.
Definition: PSD.cpp:1557
struct Parameters_structure::GridElement radialDiffusionGrid_L
grid elements for local and radial diffusion, and corresponding types of elements ...
bool ActivateAndScale(double time, double Kp)
void SourcesAndLosses(Parameters_structure parameters, GridElement &L, GridElement &pc, GridElement &alpha, AdditionalSourcesAndLosses &SL, double dt, double Lpp, double tau, double tauLpp, double Kp)
Definition: PSD.cpp:319
Calculates any additional sources and losses due to magnetopause.
struct Parameters_structure::Interpolation interpolation
Interpolation parameters.
Matrix2D< T > ySlice(int p_y) const
Definition: Matrix.cpp:1373
bool useLossCone
specify whether using loss cone or not
Definition: Parameters.h:128
#define __VERB_VERSION_NUMBER__
Number of current VERB code version.
Definition: Main.cpp:124
void Output1DValues(ofstream &output1D, DiffusionCoefficientsGroup &Daa, DiffusionCoefficientsGroup &Dpcpc, DiffusionCoefficientsGroup &Dpca, DiffusionCoefficientsGroup &DaaLpp, DiffusionCoefficientsGroup &DpcpcLpp, DiffusionCoefficientsGroup &DpcaLpp, double &time, Parameters_structure &parameters, int iteration)
Function used in main() to print values of diffusion coefficients in output file. ...
void Diffusion_alpha(double dt, double Lpp, DiffusionCoefficient &Daa, DiffusionCoefficient &DaaLpp, GridElement &L, GridElement &pc, GridElement &alpha, Matrix3D< double > Jacobian, Matrix2D< double > alpha_lowerBoundaryCondition, Matrix2D< double > alpha_upperBoundaryCondition, string alpha_lowerBoundaryCondition_calculationType, string alpha_upperBoundaryCondition_calculationType)
Definition: PSD.cpp:599
void writeToFile(string filename)
Save matrix to a file.
Definition: Matrix.cpp:1164
Holds diffusion coefficient matrix and routines to load and calculate it.
bool useAlphaDiffusion
Using diffusions flags.
Definition: Parameters.h:101
void Diffusion_pc_alpha(double dt, double Lpp, DiffusionCoefficient &Dpcpc, DiffusionCoefficient &DpcpcLpp, DiffusionCoefficient &Daa, DiffusionCoefficient &DaaLpp, DiffusionCoefficient &Dpca, DiffusionCoefficient &DpcaLpp, GridElement &L, GridElement &pc, GridElement &alpha, Matrix3D< double > Jacobian, Matrix2D< double > pc_lowerBoundaryCondition, Matrix2D< double > pc_upperBoundaryCondition, Matrix2D< double > alpha_lowerBoundaryCondition, Matrix2D< double > alpha_upperBoundaryCondition, string pc_lowerBoundaryCondition_calculationType, string pc_upperBoundaryCondition_calculationType, string alpha_lowerBoundaryCondition_calculationType, string alpha_upperBoundaryCondition_calculationType)
Definition: PSD.cpp:1053
string calculationType
Type: constant function/constant derivative.
string what()
Definition: error.h:117
void Create_Initial_PSD(Parameters_structure::PSD parameters, Grid &grid, BoundaryCondition L_UpperBoundaryCondition)
Create inital PSD (Steady State)
Definition: PSD.cpp:47
void Output_without_grid(double time)
Definition: PSD.cpp:1720
void set_output_lvl(int new_outputLvl)
Set the verbose level until redefine.
Definition: Output.cpp:40
struct Parameters_structure::General_Output_parameters general_Output_parameters
instance of the struct
Matrix2D< double > arr
2D matrix of boundary conditions
DiffusionCoefficient CurrentDxx
flag, indicated that the initialization was passed
Matrix1D< double > Bf
Boundary flux array.
Definition: Parameters.h:110
static const double exp
exponential
void echo(int msglvl, const char *format,...)
Basic function! Write the message to the file.
Definition: Output.cpp:44
void open_log_file(string filename)
Open file stream for log file.
Definition: Output.cpp:23
void AllocateMemory(int size_x, int size_y, int size_z)
Definition: Matrix.cpp:841
void flush_log_file()
flush file stream for log file
Definition: Output.cpp:36
string radialDiffusionGrid_filename
filename, if grid load from file
Definition: Parameters.h:145
struct Parameters_structure::SL_structure SL
instance of the struct
void Create_Grid(Parameters_structure::GridElement parameters_L, Parameters_structure::GridElement parameters_pc, Parameters_structure::GridElement parameters_alpha, Parameters_structure::GridElement parameters_epc, string grid_filename, string gridType, Grid SecondGrid=Grid())
Definition: Grid.cpp:79
int main(int argc, char *argv[])
Main code.
Definition: Main.cpp:195
int iterStep
Number of iterations between outputs.
Definition: Parameters.h:135
string approximationMethod
Approximation method. Check StrToVal(string input, ApproximationMethods &place) for known values...
Definition: Parameters.h:216
static const double pi
π
void MakeBoundaryCondition(Parameters_structure::BoundaryCondition parameters, Matrix2D< double > psd2DSlice, Matrix2D< double > gridElement1, Matrix2D< double > gridElement2)
void DiffusionMixTermExplicit(double dt, double Lpp, DiffusionCoefficient &Dpca, DiffusionCoefficient &DpcaLpp, GridElement &L, GridElement &pc, GridElement &alpha, Matrix3D< double > Jacobian, Matrix2D< double > pc_lowerBoundaryCondition, Matrix2D< double > pc_upperBoundaryCondition, Matrix2D< double > alpha_lowerBoundaryCondition, Matrix2D< double > alpha_upperBoundaryCondition, string pc_lowerBoundaryCondition_calculationType, string pc_upperBoundaryCondition_calculationType, string alpha_lowerBoundaryCondition_calculationType, string alpha_upperBoundaryCondition_calculationType)
Definition: PSD.cpp:414
Matrix1D< double > Lpp
Lpp array.
Definition: Parameters.h:114
string DLLType
DLL type (which method to use to calculate). Check StrToVal(string input, DLLTypes &place) for known ...
Definition: Parameters.h:105
string type
Type. Check StrToVal(string input, BoundaryConditionTypes &place) for known values.
Definition: Parameters.h:173
void Diffusion_pc(double dt, double Lpp, DiffusionCoefficient &Dpcpc, DiffusionCoefficient &DpcpcLpp, GridElement &L, GridElement &pc, GridElement &alpha, Matrix3D< double > Jacobian, Matrix2D< double > pc_lowerBoundaryCondition, Matrix2D< double > pc_upperBoundaryCondition, string pc_lowerBoundaryCondition_calculationType, string pc_upperBoundaryCondition_calculationType)
Definition: PSD.cpp:747
Main parameters structure that holds smaller structures for individual parameters.
Definition: Parameters.h:94
bool useEnergyDiffusion
Using diffusions flags.
Definition: Parameters.h:102
int size
size of grid element
Definition: Grid.h:39
int totalIterationsNumber
Total number of iterations.
Definition: Parameters.h:99
Matrix3D< double > Jacobian
Definition: Grid.h:84
void close_log_file()
Close file stream for log file.
Definition: Output.cpp:32
double bounce_time_new(double L, double pc, double alpha)
bounce time
void Create_SL(Parameters_structure::SL_structure &SL, Grid &grid, int numberOfIterations, double timeStep)
Matrix3D< double > arr
array of diffusion coefficients
bool NoNegative
Artificially no negative PSD.
Definition: Parameters.h:126
bool Load_parameters(string filename)
Definition: Parameters.cpp:72
void CreateAllDiffusionCoefficients(DiffusionCoefficient &DLL, DiffusionCoefficientsGroup &Daa, DiffusionCoefficientsGroup &Dpcpc, DiffusionCoefficientsGroup &Dpca, DiffusionCoefficientsGroup &DaaLpp, DiffusionCoefficientsGroup &DpcpcLpp, DiffusionCoefficientsGroup &DpcaLpp, Parameters_structure &parameters, Grid &radialDiffusionGrid, Grid &localDiffusionsGrid)
string radialDiffusionGrid_type
Local diffusions grid type.
Definition: Parameters.h:142
void Output1DHeaders(ofstream &output1D, DiffusionCoefficientsGroup &Daa, DiffusionCoefficientsGroup &Dpcpc, DiffusionCoefficientsGroup &Dpca, DiffusionCoefficientsGroup &DaaLpp, DiffusionCoefficientsGroup &DpcpcLpp, DiffusionCoefficientsGroup &DpcaLpp)
Function used in main() to print header in output file. Header needs to understand the output file...
void Diffusion_pc_alpha_KC(double dt, double Lpp, DiffusionCoefficient &Dpcpc, DiffusionCoefficient &DpcpcLpp, DiffusionCoefficient &Daa, DiffusionCoefficient &DaaLpp, DiffusionCoefficient &Dpca, DiffusionCoefficient &DpcaLpp, GridElement &L, GridElement &pc, GridElement &alpha, Matrix3D< double > Jacobian, Matrix2D< double > pc_lowerBoundaryCondition, Matrix2D< double > pc_upperBoundaryCondition, Matrix2D< double > alpha_lowerBoundaryCondition, Matrix2D< double > alpha_upperBoundaryCondition, string pc_lowerBoundaryCondition_calculationType, string pc_upperBoundaryCondition_calculationType, string alpha_lowerBoundaryCondition_calculationType, string alpha_upperBoundaryCondition_calculationType)
Definition: PSD.cpp:1269
bool outputModelMatrix
Output model matrix.
Definition: Parameters.h:125
Matrix3D< double > arr
array of PSD values
Definition: PSD.h:53
Phase Space Density class.
Definition: PSD.h:48
GridElement L
grid elements to be used
Definition: Grid.h:59
double timeStep
Time step, hours.
Definition: Parameters.h:98
double alc(double L)
Loss cone calculations.
Holds upper and lower boundary conditions.
Matrix2D< T > xSlice(int p_x) const
Definition: Matrix.cpp:1357
Parameters_structure parameters
Parameters structure, with all parameters from the parameters.ini file. The default parameters define...
Definition: Main.cpp:187
GridElement epc
grid elements to be used
Definition: Grid.h:59
Error message - stack of single_errors.
Definition: error.h:54
string localDiffusionsGrid_filename
filename, if grid load from file
Definition: Parameters.h:146
Matrix1D< double > tauLpp
Lifetime inside of the plasmasphere.
Definition: Parameters.h:123
void Output(string filename)
Definition: Grid.cpp:371
bool check_time(int iter, int d_iter)
Small routine for checking if it is right iteration for the next output.
string fileName1D
1d-output file name.
Definition: Parameters.h:138
string localDiffusionsGrid_type
Radial diffusions grid type.
Definition: Parameters.h:142
Matrix1D< double > tau
Lifetime out of the plasmasphere.
Definition: Parameters.h:122
T max()
Definition: Matrix.cpp:1305
Computational grid composed of 3 different GridElement.
Definition: Grid.h:53
void writeToFile(string filename)
Definition: Matrix.cpp:2000
static const double zero_f
Minimum value for PSD.
CalculationMatrix matr_C
calculation matrices A, B, and C
Definition: PSD.h:167
CalculationMatrix matr_B
calculation matrices A, B, and C
Definition: PSD.h:167
void Update(int iteration, Matrix2D< double > PSD_2D_Slice, double time=-1, double dt=-1)
bool useRadialDiffusion
Using diffusions flags.
Definition: Parameters.h:100
void LoadBoundaryCondition(Parameters_structure::BoundaryCondition parameters, Matrix2D< double > gridElement1, Matrix2D< double > gridElement2)
struct Parameters_structure::PSD psdRadialDiffusion
PSD parameters for radial diffusion.
int outputLvl
Detalization level of screen output.
Definition: Parameters.h:96
void MakeDLL(double Kp)
Making DLL.
Matrix2D< T > zSlice(int p_z) const
Definition: Matrix.cpp:1389
struct Parameters_structure::PSD psdLocalDiffusions
PSD parameters for local diffusions.
GridElement pc
grid elements to be used
Definition: Grid.h:59
CalculationMatrix matr_A
calculation matrices A, B, and C
Definition: PSD.h:167