VERB_code_2.2  2
 All Classes Namespaces Files Functions Variables Typedefs Macros Pages
Parameters.cpp
Go to the documentation of this file.
1 /**
2  * \file Parameters.cpp
3  * \brief Defenition of the struct Parameters_structure, struct DxxParameters_structure and supporting function
4  * Loads parameters from file into map, then use parameters-structure related functions to get parameters values from the map and store them into the parameters structures.
5  *
6  * \author Developed by Dmitry Subbotin under supervision of the PI Yuri Shprits
7  */
8 
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 
17 using namespace std;
18 
19 
20 /**
21  * Read each parameter from parameters file
22  */
23 template <typename T>
24 bool ReadFromFile(T &variable, ifstream &parameters_file, string variable_name, string default_value = "") {
25 
26  Output::echo(0, "%s = ", variable_name.c_str());
27 
28  // Reset read pointer to beginning of file
29  parameters_file.clear();
30  parameters_file.seekg(0, ios::beg);
31 
32  string word1="", word2="";
33  string string_value;
34  if (!parameters_file.eof()) {
35  while (parameters_file) {
36  word2 = word1;
37  // read one word
38  parameters_file >> word1;
39  // if last word (word2 equals previous word1) is "variableName" then read next word for the variable value
40  if (word1 == "=") {
41  if (word2 == variable_name) {
42  parameters_file >> string_value;
43  break;
44  }
45  }
46  }
47  }
48 
49  ostringstream ostr; // temporary
50  if (string_value != "") {
51  StrToVal(string_value, variable);
52  // output whatever was loaded
53  ostr << variable;
54  Output::echo(0, "%s\n", ostr.str().c_str());
55  return true;
56  } else if (default_value != "") {
57  StrToVal(default_value, variable);
58  // output whatever was loaded
59  ostr << variable;
60  Output::echo(0, "%s (default value)\n", ostr.str().c_str());
61  return true;
62  } else {
63  throw error_msg("ERROR", "Parameter value %s was not found.\n", variable_name.c_str());
64  }
65 }
66 
67 
68 /**
69  * Loads parameters from file
70  */
72 
73  ifstream parameters_file;
74  parameters_file.open(filename.c_str());
75  if (parameters_file == NULL) {
76  throw error_msg("CONST_FILE_OPEN_ERR", "Error reading %s for parameters.", filename.c_str());
77  }
78 
79  ReadFromFile(outputLvl, parameters_file, "outputLvl", "0");
80 
81  // look into map for value for "nDays" and store it into nDays
82  ReadFromFile(nDays, parameters_file, "nDays");
83  // look into map for value for "timeStep" and store it into timeStep
84  ReadFromFile(timeStep, parameters_file, "timeStep");
85  // convert hours into days
86  timeStep = timeStep/24.0;
87  // calculating total number of iterations
88  totalIterationsNumber = int(nDays/timeStep)+1;
89 
90  ReadFromFile(useRadialDiffusion, parameters_file, "useRadialDiffusion", "true");
91  ReadFromFile(useAlphaDiffusion, parameters_file, "useAlphaDiffusion", "true");
92  ReadFromFile(useEnergyDiffusion, parameters_file, "useEnergyDiffusion", "true");
93  ReadFromFile(useEnergyAlphaMixedTerms, parameters_file, "usePcAlphaMixedTerms", "true");
94 
95  if (useRadialDiffusion) {
96  // if 4-th parameter exists - it is default value which is used if we didn't find any value inside the map
97  ReadFromFile(DLLType, parameters_file, "DLLType", "DLLT_B");
98  }
99 
100  Kp.AllocateMemory(totalIterationsNumber);
101  ReadFromFile(useKp, parameters_file, "useKp", string("file"));
102  if (useKp == "file") {
103  ReadFromFile(fileKp, parameters_file, "fileKp");
104  load_1d(Kp, fileKp, timeStep);
105  } else {
106  ReadFromFile(constKp, parameters_file, "constKp");
107  Kp = constKp;
108  }
109 
110  ReadFromFile(useBf, parameters_file, "useBf", string("file"));
111  Bf.AllocateMemory(totalIterationsNumber);
112  // check for useBf parameter
113  if (useBf == "file") {
114  // we need to read values from file. Getting filename from parameters
115  ReadFromFile(fileBf, parameters_file, "fileBf");
116  // loading that 1d file and interpolating in on out time-axis
117  load_1d(Bf, fileBf, timeStep);
118  } else if (useBf == "constant") {
119  // if Bf is constant - using constBf parameter from map
120  ReadFromFile(constBf, parameters_file, "constBf");
121  Bf = constBf;
122  } else {
123  Bf = 1;
124  }
125  for (int i = 0; i < totalIterationsNumber; i++) {
126  if (Bf[i] <= 0) {
127  throw error_msg("BF_ERROR", "Negative or zero Bf[time[%d]=%f]=%f after loading and interpolation.\n", i, timeStep*i, Bf[i]);
128  }
129  }
130 
131 
132  ReadFromFile(useLpp, parameters_file, "useLpp", string("No"));
133  Lpp.AllocateMemory(totalIterationsNumber);
134  if (useLpp == "file") {
135  ReadFromFile(fileLpp, parameters_file, "fileLpp", string(""));
136  load_1d(Lpp, fileLpp, timeStep);
137  } else if (useLpp == "constant") {
138  ReadFromFile(constLpp, parameters_file, "constLpp", "0.0");
139  Lpp = constLpp;
140  } else if (useLpp == "calculate" || str2bool(useLpp) == true) {
141  // using Kp24 - the max Kp index past 24 hours for calculate the location of the plasma pause
142  double Kp24 = 0;
143  int iteration, back_it;
144  for (iteration = 0; iteration < totalIterationsNumber; iteration++) {
145  for (back_it = iteration; (fabs(timeStep * back_it - timeStep * iteration) <= 1.0) && (back_it >= 0); back_it--) {
146  Kp24 = (Kp24 > Kp[back_it]) ? Kp24 : Kp[back_it];
147  }
148  Lpp[iteration] = (5.6 - 0.46*Kp24);
149  Kp24 = 0;
150  }
151  } else {
152  Lpp = 0;
153  }
154 
155  //string usetau, usetauLpp;
156 
157  ReadFromFile(usetau, parameters_file, "usetau", string("constant"));
158  tau.AllocateMemory(totalIterationsNumber);
159  if (usetau == "x/Kp") {
160  double tmp_tau;
161  ReadFromFile(tmp_tau, parameters_file, "tau");
162  for (int i = 0; i < totalIterationsNumber; i++) {
163  tau[i] = tmp_tau/max(Kp[i], 0.01);
164  }
165  // tau = Kp/Kp/Kp * tmp_tau;
166  } if (str2bool(usetau) == true || usetau == "constant") {
167  double tmp_tau;
168  ReadFromFile(tmp_tau, parameters_file, "tau");
169  tau = tmp_tau;
170  } else {
171  tau = 1e99;
172  }
173  ReadFromFile(usetauLpp, parameters_file, "usetauLpp", string("constant"));
174  tauLpp.AllocateMemory(totalIterationsNumber);
175  if (str2bool(usetauLpp) == true || usetauLpp == "constant") {
176  double tmp_tauLpp;
177  ReadFromFile(tmp_tauLpp, parameters_file, "tauLpp");
178  tauLpp = tmp_tauLpp;
179  } else {
180  tauLpp = 1e99;
181  }
182 
183  ReadFromFile(outputModelMatrix, parameters_file, "outputModelMatrix", "false");
184  ReadFromFile(NoNegative, parameters_file, "NoNegative", "false");
185  ReadFromFile(useLossCone, parameters_file, "useLossCone", "true");
186 
187  // Read output parameters from map by calling function, which read output parameters from the map
188  ReadFromFile(general_Output_parameters.timeStep, parameters_file, "general_Output_parameters.timeStep", "1");
189  ReadFromFile(general_Output_parameters.iterStep, parameters_file, "general_Output_parameters.iterStep", "1");
190  // change dt from hours into days
191  general_Output_parameters.timeStep = general_Output_parameters.timeStep/24.0;
192  ReadFromFile(general_Output_parameters.logFileName, parameters_file, "general_Output_parameters.logFileName", "logfile.log");
193  ReadFromFile(general_Output_parameters.folderName, parameters_file, "general_Output_parameters.folderName", "./");
194  ReadFromFile(general_Output_parameters.fileName1D, parameters_file, "general_Output_parameters.fileName1D", "Out1D.plt");
195 
196  // Grid
197  ReadFromFile(radialDiffusionGrid_type, parameters_file, "radialDiffusionGrid_type");
198  ReadFromFile(localDiffusionsGrid_type, parameters_file, "localDiffusionsGrid_type");
199 
200  if (radialDiffusionGrid_type == "GT_FILE")
201  ReadFromFile(radialDiffusionGrid_filename, parameters_file, "radialDiffusionGrid_filename");
202  if (localDiffusionsGrid_type == "GT_FILE")
203  ReadFromFile(localDiffusionsGrid_filename, parameters_file, "localDiffusionsGrid_filename");
204 
205  ReadFromFile(localDiffusionsGrid_L.name, parameters_file, "localDiffusionsGrid_L.name", string("L"));
206  ReadFromFile(localDiffusionsGrid_L.size, parameters_file, "localDiffusionsGrid_L.size");
207  ReadFromFile(localDiffusionsGrid_L.useLogScale, parameters_file, "localDiffusionsGrid_L.useLogScale");
208  ReadFromFile(localDiffusionsGrid_L.min, parameters_file, "localDiffusionsGrid_L.min");
209  ReadFromFile(localDiffusionsGrid_L.max, parameters_file, "localDiffusionsGrid_L.max");
210 
211 
212  ReadFromFile(localDiffusionsGrid_epc.name, parameters_file, "localDiffusionsGrid_epc.name", string("Energy, MeV"));
213  ReadFromFile(localDiffusionsGrid_epc.size, parameters_file, "localDiffusionsGrid_epc.size");
214  ReadFromFile(localDiffusionsGrid_epc.useLogScale, parameters_file, "localDiffusionsGrid_epc.useLogScale");
215  ReadFromFile(localDiffusionsGrid_epc.min, parameters_file, "localDiffusionsGrid_epc.min");
216  ReadFromFile(localDiffusionsGrid_epc.max, parameters_file, "localDiffusionsGrid_epc.max");
217 
218  localDiffusionsGrid_pc = localDiffusionsGrid_epc;
219  localDiffusionsGrid_pc.min = VF::pfunc(localDiffusionsGrid_epc.min);
220  localDiffusionsGrid_pc.max = VF::pfunc(localDiffusionsGrid_epc.max);
221 
222  ReadFromFile(localDiffusionsGrid_alpha.name, parameters_file, "localDiffusionsGrid_alpha.name", string("Pitch angle, deg"));
223  ReadFromFile(localDiffusionsGrid_alpha.size, parameters_file, "localDiffusionsGrid_alpha.size");
224  ReadFromFile(localDiffusionsGrid_alpha.useLogScale, parameters_file, "localDiffusionsGrid_alpha.useLogScale");
225  ReadFromFile(localDiffusionsGrid_alpha.min, parameters_file, "localDiffusionsGrid_alpha.min");
226  ReadFromFile(localDiffusionsGrid_alpha.max, parameters_file, "localDiffusionsGrid_alpha.max");
227 
228  // change alpha from graduss to radians
229  localDiffusionsGrid_alpha.min = localDiffusionsGrid_alpha.min*VC::pi/180;
230  localDiffusionsGrid_alpha.max = localDiffusionsGrid_alpha.max*VC::pi/180;
231 
232 
233  // copy grid parameters from the grid for the local proceses to the grid for the radial diffusion
234  // cause most of them are the same.
235  radialDiffusionGrid_L = localDiffusionsGrid_L;
236  radialDiffusionGrid_epc = localDiffusionsGrid_epc;
237  radialDiffusionGrid_alpha = localDiffusionsGrid_alpha;
238  radialDiffusionGrid_pc = localDiffusionsGrid_pc;
239 
240  // It's repeated 6 times, cause Yuri doesn't understand it the other way
241  ReadFromFile(L_LowerBoundaryCondition.type, parameters_file, "L_LowerBoundaryCondition.type", "BCT_CONSTANT_PSD");
242  if (L_LowerBoundaryCondition.type == "BCT_FILE"
243  || L_LowerBoundaryCondition.type == "BCT_FILE_GRID"
244  || L_LowerBoundaryCondition.type == "BCT_MULTIPLE_FILES") {
245  ReadFromFile(L_LowerBoundaryCondition.filename, parameters_file, "L_LowerBoundaryCondition.filename");
246  } else {
247  ReadFromFile(L_LowerBoundaryCondition.value, parameters_file, "L_LowerBoundaryCondition.value", "1e-9");
248  }
249 
250  ReadFromFile(L_UpperBoundaryCondition.type, parameters_file, "L_UpperBoundaryCondition.type", "BCT_CONSTANT_PSD");
251  if (L_UpperBoundaryCondition.type == "BCT_FILE"
252  || L_UpperBoundaryCondition.type == "BCT_FILE_GRID"
253  || L_UpperBoundaryCondition.type == "BCT_MULTIPLE_FILES") {
254  ReadFromFile(L_UpperBoundaryCondition.filename, parameters_file, "L_UpperBoundaryCondition.filename");
255  } else {
256  ReadFromFile(L_UpperBoundaryCondition.value, parameters_file, "L_UpperBoundaryCondition.value", "1e-9");
257  }
258 
259  ReadFromFile(pc_LowerBoundaryCondition.type, parameters_file, "pc_LowerBoundaryCondition.type", "BCT_CONSTANT_PSD");
260  if (pc_LowerBoundaryCondition.type == "BCT_FILE"
261  || pc_LowerBoundaryCondition.type == "BCT_FILE_GRID"
262  || pc_LowerBoundaryCondition.type == "BCT_MULTIPLE_FILES") {
263  ReadFromFile(pc_LowerBoundaryCondition.filename, parameters_file, "pc_LowerBoundaryCondition.filename");
264  } else {
265  ReadFromFile(pc_LowerBoundaryCondition.value, parameters_file, "pc_LowerBoundaryCondition.value", "1e-9");
266  }
267 
268  ReadFromFile(pc_UpperBoundaryCondition.type, parameters_file, "pc_UpperBoundaryCondition.type", "BCT_CONSTANT_PSD");
269  if (pc_UpperBoundaryCondition.type == "BCT_FILE"
270  || pc_UpperBoundaryCondition.type == "BCT_FILE_GRID"
271  || pc_UpperBoundaryCondition.type == "BCT_MULTIPLE_FILES") {
272  ReadFromFile(pc_UpperBoundaryCondition.filename, parameters_file, "pc_UpperBoundaryCondition.filename");
273  } else {
274  ReadFromFile(pc_UpperBoundaryCondition.value, parameters_file, "pc_UpperBoundaryCondition.value", "1e-9");
275  }
276 
277  ReadFromFile(alpha_LowerBoundaryCondition.type, parameters_file, "alpha_LowerBoundaryCondition.type", "BCT_CONSTANT_PSD");
278  if (alpha_LowerBoundaryCondition.type == "BCT_FILE"
279  || alpha_LowerBoundaryCondition.type == "BCT_FILE_GRID"
280  || alpha_LowerBoundaryCondition.type == "BCT_MULTIPLE_FILES") {
281  ReadFromFile(alpha_LowerBoundaryCondition.filename, parameters_file, "alpha_LowerBoundaryCondition.filename");
282  } else {
283  ReadFromFile(alpha_LowerBoundaryCondition.value, parameters_file, "alpha_LowerBoundaryCondition.value", "1e-9");
284  }
285 
286  ReadFromFile(alpha_UpperBoundaryCondition.type, parameters_file, "alpha_UpperBoundaryCondition.type", "BCT_CONSTANT_PSD");
287  if (alpha_UpperBoundaryCondition.type == "BCT_FILE"
288  || alpha_UpperBoundaryCondition.type == "BCT_FILE_GRID"
289  || alpha_UpperBoundaryCondition.type == "BCT_MULTIPLE_FILES") {
290  ReadFromFile(alpha_UpperBoundaryCondition.filename, parameters_file, "alpha_UpperBoundaryCondition.filename");
291  } else {
292  ReadFromFile(alpha_UpperBoundaryCondition.value, parameters_file, "alpha_UpperBoundaryCondition.value", "1e-9");
293  }
294 
295  ReadFromFile(psdRadialDiffusion.initial_PSD_Kp0, parameters_file, "psdRadialDiffusion.initial_PSD_Kp0", "0.0");
296  // write down first Kp0 for initializing psd
297  if (psdRadialDiffusion.initial_PSD_Kp0 <= 1e-99) {
298  psdRadialDiffusion.initial_PSD_Kp0 = Kp[0];
299  Output::echo(0, "psdRadialDiffusion.initial_PSD_Kp0 set to %f.\n", psdRadialDiffusion.initial_PSD_Kp0);
300  }
301 
302  // Radial diffusion parameters
303  ReadFromFile(psdRadialDiffusion.initial_PSD_Type, parameters_file, "psdRadialDiffusion.initial_PSD_Type", "IPSDT_STEADY_STATE");
304  if (psdRadialDiffusion.initial_PSD_Type == "IPSDT_STEADY_STATE"
305  || psdRadialDiffusion.initial_PSD_Type == "IPSDT_STEADY_STATE_FROM_BOUNDARY") {
306  ReadFromFile(psdRadialDiffusion.initial_PSD_tauSteadyState, parameters_file, "psdRadialDiffusion.initial_PSD_tauSteadyState", "0.0");
307  if (psdRadialDiffusion.initial_PSD_tauSteadyState <= 1.e-98) {
308  psdRadialDiffusion.initial_PSD_tauSteadyState = 4.0/psdRadialDiffusion.initial_PSD_Kp0;
309  Output::echo(0, "psdRadialDiffusion.initial_PSD_tauSteadyState set to %f.\n", psdRadialDiffusion.initial_PSD_tauSteadyState);
310  }
311 
312  ///////////////////////////////// kckim ////////////////////////////////////////
313  } else if (psdRadialDiffusion.initial_PSD_Type == "IPSDT_STEADY_STATE_TWO_ZONE") {
314  // Lpp for two zone steady state
315  ReadFromFile(psdRadialDiffusion.initial_PSD_tauSteadyState, parameters_file, "psdRadialDiffusion.initial_PSD_tauSteadyState", "0.0");
316  if (psdRadialDiffusion.initial_PSD_tauSteadyState <= 1.e-98) {
317  psdRadialDiffusion.initial_PSD_tauSteadyState = 4.0/psdRadialDiffusion.initial_PSD_Kp0;
318  Output::echo(0, "psdRadialDiffusion.initial_PSD_tauSteadyState set to %f.\n", psdRadialDiffusion.initial_PSD_tauSteadyState);
319  }
320  // tauLpp for two zone steady state
321  ReadFromFile(psdRadialDiffusion.usetauLpp_SteadyState, parameters_file, "psdRadialDiffusion.usetauLpp_SteadyState", string("combined"));
322  if (psdRadialDiffusion.usetauLpp_SteadyState == "constant") {
323  ReadFromFile(psdRadialDiffusion.initial_PSD_tauLppSteadyState, parameters_file, "psdRadialDiffusion.initial_PSD_tauLppSteadyState");
324  }
325  ////////////////////////////////////////////////////////////////////////////////
326  } else if (psdRadialDiffusion.initial_PSD_Type == "IPSDT_FILE"
327  || psdRadialDiffusion.initial_PSD_Type == "IPSDT_FILE_GRID"
328  || psdRadialDiffusion.initial_PSD_Type == "IPSDT_ORBIT_FLUX_2D"){
329  //////////////////////////// kckim ////////////////////////////////////
330  //|| psdRadialDiffusion.initial_PSD_Type == "IPSDT_FLUX_2D_MAWELLIAN"){
331  ///////////////////////////////////////////////////////////////////////
332  ReadFromFile(psdRadialDiffusion.initial_PSD_fileName, parameters_file, "psdRadialDiffusion.initial_PSD_fileName");
333  }
334 
335 
336  ReadFromFile(psdRadialDiffusion.initial_PSD_some_constant_value, parameters_file, "psdRadialDiffusion.initial_PSD_some_constant_value", "1e-99"); // value_psd
337  ReadFromFile(psdRadialDiffusion.initial_PSD_J_L7_function, parameters_file, "psdRadialDiffusion.initial_PSD_J_L7_function", "J_L7_corrected");
338 
339  ReadFromFile(psdRadialDiffusion.initial_PSD_outer_psd, parameters_file, "psdRadialDiffusion.initial_PSD_outer_psd", "1.0");
340  ReadFromFile(psdRadialDiffusion.initial_PSD_inner_psd, parameters_file, "psdRadialDiffusion.initial_PSD_inner_psd", "0.0");
341 
342  ReadFromFile(psdRadialDiffusion.output_PSD_folderName, parameters_file, "psdRadialDiffusion.output_PSD_folderName", general_Output_parameters.folderName );
343  ReadFromFile(psdRadialDiffusion.output_PSD_fileName4D, parameters_file, "psdRadialDiffusion.output_PSD_fileName4D", "OutPSD_rad.dat");
344 
345  ReadFromFile(psdRadialDiffusion.output_PSD_timeStep, parameters_file, "psdRadialDiffusion.output_PSD_timeStep", "1");
346  // convert hours into days
347  psdRadialDiffusion.output_PSD_timeStep = psdRadialDiffusion.output_PSD_timeStep/24.0;
348 
349  psdRadialDiffusion.approximationMethod = "AM_Split_C";
350  psdRadialDiffusion.solutionMethod = "SM_Tridiag";
351 
352  // Local diffusions parameters
353  psdLocalDiffusions.initial_PSD_Type = "IPSDT_CONSTANT";
354  psdLocalDiffusions.initial_PSD_some_constant_value = psdRadialDiffusion.initial_PSD_some_constant_value; // value_psd
355 
356  ReadFromFile(psdLocalDiffusions.output_PSD_folderName, parameters_file, "psdLocalDiffusions.output_PSD_folderName", general_Output_parameters.folderName);
357  ReadFromFile(psdLocalDiffusions.output_PSD_fileName4D, parameters_file, "psdLocalDiffusions.output_PSD_fileName4D", string("OutPSD.dat"));
358 
359  ReadFromFile(psdLocalDiffusions.output_PSD_timeStep, parameters_file, "psdLocalDiffusions.output_PSD_timeStep", "1");
360  // convert hours into days
361  psdLocalDiffusions.output_PSD_timeStep = psdLocalDiffusions.output_PSD_timeStep/24.0;
362 
363  ReadFromFile(psdLocalDiffusions.approximationMethod, parameters_file, "psdLocalDiffusions.approximationMethod", "AM_Split_C");
364  if (psdLocalDiffusions.approximationMethod == "AM_Split_C" || psdLocalDiffusions.approximationMethod == "AM_Split_LR") {
365  ReadFromFile(psdLocalDiffusions.solutionMethod, parameters_file, "psdLocalDiffusions.solutionMethod", "SM_Tridiag");
366  } else {
367  ReadFromFile(psdLocalDiffusions.solutionMethod, parameters_file, "psdLocalDiffusions.solutionMethod", "SM_Relaxation");
368 
369  if (psdLocalDiffusions.solutionMethod == "SM_GMRES" || psdLocalDiffusions.solutionMethod == "SM_GMRES2") {
370  ReadFromFile(psdLocalDiffusions.GMRES_parameters.SOL_maxiter, parameters_file, "psdLocalDiffusions.GMRES_parameters.SOL_maxiter", "100");
371  ReadFromFile(psdLocalDiffusions.GMRES_parameters.SOL_i_max, parameters_file, "psdLocalDiffusions.GMRES_parameters.SOL_i_max", "25");
372  ReadFromFile(psdLocalDiffusions.GMRES_parameters.SOL_max_iter_err, parameters_file, "psdLocalDiffusions.GMRES_parameters.SOL_max_iter_err", "1e-1");
373  ReadFromFile(psdLocalDiffusions.GMRES_parameters.preconditioner_type, parameters_file, "psdLocalDiffusions.GMRES_parameters.preconditioner_type", "none");
374  ReadFromFile(psdLocalDiffusions.GMRES_parameters.use_normalization, parameters_file, "psdLocalDiffusions.GMRES_parameters.use_normalization", "Yes");
375  }
376  }
377 
378  // Sources and losses parameters
379  ReadFromFile(SL.SL_L_top, parameters_file, "SL.SL_L_top", "false");
380  if (SL.SL_L_top) {
381  ReadFromFile(SL.SL_L_top_filename, parameters_file, "SL.SL_L_top_filename");
382  }
383  ReadFromFile(SL.SL_E_min, parameters_file, "SL.SL_E_min", "false");
384  if (SL.SL_E_min) {
385  ReadFromFile(SL.SL_E_min_filename, parameters_file, "SL.SL_E_min_filename");
386  }
387 
388  // Interpolation parameters
389  ReadFromFile(interpolation.type, parameters_file, "interpolation.type", "IT_SPLINE");
390  ReadFromFile(interpolation.useLog, parameters_file, "interpolation.useLog", "Log_10");
391  ReadFromFile(interpolation.linearSplineCoef, parameters_file, "interpolation.linearSplineCoef", "1");
392  ReadFromFile(interpolation.maxSecondDerivative, parameters_file, "interpolation.maxSecondDerivative", "1e2");
393 
394  //ReadFromFile(DxxParameters_structureList, parameters_file, "CoeffFileName");
395 
396  // Checks
397 
398  // Check outer L boundry condition creation parameters
399  if (psdRadialDiffusion.initial_PSD_Type == "IPSDT_STEADY_STATE_FROM_BOUNDARY"
400  && (L_UpperBoundaryCondition.type != "BCT_FILE"
401  && L_UpperBoundaryCondition.type != "BCT_FILE_GRID"
402  && L_UpperBoundaryCondition.type != "BCT_MULTIPLE_FILES")
403  ) 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.");
404 
405  if ((psdLocalDiffusions.approximationMethod == "AM_Block_LR" || psdLocalDiffusions.approximationMethod == "AM_Block_C") && !(localDiffusionsGrid_pc.size > 1 && localDiffusionsGrid_alpha.size > 1)) {
406  throw error_msg("BLOCK_WITH_1D", "Block method can not be used with 1D");
407  }
408 
409  if (localDiffusionsGrid_pc.size < 3 && useEnergyDiffusion)
410  throw error_msg("ERROR", "Not enough grid points for energy diffusion (localDiffusionsGrid.pc.size = %d and useEnergyDiffusion = Yes)", localDiffusionsGrid_pc.size);
411 
412  if (localDiffusionsGrid_alpha.size < 3 && useAlphaDiffusion)
413  throw error_msg("ERROR", "Not enough grid points for pitch angle diffusion (localDiffusionsGrid.alpha.size = %d and useAlphaDiffusion = Yes)", localDiffusionsGrid_alpha.size);
414 
415  if (localDiffusionsGrid_L.size < 3 && useRadialDiffusion)
416  throw error_msg("ERROR", "Not enough grid points for radial diffusion (localDiffusionsGrid.L.size = %d and useRadialDiffusion = Yes)", localDiffusionsGrid_L.size);
417 
418  // Read diffusion coefficients
419  parameters_file.seekg(0,ios::beg);
420  string word1="", word2="";
421  string string_value, dxx_parameters_file_name;
422  if (parameters_file != NULL) {
423  while (parameters_file) {
424  word2 = word1;
425  // read one word
426  parameters_file >> word1;
427  // if last word (word2 equals previous word1) is "variableName" then read next word for the variable value
428  if (word1 == "=") {
429  if (word2 == "CoeffFileName") {
430  parameters_file >> dxx_parameters_file_name;
431 
432  // add a value to the vector of diffusion coefficients parameters structures (well, a structure for a diffusion coefficient)
433  DxxParametersList.push_back(DxxParameters_structure());
434  // Read that diffusion coefficient parameters from it's file
435  // ReadFromFile(DxxParameters_structureList.back());
436 
437  // Load dxx parameters into the just added dxx parameter structure
438  DxxParametersList.back().Load_dxx_parameters(dxx_parameters_file_name);
439 
440  }
441  }
442  }
443  }
444 
445  return true;
446 }
447 
448 
449 /**
450  * Loads parameters from file
451  */
452 bool DxxParameters_structure::Load_dxx_parameters(string dxx_parameters_file_name) {
453  ifstream dxx_parameters_file;
454  dxx_parameters_file.open(dxx_parameters_file_name.c_str());
455  if (dxx_parameters_file == NULL) {
456  throw error_msg("DXX_PARAM_FILE_OPEN_ERR", "Error reading diffusion coefficients parameters file %s.", dxx_parameters_file_name.c_str());
457  }
458 
459  ReadFromFile(DxxType, dxx_parameters_file, "DxxType");
460  ReadFromFile(DxxName, dxx_parameters_file, "DxxType");
461 
462  ReadFromFile(waveType, dxx_parameters_file, "waveType");
463  ReadFromFile(waveName, dxx_parameters_file, "waveType");
464 
465  ReadFromFile(time_start, dxx_parameters_file, "time_start");
466  ReadFromFile(time_end, dxx_parameters_file, "time_end");
467 
468  ReadFromFile(useScale, dxx_parameters_file, "useScale", "false");
469  if (useScale) {
470  ReadFromFile(DxxKp , dxx_parameters_file, "DxxKp");
471  }
472 
473  ReadFromFile(loadOrCalculate, dxx_parameters_file, "loadOrCalculate");
474 
475  ReadFromFile(filetype, dxx_parameters_file, "filetype");
476  ReadFromFile(filename, dxx_parameters_file, "filename");
477 
478  ReadFromFile(multiplicator, dxx_parameters_file, "multiplicator", "1.0");
479  ReadFromFile(MLT_averaging, dxx_parameters_file, "MLT_averaging", "1.0");
480 
481  if (loadOrCalculate == "LOC_CALCULATE" || loadOrCalculate == "LOC_LOADORCALCULATE") {
482  ReadFromFile(numberDensity, dxx_parameters_file, "numberDensity");
483  ReadFromFile(Omega_mType, dxx_parameters_file, "Omega_mType");
484  ReadFromFile(Omega_m, dxx_parameters_file, "Omega_m");
485  ReadFromFile(d_omega, dxx_parameters_file, "d_omega");
486  ReadFromFile(omega_lc, dxx_parameters_file, "omega_lc");
487  ReadFromFile(omega_uc, dxx_parameters_file, "omega_uc");
488  ReadFromFile(eta1, dxx_parameters_file, "eta1");
489  ReadFromFile(eta2, dxx_parameters_file, "eta2");
490  ReadFromFile(eta3, dxx_parameters_file, "eta3");
491  ReadFromFile(Bw, dxx_parameters_file, "Bw");
492  ReadFromFile(BwFromLambda, dxx_parameters_file, "BwFromLambda");
493 
494  ReadFromFile(nint, dxx_parameters_file, "nint");
495  ReadFromFile(mirror_point_coeff, dxx_parameters_file, "mirror_point_coeff", "0.999");
496 
497  ReadFromFile(lam_min, dxx_parameters_file, "lam_min", "0.0");
498  ReadFromFile(lam_max, dxx_parameters_file, "lam_max");
499  ReadFromFile(nu, dxx_parameters_file, "nu", "0.0");
500  ReadFromFile(s, dxx_parameters_file, "sign");
501  if (numberDensity == "ND_CONSTANT") {
502  ReadFromFile(f, dxx_parameters_file, "f_ratio");
503  }
504  ReadFromFile(particle, dxx_parameters_file, "particle");
505 
506  double sigma = (omega_uc - omega_lc) / (2.0 * d_omega);
507  Output::echo(0, "sigma = %f\n", sigma);
508  nu = erf(sigma)*sqrt(VC::pi);
509  Output::echo(0, "nu = %f\n", nu);
510 
511  if (Bw <= 0 && BwFromLambda == false) {
512  throw error_msg("DXX_ERROR", "BwFromLambda is false and Bw = %f\n", Bw);
513  } else if (BwFromLambda == true && Bw <= 0) {
514  Bw = 1;
515  }
516 
517  }
518 
519  return true;
520 }
521 
522 
523 /**
524  * Converting string to double.
525  */
526 void StrToVal(string input, double &place) {
527  istringstream(input) >> place;
528 // Output::echo("%e", place);
529 }
530 
531 /**
532  * Converting string to int.
533  */
534 void StrToVal(string input, int &place) {
535  istringstream(input) >> place;
536 // Output::echo(place);
537 }
538 
539 /**
540  * Converting string to string. Function StrToVal is used in template, so we need thad function to make template function works in case of string.
541  */
542 void StrToVal(string input, string &place) {
543  place = input;
544 // Output::echo(place);
545 }
546 
547 /**
548  * Converting string to bool.
549  */
550 void StrToVal(string input, bool &place) {
551  place = str2bool(input);;
552 // Output::echo(place);
553 }
554 
555 
556 /**
557  * Converting string to boolean.
558  */
559 bool str2bool(string str) {
560  if (str == "Yes" || str == "yes" || str == "1" || str == "True" || str == "true") return true;
561  else return false;
562 }
563 
564 /**
565  * Converting boolean to string.
566  */
567 string bool2str(bool b) {
568  return (b==true)?"Yes":"No";
569 }
570 
571 
572 /** Read 1d matrix_array from txt-file and interpolate to out time-axis.
573  */
574 void load_1d(Matrix1D<double> &var, string filename, double dt, int var_size) {
575  // var_size is a flag, if it is equal to zero, then we assume memory 1d-matrix allocated
576  // if it is equal to some value, then we need to allocate the memory
577  if (var_size != 0) {
578  // allocate memory
579  var.AllocateMemory(var_size);
580  } else {
581  // if var_size==0, assume memory already allocated
582  var_size = var.size_x;
583  }
584 
585  ifstream input_1d(filename.c_str());
586  if (input_1d != NULL) {
587  vector <double> time1d, value1d;
588  double time_tmp, value_tmp;
589  while (!input_1d.eof()) {
590  // reading time values and function values from the file into vectors
591  input_1d >> time_tmp >> value_tmp;
592  time1d.push_back(time_tmp);
593  value1d.push_back(value_tmp);
594  }
595  double time;
596  double a;
597  // interpolating loaded values into out time-axis
598  for (int it = 0; it < var.size_x; it++) {
599  time = dt*it;
600 
601  unsigned int i = 0;
602  while (((++i) < (time1d.size()-1)) && time > time1d[i]) {
603  // Not enough Kp error here
604  };
605 
606  a = (time - time1d[i-1]) / (time1d[i] - time1d[i-1]);
607  var[it] = value1d[i-1] + a*(value1d[i] - value1d[i-1]);
608  }
609  } else {
610  throw error_msg("LOAD_1D_ERROR", "Load 1d file error: file not found %f\n", filename.c_str());
611  }
612 }