00001
00062
00063
00064
00065
00066 #define _CRT_SECURE_NO_DEPRECATE
00067
00068
00069
00070 #define VERB_VERSION_NUMBER 2.00
00071
00072
00073 #include <valarray>
00074 #include <math.h>
00075 #include <ctime>
00076 #include <time.h>
00077
00078
00079 #include "../Logging/Output.h"
00080
00081
00082 #include "../VariousFunctions/variousFunctions.h"
00083
00084
00085 #include "../Diffusion/PSD.h"
00086
00087
00088 #include "../Grid/Grid.h"
00089
00090
00091 #include "../Grid/BoundaryConditions.h"
00092
00093
00094 #include "../Grid/SourcesAndLosses.h"
00095
00096
00097 #include "../Parameters/Parameters.h"
00098
00099
00101 bool check_time(double time, double interval) {
00102 if (fabs(fmod(double(time + 1e-11), interval)) <= 1e-4) {
00103 return true;
00104 }
00105 return false;
00106 }
00107
00109 bool check_time(int iter, int d_iter) {
00110 if ((iter % d_iter) == 0) {
00111 return true;
00112 }
00113 return false;
00114 }
00115
00116
00117
00118 ParamStructure parameters;
00119
00124 int main(int argc, char* argv[]) {
00125
00126 try {
00127
00128 std::clock_t zero_time, start_time, current_time;
00129 zero_time = clock();
00130 time_t now;
00131
00132
00133
00134 Output::open_log_file("logfile.log");
00135
00136 time(&now);
00137 Output::echo(0, "VERB code VERB_VERSION_NUMBER\n");
00138 Output::echo(0, "Starting at %s", ctime(&now) );
00139
00140
00141 if (argc > 1) {
00142
00143
00144 parameters.Load_parameters(argv[1]);
00145 } else {
00146
00147 parameters.Load_parameters("Parameters.ini");
00148 }
00149
00150
00151 Output::set_output_lvl(parameters.outputLvl);
00152
00153
00154 Output::echo(1, "Creating radial diffusion gird... ");
00155 start_time = std::clock();
00156 Grid radialDiffusionGrid(parameters.radialDiffusionGrid_L,
00157 parameters.radialDiffusionGrid_pc,
00158 parameters.radialDiffusionGrid_alpha,
00159 parameters.radialDiffusionGrid_epc,
00160 parameters.radialDiffusionGrid_filename,
00161 parameters.radialDiffusionGrid_type);
00162 current_time = std::clock();
00163 Output::echo(1, "done (%.2f sec)\n", ( current_time - start_time ) / (double)CLOCKS_PER_SEC );
00164 Output::echo(1, "Creating local diffusions gird... ");
00165 Grid localDiffusionsGrid(parameters.localDiffusionsGrid_L,
00166 parameters.localDiffusionsGrid_pc,
00167 parameters.localDiffusionsGrid_alpha,
00168 parameters.localDiffusionsGrid_epc,
00169 parameters.localDiffusionsGrid_filename,
00170 parameters.localDiffusionsGrid_type, radialDiffusionGrid);
00171 Output::echo(1, "done!\n");
00172
00173
00174 localDiffusionsGrid.Output((parameters.general_Output_parameters.folderName + "perp_grid.plt").c_str());
00175 radialDiffusionGrid.Output((parameters.general_Output_parameters.folderName + "rad_grid.plt").c_str());
00176
00177
00178 BoundaryCondition
00179 L_LowerBoundaryCondition(
00180
00181 1,
00182 radialDiffusionGrid.pc.size,
00183 radialDiffusionGrid.alpha.size,
00184 parameters.L_LowerBoundaryCondition),
00185 L_UpperBoundaryCondition(
00186
00187 1,
00188 radialDiffusionGrid.pc.size,
00189 radialDiffusionGrid.alpha.size,
00190 parameters.L_UpperBoundaryCondition),
00191 pc_LowerBoundaryCondition(
00192
00193 1,
00194 localDiffusionsGrid.L.size,
00195 localDiffusionsGrid.alpha.size,
00196 parameters.pc_LowerBoundaryCondition),
00197 pc_UpperBoundaryCondition(
00198
00199 1,
00200 localDiffusionsGrid.L.size,
00201 localDiffusionsGrid.alpha.size,
00202 parameters.pc_UpperBoundaryCondition),
00203 alpha_LowerBoundaryCondition(
00204
00205 1,
00206 localDiffusionsGrid.L.size,
00207 localDiffusionsGrid.pc.size,
00208 parameters.alpha_LowerBoundaryCondition),
00209 alpha_UpperBoundaryCondition(
00210
00211 1,
00212 localDiffusionsGrid.L.size,
00213 localDiffusionsGrid.pc.size,
00214 parameters.alpha_UpperBoundaryCondition);
00215
00216
00217 L_UpperBoundaryCondition.LoadBoundaryCondition(
00218 radialDiffusionGrid.epc.xSlice(radialDiffusionGrid.L.size-1),
00219 radialDiffusionGrid.alpha.xSlice(radialDiffusionGrid.L.size-1));
00220 pc_LowerBoundaryCondition.LoadBoundaryCondition(
00221 localDiffusionsGrid.L.ySlice(0),
00222 localDiffusionsGrid.alpha.ySlice(0));
00223
00224
00225 Output::echo(1, "Radial diffusion initialization... ");
00226 PSD psdRadialDiffusion(parameters.psdRadialDiffusion, radialDiffusionGrid, L_UpperBoundaryCondition);
00227 Output::echo(1, "done!\n");
00228
00229
00230
00231 Output::echo(1, "Making boundary conditions... ");
00232
00233 L_LowerBoundaryCondition.MakeBoundaryCondition(
00234 psdRadialDiffusion.xSlice(0),
00235 radialDiffusionGrid.epc.xSlice(0),
00236 radialDiffusionGrid.alpha.xSlice(0));
00237 L_UpperBoundaryCondition.MakeBoundaryCondition(
00238 psdRadialDiffusion.xSlice(radialDiffusionGrid.L.size-1),
00239 radialDiffusionGrid.epc.xSlice(radialDiffusionGrid.L.size-1),
00240 radialDiffusionGrid.alpha.xSlice(radialDiffusionGrid.L.size-1));
00241
00242
00243 pc_LowerBoundaryCondition.MakeBoundaryCondition(
00244 psdRadialDiffusion.ySlice(0),
00245 localDiffusionsGrid.L.ySlice(0),
00246 localDiffusionsGrid.alpha.ySlice(0));
00247 pc_UpperBoundaryCondition.MakeBoundaryCondition(
00248 psdRadialDiffusion.ySlice(radialDiffusionGrid.pc.size-1),
00249 localDiffusionsGrid.L.ySlice(radialDiffusionGrid.pc.size-1),
00250 localDiffusionsGrid.alpha.ySlice(radialDiffusionGrid.pc.size-1));
00251
00252 alpha_LowerBoundaryCondition.MakeBoundaryCondition(
00253 psdRadialDiffusion.zSlice(0),
00254 localDiffusionsGrid.L.zSlice(0),
00255 localDiffusionsGrid.epc.zSlice(0));
00256 alpha_UpperBoundaryCondition.MakeBoundaryCondition(
00257 psdRadialDiffusion.zSlice(radialDiffusionGrid.alpha.size-1),
00258 localDiffusionsGrid.L.zSlice(radialDiffusionGrid.alpha.size-1),
00259 localDiffusionsGrid.epc.zSlice(radialDiffusionGrid.alpha.size-1));
00260 Output::echo(1, "done!\n");
00261
00262 Output::echo(1, "Local diffusions initialization... ");
00263 PSD psdLocalDiffusions(parameters.psdLocalDiffusions, localDiffusionsGrid);
00264
00265 psdLocalDiffusions.Interpolate(psdRadialDiffusion, parameters.interpolation, radialDiffusionGrid, localDiffusionsGrid, pc_LowerBoundaryCondition.xSlice(0), pc_UpperBoundaryCondition.xSlice(0));
00266 Output::echo(1, "done!\n");
00267
00269 SourcesAndLosses SL;
00270 if (parameters.sourcesAndLosses.SL_L_top || parameters.sourcesAndLosses.SL_E_min) {
00271 SL.Initialize(parameters.sourcesAndLosses, localDiffusionsGrid, parameters.totalIterationsNumber, parameters.timeStep);
00272 } else {
00273 SL.Initialize(parameters.sourcesAndLosses, localDiffusionsGrid, 1, parameters.timeStep);
00274 }
00275
00276
00277
00278 psdLocalDiffusions.Output_without_grid(0);
00279 psdRadialDiffusion.Output_without_grid(0);
00280
00281
00282 DiffusionCoefficient DLL(radialDiffusionGrid);
00283 DiffusionCoefficientsGroup Dpcpc("DCT_Dpp", localDiffusionsGrid);
00284 DiffusionCoefficientsGroup DpcpcLpp("DCT_DppLpp", localDiffusionsGrid);
00285 DiffusionCoefficientsGroup Daa("DCT_Daa", localDiffusionsGrid);
00286 DiffusionCoefficientsGroup DaaLpp("DCT_DaaLpp", localDiffusionsGrid);
00287
00288 DiffusionCoefficientsGroup Dpca("DCT_Dpa", localDiffusionsGrid);
00289 DiffusionCoefficientsGroup DpcaLpp("DCT_DpaLpp", localDiffusionsGrid);
00290
00291
00292 if (parameters.useAlphaDifusion) {
00293 Daa.Get(localDiffusionsGrid, parameters.DxxParamStructureList);
00294 DaaLpp.Get(localDiffusionsGrid, parameters.DxxParamStructureList);
00295 } else {
00296 Daa = 0;
00297 DaaLpp = 0;
00298 }
00299
00300 if (parameters.useEnergyDifusion) {
00301
00302
00303 Dpcpc.Get(localDiffusionsGrid, parameters.DxxParamStructureList, "DCT_Dpp");
00304 Dpcpc.Get(localDiffusionsGrid, parameters.DxxParamStructureList, "DCT_Dpcpc");
00305 DpcpcLpp.Get(localDiffusionsGrid, parameters.DxxParamStructureList, "DCT_DppLpp");
00306 DpcpcLpp.Get(localDiffusionsGrid, parameters.DxxParamStructureList, "DCT_DpcpcLpp");
00307 } else {
00308 Dpcpc = 0;
00309 DpcpcLpp = 0;
00310 }
00311
00312
00313 if (parameters.useRadialDifusion) {
00314 DLL.MakeDLL(radialDiffusionGrid.L, radialDiffusionGrid.pc, radialDiffusionGrid.alpha, parameters.Kp[0], parameters.DLLType);
00315
00316 } else {
00317 DLL = 0;
00318
00319 }
00320
00321 if (parameters.useAlphaDifusion && parameters.useEnergyDifusion) {
00322 Dpca.Get(localDiffusionsGrid, parameters.DxxParamStructureList, "DCT_Dpa");
00323 Dpca.Get(localDiffusionsGrid, parameters.DxxParamStructureList, "DCT_Dpca");
00324 DpcaLpp.Get(localDiffusionsGrid, parameters.DxxParamStructureList, "DCT_DpaLpp");
00325 DpcaLpp.Get(localDiffusionsGrid, parameters.DxxParamStructureList, "DCT_DpcaLpp");
00326 } else {
00327 Dpca = 0;
00328 DpcaLpp = 0;
00329 }
00330
00331
00332 ofstream output1D((parameters.general_Output_parameters.folderName + parameters.general_Output_parameters.fileName1D).c_str());
00333
00334
00335 output1D << "Variables = \"time\", \"Kp\", \"Boundary_fluxes\", \"Lpp\", \"tau\" ";
00336
00337
00340 unsigned int Dxx_it;
00341 for (Dxx_it = 0; Dxx_it < Daa.DxxList.size(); Dxx_it++) {
00342 output1D << "\"" << Daa.DxxList[Dxx_it].name << "\"" << ", ";
00343 }
00344 for (Dxx_it = 0; Dxx_it < Dpcpc.DxxList.size(); Dxx_it++) {
00345 output1D << "\"" << Dpcpc.DxxList[Dxx_it].name << "\"" << ", ";
00346 }
00347 for (Dxx_it = 0; Dxx_it < DaaLpp.DxxList.size(); Dxx_it++) {
00348 output1D << "\"" << DaaLpp.DxxList[Dxx_it].name << "\"" << ", ";
00349 }
00350 for (Dxx_it = 0; Dxx_it < DpcpcLpp.DxxList.size(); Dxx_it++) {
00351 output1D << "\"" << DpcpcLpp.DxxList[Dxx_it].name << "\"" << ", ";
00352 }
00353 output1D << endl;
00354
00355
00356 output1D << "ZONE T=\"1d-output\" " << endl;
00357 output1D << 0 << "\t" << parameters.Kp[0] << "\t" << parameters.Bf[0] << "\t" << parameters.Lpp[0] << "\t" << parameters.tau[0] << "\t";
00358 for (Dxx_it = 0; Dxx_it < Daa.DxxList.size(); Dxx_it++) {
00359 output1D << Daa.DxxList[Dxx_it].Scale(parameters.Kp[0]) << "\t";
00360 }
00361 for (Dxx_it = 0; Dxx_it < Dpcpc.DxxList.size(); Dxx_it++) {
00362 output1D << Dpcpc.DxxList[Dxx_it].Scale(parameters.Kp[0]) << "\t";
00363 }
00364 for (Dxx_it = 0; Dxx_it < DaaLpp.DxxList.size(); Dxx_it++) {
00365 output1D << DaaLpp.DxxList[Dxx_it].Scale(parameters.Kp[0]) << "\t";
00366 }
00367 for (Dxx_it = 0; Dxx_it < DpcpcLpp.DxxList.size(); Dxx_it++) {
00368 output1D << DpcpcLpp.DxxList[Dxx_it].Scale(parameters.Kp[0]) << "\t";
00369 }
00370 output1D << endl;
00371
00372
00373
00375
00376 double time = 0;
00377 int il, im, ia;
00378 Output::echo(0, "Starting...\n");
00379 for (int iteration = 1; iteration < parameters.totalIterationsNumber; iteration++) {
00380
00381 time = parameters.timeStep * iteration;
00382
00383
00384 Output::echo(1, "Step %d of %d. Time = %d days, %.0f hours.\n", iteration, parameters.totalIterationsNumber, int(time), (time-int(time))*24);
00385
00386
00388 if (radialDiffusionGrid.L.size > 1) {
00389
00390
00391 Output::echo(3, "Updating L-boundaries... ");
00392 L_LowerBoundaryCondition.Update(iteration, psdRadialDiffusion.xSlice(0));
00393 L_UpperBoundaryCondition.Update(iteration, psdRadialDiffusion.xSlice(radialDiffusionGrid.L.size-1));
00394 Output::echo(3, "done.\n");
00395
00396
00397 if (parameters.useRadialDifusion) {
00398 Output::echo(3, "Radial diffusion (%e)\n", DLL.max());
00399 Output::echo(3, "Making DLL... ");
00400
00401 DLL.MakeDLL(radialDiffusionGrid.L, radialDiffusionGrid.pc, radialDiffusionGrid.alpha, parameters.Kp[iteration], parameters.DLLType);
00402 Output::echo(3, "done!\n");
00403
00404 }
00405
00406
00407 psdRadialDiffusion.Diffusion_L(parameters.timeStep, parameters.Lpp[iteration],
00408 DLL,
00409 radialDiffusionGrid.L, radialDiffusionGrid.pc, radialDiffusionGrid.alpha, radialDiffusionGrid.Jacobian,
00410
00411
00412 L_LowerBoundaryCondition.xSlice(0),
00413 L_UpperBoundaryCondition.xSlice(0) * parameters.Bf[iteration],
00414 L_LowerBoundaryCondition.calculationType,
00415 L_UpperBoundaryCondition.calculationType,
00416 parameters.tau[iteration],
00417 parameters.tauLpp[iteration]);
00418
00419
00420 Output::echo(3, "Interpolation radial->local... ");
00421
00422
00423
00424
00425 psdLocalDiffusions.Interpolate(psdRadialDiffusion, parameters.interpolation, radialDiffusionGrid, localDiffusionsGrid, pc_LowerBoundaryCondition.xSlice(0), pc_UpperBoundaryCondition.xSlice(0));
00426 Output::echo(3, "done!\n");
00427
00428
00429 }
00430
00431
00432 if (parameters.psdLocalDiffusions.approximationMethod == "AM_Block_C" || parameters.psdLocalDiffusions.approximationMethod == "AM_Block_LR") {
00433
00434
00435
00436 if (localDiffusionsGrid.pc.size > 1 && localDiffusionsGrid.alpha.size > 1) {
00437
00438
00439 pc_LowerBoundaryCondition.Update(iteration, psdLocalDiffusions.ySlice(0));
00440 pc_UpperBoundaryCondition.Update(iteration, psdLocalDiffusions.ySlice(radialDiffusionGrid.pc.size-1));
00441
00442 alpha_LowerBoundaryCondition.Update(iteration, psdLocalDiffusions.zSlice(0));
00443 alpha_UpperBoundaryCondition.Update(iteration, psdLocalDiffusions.zSlice(radialDiffusionGrid.alpha.size-1));
00444
00445 if (parameters.useEnergyDifusion) {
00446
00447
00448
00449 Dpcpc.ActivateAndScale(time, parameters.Kp[iteration]);
00450 DpcpcLpp.ActivateAndScale(time, parameters.Kp[iteration]);
00451 Output::echo(3, "Momentum diffusion (%e/%e)... \n", Dpcpc.max(), DpcpcLpp.max());
00452 }
00453 if (parameters.useAlphaDifusion) {
00454
00455
00456
00457 Daa.ActivateAndScale(time, parameters.Kp[iteration]);
00458 DaaLpp.ActivateAndScale(time, parameters.Kp[iteration]);
00459 Output::echo(3, "Pitch-angle diffusion (%e/%e)... \n", Daa.max(), DaaLpp.max());
00460 }
00461 if (parameters.useAlphaDifusion && parameters.useEnergyDifusion) {
00462 Dpca.ActivateAndScale(time, parameters.Kp[iteration]);
00463 DpcaLpp.ActivateAndScale(time, parameters.Kp[iteration]);
00464 Output::echo(3, "Mixed term (implicit) (%e/%e)... \n", Dpca.max(), DpcaLpp.max());
00465 }
00466
00467
00468
00469
00470
00471
00472
00473
00474
00475
00476
00477
00478
00479
00480
00481 psdLocalDiffusions.Diffusion_pc_alpha(parameters.timeStep,
00482 parameters.Lpp[iteration],
00483 Dpcpc, DpcpcLpp,
00484 Daa, DaaLpp,
00485 Dpca, DpcaLpp,
00486 localDiffusionsGrid.L, localDiffusionsGrid.pc, localDiffusionsGrid.alpha, localDiffusionsGrid.Jacobian,
00487
00488
00489
00490
00491 pc_LowerBoundaryCondition.xSlice(0),
00492 pc_UpperBoundaryCondition.xSlice(0),
00493 alpha_LowerBoundaryCondition.xSlice(0),
00494 alpha_UpperBoundaryCondition.xSlice(0),
00495 pc_LowerBoundaryCondition.calculationType,
00496 pc_UpperBoundaryCondition.calculationType,
00497 alpha_LowerBoundaryCondition.calculationType,
00498 alpha_UpperBoundaryCondition.calculationType);
00499
00500
00501 if (iteration == 1 && parameters.outputModelMatrix) {
00502
00503 }
00504 if (parameters.useAlphaDifusion || parameters.useEnergyDifusion) Output::echo(3, "done!\n");
00505 } else {
00506 throw error_msg("BLOCK_WITH_1D", "Block method can not be used with 1D");
00507 }
00508 } else if (parameters.psdLocalDiffusions.approximationMethod == "AM_Split_C"
00509 || parameters.psdLocalDiffusions.approximationMethod == "AM_Split_LR"
00510 || parameters.psdLocalDiffusions.approximationMethod == "AM_Split_C_oldway"
00511 || parameters.psdLocalDiffusions.approximationMethod == "AM_Split_LR_oldway") {
00512
00513
00515 if (localDiffusionsGrid.pc.size > 1) {
00516
00517
00518 pc_LowerBoundaryCondition.Update(iteration, psdLocalDiffusions.ySlice(0));
00519 pc_UpperBoundaryCondition.Update(iteration, psdLocalDiffusions.ySlice(radialDiffusionGrid.pc.size-1));
00520
00521
00522 if (parameters.useEnergyDifusion) {
00523
00524
00525
00526
00527 if (Dpcpc.ActivateAndScale(time, parameters.Kp[iteration])) Dpcpc.writeToFile((parameters.general_Output_parameters.folderName + Dpcpc.name + "_" + VF::dtostr(time) + ".plt").c_str(), localDiffusionsGrid.L, localDiffusionsGrid.epc, localDiffusionsGrid.alpha);
00528 if (DpcpcLpp.ActivateAndScale(time, parameters.Kp[iteration])) DpcpcLpp.writeToFile((parameters.general_Output_parameters.folderName + DpcpcLpp.name + "_" + VF::dtostr(time) + ".plt").c_str(), localDiffusionsGrid.L, localDiffusionsGrid.epc, localDiffusionsGrid.alpha);
00529 Output::echo(3, "Momentum diffusion (%e/%e)... ", Dpcpc.max(), DpcpcLpp.max());
00530 }
00531
00532
00533
00534
00535
00536
00537
00538
00539
00540
00541 psdLocalDiffusions.Diffusion_pc(parameters.timeStep,
00542 parameters.Lpp[iteration],
00543 Dpcpc, DpcpcLpp,
00544 localDiffusionsGrid.L, localDiffusionsGrid.pc, localDiffusionsGrid.alpha, localDiffusionsGrid.Jacobian,
00545
00546
00547 pc_LowerBoundaryCondition.xSlice(0),
00548 pc_UpperBoundaryCondition.xSlice(0),
00549 pc_LowerBoundaryCondition.calculationType,
00550 pc_UpperBoundaryCondition.calculationType);
00551
00552 if (iteration == 1 && parameters.outputModelMatrix) {
00553
00554 }
00555 if (parameters.useEnergyDifusion) Output::echo(3, "done!\n");
00556 }
00558
00560 if (localDiffusionsGrid.alpha.size > 1) {
00561
00562
00563 alpha_LowerBoundaryCondition.Update(iteration, psdLocalDiffusions.zSlice(0));
00564 alpha_UpperBoundaryCondition.Update(iteration, psdLocalDiffusions.zSlice(radialDiffusionGrid.alpha.size-1));
00565
00566
00567 if (parameters.useAlphaDifusion) {
00568
00569
00570
00571
00572 if (Daa.ActivateAndScale(time, parameters.Kp[iteration]))
00573 Daa.writeToFile((parameters.general_Output_parameters.folderName + Daa.name + "_" + VF::dtostr(time) + ".plt").c_str(), localDiffusionsGrid.L, localDiffusionsGrid.epc, localDiffusionsGrid.alpha);
00574 if (DaaLpp.ActivateAndScale(time, parameters.Kp[iteration]))
00575 DaaLpp.writeToFile((parameters.general_Output_parameters.folderName + DaaLpp.name + "_" + VF::dtostr(time) + ".plt").c_str(), localDiffusionsGrid.L, localDiffusionsGrid.epc, localDiffusionsGrid.alpha);
00576 Output::echo(3, "Pitch-angle diffusion (%e/%e)... ", Daa.max(), DaaLpp.max());
00577 }
00578
00579
00580 psdLocalDiffusions.Diffusion_alpha(parameters.timeStep,
00581 parameters.Lpp[iteration],
00582 Daa, DaaLpp,
00583 localDiffusionsGrid.L, localDiffusionsGrid.pc, localDiffusionsGrid.alpha, localDiffusionsGrid.Jacobian,
00584
00585
00586 alpha_LowerBoundaryCondition.xSlice(0),
00587 alpha_UpperBoundaryCondition.xSlice(0),
00588 alpha_LowerBoundaryCondition.calculationType,
00589 alpha_UpperBoundaryCondition.calculationType);
00590
00591
00592
00593 if (iteration == 1 && parameters.outputModelMatrix) {
00594
00595 }
00596 if (parameters.useAlphaDifusion) Output::echo(3, "done!\n");
00597 }
00598
00599
00600 if (parameters.useEnergyAlphaMixedTerms && localDiffusionsGrid.pc.size > 1 && localDiffusionsGrid.alpha.size > 1) {
00601
00602 if (parameters.useEnergyAlphaMixedTerms) {
00603 if (Dpca.ActivateAndScale(time, parameters.Kp[iteration]))
00604 Dpca.writeToFile((parameters.general_Output_parameters.folderName + Dpca.name + "_" + VF::dtostr(time) + ".plt").c_str(), localDiffusionsGrid.L, localDiffusionsGrid.epc, localDiffusionsGrid.alpha);
00605 if (DpcaLpp.ActivateAndScale(time, parameters.Kp[iteration]))
00606 DpcaLpp.writeToFile((parameters.general_Output_parameters.folderName + DpcaLpp.name + "_" + VF::dtostr(time) + ".plt").c_str(), localDiffusionsGrid.L, localDiffusionsGrid.epc, localDiffusionsGrid.alpha);
00607 Output::echo(3, "Mixed term (Energy-Alpha, explicit) (%e/%e)... \n", Dpca.max(), DpcaLpp.max());
00608 }
00609
00610 psdLocalDiffusions.DiffusionMixTermExplicit(parameters.timeStep,
00611 parameters.Lpp[iteration],
00612 Dpca, DpcaLpp,
00613 localDiffusionsGrid.L, localDiffusionsGrid.pc, localDiffusionsGrid.alpha, localDiffusionsGrid.Jacobian,
00614
00615
00616
00617
00618 pc_LowerBoundaryCondition.xSlice(0),
00619 pc_UpperBoundaryCondition.xSlice(0),
00620 alpha_LowerBoundaryCondition.xSlice(0),
00621 alpha_UpperBoundaryCondition.xSlice(0),
00622 pc_LowerBoundaryCondition.calculationType,
00623 pc_UpperBoundaryCondition.calculationType,
00624 alpha_LowerBoundaryCondition.calculationType,
00625 alpha_UpperBoundaryCondition.calculationType);
00626
00627 if (parameters.useAlphaDifusion) Output::echo(3, "done!\n");
00628 }
00629 } else {
00630 throw error_msg("APPROX_ERROR", "Approximation unknown");
00631 }
00632
00634
00635
00636 psdLocalDiffusions.SourcesAndLosses(
00637 localDiffusionsGrid.L, localDiffusionsGrid.pc, localDiffusionsGrid.alpha,
00638
00639 SL[0],
00640 parameters.timeStep,
00641 parameters.Lpp[iteration],
00642 parameters.tau[iteration],
00643 parameters.tauLpp[iteration]);
00644
00645
00646 if (parameters.NoNegative)
00647 for (il = 0; il < localDiffusionsGrid.L.size; il++) {
00648 for (im = 0; im < localDiffusionsGrid.pc.size; im++) {
00649 for (ia = 0; ia < localDiffusionsGrid.alpha.size; ia++) {
00650 if (psdLocalDiffusions[il][im][ia] < VC::zero_f) {
00651
00652 psdLocalDiffusions[il][im][ia] = VC::zero_f;
00653 }
00654 }
00655 }
00656 }
00657
00658 if (radialDiffusionGrid.L.size > 1) {
00659
00660 Output::echo(3, "Interpolation local->radial... ");
00661
00662
00663
00664
00665 psdRadialDiffusion.Interpolate(psdLocalDiffusions, parameters.interpolation, localDiffusionsGrid, radialDiffusionGrid, pc_LowerBoundaryCondition.xSlice(0), pc_UpperBoundaryCondition.xSlice(0));
00666 Output::echo(3, "done!\n");
00667 }
00668
00669
00670
00671 if (check_time(iteration, parameters.general_Output_parameters.iterStep)) {
00672 Output::echo(1, "Step %d of %d. Time = %d days, %.0f hours. ", iteration, parameters.totalIterationsNumber, int(time), (time-int(time))*24);
00673 Output::echo(1, "Writing data files.\n");
00674
00675 output1D << time << "\t" << parameters.Kp[iteration] << "\t" << parameters.Bf[iteration] << "\t" << parameters.Lpp[iteration] << "\t" << parameters.tau[iteration] << "\t" ;
00676 for (Dxx_it = 0; Dxx_it < Daa.DxxList.size(); Dxx_it++) {
00677 output1D << Daa.DxxList[Dxx_it].Scale(parameters.Kp[iteration]) << "\t";
00678 }
00679 for (Dxx_it = 0; Dxx_it < Dpcpc.DxxList.size(); Dxx_it++) {
00680 output1D << Dpcpc.DxxList[Dxx_it].Scale(parameters.Kp[iteration]) << "\t";
00681 }
00682 for (Dxx_it = 0; Dxx_it < DaaLpp.DxxList.size(); Dxx_it++) {
00683 output1D << DaaLpp.DxxList[Dxx_it].Scale(parameters.Kp[iteration]) << "\t";
00684
00685 }
00686 for (Dxx_it = 0; Dxx_it < DpcpcLpp.DxxList.size(); Dxx_it++) {
00687 output1D << DpcpcLpp.DxxList[Dxx_it].Scale(parameters.Kp[iteration]) << "\t";
00688 }
00689 output1D << endl;
00690
00691 psdLocalDiffusions.Output_without_grid(time);
00692 if (radialDiffusionGrid.L.size > 1) psdRadialDiffusion.Output_without_grid(time);
00693 }
00694
00695 }
00696
00697
00698 std::time(&now);
00699 current_time = std::clock();
00700 Output::echo(0, "End of simulation %s", ctime(&now));
00701 Output::echo(0, "Total time %0.2f sec", ( current_time - zero_time ) / (double)CLOCKS_PER_SEC );
00702
00703 Output::close_log_file();
00704 return 0;
00705
00706 } catch (error_msg err) {
00707 Output::echo(0, "\nCRITICAL ERROR:\n%s\n", err.what().c_str());
00708 Output::close_log_file();
00709 return 1;
00710 } catch (exception exp) {
00711 Output::echo(0, "\nCRITICAL ERROR:\n%s\n", exp.what());
00712 Output::close_log_file();
00713 return 1;
00714 }
00715
00716
00717
00718
00719
00720 Output::close_log_file();
00721 return 0;
00722 }