104 #define _CRT_SECURE_NO_DEPRECATE
105 #define __VERB_VERSION_NUMBER__ "2.03.000"
118 #include "../Logging/Output.h"
121 #include "../VariousFunctions/variousFunctions.h"
124 #include "../Diffusion/PSD.h"
127 #include "../Grid/Grid.h"
130 #include "../Grid/BoundaryConditions.h"
133 #include "../Grid/AdditionalSourcesAndLosses.h"
136 #include "../Parameters/Parameters.h"
140 #if defined(_WIN32) || defined(_WIN64)
141 #define _CRTDBG_MAP_ALLOC
171 int main(
int argc,
char* argv[]) {
188 radialDiffusionGrid, localDiffusionsGrid;
192 psdRadialDiffusion, psdLocalDiffusions;
200 L_LowerBoundaryCondition, L_UpperBoundaryCondition,
201 pc_LowerBoundaryCondition, pc_UpperBoundaryCondition,
202 alpha_LowerBoundaryCondition, alpha_UpperBoundaryCondition;
217 double simulation_time = 0;
225 std::clock_t zero_time, start_time, current_time;
227 time_t now = time(NULL);
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;
242 Output::echo(0,
"Compiled time: %s %s\n", __DATE__, __TIME__);
377 psdLocalDiffusions.
Interpolate(psdRadialDiffusion, parameters.
interpolation, radialDiffusionGrid, localDiffusionsGrid, pc_LowerBoundaryCondition.
arr, pc_UpperBoundaryCondition.
arr);
407 Output1DHeaders(output1D, Daa, Dpcpc, Dpca, DaaLpp, DpcpcLpp, DpcaLpp);
409 Output1DValues(output1D, Daa, Dpcpc, Dpca, DaaLpp, DpcpcLpp, DpcaLpp, simulation_time, parameters, 0);
425 simulation_time = parameters.
timeStep * iteration;
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);
433 if (radialDiffusionGrid.
L.
size > 1) {
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);
443 Output::echo(3,
"Updating radial diffusion coefficient...\n");
444 DLL.
MakeDLL(radialDiffusionGrid.
L, radialDiffusionGrid.
pc, radialDiffusionGrid.
alpha, parameters.
Kp[iteration], parameters.
DLLType);
459 radialDiffusionGrid.
L, radialDiffusionGrid.
pc, radialDiffusionGrid.
alpha, radialDiffusionGrid.
Jacobian,
460 L_LowerBoundaryCondition.
arr,
461 L_UpperBoundaryCondition.
arr * parameters.
Bf[iteration],
464 parameters.
tau[iteration],
465 parameters.
tauLpp[iteration]);
486 psdLocalDiffusions.
Interpolate(psdRadialDiffusion, parameters.
interpolation, radialDiffusionGrid, localDiffusionsGrid, pc_LowerBoundaryCondition.
arr, psdRadialDiffusion.
arr.
ySlice(radialDiffusionGrid.
pc.
size-1));
496 if (localDiffusionsGrid.
pc.
size > 1) {
498 pc_LowerBoundaryCondition.
Update(iteration, psdLocalDiffusions.
arr.
ySlice(0));
499 pc_UpperBoundaryCondition.
Update(iteration, psdLocalDiffusions.
arr.
ySlice(radialDiffusionGrid.
pc.
size-1));
502 if (localDiffusionsGrid.
alpha.
size > 1) {
504 alpha_LowerBoundaryCondition.
Update(iteration, psdLocalDiffusions.
arr.
zSlice(0));
512 Output::echo(3,
"Updating momentum diffusion coefficients...\n");
518 Output::echo(3,
"Updating pitch-angle diffusion coefficients...\n");
524 Output::echo(3,
"Updating mixed diffusion coefficients... \n");
538 if (localDiffusionsGrid.
pc.
size > 1 && localDiffusionsGrid.
alpha.
size > 1) {
539 Output::echo(3,
"Calculating local diffusions using BLOCK method... \n");
543 parameters.
Lpp[iteration],
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,
574 Output::echo(3,
"Calculating energy diffusions using SPLIT method... \n");
576 parameters.
Lpp[iteration],
578 localDiffusionsGrid.
L, localDiffusionsGrid.
pc, localDiffusionsGrid.
alpha, localDiffusionsGrid.
Jacobian,
579 pc_LowerBoundaryCondition.
arr,
580 pc_UpperBoundaryCondition.
arr,
591 Output::echo(3,
"Calculating pitch angle diffusions using SPLIT method... \n");
593 parameters.
Lpp[iteration],
595 localDiffusionsGrid.
L, localDiffusionsGrid.
pc, localDiffusionsGrid.
alpha, localDiffusionsGrid.
Jacobian,
596 alpha_LowerBoundaryCondition.
arr,
597 alpha_UpperBoundaryCondition.
arr,
608 Output::echo(3,
"Calculating mixed terms EXPLICITLY ... \n");
610 parameters.
Lpp[iteration],
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,
627 if (localDiffusionsGrid.
pc.
size > 1 && localDiffusionsGrid.
alpha.
size > 1) {
628 Output::echo(3,
"Calculating local diffusions using Kyung-Chan's approximation... \n");
631 parameters.
Lpp[iteration],
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,
654 throw error_msg(
"APPROX_ERROR",
"Approximation unknown");
664 localDiffusionsGrid.
L, localDiffusionsGrid.
pc, localDiffusionsGrid.
alpha,
667 parameters.
Lpp[iteration],
668 parameters.
tau[iteration],
669 parameters.
tauLpp[iteration], parameters.
Kp[iteration]);
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++)
697 double LaaLR, LaaRL, LmmLR, LmmRL;
698 double LmaLL, LmaLR, LmaRL, LmaRR;
699 double LamLL, LamLR, LamRL, LamRR;
703 double maxerror = 0, aveerror = 0;
705 for (il = 0; il < localDiffusionsGrid.
L.
size; il++) {
706 for (im = 0; im < localDiffusionsGrid.
pc.
size; im++) {
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]);
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]);
716 for (ia = 0; ia < localDiffusionsGrid.
alpha.
size; ia++) {
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]);
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]);
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++) {
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]);
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]);
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]);
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]);
777 error[il][im][ia] = psdLocalDiffusions.
arr[il][im][ia] - psdPrevious[il][im][ia]
778 - (LaaLR + LaaRL) / 2 * parameters.
timeStep
779 - (LmmLR + LmmRL) / 2 * parameters.
timeStep
780 - (LmaRR + LmaRL + LmaLR + LmaLL) / 4 * parameters.
timeStep
781 - (LamRR + LamRL + LamLR + LamLL) / 4 * parameters.
timeStep;
785 error[il][im][ia] += psdLocalDiffusions.
arr[il][im][ia]/taulc * parameters.
timeStep;
788 aveerror += fabs(error[il][im][ia]);
789 maxerror = (fabs(maxerror) > fabs(error[il][im][ia]))?maxerror:error[il][im][ia];
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);
796 stringstream filename;
797 filename <<
"./Debug/Calc_error_" << iteration <<
".plt";
801 psdPrevious = psdLocalDiffusions.
arr;
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);
816 Output1DValues(output1D, Daa, Dpcpc, Dpca, DaaLpp, DpcpcLpp, DpcaLpp, simulation_time, parameters, iteration);
831 if (radialDiffusionGrid.
L.
size > 1) {
845 psdRadialDiffusion.
Interpolate(psdLocalDiffusions, parameters.
interpolation, localDiffusionsGrid, radialDiffusionGrid, pc_LowerBoundaryCondition.
arr, psdLocalDiffusions.
arr.
ySlice(radialDiffusionGrid.
pc.
size-1));
857 current_time = std::clock();
859 Output::echo(0,
"Total time %0.2f sec\n", ( current_time - zero_time ) / (
double)CLOCKS_PER_SEC );
876 }
catch (exception
exp) {