VERB_code_2.3
Parameters.cpp
Go to the documentation of this file.
1 
9 #include "Parameters.h"
10 
11 #include "../VariousFunctions/variousFunctions.h"
12 #include "../Exceptions/error.h"
13 #include "../Logging/Output.h"
14 #include "../VariousFunctions/erf.h"
15 #include <iostream>
16 #include <iomanip> // for setpersidion
17 
18 using namespace std;
19 
20 
24 template <typename T>
25 bool ReadFromFile(T &variable, ifstream &parameters_file, string variable_name, string default_value = "") {
26 
27  Output::echo(0, "%s = ", variable_name.c_str());
28 
29  // Reset read pointer to beginning of file
30  parameters_file.clear();
31  parameters_file.seekg(0, ios::beg);
32 
33  string word1="", word2="";
34  string string_value;
35  if (!parameters_file.eof()) {
36  while (parameters_file) {
37  word2 = word1;
38  // read one word
39  parameters_file >> word1;
40  // if last word (word2 equals previous word1) is "variableName" then read next word for the variable value
41  if (word1 == "=") {
42  if (word2 == variable_name) {
43  parameters_file >> string_value;
44  break;
45  }
46  }
47  }
48  }
49 
50  ostringstream ostr; // temporary
51  if (string_value != "") {
52  StrToVal(string_value, variable);
53  // output whatever was loaded
54  ostr << variable;
55  Output::echo(0, "%s\n", ostr.str().c_str());
56  return true;
57  } else if (default_value != "") {
58  StrToVal(default_value, variable);
59  // output whatever was loaded
60  ostr << variable;
61  Output::echo(0, "%s (default value)\n", ostr.str().c_str());
62  return true;
63  } else {
64  throw error_msg("ERROR", "Parameter value %s was not found.\n", variable_name.c_str());
65  }
66 }
67 
68 
73 
74  ifstream parameters_file;
75  parameters_file.open(filename.c_str());
76  if (parameters_file == NULL) {
77  throw error_msg("CONST_FILE_OPEN_ERR", "Error reading %s for parameters.", filename.c_str());
78  }
79 
80  ReadFromFile(outputLvl, parameters_file, "outputLvl", "0");
81 
82  // look into map for value for "nDays" and store it into nDays
83  ReadFromFile(nDays, parameters_file, "nDays");
84  // look into map for value for "timeStep" and store it into timeStep
85  ReadFromFile(timeStep, parameters_file, "timeStep");
86  // convert hours into days
87  timeStep = timeStep/24.0;
88  // calculating total number of iterations
89  totalIterationsNumber = int(nDays/timeStep)+1;
90 
91  ReadFromFile(useRadialDiffusion, parameters_file, "useRadialDiffusion", "true");
92  ReadFromFile(useAlphaDiffusion, parameters_file, "useAlphaDiffusion", "true");
93  ReadFromFile(useEnergyDiffusion, parameters_file, "useEnergyDiffusion", "true");
94  ReadFromFile(useEnergyAlphaMixedTerms, parameters_file, "usePcAlphaMixedTerms", "true");
95 
96  if (useRadialDiffusion) {
97  // if 4-th parameter exists - it is default value which is used if we didn't find any value inside the map
98  ReadFromFile(DLLType, parameters_file, "DLLType", "DLLT_B");
99  }
100 
101  Kp.AllocateMemory(totalIterationsNumber);
102  ReadFromFile(useKp, parameters_file, "useKp", string("file"));
103  if (useKp == "file") {
104  ReadFromFile(fileKp, parameters_file, "fileKp");
105  load_1d(Kp, fileKp, timeStep);
106  } else {
107  ReadFromFile(constKp, parameters_file, "constKp");
108  Kp = constKp;
109  }
110 
111  ReadFromFile(useBf, parameters_file, "useBf", string("file"));
112  Bf.AllocateMemory(totalIterationsNumber);
113  // check for useBf parameter
114  if (useBf == "file") {
115  // we need to read values from file. Getting filename from parameters
116  ReadFromFile(fileBf, parameters_file, "fileBf");
117  // loading that 1d file and interpolating in on out time-axis
118  load_1d(Bf, fileBf, timeStep);
119  } else if (useBf == "constant") {
120  // if Bf is constant - using constBf parameter from map
121  ReadFromFile(constBf, parameters_file, "constBf");
122  Bf = constBf;
123  } else {
124  Bf = 1;
125  }
126  for (int i = 0; i < totalIterationsNumber; i++) {
127  if (Bf[i] <= 0) {
128  throw error_msg("BF_ERROR", "Negative or zero Bf[time[%d]=%f]=%f after loading and interpolation.\n", i, timeStep*i, Bf[i]);
129  }
130  }
131 
132 
133  ReadFromFile(useLpp, parameters_file, "useLpp", string("No"));
134  Lpp.AllocateMemory(totalIterationsNumber);
135  if (useLpp == "file") {
136  ReadFromFile(fileLpp, parameters_file, "fileLpp", string(""));
137  load_1d(Lpp, fileLpp, timeStep);
138  } else if (useLpp == "constant") {
139  ReadFromFile(constLpp, parameters_file, "constLpp", "0.0");
140  Lpp = constLpp;
141  } else if (useLpp == "calculate" || str2bool(useLpp) == true) {
142  // using Kp24 - the max Kp index past 24 hours for calculate the location of the plasma pause
143  double Kp24 = 0;
144  int iteration, back_it;
145  for (iteration = 0; iteration < totalIterationsNumber; iteration++) {
146  for (back_it = iteration; (fabs(timeStep * back_it - timeStep * iteration) <= 1.0) && (back_it >= 0); back_it--) {
147  Kp24 = (Kp24 > Kp[back_it]) ? Kp24 : Kp[back_it];
148  }
149  Lpp[iteration] = (5.6 - 0.46*Kp24);
150  Kp24 = 0;
151  }
152  } else {
153  Lpp = 0;
154  }
155 
156  //string usetau, usetauLpp;
157 
158  ReadFromFile(usetau, parameters_file, "usetau", string("constant"));
159  tau.AllocateMemory(totalIterationsNumber);
160  if (usetau == "x/Kp") {
161  double tmp_tau;
162  ReadFromFile(tmp_tau, parameters_file, "tau");
163  for (int i = 0; i < totalIterationsNumber; i++) {
164  tau[i] = tmp_tau/max(Kp[i], 0.01);
165  }
166  // tau = Kp/Kp/Kp * tmp_tau;
167  } if (str2bool(usetau) == true || usetau == "constant") {
168  double tmp_tau;
169  ReadFromFile(tmp_tau, parameters_file, "tau");
170  tau = tmp_tau;
171  } else {
172  tau = 1e99;
173  }
174  ReadFromFile(usetauLpp, parameters_file, "usetauLpp", string("constant"));
175  tauLpp.AllocateMemory(totalIterationsNumber);
176  if (str2bool(usetauLpp) == true || usetauLpp == "constant") {
177  double tmp_tauLpp;
178  ReadFromFile(tmp_tauLpp, parameters_file, "tauLpp");
179  tauLpp = tmp_tauLpp;
180  } else {
181  tauLpp = 1e99;
182  }
183 
184  ReadFromFile(outputModelMatrix, parameters_file, "outputModelMatrix", "false");
185  ReadFromFile(NoNegative, parameters_file, "NoNegative", "false");
186  ReadFromFile(useLossCone, parameters_file, "useLossCone", "true");
187 
188  // Read output parameters from map by calling function, which read output parameters from the map
189  ReadFromFile(general_Output_parameters.timeStep, parameters_file, "general_Output_parameters.timeStep", "1");
190  ReadFromFile(general_Output_parameters.iterStep, parameters_file, "general_Output_parameters.iterStep", "1");
191  // change dt from hours into days
192  general_Output_parameters.timeStep = general_Output_parameters.timeStep/24.0;
193  ReadFromFile(general_Output_parameters.logFileName, parameters_file, "general_Output_parameters.logFileName", "logfile.log");
194  ReadFromFile(general_Output_parameters.folderName, parameters_file, "general_Output_parameters.folderName", "./");
195  ReadFromFile(general_Output_parameters.fileName1D, parameters_file, "general_Output_parameters.fileName1D", "Out1D.plt");
196 
197  // Grid
198  ReadFromFile(radialDiffusionGrid_type, parameters_file, "radialDiffusionGrid_type");
199  ReadFromFile(localDiffusionsGrid_type, parameters_file, "localDiffusionsGrid_type");
200 
201  if (radialDiffusionGrid_type == "GT_FILE")
202  ReadFromFile(radialDiffusionGrid_filename, parameters_file, "radialDiffusionGrid_filename");
203  if (localDiffusionsGrid_type == "GT_FILE")
204  ReadFromFile(localDiffusionsGrid_filename, parameters_file, "localDiffusionsGrid_filename");
205 
206  ReadFromFile(localDiffusionsGrid_L.name, parameters_file, "localDiffusionsGrid_L.name", string("L"));
207  ReadFromFile(localDiffusionsGrid_L.size, parameters_file, "localDiffusionsGrid_L.size");
208  ReadFromFile(localDiffusionsGrid_L.useLogScale, parameters_file, "localDiffusionsGrid_L.useLogScale");
209  ReadFromFile(localDiffusionsGrid_L.min, parameters_file, "localDiffusionsGrid_L.min");
210  ReadFromFile(localDiffusionsGrid_L.max, parameters_file, "localDiffusionsGrid_L.max");
211 
212 
213  ReadFromFile(localDiffusionsGrid_epc.name, parameters_file, "localDiffusionsGrid_epc.name", string("Energy, MeV"));
214  ReadFromFile(localDiffusionsGrid_epc.size, parameters_file, "localDiffusionsGrid_epc.size");
215  ReadFromFile(localDiffusionsGrid_epc.useLogScale, parameters_file, "localDiffusionsGrid_epc.useLogScale");
216  ReadFromFile(localDiffusionsGrid_epc.min, parameters_file, "localDiffusionsGrid_epc.min");
217  ReadFromFile(localDiffusionsGrid_epc.max, parameters_file, "localDiffusionsGrid_epc.max");
218 
219  localDiffusionsGrid_pc = localDiffusionsGrid_epc;
220  localDiffusionsGrid_pc.min = VF::pfunc(localDiffusionsGrid_epc.min);
221  localDiffusionsGrid_pc.max = VF::pfunc(localDiffusionsGrid_epc.max);
222 
223  ReadFromFile(localDiffusionsGrid_alpha.name, parameters_file, "localDiffusionsGrid_alpha.name", string("Pitch angle, deg"));
224  ReadFromFile(localDiffusionsGrid_alpha.size, parameters_file, "localDiffusionsGrid_alpha.size");
225  ReadFromFile(localDiffusionsGrid_alpha.useLogScale, parameters_file, "localDiffusionsGrid_alpha.useLogScale");
226  ReadFromFile(localDiffusionsGrid_alpha.min, parameters_file, "localDiffusionsGrid_alpha.min");
227  ReadFromFile(localDiffusionsGrid_alpha.max, parameters_file, "localDiffusionsGrid_alpha.max");
228 
229  // change alpha from graduss to radians
230  localDiffusionsGrid_alpha.min = localDiffusionsGrid_alpha.min*VC::pi/180;
231  localDiffusionsGrid_alpha.max = localDiffusionsGrid_alpha.max*VC::pi/180;
232 
233 
234  // copy grid parameters from the grid for the local proceses to the grid for the radial diffusion
235  // cause most of them are the same.
236  radialDiffusionGrid_L = localDiffusionsGrid_L;
237  radialDiffusionGrid_epc = localDiffusionsGrid_epc;
238  radialDiffusionGrid_alpha = localDiffusionsGrid_alpha;
239  radialDiffusionGrid_pc = localDiffusionsGrid_pc;
240 
241  // It's repeated 6 times, cause Yuri doesn't understand it the other way
242  ReadFromFile(L_LowerBoundaryCondition.type, parameters_file, "L_LowerBoundaryCondition.type", "BCT_CONSTANT_PSD");
243  if (L_LowerBoundaryCondition.type == "BCT_FILE"
244  || L_LowerBoundaryCondition.type == "BCT_FILE_GRID"
245  || L_LowerBoundaryCondition.type == "BCT_MULTIPLE_FILES") {
246  ReadFromFile(L_LowerBoundaryCondition.filename, parameters_file, "L_LowerBoundaryCondition.filename");
247  } else {
248  ReadFromFile(L_LowerBoundaryCondition.value, parameters_file, "L_LowerBoundaryCondition.value", "1e-9");
249  }
250 
251  ReadFromFile(L_UpperBoundaryCondition.type, parameters_file, "L_UpperBoundaryCondition.type", "BCT_CONSTANT_PSD");
252  if (L_UpperBoundaryCondition.type == "BCT_FILE"
253  || L_UpperBoundaryCondition.type == "BCT_FILE_GRID"
254  || L_UpperBoundaryCondition.type == "BCT_MULTIPLE_FILES") {
255  ReadFromFile(L_UpperBoundaryCondition.filename, parameters_file, "L_UpperBoundaryCondition.filename");
256  } else {
257  ReadFromFile(L_UpperBoundaryCondition.value, parameters_file, "L_UpperBoundaryCondition.value", "1e-9");
258  }
259 
260  ReadFromFile(pc_LowerBoundaryCondition.type, parameters_file, "pc_LowerBoundaryCondition.type", "BCT_CONSTANT_PSD");
261  if (pc_LowerBoundaryCondition.type == "BCT_FILE"
262  || pc_LowerBoundaryCondition.type == "BCT_FILE_GRID"
263  || pc_LowerBoundaryCondition.type == "BCT_MULTIPLE_FILES") {
264  ReadFromFile(pc_LowerBoundaryCondition.filename, parameters_file, "pc_LowerBoundaryCondition.filename");
265  } else {
266  ReadFromFile(pc_LowerBoundaryCondition.value, parameters_file, "pc_LowerBoundaryCondition.value", "1e-9");
267  }
268 
269  ReadFromFile(pc_UpperBoundaryCondition.type, parameters_file, "pc_UpperBoundaryCondition.type", "BCT_CONSTANT_PSD");
270  if (pc_UpperBoundaryCondition.type == "BCT_FILE"
271  || pc_UpperBoundaryCondition.type == "BCT_FILE_GRID"
272  || pc_UpperBoundaryCondition.type == "BCT_MULTIPLE_FILES") {
273  ReadFromFile(pc_UpperBoundaryCondition.filename, parameters_file, "pc_UpperBoundaryCondition.filename");
274  } else {
275  ReadFromFile(pc_UpperBoundaryCondition.value, parameters_file, "pc_UpperBoundaryCondition.value", "1e-9");
276  }
277 
278  ReadFromFile(alpha_LowerBoundaryCondition.type, parameters_file, "alpha_LowerBoundaryCondition.type", "BCT_CONSTANT_PSD");
279  if (alpha_LowerBoundaryCondition.type == "BCT_FILE"
280  || alpha_LowerBoundaryCondition.type == "BCT_FILE_GRID"
281  || alpha_LowerBoundaryCondition.type == "BCT_MULTIPLE_FILES") {
282  ReadFromFile(alpha_LowerBoundaryCondition.filename, parameters_file, "alpha_LowerBoundaryCondition.filename");
283  } else {
284  ReadFromFile(alpha_LowerBoundaryCondition.value, parameters_file, "alpha_LowerBoundaryCondition.value", "1e-9");
285  }
286 
287  ReadFromFile(alpha_UpperBoundaryCondition.type, parameters_file, "alpha_UpperBoundaryCondition.type", "BCT_CONSTANT_PSD");
288  if (alpha_UpperBoundaryCondition.type == "BCT_FILE"
289  || alpha_UpperBoundaryCondition.type == "BCT_FILE_GRID"
290  || alpha_UpperBoundaryCondition.type == "BCT_MULTIPLE_FILES") {
291  ReadFromFile(alpha_UpperBoundaryCondition.filename, parameters_file, "alpha_UpperBoundaryCondition.filename");
292  } else {
293  ReadFromFile(alpha_UpperBoundaryCondition.value, parameters_file, "alpha_UpperBoundaryCondition.value", "1e-9");
294  }
295 
296  ReadFromFile(psdRadialDiffusion.initial_PSD_Kp0, parameters_file, "psdRadialDiffusion.initial_PSD_Kp0", "0.0");
297  // write down first Kp0 for initializing psd
298  if (psdRadialDiffusion.initial_PSD_Kp0 <= 1e-99) {
299  psdRadialDiffusion.initial_PSD_Kp0 = Kp[0];
300  Output::echo(0, "psdRadialDiffusion.initial_PSD_Kp0 set to %f.\n", psdRadialDiffusion.initial_PSD_Kp0);
301  }
302 
303  // Radial diffusion parameters
304  ReadFromFile(psdRadialDiffusion.initial_PSD_Type, parameters_file, "psdRadialDiffusion.initial_PSD_Type", "IPSDT_STEADY_STATE");
305  if (psdRadialDiffusion.initial_PSD_Type == "IPSDT_STEADY_STATE" || psdRadialDiffusion.initial_PSD_Type == "IPSDT_STEADY_STATE_L7"
306  || psdRadialDiffusion.initial_PSD_Type == "IPSDT_STEADY_STATE_FROM_BOUNDARY") {
307  ReadFromFile(psdRadialDiffusion.initial_PSD_tauSteadyState, parameters_file, "psdRadialDiffusion.initial_PSD_tauSteadyState", "0.0");
308  if (psdRadialDiffusion.initial_PSD_tauSteadyState <= 1.e-98) {
309  psdRadialDiffusion.initial_PSD_tauSteadyState = 4.0/psdRadialDiffusion.initial_PSD_Kp0;
310  Output::echo(0, "psdRadialDiffusion.initial_PSD_tauSteadyState set to %f.\n", psdRadialDiffusion.initial_PSD_tauSteadyState);
311  }
312 
314  } else if (psdRadialDiffusion.initial_PSD_Type == "IPSDT_STEADY_STATE_TWO_ZONE") {
315  // Lpp for two zone steady state
316  ReadFromFile(psdRadialDiffusion.initial_PSD_tauSteadyState, parameters_file, "psdRadialDiffusion.initial_PSD_tauSteadyState", "0.0");
317  if (psdRadialDiffusion.initial_PSD_tauSteadyState <= 1.e-98) {
318  psdRadialDiffusion.initial_PSD_tauSteadyState = 4.0/psdRadialDiffusion.initial_PSD_Kp0;
319  Output::echo(0, "psdRadialDiffusion.initial_PSD_tauSteadyState set to %f.\n", psdRadialDiffusion.initial_PSD_tauSteadyState);
320  }
321  // tauLpp for two zone steady state
322  ReadFromFile(psdRadialDiffusion.usetauLpp_SteadyState, parameters_file, "psdRadialDiffusion.usetauLpp_SteadyState", string("combined"));
323  if (psdRadialDiffusion.usetauLpp_SteadyState == "constant") {
324  ReadFromFile(psdRadialDiffusion.initial_PSD_tauLppSteadyState, parameters_file, "psdRadialDiffusion.initial_PSD_tauLppSteadyState");
325  }
327  } else if (psdRadialDiffusion.initial_PSD_Type == "IPSDT_FILE"
328  || psdRadialDiffusion.initial_PSD_Type == "IPSDT_FILE_GRID"
329  || psdRadialDiffusion.initial_PSD_Type == "IPSDT_ORBIT_FLUX_2D"){
331  //|| psdRadialDiffusion.initial_PSD_Type == "IPSDT_FLUX_2D_MAWELLIAN"){
333  ReadFromFile(psdRadialDiffusion.initial_PSD_fileName, parameters_file, "psdRadialDiffusion.initial_PSD_fileName");
334  }
335 
336 
337  ReadFromFile(psdRadialDiffusion.initial_PSD_some_constant_value, parameters_file, "psdRadialDiffusion.initial_PSD_some_constant_value", "1e-99"); // value_psd
338  ReadFromFile(psdRadialDiffusion.initial_PSD_J_L7_function, parameters_file, "psdRadialDiffusion.initial_PSD_J_L7_function", "J_L7_corrected");
339 
340  ReadFromFile(psdRadialDiffusion.initial_PSD_outer_psd, parameters_file, "psdRadialDiffusion.initial_PSD_outer_psd", "1.0");
341  ReadFromFile(psdRadialDiffusion.initial_PSD_inner_psd, parameters_file, "psdRadialDiffusion.initial_PSD_inner_psd", "0.0");
342 
343  ReadFromFile(psdRadialDiffusion.output_PSD_folderName, parameters_file, "psdRadialDiffusion.output_PSD_folderName", general_Output_parameters.folderName );
344  ReadFromFile(psdRadialDiffusion.output_PSD_fileName4D, parameters_file, "psdRadialDiffusion.output_PSD_fileName4D", "OutPSD_rad.dat");
345 
346  ReadFromFile(psdRadialDiffusion.output_PSD_timeStep, parameters_file, "psdRadialDiffusion.output_PSD_timeStep", "1");
347  // convert hours into days
348  psdRadialDiffusion.output_PSD_timeStep = psdRadialDiffusion.output_PSD_timeStep/24.0;
349 
350  psdRadialDiffusion.approximationMethod = "AM_Split_C";
351  psdRadialDiffusion.solutionMethod = "SM_Tridiag";
352 
353  // Local diffusions parameters
354  psdLocalDiffusions.initial_PSD_Type = "IPSDT_CONSTANT";
355  psdLocalDiffusions.initial_PSD_some_constant_value = psdRadialDiffusion.initial_PSD_some_constant_value; // value_psd
356 
357  ReadFromFile(psdLocalDiffusions.output_PSD_folderName, parameters_file, "psdLocalDiffusions.output_PSD_folderName", general_Output_parameters.folderName);
358  ReadFromFile(psdLocalDiffusions.output_PSD_fileName4D, parameters_file, "psdLocalDiffusions.output_PSD_fileName4D", string("OutPSD.dat"));
359 
360  ReadFromFile(psdLocalDiffusions.output_PSD_timeStep, parameters_file, "psdLocalDiffusions.output_PSD_timeStep", "1");
361  // convert hours into days
362  psdLocalDiffusions.output_PSD_timeStep = psdLocalDiffusions.output_PSD_timeStep/24.0;
363 
364  ReadFromFile(psdLocalDiffusions.approximationMethod, parameters_file, "psdLocalDiffusions.approximationMethod", "AM_Split_C");
365  if (psdLocalDiffusions.approximationMethod == "AM_Split_C" || psdLocalDiffusions.approximationMethod == "AM_Split_LR") {
366  ReadFromFile(psdLocalDiffusions.solutionMethod, parameters_file, "psdLocalDiffusions.solutionMethod", "SM_Tridiag");
367  } else {
368  ReadFromFile(psdLocalDiffusions.solutionMethod, parameters_file, "psdLocalDiffusions.solutionMethod", "SM_Relaxation");
369 
370  if (psdLocalDiffusions.solutionMethod == "SM_GMRES" || psdLocalDiffusions.solutionMethod == "SM_GMRES2") {
371  ReadFromFile(psdLocalDiffusions.GMRES_parameters.SOL_maxiter, parameters_file, "psdLocalDiffusions.GMRES_parameters.SOL_maxiter", "100");
372  ReadFromFile(psdLocalDiffusions.GMRES_parameters.SOL_i_max, parameters_file, "psdLocalDiffusions.GMRES_parameters.SOL_i_max", "25");
373  ReadFromFile(psdLocalDiffusions.GMRES_parameters.SOL_max_iter_err, parameters_file, "psdLocalDiffusions.GMRES_parameters.SOL_max_iter_err", "1e-1");
374  ReadFromFile(psdLocalDiffusions.GMRES_parameters.preconditioner_type, parameters_file, "psdLocalDiffusions.GMRES_parameters.preconditioner_type", "none");
375  ReadFromFile(psdLocalDiffusions.GMRES_parameters.use_normalization, parameters_file, "psdLocalDiffusions.GMRES_parameters.use_normalization", "Yes");
376  }
377  }
378 
379  // Sources and losses parameters
380  ReadFromFile(SL.SL_L_top, parameters_file, "SL.SL_L_top", "false");
381  if (SL.SL_L_top) {
382  ReadFromFile(SL.SL_L_top_filename, parameters_file, "SL.SL_L_top_filename");
383  }
384  ReadFromFile(SL.SL_E_min, parameters_file, "SL.SL_E_min", "false");
385  if (SL.SL_E_min) {
386  ReadFromFile(SL.SL_E_min_filename, parameters_file, "SL.SL_E_min_filename");
387  }
388 
389  // Interpolation parameters
390  ReadFromFile(interpolation.type, parameters_file, "interpolation.type", "IT_SPLINE");
391  ReadFromFile(interpolation.useLog, parameters_file, "interpolation.useLog", "Log_10");
392  ReadFromFile(interpolation.linearSplineCoef, parameters_file, "interpolation.linearSplineCoef", "1");
393  ReadFromFile(interpolation.maxSecondDerivative, parameters_file, "interpolation.maxSecondDerivative", "1e2");
394 
395  //ReadFromFile(DxxParameters_structureList, parameters_file, "CoeffFileName");
396 
397  // Checks
398 
399  // Check outer L boundry condition creation parameters
400  if (psdRadialDiffusion.initial_PSD_Type == "IPSDT_STEADY_STATE_FROM_BOUNDARY"
401  && (L_UpperBoundaryCondition.type != "BCT_FILE"
402  && L_UpperBoundaryCondition.type != "BCT_FILE_GRID"
403  && L_UpperBoundaryCondition.type != "BCT_MULTIPLE_FILES")
404  ) throw error_msg("PARAMETERS_ERROR", "Outer L boundary can not be used for PSD steady state, cause it is set to be created from PSD itself.");
405 
406  if ((psdLocalDiffusions.approximationMethod == "AM_Block_LR" || psdLocalDiffusions.approximationMethod == "AM_Block_C") && !(localDiffusionsGrid_pc.size > 1 && localDiffusionsGrid_alpha.size > 1)) {
407  throw error_msg("BLOCK_WITH_1D", "Block method can not be used with 1D");
408  }
409 
410  if (localDiffusionsGrid_pc.size < 3 && useEnergyDiffusion)
411  throw error_msg("ERROR", "Not enough grid points for energy diffusion (localDiffusionsGrid.pc.size = %d and useEnergyDiffusion = Yes)", localDiffusionsGrid_pc.size);
412 
413  if (localDiffusionsGrid_alpha.size < 3 && useAlphaDiffusion)
414  throw error_msg("ERROR", "Not enough grid points for pitch angle diffusion (localDiffusionsGrid.alpha.size = %d and useAlphaDiffusion = Yes)", localDiffusionsGrid_alpha.size);
415 
416  if (localDiffusionsGrid_L.size < 3 && useRadialDiffusion)
417  throw error_msg("ERROR", "Not enough grid points for radial diffusion (localDiffusionsGrid.L.size = %d and useRadialDiffusion = Yes)", localDiffusionsGrid_L.size);
418 
419  // Read diffusion coefficients
420  parameters_file.seekg(0,ios::beg);
421  string word1="", word2="";
422  string string_value, dxx_parameters_file_name;
423  if (parameters_file != NULL) {
424  while (parameters_file) {
425  word2 = word1;
426  // read one word
427  parameters_file >> word1;
428  // if last word (word2 equals previous word1) is "variableName" then read next word for the variable value
429  if (word1 == "=") {
430  if (word2 == "CoeffFileName") {
431  parameters_file >> dxx_parameters_file_name;
432 
433  // add a value to the vector of diffusion coefficients parameters structures (well, a structure for a diffusion coefficient)
434  DxxParametersList.push_back(DxxParameters_structure());
435  // Read that diffusion coefficient parameters from it's file
436  // ReadFromFile(DxxParameters_structureList.back());
437 
438  // Load dxx parameters into the just added dxx parameter structure
439  DxxParametersList.back().Load_dxx_parameters(dxx_parameters_file_name);
440 
441  }
442  }
443  }
444  }
445 
446 
447  ReadFromFile(SL.useLmp, parameters_file, "sourcesAndLosses.useLmp", "false");
448  if (SL.useLmp) {
449  ReadFromFile(SL.fileLmp, parameters_file, "sourcesAndLosses.fileLmp");
450  }
451 
452  if (SL.useLmp == true){
453 
454  // check localDiffusionsGrid_alpha.size
455  if (localDiffusionsGrid_alpha.size < 1 ){
456  throw error_msg("ERROR", "Not enough alpha grid points for magnitopause shadowing (localDiffusionsGrid_alpha.size = %d)", localDiffusionsGrid_alpha.size);
457  }
458 
459  // load file Matrix2D Lmp
460  SL.Lmp.AllocateMemory(totalIterationsNumber, localDiffusionsGrid_alpha.size);
461  load_2d(SL.Lmp, SL.fileLmp, timeStep);
462 
463  // // // load_2d(); // // //
464 
465 
466  }
467 
468 
469  return true;
470 }
471 
472 
476 bool DxxParameters_structure::Load_dxx_parameters(string dxx_parameters_file_name) {
477  ifstream dxx_parameters_file;
478  dxx_parameters_file.open(dxx_parameters_file_name.c_str());
479  if (dxx_parameters_file == NULL) {
480  throw error_msg("DXX_PARAM_FILE_OPEN_ERR", "Error reading diffusion coefficients parameters file %s.", dxx_parameters_file_name.c_str());
481  }
482 
483  ReadFromFile(DxxType, dxx_parameters_file, "DxxType");
484  ReadFromFile(DxxName, dxx_parameters_file, "DxxType");
485 
486  ReadFromFile(waveType, dxx_parameters_file, "waveType");
487  ReadFromFile(waveName, dxx_parameters_file, "waveType");
488 
489  ReadFromFile(time_start, dxx_parameters_file, "time_start");
490  ReadFromFile(time_end, dxx_parameters_file, "time_end");
491 
492  ReadFromFile(useScale, dxx_parameters_file, "useScale", "false");
493  if (useScale) {
494  ReadFromFile(DxxKp , dxx_parameters_file, "DxxKp");
495  }
496 
497  ReadFromFile(loadOrCalculate, dxx_parameters_file, "loadOrCalculate");
498 
499  ReadFromFile(filetype, dxx_parameters_file, "filetype");
500  ReadFromFile(filename, dxx_parameters_file, "filename");
501 
502  ReadFromFile(multiplicator, dxx_parameters_file, "multiplicator", "1.0");
503  ReadFromFile(MLT_averaging, dxx_parameters_file, "MLT_averaging", "1.0");
504 
505  if (loadOrCalculate == "LOC_CALCULATE" || loadOrCalculate == "LOC_LOADORCALCULATE") {
506  ReadFromFile(numberDensity, dxx_parameters_file, "numberDensity");
507  ReadFromFile(Omega_mType, dxx_parameters_file, "Omega_mType");
508  ReadFromFile(Omega_m, dxx_parameters_file, "Omega_m");
509  ReadFromFile(d_omega, dxx_parameters_file, "d_omega");
510  ReadFromFile(omega_lc, dxx_parameters_file, "omega_lc");
511  ReadFromFile(omega_uc, dxx_parameters_file, "omega_uc");
512  ReadFromFile(eta1, dxx_parameters_file, "eta1");
513  ReadFromFile(eta2, dxx_parameters_file, "eta2");
514  ReadFromFile(eta3, dxx_parameters_file, "eta3");
515  ReadFromFile(Bw, dxx_parameters_file, "Bw");
516  ReadFromFile(BwFromLambda, dxx_parameters_file, "BwFromLambda");
517 
518  ReadFromFile(nint, dxx_parameters_file, "nint");
519  ReadFromFile(mirror_point_coeff, dxx_parameters_file, "mirror_point_coeff", "0.999");
520 
521  ReadFromFile(lam_min, dxx_parameters_file, "lam_min", "0.0");
522  ReadFromFile(lam_max, dxx_parameters_file, "lam_max");
523  ReadFromFile(nu, dxx_parameters_file, "nu", "0.0");
524  ReadFromFile(s, dxx_parameters_file, "sign");
525  if (numberDensity == "ND_CONSTANT") {
526  ReadFromFile(f, dxx_parameters_file, "f_ratio");
527  }
528  ReadFromFile(particle, dxx_parameters_file, "particle");
529 
530  double sigma = (omega_uc - omega_lc) / (2.0 * d_omega);
531  Output::echo(0, "sigma = %f\n", sigma);
532  nu = erf(sigma)*sqrt(VC::pi);
533  Output::echo(0, "nu = %f\n", nu);
534 
535  if (Bw <= 0 && BwFromLambda == false) {
536  throw error_msg("DXX_ERROR", "BwFromLambda is false and Bw = %f\n", Bw);
537  } else if (BwFromLambda == true && Bw <= 0) {
538  Bw = 1;
539  }
540 
541  }
542 
543  return true;
544 }
545 
546 
550 void StrToVal(string input, double &place) {
551  istringstream(input) >> place;
552 // Output::echo("%e", place);
553 }
554 
558 void StrToVal(string input, int &place) {
559  istringstream(input) >> place;
560 // Output::echo(place);
561 }
562 
566 void StrToVal(string input, string &place) {
567  place = input;
568 // Output::echo(place);
569 }
570 
574 void StrToVal(string input, bool &place) {
575  place = str2bool(input);;
576 // Output::echo(place);
577 }
578 
579 
583 bool str2bool(string str) {
584  if (str == "Yes" || str == "yes" || str == "1" || str == "True" || str == "true") return true;
585  else return false;
586 }
587 
591 string bool2str(bool b) {
592  return (b==true)?"Yes":"No";
593 }
594 
595 
598 void load_1d(Matrix1D<double> &var, string filename, double dt, int var_size) {
599  // var_size is a flag, if it is equal to zero, then we assume memory 1d-matrix allocated
600  // if it is equal to some value, then we need to allocate the memory
601  if (var_size != 0) {
602  // allocate memory
603  var.AllocateMemory(var_size);
604  } else {
605  // if var_size==0, assume memory already allocated
606  var_size = var.size_x;
607  }
608 
609  ifstream input_1d(filename.c_str());
610  if (input_1d != NULL) {
611  vector <double> time1d, value1d;
612  double time_tmp, value_tmp;
613  while (!input_1d.eof()) {
614  // reading time values and function values from the file into vectors
615  input_1d >> time_tmp >> value_tmp;
616  time1d.push_back(time_tmp);
617  value1d.push_back(value_tmp);
618  }
619  double time;
620  double a;
621  // interpolating loaded values into out time-axis
622  for (int it = 0; it < var.size_x; it++) {
623  time = dt*it;
624 
625  unsigned int i = 0;
626  while (((++i) < (time1d.size()-1)) && time > time1d[i]) {
627  // Not enough Kp error here
628  };
629 
630  a = (time - time1d[i-1]) / (time1d[i] - time1d[i-1]);
631  var[it] = value1d[i-1] + a*(value1d[i] - value1d[i-1]);
632  }
633  } else {
634  throw error_msg("LOAD_1D_ERROR", "Load 1d file error: file not found %f\n", filename.c_str());
635  }
636 }
637 
638 
642 void load_2d(Matrix2D<double> &var, string filename, double dt, int var_size_x, int var_size_y) {
643 
644 
645  // vector of time value
646  vector <double> time1d;
647 
648 
649  // vectors of 2d values
650  T_2d_double_vector value2d;
651 
652  // var_size is a flag, if it is equal to zero, then we assume memory 1d-matrix allocated
653  // if it is equal to some value, then we need to allocate the memory
654  if (var_size_x != 0 & var_size_y != 0) {
655  // allocate memory
656  var.AllocateMemory(var_size_x,var_size_y);
657  } /* else {
658  // if var_size==0, assume memory already allocated
659  var_size = var.size_x * var.size_y;
660  }*/
661 
662  ifstream input_2d(filename.c_str());
663 
664  if (input_2d != NULL) {
665 
666  double time_tmp, value_tmp;
667  int it = 0;
668  while (input_2d >> time_tmp) { // reading time values and function values from the file into vectors
669  // This method of reading prevent emty lines reading
670  time1d.push_back(time_tmp);
671 
672  for(int f=0; f < var.size_y; f++){
673  input_2d >> value_tmp;
674 
675  // Add value value_tmp to value2d (vector of changing size)
676  T_2d_double_vector_add(value2d, it, f) = value_tmp;
677  }
678  it++;
679 
680  }
681  double time;
682  double a;
683  // interpolating loaded values into out time-axis
684  for (it = 0; it < var.size_x; it++) {
685  time = dt*it;
686 
687  unsigned int i = 0;
688  while (((++i) < (time1d.size()-1)) && time > time1d[i]) {
689  // Not enough Lmp error here (stupid check)
690  };
691 
692 
693  #ifdef DEBUG_MODE
694  cout << "\n" << setprecision(2) << setw(5) << time << ' ';
695  #endif
696 
697 
698  for (int ia = 0; ia < var.size_y; ia++)
699  {
700  a = (time - time1d[i-1]) / (time1d[i] - time1d[i-1]);
701  var[it][ia] = value2d[i-1][ia] + a*(value2d[i][ia] - value2d[i-1][ia]);
702 
703  #ifdef DEBUG_MODE
704  if (var[it][ia] > 6) cout << 0 << ' ';
705  else cout << 1 << ' ';
706  #endif
707 
708  }
709  }
710  } else {
711  throw error_msg("LOAD_2D_ERROR", "Load 2d file error: file not found %f\n", filename.c_str());
712  }
713 }
714 
715 
723  double& T_2d_double_vector_add(T_2d_double_vector &v, size_t x, size_t y)
724  {
725  if(v.size() <= x) v.resize(x + 1);
726  if(v[x].size() <= y) v[x].resize(y + 1, -1);
727  return v[x][y];
728  }
729 
730 
int size_y
size x, size y
Definition: Matrix.h:141
void load_1d(Matrix1D< double > &var, string filename, double dt, int var_size)
Definition: Parameters.cpp:598
double max(double v1, double v2)
Return maximum.
bool str2bool(string str)
Definition: Parameters.cpp:583
void AllocateMemory(int size_x, int size_y)
Definition: Matrix.cpp:524
double erf(double x)
Returns the error function erf(x).
Definition: erf.cpp:137
string bool2str(bool b)
Definition: Parameters.cpp:591
int size_x
size x, size y
Definition: Matrix.h:141
int outputLvl
verbose level defined in namespase Output
Definition: Output.cpp:21
void echo(int msglvl, const char *format,...)
Basic function! Write the message to the file.
Definition: Output.cpp:44
Struct that holds various parameters to be used for Dxx.
Definition: Parameters.h:27
void AllocateMemory(int x_size)
Definition: Matrix.cpp:173
vector< T_double_vector > T_2d_double_vector
typedef vector of vector of doubles
Definition: Parameters.h:286
double pfunc(double K)
Computation of moumentum from Kinetic energy .
bool Load_dxx_parameters(string filename)
Definition: Parameters.cpp:476
static const double pi
π
bool Load_parameters(string filename)
Definition: Parameters.cpp:72
bool ReadFromFile(T &variable, ifstream &parameters_file, string variable_name, string default_value="")
Definition: Parameters.cpp:25
void load_2d(Matrix2D< double > &var, string filename, double dt, int var_size_x, int var_size_y)
Definition: Parameters.cpp:642
void StrToVal(string input, double &place)
Definition: Parameters.cpp:550
double T(double alpha)
Function for bounce time.
Error message - stack of single_errors.
Definition: error.h:54
double & T_2d_double_vector_add(T_2d_double_vector &v, size_t x, size_t y)
Definition: Parameters.cpp:723
int size_x
size x
Definition: Matrix.h:46