VERB_code_2.2  2
 All Classes Namespaces Files Functions Variables Typedefs Macros Pages
Grid.cpp
Go to the documentation of this file.
1 /**
2  * \file Grid.cpp
3  * Grids, grid elements and boundary conditions source.
4  *
5  * \author Developed under supervision of the PI Yuri Shprits
6  */
7 
8 #include <iostream>
9 #include <iomanip>
10 
11 #include "Grid.h"
12 #include <math.h>
13 #include "../VariousFunctions/variousFunctions.h"
14 #include "../Exceptions/error.h"
15 
16 using namespace std;
17 
18 #define maxERR 1.e-6
19 
20 #if defined(_WIN32) || defined(_WIN64)
21 #define strncasecmp _strnicmp
22 #define strcasecmp _stricmp
23 #endif
24 
25 /**
26  * Convert agr to epc. Arg should be pc.
27  *
28  * \todo Wired function, should be removed.
29  *
30  * \param arg - should be pc.
31  */
33  arr.AllocateMemory(pc.arr.size_x, pc.arr.size_y, pc.arr.size_z);
34  int x, y, z;
35  for (x = 0; x < pc.arr.size_x; x++)
36  for (y = 0; y < pc.arr.size_y; y++)
37  for (z = 0; z < pc.arr.size_z; z++)
38  arr[x][y][z] = VF::Kfunc(pc.arr[x][y][z]);
39 }
40 
41 /**
42  * Function return the value for the specified point on regular grid.
43  *
44  * Very usefull function, allows to avoid a lot of typos in creation of regular grids everywhere.
45  *
46  * Just don't have to repeat regular grid creation code a lot of times.
47  * The function create logarithmic grid if useLogScale == true as: 10^(log10(min) + x * (log10(max) - log10(min))/(size-1);
48  * if useLogScale == false, then grid is: min + x * (max - min)/(size-1).
49  * If grid size == 1 the function returns max.
50  *
51  * PS: size, min and max stored in the gridElement already.
52  * PPS: useLogScale is also there.
53  *
54  * \param iteratorL - index on L-grid.
55  * \param iteratorPc - index on pc-grid.
56  * \param iteratorAlpha - index on alpha-grid.
57  * \param gridElementDirection - index on the grid, which we are making.
58  */
59 void GridElement::SetRegularGridValue(int iteratorL, int iteratorPc, int iteratorAlpha, int gridElementDirection) {
60  if (GridElement_parameters.size > 1)
61  // in case of 3d
62  if (GridElement_parameters.useLogScale) // in case of log scale
63  arr[iteratorL][iteratorPc][iteratorAlpha] = pow(10, log10(GridElement_parameters.min) + gridElementDirection*(log10(GridElement_parameters.max) - log10(GridElement_parameters.min))/(GridElement_parameters.size - 1));
64  else // in case of regular scale
65  arr[iteratorL][iteratorPc][iteratorAlpha] = GridElement_parameters.min + gridElementDirection*(GridElement_parameters.max - GridElement_parameters.min)/(GridElement_parameters.size - 1);
66  else // in case of 2d
67  arr[iteratorL][iteratorPc][iteratorAlpha] = GridElement_parameters.max;
68 }
69 
70 
73  Parameters_structure::GridElement parameters_alpha,
74  Parameters_structure::GridElement parameters_epc,
75  string grid_filename, string gridType,
76  Grid SecondGrid) {
77 
78  //if (!L.initialized) {
79  //throw error_msg("GRID_INITIALIZATION_ERROR", "Grid was not initialized");
80 
81  // allocating memory for all member classes
82  L.arr.AllocateMemory(parameters_L.size, parameters_pc.size, parameters_alpha.size);
83  pc.arr.AllocateMemory(parameters_L.size, parameters_pc.size, parameters_alpha.size);
84  alpha.arr.AllocateMemory(parameters_L.size, parameters_pc.size, parameters_alpha.size);
85  epc.arr.AllocateMemory(parameters_L.size, parameters_pc.size, parameters_alpha.size);
86  Jacobian.AllocateMemory(parameters_L.size, parameters_pc.size, parameters_alpha.size);
87 
88  //L.Initialize(parameters_L);
89  L.GridElement_parameters = parameters_L;
90  L.arr.name = parameters_L.name;
91  L.size = parameters_L.size;
92 
93  //pc.Initialize(parameters_pc);
94  pc.GridElement_parameters = parameters_pc;
95  pc.arr.name = parameters_pc.name;
96  pc.size = parameters_pc.size;
97 
98  //alpha.Initialize(parameters_alpha);
99  alpha.GridElement_parameters = parameters_alpha;
100  alpha.arr.name = parameters_alpha.name;
101  alpha.size = parameters_alpha.size;
102 
103  //epc.Initialize(parameters_epc);
104  epc.GridElement_parameters = parameters_epc;
105  epc.arr.name = parameters_epc.name;
106  epc.size = parameters_epc.size;
107 
108  type = gridType; // save the grid type for future use
109 
110  //LSize = L.size;
111  //pcSize = pc.size;
112  //alphaSize = alpha.size;
113  //epcSize = epc.size;
114 
115  //}
116 
117  // iterators, cold be useful later
118  int il, im, ia;
119  // temparary matrix_arrays
120  Matrix3D<double> mu_arr(L.size, pc.size, alpha.size);
121  // tempporary variables
122  double mu, Jc, mu_check, Jc_check;
123  if (gridType == "GT_FILE") {
124  int il, im, ia;
125  string inBuf;
126  // making stresm
127  ifstream input(grid_filename.c_str());
128  if (input != NULL && !input.eof()) {
129  // Skipping header.
130  input >> inBuf;
131  while (strcasecmp(inBuf.c_str(), "ZONE") != 0 && !input.eof() ) {
132  input >> inBuf;
133  }
134  // read to the end of the line with 'zone'
135  input.ignore(9999, '\n');
136  for (il = 0; il < L.size; il++) {
137  for (im = 0; im < epc.size; im++) {
138  for (ia = 0; ia < alpha.size; ia++) {
139  input >> L.arr[il][im][ia] >> epc.arr[il][im][ia] >> alpha.arr[il][im][ia] >> pc.arr[il][im][ia];
140  // skip till the end of the line
141  input.ignore(9999, '\n');
142  if (input.eof()) {
143  throw error_msg("GRID_INITIALIZATION_ERROR", "Not enough grid points in file %s - check grid resolution.", grid_filename.c_str());
144  }
145  }
146  }
147  }
148  input.ignore(1, '\n');
149  if (!input.eof()) {
150  Output::echo(0, "Warning: probably too many points in the grid file %s - check grid resolution.", grid_filename.c_str());
151  }
152 
153  // calculating Jacobian for the grid
154  // \todo Calculate Jacobian from the grid data
155  int il, im, ia;
156  for (il = 0; il < L.size; il++) {
157  for (im = 0; im < pc.size; im++) {
158  for (ia = 0; ia < alpha.size; ia++) {
159  // remove all constant - they don't matter
160  Jacobian[il][im][ia] = pow(pc.arr[il][im][ia],2) * pow(L.arr[il][im][ia],2) * VF::G(alpha.arr[il][im][ia]);
161  // Jacobian[il][im][ia] = -4.0 * pow(pc[il][im][ia],2) * pow(L[il][im][ia],2) * VF::G(alpha[il][im][ia])/VF::B(L[il][im][ia]);
162  //psd.G[q0][q1][q2] = -4.0 * pow(psd.pc[q0][q1][q2],2) * pow(psd.L[q0][q1][q2],2) * G(psd.alpha[q0][q1][q2])/B(psd.L[q0][q1][q2]);
163  }
164  }
165  }
166 
167  } else {
168  throw error_msg("GRID_INITIALIZATION_ERROR", "Grid file %s not found.", grid_filename.c_str());
169  }
170  } else if (gridType == "GT_L_PC_ALPHA") {
171  // making grid
172  for (il = 0; il < L.size ; il++){
173  for (im = 0; im < pc.size; im++){
174  for (ia = 0; ia < alpha.size; ia++){
175  // You see, instead of that multiline formula for grid, appears everywherem we have just one.
176  L.SetRegularGridValue(il, im, ia, il);
177  pc.SetRegularGridValue(il, im, ia, im);
178  alpha.SetRegularGridValue(il, im, ia, ia);
179  }
180  }
181  }
182  //epc = pc.Kfunc();
183  epc.Kfunc_equal(pc);
184 
185  // calculating Jacobian for the grid
186  for (il = 0; il < L.size; il++) {
187  for (im = 0; im < pc.size; im++) {
188  for (ia = 0; ia < alpha.size; ia++) {
189  // remove all constant - they don't matter
190  // !!!!!!!!!!!!!Jacobian[il][im][ia] = pow(pc.arr[il][im][ia],2) * pow(L.arr[il][im][ia],2) * VF::G(alpha.arr[il][im][ia]);
191  // Jacobian[il][im][ia] = -4.0 * pow(pc[il][im][ia],2) * pow(L[il][im][ia],2) * VF::G(alpha[il][im][ia])/VF::B(L[il][im][ia]);
192  Jacobian[il][im][ia] = pow(pc.arr[il][im][ia],2) * VF::G(alpha.arr[il][im][ia]);
193 
194  }
195  }
196  }
197  //Jacobian = 1;
198  //Output::echo(0, "Jacobian == 1 !!!!!!!!!!!!!!!!\n");
199 
200 
201  } else if (gridType == "GT_L_MU_J") {
202  // Make L regular grid
203  for (il = 0; il < L.size ; il++) {
204  for (im = 0; im < pc.size; im++) {
205  for (ia = 0; ia < alpha.size; ia++) {
206  // You see, instead of that multiline formula for grid, appears everywherem we have just one.
207  L.SetRegularGridValue(il, im, ia, il);
208  }
209  }
210  }
211  // Making a 2-D grid on the top level at L=Lmax
212  il = L.size-1;
213  for (im = 0; im < pc.size; im++){
214  for (ia = 0; ia < alpha.size; ia++) {
215  // You see, instead of that multiline formula for grid, appears everywherem we have just one.
216  pc.SetRegularGridValue(il, im, ia, im);
217  alpha.SetRegularGridValue(il, im, ia, ia);
218  }
219  }
220  // Making a grid which conserves 1st J1 and second J2 adiabatic invariant
221  if (L.size > 1) {
222  for (im = 0; im < pc.size; im++) {
223  for (ia = 0; ia < alpha.size; ia++) {
224  mu = VF::mu_calc(L.arr[L.size-1][im][ia], pc.arr[L.size-1][im][ia], alpha.arr[L.size-1][im][ia]);
225  Jc = VF::Jc_calc(L.arr[L.size-1][im][ia], pc.arr[L.size-1][im][ia], alpha.arr[L.size-1][im][ia]);
226  for (il = 0; il < L.size-1 ; il++) {
227  if (alpha.arr[L.size-1][im][ia] < VC::pi/2)
228  alpha.arr[il][im][ia] = find_alpha(Jc * sqrt(L.arr[il][im][ia]) / sqrt(8.0 * VF::B(1) * VC::mc2 * mu), alpha.GridElement_parameters.min, VC::pi/2);
229  else if (alpha.arr[L.size-1][im][ia] > VC::pi/2)
230  alpha.arr[il][im][ia] = find_alpha(Jc * sqrt(L.arr[il][im][ia]) / sqrt(8.0 * VF::B(1) * VC::mc2 * mu), VC::pi/2, alpha.GridElement_parameters.max);
231  else
232  alpha.arr[il][im][ia] = VC::pi/2;
233  pc.arr[il][im][ia] = sqrt(2.0 * mu * VC::mc2 * VF::B(L.arr[il][im][ia])) / sin(alpha.arr[il][im][ia]);
234  mu_check = VF::mu_calc(L.arr[il][im][ia], pc.arr[il][im][ia], alpha.arr[il][im][ia]);
235  Jc_check = VF::Jc_calc(L.arr[il][im][ia], pc.arr[il][im][ia], alpha.arr[il][im][ia]);
236  if (fabs(mu - mu_check) > maxERR || fabs(Jc - Jc_check) > maxERR) {
237  throw error_msg("L_MU_J_GRID_ERROR", "Grid computation error (L, pc, a) = (%f, %f, %f) :\n mu=%e, Jc=%e\n mu_chk=%e, Jc_chk=%e\n mu_err=%e, Jc_err=%e.\n", L.arr[il][im][ia], pc.arr[il][im][ia], alpha.arr[il][im][ia], mu, Jc, mu_check, Jc_check, fabs(mu - mu_check), fabs(Jc - Jc_check));
238  }
239  }
240  }
241  }
242  }
243  epc.Kfunc_equal(pc);
244 
245  // calculating Jacobian for the grid
246  for (il = 0; il < L.size; il++) {
247  for (im = 0; im < pc.size; im++) {
248  for (ia = 0; ia < alpha.size; ia++) {
249  // remove all constant - they don't matter
250 
251 // Jacobian[il][im][ia] = pow(pc[il][im][ia],2) * pow(L[il][im][ia],-2) * VF::G(alpha[il][im][ia]);;
252  Jacobian[il][im][ia] = pow(L.arr[il][im][ia],-2); // <--
253 // Jacobian[il][im][ia] = 1;
254 // Jacobian[il][im][ia] = pow(L[il][im][ia],2);
255 
256  // Jacobian[il][im][ia] = -4.0 * pow(pc[il][im][ia],2) * pow(L[il][im][ia],2) * VF::G(alpha[il][im][ia])/VF::B(L[il][im][ia]);
257  //psd.G[q0][q1][q2] = -4.0 * pow(psd.pc[q0][q1][q2],2) * pow(psd.L[q0][q1][q2],2) * G(psd.alpha[q0][q1][q2])/B(psd.L[q0][q1][q2]);
258  }
259  }
260  }
261 
262  } else if (gridType == "GT_L_MU_ALPHA") {
263  // Make L regular grid
264  for (il = 0; il < L.size ; il++) {
265  for (im = 0; im < pc.size; im++) {
266  for (ia = 0; ia < alpha.size; ia++) {
267  // You see, instead of that multiline formula for grid, appears everywherem we have just one.
268  L.SetRegularGridValue(il, im, ia, il);
269  }
270  }
271  }
272  // Making a 2-D grid on the top level at L=Lmax
273  il = L.size-1;
274  for (im = 0; im < pc.size; im++){
275  for (ia = 0; ia < alpha.size; ia++) {
276  pc.SetRegularGridValue(il, im, ia, im);
277  alpha.SetRegularGridValue(il, im, ia, ia);
278 // mu_arr[il][im][ia] = VF::pc2mu(L[il][im][ia], pc[il][im][ia], alpha[il][im][ia]);
279  }
280  }
281  for (im = 0; im < pc.size; im++) {
282  for (ia = 0; ia < alpha.size; ia++) {
283  for (il = 0; il < L.size-1 ; il++) {
284  alpha.SetRegularGridValue(il, im, ia, ia);
285  // mu on the top
286  mu = VF::pc2mu(L.arr[L.size-1][im][ia], pc.arr[L.size-1][im][ia], alpha.arr[L.size-1][im][ia]);
287  pc.arr[il][im][ia] = VF::mu2pc(L.arr[il][im][ia], mu, alpha.arr[il][im][ia]);
288  }
289  }
290  }
291  epc.Kfunc_equal(pc);
292 
293  // calculating Jacobian for the grid
294  for (il = 0; il < L.size; il++) {
295  for (im = 0; im < pc.size; im++) {
296  for (ia = 0; ia < alpha.size; ia++) {
297  // remove all constant - they don't matter
298  // L^-2?? If it's on (L, pc, alpha) grid, I think it should be L^2. It does not matter cause it cancels anyway, but still.
299  Jacobian[il][im][ia] = pow(pc.arr[il][im][ia],2) * pow(L.arr[il][im][ia],-2) * VF::G(alpha.arr[il][im][ia]);
300  // Jacobian[il][im][ia] = -4.0 * pow(pc[il][im][ia],2) * pow(L[il][im][ia],2) * VF::G(alpha[il][im][ia])/VF::B(L[il][im][ia]);
301  //psd.G[q0][q1][q2] = -4.0 * pow(psd.pc[q0][q1][q2],2) * pow(psd.L[q0][q1][q2],2) * G(psd.alpha[q0][q1][q2])/B(psd.L[q0][q1][q2]);
302  }
303  }
304  }
305 
306  } else if (gridType == "GT_PERP_TO_L_MU_J") {
307  if (!SecondGrid.Grid_initialized) {
308  throw error_msg("GRID_INIT_ERR", "GT_PERP_TO_L_MU_J with empty second grid");
309  }
310  L = SecondGrid.L;
311  alpha = SecondGrid.alpha;
312 
313  // make new pc
314  double pc_local_max, pc_local_min;
315  for (il = 0; il < L.size ; il++){
316  // maximum from the upper ENERGY boundary
317  pc_local_max = SecondGrid.pc.arr[il][SecondGrid.pc.size-1][alpha.size-1];
318  for (ia = 0; ia < alpha.size; ia++) if (pc_local_max < SecondGrid.pc.arr[il][SecondGrid.pc.size-1][ia]) pc_local_max = SecondGrid.pc.arr[il][SecondGrid.pc.size-1][ia];
319  // MAXIMUM from the LOWER energy boundary
320  pc_local_min = SecondGrid.pc.arr[il][0][alpha.size-1];
321  for (ia = 0; ia < alpha.size; ia++) if (pc_local_min < SecondGrid.pc.arr[il][0][ia]) pc_local_min = SecondGrid.pc.arr[il][0][ia];
322  for (im = 0; im < pc.size; im++){
323  for (ia = 0; ia < alpha.size; ia++){
324  if (pc.size > 1) {
325  if (!pc.GridElement_parameters.useLogScale) pc.arr[il][im][ia] = pc_local_min + im*(pc_local_max - pc_local_min)/(pc.size - 1);
326  else pc.arr[il][im][ia] = pow(10, log10(pc_local_min) + im*(log10(pc_local_max) - log10(pc_local_min))/(pc.size - 1));
327  } else pc.arr[il][im][ia] = pc.GridElement_parameters.max;
328  }
329  }
330  }
331 
332  epc.Kfunc_equal(pc);
333 
334  // calculating Jacobian for the grid
335  for (il = 0; il < L.size; il++) {
336  for (im = 0; im < pc.size; im++) {
337  for (ia = 0; ia < alpha.size; ia++) {
338  // remove all constant - they don't matter
339  Jacobian[il][im][ia] = pow(pc.arr[il][im][ia],2) * pow(L.arr[il][im][ia],2) * VF::G(alpha.arr[il][im][ia]);
340  // Jacobian[il][im][ia] = -4.0 * pow(pc[il][im][ia],2) * pow(L[il][im][ia],2) * VF::G(alpha[il][im][ia])/VF::B(L[il][im][ia]);
341  //psd.G[q0][q1][q2] = -4.0 * pow(psd.pc[q0][q1][q2],2) * pow(psd.L[q0][q1][q2],2) * G(psd.alpha[q0][q1][q2])/B(psd.L[q0][q1][q2]);
342  }
343  }
344  }
345 
346  } else {
347  throw error_msg("GRID_UNKNOWN_TYPE", "Unknown grid type: %s", gridType.c_str());
348  }
349 
350  Jacobian.writeToFile("./Debug/Jacobian.plt");
351 
352 
353  Grid_initialized = true;
354 
355 }
356 
357 
358 /**
359  * Output grid to a file.
360  *
361  * \param filename - file name
362  */
363 void Grid::Output(string filename) {
364  ofstream output(filename.c_str());
365  int il, im, ia;
366 
367  // output header of the file
368  output << "VARIABLES = ";
369  output << "\"L, Earth radius\", \"Energy, MeV\", \"Pitch-angle, deg\", \"pc\" " << endl;
370 
371  // output values...
372  output << "ZONE T=\"Grid\" " << " I=" << alpha.size << ", J=" << epc.size << ", K=" << L.size << endl;
373  output.setf(ios_base::scientific, ios_base::floatfield);
374  for (il = 0; il < L.size; il++) {
375  for (im = 0; im < epc.size; im++) {
376  for (ia = 0; ia < alpha.size; ia++) {
377  output << L.arr[il][im][ia] << "\t" << epc.arr[il][im][ia] << "\t" << alpha.arr[il][im][ia] << "\t" << VF::pfunc(epc.arr[il][im][ia]) << endl;
378  }
379  }
380  }
381 }
382 
383 /**
384  * Function finds alpha-root of equation Y(alpha)/sin(alpha) = RHS.. Used in grid-bilding.
385  *
386  * \param RHS - right hand side
387  * \param alpha_min - min alpha value to look for root
388  * \param alpha_max - max alpha value to look for root
389  * \param ERR - max error (min accuracy)
390  * \param max_it - max number of iterations
391  * \param it - current iteration number (for recurcion)
392  */
393 
394 
395 
396 double find_alpha(double RHS, double alpha_min, double alpha_max, double ERR, int max_it, int it) {
397 
398  double alpha_root;
399  do {
400  it++;
401 
402  alpha_root = (alpha_min + alpha_max) / 2;
403  double y0 = VF::Y(alpha_min)/sin(alpha_min) - RHS, yy = VF::Y(alpha_root)/sin(alpha_root) - RHS;
404 
405  if (yy * y0 < 0) alpha_max = alpha_root;
406  else if (yy * y0 > 0) alpha_min = alpha_root;
407  else return alpha_root;
408 
409  } while (fabs(alpha_max - alpha_min) > ERR);
410 
411  return alpha_root;
412 
413 }