ARTS 2.5.0 (git: 9ee3ac6c)
m_batch.cc
Go to the documentation of this file.
1/* Copyright (C) 2004-2012
2 Patrick Eriksson <Patrick.Eriksson@chalmers.se>
3 Stefan Buehler <sbuehler@ltu.se>
4
5 This program is free software; you can redistribute it and/or modify it
6 under the terms of the GNU General Public License as published by the
7 Free Software Foundation; either version 2, or (at your option) any
8 later version.
9
10 This program is distributed in the hope that it will be useful,
11 but WITHOUT ANY WARRANTY; without even the implied warranty of
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 GNU General Public License for more details.
14
15 You should have received a copy of the GNU General Public License
16 along with this program; if not, write to the Free Software
17 Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
18 USA. */
19
20/*===========================================================================
21 === File description
22 ===========================================================================*/
23
35/*===========================================================================
36 === External declarations
37 ===========================================================================*/
38
39#include <cmath>
40using namespace std;
41
42#include "arts.h"
43#include "arts_omp.h"
44#include "auto_md.h"
45#include "math_funcs.h"
46#include "physics_funcs.h"
47#include "rte.h"
48#include "xml_io.h"
49
50extern const Numeric PI;
51extern const Numeric DEG2RAD;
52extern const Numeric RAD2DEG;
53extern const Index GFIELD3_P_GRID;
54
55/*===========================================================================
56 === The functions (in alphabetical order)
57 ===========================================================================*/
58
59/* Workspace method: Doxygen documentation will be auto-generated
60
61 2008-07-21 Stefan Buehler */
63 // WS Input:
64 const Agenda& forloop_agenda,
65 // Control Parameters:
66 const Index& start,
67 const Index& stop,
68 const Index& step,
69 const Verbosity& verbosity) {
71
72 for (Index i = start; i <= stop; i += step) {
73 out1 << " Executing for loop body, index: " << i << "\n";
74 forloop_agendaExecute(ws, i, forloop_agenda);
75 }
76}
77
78/* Workspace method: Doxygen documentation will be auto-generated */
80 // WS Output:
81 ArrayOfVector& ybatch,
82 ArrayOfArrayOfVector& ybatch_aux,
83 ArrayOfMatrix& ybatch_jacobians,
84 // WS Input:
85 const Index& ybatch_start,
86 const Index& ybatch_n,
87 const Agenda& ybatch_calc_agenda,
88 // Control Parameters:
89 const Index& robust,
90 const Verbosity& verbosity) {
92
93 Index first_ybatch_index = 0;
94
95 ArrayOfString fail_msg;
96 bool do_abort = false;
97
98 // We allow a start index ybatch_start that is different from 0. We
99 // will calculate ybatch_n jobs starting at the start
100 // index. Internally, we count from zero, which is the right
101 // counting for the output array ybatch. When we call
102 // ybatch_calc_agenda, we add ybatch_start to the internal index
103 // count.
104
105 // We create a counter, so that we can generate nice output about
106 // how many jobs are already done. (All parallel threads later will
107 // increment this, so that we really get an accurate total count!)
108 Index job_counter = 0;
109
110 // Resize the output arrays:
111 ybatch.resize(ybatch_n);
112 ybatch_aux.resize(ybatch_n);
113 ybatch_jacobians.resize(ybatch_n);
114 for (Index i = 0; i < ybatch_n; i++) {
115 ybatch[i].resize(0);
116 ybatch_aux[i].resize(0);
117 ybatch_jacobians[i].resize(0, 0);
118 }
119
120 // We have to make a local copy of the Workspace and the agendas because
121 // only non-reference types can be declared firstprivate in OpenMP
122 Workspace l_ws(ws);
123 Agenda l_ybatch_calc_agenda(ybatch_calc_agenda);
124
125 // Go through the batch:
126
127 if (ybatch_n)
128#pragma omp parallel for schedule(dynamic) if (!arts_omp_in_parallel() && \
129 ybatch_n > 1) \
130 firstprivate(l_ws, l_ybatch_calc_agenda)
131 for (Index ybatch_index = first_ybatch_index; ybatch_index < ybatch_n;
132 ybatch_index++) {
133 Index l_job_counter; // Thread-local copy of job counter.
134
135 if (do_abort) continue;
136#pragma omp critical(ybatchCalc_job_counter)
137 { l_job_counter = ++job_counter; }
138
139 {
140 ostringstream os;
141 os << " Job " << l_job_counter << " of " << ybatch_n << ", Index "
142 << ybatch_start + ybatch_index << ", Thread-Id "
143 << arts_omp_get_thread_num() << "\n";
144 out2 << os.str();
145 }
146
147 try {
148 Vector y;
149 ArrayOfVector y_aux;
150 Matrix jacobian;
151
153 y,
154 y_aux,
155 jacobian,
156 ybatch_start + ybatch_index,
157 l_ybatch_calc_agenda);
158
159 if (y.nelem()) {
160#pragma omp critical(ybatchCalc_assign_y)
161 ybatch[ybatch_index] = y;
162#pragma omp critical(ybatchCalc_assign_y_aux)
163 ybatch_aux[ybatch_index] = y_aux;
164
165 // Dimensions of Jacobian:
166 const Index Knr = jacobian.nrows();
167 const Index Knc = jacobian.ncols();
168
169 if (Knr != 0 || Knc != 0) {
170 if (Knr != y.nelem()) {
171 ostringstream os;
172 os << "First dimension of Jacobian must have same length as the measurement *y*.\n"
173 << "Length of *y*: " << y.nelem() << "\n"
174 << "Dimensions of *jacobian*: (" << Knr << ", " << Knc
175 << ")\n";
176 // A mismatch of the Jacobian dimension is a fatal error
177 // and should result in program termination. By setting abort
178 // to true, this will result in a runtime error in the catch
179 // block even if robust == 1
180#pragma omp critical(ybatchCalc_setabort)
181 do_abort = true;
182
183 throw runtime_error(os.str());
184 }
185
186 ybatch_jacobians[ybatch_index] = jacobian;
187
188 // After creation, all individual Jacobi matrices in the array will be
189 // empty (size zero). No need for explicit initialization.
190 }
191 }
192 } catch (const std::exception& e) {
193 if (robust && !do_abort) {
194 ostringstream os;
195 os << "WARNING! Job at ybatch_index " << ybatch_start + ybatch_index
196 << " failed.\n"
197 << "y Vector in output variable ybatch will be empty for this job.\n"
198 << "The runtime error produced was:\n"
199 << e.what() << "\n";
200 out0 << os.str();
201 } else {
202 // The user wants the batch job to fail if one of the
203 // jobs goes wrong.
204#pragma omp critical(ybatchCalc_setabort)
205 do_abort = true;
206
207 ostringstream os;
208 os << " Job at ybatch_index " << ybatch_start + ybatch_index
209 << " failed. Aborting...\n";
210 out1 << os.str();
211 }
212 ostringstream os;
213 os << "Run-time error at ybatch_index " << ybatch_start + ybatch_index
214 << ": \n"
215 << e.what();
216#pragma omp critical(ybatchCalc_push_fail_msg)
217 fail_msg.push_back(os.str());
218 }
219 }
220
221 if (fail_msg.nelem()) {
222 ostringstream os;
223
224 if (!do_abort) os << "\nError messages from failed batch cases:\n";
225 for (ArrayOfString::const_iterator it = fail_msg.begin();
226 it != fail_msg.end();
227 it++)
228 os << *it << '\n';
229
230 if (do_abort)
231 throw runtime_error(os.str());
232 else
233 out0 << os.str();
234 }
235}
236
237/* Workspace method: Doxygen documentation will be auto-generated */
239 //Output
240 ArrayOfVector& ybatch,
241 //Input
242 const ArrayOfArrayOfSpeciesTag& abs_species,
243 const Agenda& met_profile_calc_agenda,
244 const Vector& f_grid,
245 const Matrix& met_amsu_data,
246 const Matrix& sensor_pos,
247 const Vector& refellipsoid,
248 const Vector& lat_grid,
249 const Vector& lon_grid,
250 const Index& atmosphere_dim,
251 const ArrayOfArrayOfSingleScatteringData& scat_data,
252 //Keyword
253 const Index& nelem_p_grid,
254 const String& met_profile_path,
255 const String& met_profile_pnd_path,
256 const Verbosity& verbosity) {
257 GriddedField3 t_field_raw;
258 GriddedField3 z_field_raw;
259 ArrayOfGriddedField3 vmr_field_raw;
260 ArrayOfGriddedField3 pnd_field_raw;
261 Vector p_grid;
262 Matrix sensor_los;
263 Index cloudbox_on;
264 ArrayOfIndex cloudbox_limits;
265 Matrix z_surface;
266 Vector y;
267 Index no_profiles = met_amsu_data.nrows();
268
269 // *vmr_field_raw* is an ArrayOfArrayOfTensor3 where the first array
270 //holds the gaseous species.
271 //Resize *vmr_field_raw* according to the number of gaseous species
272 //elements
273 vmr_field_raw.resize(abs_species.nelem());
274
275 //pnd_field_raw is an ArrayOfArrayOfTensor3 where the first array
276 //holds the scattering elements.
277 // Number of scattering elements:
278 const Index N_se = TotalNumberOfElements(scat_data);
279
280 pnd_field_raw.resize(N_se);
281
282 // The satellite zenith angle is read in from the amsu data
283 // and converted to arts sensor_los
284 ConstVectorView sat_za_from_data = met_amsu_data(Range(joker), 3);
285
286 sensor_los.resize(1, 1);
287
288 // The lat and lon are extracted to get the proper file names of
289 // profiles
290 ConstVectorView lat = met_amsu_data(Range(joker), 0);
291 ConstVectorView lon = met_amsu_data(Range(joker), 1);
292
293 z_surface.resize(1, 1);
294
295 // The spectra .
296 y.resize(f_grid.nelem());
297
298 // The batch spectra.
299 ybatch.resize(no_profiles);
300
301 // Loop over the number of profiles.
302 for (Index i = 0; i < no_profiles; ++i) {
303 ostringstream lat_os, lon_os;
304
305 Index lat_prec = 3;
306 if (lat[i] < 0) lat_prec--;
307 if (abs(lat[i]) >= 10) {
308 lat_prec--;
309 if (abs(lat[i]) >= 100) lat_prec--;
310 }
311
312 lat_os.setf(ios::showpoint | ios::fixed);
313 lat_os << setprecision((int)lat_prec) << lat[i];
314
315 Index lon_prec = 4;
316 if (lon[i] < 0) lon_prec--;
317 if (abs(lon[i]) >= 10) {
318 lon_prec--;
319 if (abs(lon[i]) >= 100) lon_prec--;
320 }
321 lon_os.setf(ios::showpoint | ios::fixed);
322 lon_os << setprecision((int)lon_prec) << lon[i];
323
324 sensor_los(0, 0) =
325 180.0 - (asin(refellipsoid[0] * sin(sat_za_from_data[i] * DEG2RAD) /
326 sensor_pos(0, 0))) *
327 RAD2DEG;
328
329 //Reads the t_field_raw from file
330 xml_read_from_file(met_profile_path + "profile.lat_" + lat_os.str() +
331 ".lon_" + lon_os.str() + ".t.xml",
332 t_field_raw,
333 verbosity);
334
335 //Reads the z_field_raw from file
336 xml_read_from_file(met_profile_path + "profile.lat_" + lat_os.str() +
337 ".lon_" + lon_os.str() + ".z.xml",
338 z_field_raw,
339 verbosity);
340
341 //Reads the humidity from file - it is only an ArrayofTensor3
342 xml_read_from_file(met_profile_path + "profile.lat_" + lat_os.str() +
343 ".lon_" + lon_os.str() + ".H2O.xml",
344 vmr_field_raw[0],
345 verbosity);
346
347 //Reads the pnd_field_raw for one scattering element
348 //xml_read_from_file("/rinax/storage/users/rekha/uk_data/profiles/new_obs/newest_forecastfields/reff100/profile.lat_"+lat_os.str()+".lon_"+lon_os.str() + ".pnd100.xml", pnd_field_raw[0]);
349
350 //xml_read_from_file(met_profile_pnd_path +"reff100_newformat/profile.lat_"+lat_os.str()+".lon_"+lon_os.str() + ".pnd100.xml", pnd_field_raw[0]);
351
352 xml_read_from_file(met_profile_pnd_path + "lwc_reff15/profile.lat_" +
353 lat_os.str() + ".lon_" + lon_os.str() + ".pnd15.xml",
354 pnd_field_raw[0],
355 verbosity);
356 //Write the profile number into a file.
357 // xml_write_to_file("profile_number.xml", i);
358
359 // Set z_surface from lowest level of z_field
360 z_surface(0, 0) = z_field_raw.data(0, 0, 0);
361
362 /* The vmr_field_raw is an ArrayofArrayofTensor3 where the outer
363 array is for species.
364
365 The oxygen and nitrogen VMRs are set to constant values of 0.209
366 and 0.782, respectively and are used along with humidity field
367 to generate *vmr_field_raw*.*/
368
369 /*The second element of the species. The first 3 Tensors in the
370 array are the same . They are pressure grid, latitude grid and
371 longitude grid. The third tensor which is the vmr is set to a
372 constant value of 0.782, corresponding to N2.*/
373
374 vmr_field_raw[1].resize(vmr_field_raw[0]);
375 vmr_field_raw[1].copy_grids(vmr_field_raw[0]);
376 vmr_field_raw[1] = 0.782; //vmr of N2
377
378 /*the third element of the species. the first 3 Tensors in the
379 array are the same . They are pressure grid, latitude grid and
380 longitude grid. The third tensor which is the vmr is set to a
381 constant value of 0.209, corresponding to O2.*/
382 vmr_field_raw[2].resize(vmr_field_raw[0]);
383 vmr_field_raw[2].copy_grids(vmr_field_raw[0]);
384 vmr_field_raw[2] = 0.209; //vmr of O2
385
386 const Vector& tfr_p_grid =
387 t_field_raw.get_numeric_grid(GFIELD3_P_GRID);
388 // N_p is the number of elements in the pressure grid
389 Index N_p = tfr_p_grid.nelem();
390
391 //Making a p_grid with the first and last element taken from the profile.
393 p_grid, nelem_p_grid, tfr_p_grid[0], tfr_p_grid[N_p - 1], verbosity);
394
395 /*To set the cloudbox limits, the lower and upper cloudbox limits
396 are to be set. The lower cloudbox limit is set to the lowest
397 pressure level. The upper level is the highest level where the
398 ice water content is non-zero.*/
399 Numeric cl_grid_min, cl_grid_max;
400
401 //Lower limit = lowest pressure level of the original grid.
402 //Could it be the interpolated p_grid? FIXME STR
403 cl_grid_min = tfr_p_grid[0];
404
405 // A counter for non-zero ice content
406 Index level_counter = 0;
407
408 // Loop over all pressure levels
409 for (Index ip = 0; ip < N_p; ++ip) {
410 //Checking for non-zero ice content. 0.001 is a threshold for
411 //ice water content.
412 // if((pnd_field_raw[0].data()(ip, 0, 0) > 0.001) || (pnd_field_raw[1](ip, 0, 0) > 0.001))
413 if (pnd_field_raw[0].data(ip, 0, 0) > 0.001) {
414 ++level_counter;
415 //if non-zero ice content is found, it is set to upper
416 //cloudbox limit. Moreover, we take one level higher
417 // than the upper limit because we want the upper limit
418 //to have 0 pnd.
419 cl_grid_max = tfr_p_grid[ip + 1];
420 }
421 }
422
423 //cloudbox limits have dimensions 2*atmosphere_dim
424 cloudbox_limits.resize(atmosphere_dim * 2);
425
426 //if there is no cloud in the considered profile, still we
427 //need to set the upper limit. I here set the first level
428 //for the upper cloudbox limit.
429 if (level_counter == 0) {
430 cl_grid_max = p_grid[1];
431 }
432
433 //Cloudbox is set.
434 cloudboxSetManually(cloudbox_on,
435 cloudbox_limits,
436 atmosphere_dim,
437 p_grid,
438 lat_grid,
439 lon_grid,
440 cl_grid_min,
441 cl_grid_max,
442 0,
443 0,
444 0,
445 0,
446 verbosity);
447
448 /*executing the met_profile_calc_agenda
449 Agenda communication variables are
450 Output of met_profile_calc_agenda : y
451 Input to met_profile_calc_agenda : t_field_raw,
452 z_field_raw, vmr_field_raw, pnd_field_raw, p_grid,
453 sensor_los, cloudbox_on, cloudbox_limits, z_surface, */
454
456 y,
457 t_field_raw,
458 vmr_field_raw,
459 z_field_raw,
460 pnd_field_raw,
461 p_grid,
462 sensor_los,
463 cloudbox_on,
464 cloudbox_limits,
465 z_surface,
466 met_profile_calc_agenda);
467
468 //putting in the spectra *y* for each profile, thus assigning y
469 //to the ith row of ybatch
470 ybatch[i] = y;
471
472 } // closing the loop over profile basenames
473}
474
475/* Workspace method: Doxygen documentation will be auto-generated */
477 //Output
478 ArrayOfVector& ybatch,
479 //Input
480 const ArrayOfArrayOfSpeciesTag& abs_species,
481 const Agenda& met_profile_calc_agenda,
482 const Vector& f_grid,
483 const Matrix& met_amsu_data,
484 const Matrix& sensor_pos,
485 const Vector& refellipsoid,
486 //Keyword
487 const Index& nelem_p_grid,
488 const String& met_profile_path,
489 const Verbosity& verbosity) {
490 GriddedField3 t_field_raw;
491 GriddedField3 z_field_raw;
492 ArrayOfGriddedField3 vmr_field_raw;
493 ArrayOfGriddedField3 pnd_field_raw;
494 Vector p_grid;
495 Matrix sensor_los;
496 Index cloudbox_on = 0;
497 ArrayOfIndex cloudbox_limits;
498 Matrix z_surface;
499 Vector y;
500 Index no_profiles = met_amsu_data.nrows();
501 //Index no_profiles = met_profile_basenames.nelem();
502 // The humidity data is stored as an ArrayOfTensor3 whereas
503 // vmr_field_raw is an ArrayOfArrayOfTensor3
504 GriddedField3 vmr_field_raw_h2o;
505
506 vmr_field_raw.resize(abs_species.nelem());
507
508 y.resize(f_grid.nelem());
509 ybatch.resize(no_profiles);
510
511 Vector sat_za_from_profile;
512 sat_za_from_profile = met_amsu_data(Range(joker), 3);
513 Numeric sat_za;
514
515 sensor_los.resize(1, 1);
516
517 Vector lat, lon;
518 lat = met_amsu_data(Range(joker), 0);
519 lon = met_amsu_data(Range(joker), 1);
520
521// Vector oro_height;
522// oro_height = met_amsu_data(Range(joker), 5);
523
524 z_surface.resize(1, 1);
525 for (Index i = 0; i < no_profiles; ++i) {
526 ostringstream lat_os, lon_os;
527
528 Index lat_prec = 3;
529 if (lat[i] < 0) lat_prec--;
530 if (abs(lat[i]) >= 10) {
531 lat_prec--;
532 if (abs(lat[i]) >= 100) lat_prec--;
533 }
534
535 lat_os.setf(ios::showpoint | ios::fixed);
536 lat_os << setprecision((int)lat_prec) << lat[i];
537
538 Index lon_prec = 4;
539 if (lon[i] < 0) lon_prec--;
540 if (abs(lon[i]) >= 10) {
541 lon_prec--;
542 if (abs(lon[i]) >= 100) lon_prec--;
543 }
544 lon_os.setf(ios::showpoint | ios::fixed);
545 lon_os << setprecision((int)lon_prec) << lon[i];
546 cout << lat_os.str() << endl;
547 cout << lon_os.str() << endl;
548
549 sat_za = sat_za_from_profile[i];
550
551 sensor_los(Range(joker), 0) =
552 180.0 -
553 (asin(refellipsoid[0] * sin(sat_za * PI / 180.) / sensor_pos(0, 0))) *
554 180. / PI;
555 cout << "sensor_los" << sat_za_from_profile[i] << endl;
556 cout << "sensor_los" << sat_za << endl;
557 cout << "sensor_los" << sensor_los << endl;
558 //Reads the t_field_raw from file
559
560 xml_read_from_file(met_profile_path + "profile.lat_" + lat_os.str() +
561 ".lon_" + lon_os.str() + ".t.xml",
562 t_field_raw,
563 verbosity);
564 //Reads the z_field_raw from file
565 xml_read_from_file(met_profile_path + "profile.lat_" + lat_os.str() +
566 ".lon_" + lon_os.str() + ".z.xml",
567 z_field_raw,
568 verbosity);
569
570 //Reads the humidity from file - it is only an ArrayofTensor3
571 // The vmr_field_raw is an ArrayofArrayofTensor3 where the outer
572 // array is for species
573 xml_read_from_file(met_profile_path + "profile.lat_" + lat_os.str() +
574 ".lon_" + lon_os.str() + ".H2O.xml",
575 vmr_field_raw_h2o,
576 verbosity);
577 //xml_read_from_file("/home/home01/rekha/uk/profiles/sat_vmr/profile.lat_"+lat_os.str()//+".lon_"+lon_os.str() + ".H2O_es.xml",
578 // vmr_field_raw_h2o, verbosity);
579
580 cout
581 << "--------------------------------------------------------------------------"
582 << endl;
583 cout << "The file"
584 << met_profile_path + "profile.lat_" + lat_os.str() + ".lon_" +
585 lon_os.str()
586 << "is executed now" << endl;
587 cout
588 << "--------------------------------------------------------------------------"
589 << endl;
590 xml_write_to_file("profile_number.xml", i, FILE_TYPE_ASCII, 0, verbosity);
591 // the first element of the species is water vapour.
592
593 // N_p is the number of elements in the pressure grid
594 //z_surface(0,0) = oro_height[i]+ 0.01;
595 z_surface(0, 0) = z_field_raw.data(0, 0, 0);
596 cout << "z_surface" << z_surface << endl;
597 const Vector& tfr_p_grid =
598 t_field_raw.get_numeric_grid(GFIELD3_P_GRID);
599 Index N_p = tfr_p_grid.nelem();
600
601 vmr_field_raw[0] = vmr_field_raw_h2o;
602
603 // the second element of the species. the first 3 Tensors in the
604 //array are the same . They are pressure grid, latitude grid and
605 // longitude grid. The third tensor which is the vmr is set to a
606 // constant value of 0.782.
607 vmr_field_raw[1].resize(vmr_field_raw[0]);
608 vmr_field_raw[1].copy_grids(vmr_field_raw[0]);
609 vmr_field_raw[1].data(joker, joker, joker) = 0.782;
610
611 // the second element of the species. the first 3 Tensors in the
612 //array are the same . They are pressure grid, latitude grid and
613 // longitude grid. The third tensor which is the vmr is set to a
614 // constant value of 0.209.
615 vmr_field_raw[2].resize(vmr_field_raw[0]);
616 vmr_field_raw[2].copy_grids(vmr_field_raw[0]);
617 vmr_field_raw[2].data(joker, joker, joker) = 0.209;
618
619 //xml_write_to_file(met_profile_basenames[i]+ ".N2.xml", vmr_field_raw[1]);
620 //xml_write_to_file(met_profile_basenames[i]+ ".O2.xml", vmr_field_raw[2]);
621
622 //Making a p_grid with the first and last element taken from the profile.
623 // this is because of the extrapolation problem.
624
626 p_grid, nelem_p_grid, tfr_p_grid[0], tfr_p_grid[N_p - 1], verbosity);
627 cout << "t_field_raw[0](0,0,0)" << tfr_p_grid[0] << endl;
628 cout << "t_field_raw[0](N_p -1,0,0)" << tfr_p_grid[N_p - 1] << endl;
629 xml_write_to_file("p_grid.xml", p_grid, FILE_TYPE_ASCII, 0, verbosity);
630
631 // executing the met_profile_calc_agenda
633 y,
634 t_field_raw,
635 vmr_field_raw,
636 z_field_raw,
637 pnd_field_raw,
638 p_grid,
639 sensor_los,
640 cloudbox_on,
641 cloudbox_limits,
642 z_surface,
643 met_profile_calc_agenda);
644
645 //putting in the spectra *y* for each profile
646 ybatch[i] = y;
647
648 } // closing the loop over profile basenames
649}
650
651/* Workspace method: Doxygen documentation will be auto-generated */
653 ArrayOfTensor7& dobatch_cloudbox_field,
654 ArrayOfTensor5& dobatch_radiance_field,
655 ArrayOfTensor4& dobatch_irradiance_field,
656 ArrayOfTensor5& dobatch_spectral_irradiance_field,
657 const Index& ybatch_start,
658 const Index& ybatch_n,
659 const Agenda& dobatch_calc_agenda,
660 const Index& robust,
661 const Verbosity& verbosity) {
663
664 Index first_ybatch_index = 0;
665
666 ArrayOfString fail_msg;
667 bool do_abort = false;
668
669 // We allow a start index ybatch_start that is different from 0. We
670 // will calculate ybatch_n jobs starting at the start
671 // index. Internally, we count from zero, which is the right
672 // counting for the output array ybatch. When we call
673 // ybatch_calc_agenda, we add ybatch_start to the internal index
674 // count.
675
676 // We create a counter, so that we can generate nice output about
677 // how many jobs are already done. (All parallel threads later will
678 // increment this, so that we really get an accurate total count!)
679 Index job_counter = 0;
680
681 // Resize the output arrays:
682 dobatch_cloudbox_field.resize(ybatch_n);
683 dobatch_radiance_field.resize(ybatch_n);
684 dobatch_irradiance_field.resize(ybatch_n);
685 dobatch_spectral_irradiance_field.resize(ybatch_n);
686
687 // We have to make a local copy of the Workspace and the agendas because
688 // only non-reference types can be declared firstprivate in OpenMP
689 Workspace l_ws(ws);
690 Agenda l_dobatch_calc_agenda(dobatch_calc_agenda);
691
692 // Go through the batch:
693
694 if (ybatch_n)
695#pragma omp parallel for schedule(dynamic) if (!arts_omp_in_parallel() && \
696 ybatch_n > 1) \
697 firstprivate(l_ws, l_dobatch_calc_agenda)
698 for (Index ybatch_index = first_ybatch_index; ybatch_index < ybatch_n;
699 ybatch_index++) {
700 Index l_job_counter; // Thread-local copy of job counter.
701
702 if (do_abort) continue;
703#pragma omp critical(dobatchCalc_job_counter)
704 { l_job_counter = ++job_counter; }
705
706 {
707 ostringstream os;
708 os << " Job " << l_job_counter << " of " << ybatch_n << ", Index "
709 << ybatch_start + ybatch_index << ", Thread-Id "
710 << arts_omp_get_thread_num() << "\n";
711 out2 << os.str();
712 }
713
714 try {
715 Tensor7 cloudbox_field;
716 Tensor5 radiance_field;
717 Tensor4 irradiance_field;
718 Tensor5 spectral_irradiance_field;
719
721 cloudbox_field,
722 radiance_field,
723 irradiance_field,
724 spectral_irradiance_field,
725 ybatch_start + ybatch_index,
726 l_dobatch_calc_agenda);
727
728#pragma omp critical(dobatchCalc_assign_cloudbox_field)
729 dobatch_cloudbox_field[ybatch_index] = cloudbox_field;
730#pragma omp critical(dobatchCalc_assign_radiance_field)
731 dobatch_radiance_field[ybatch_index] = radiance_field;
732#pragma omp critical(dobatchCalc_assign_irradiance_field)
733 dobatch_irradiance_field[ybatch_index] = irradiance_field;
734#pragma omp critical(dobatchCalc_assign_spectral_irradiance_field)
735 dobatch_spectral_irradiance_field[ybatch_index] =
736 spectral_irradiance_field;
737
738 } catch (const std::exception& e) {
739 if (robust && !do_abort) {
740 ostringstream os;
741 os << "WARNING! Job at ybatch_index " << ybatch_start + ybatch_index
742 << " failed.\n"
743 << "element in output variables will be empty for this job.\n"
744 << "The runtime error produced was:\n"
745 << e.what() << "\n";
746 out0 << os.str();
747 } else {
748 // The user wants the batch job to fail if one of the
749 // jobs goes wrong.
750#pragma omp critical(dobatchCalc_setabort)
751 do_abort = true;
752
753 ostringstream os;
754 os << " Job at ybatch_index " << ybatch_start + ybatch_index
755 << " failed. Aborting...\n";
756 out1 << os.str();
757 }
758 ostringstream os;
759 os << "Run-time error at ybatch_index " << ybatch_start + ybatch_index
760 << ": \n"
761 << e.what();
762#pragma omp critical(dobatchCalc_push_fail_msg)
763 fail_msg.push_back(os.str());
764 }
765 }
766
767 if (fail_msg.nelem()) {
768 ostringstream os;
769
770 if (!do_abort) os << "\nError messages from failed batch cases:\n";
771 for (ArrayOfString::const_iterator it = fail_msg.begin();
772 it != fail_msg.end();
773 it++)
774 os << *it << '\n';
775
776 if (do_abort)
777 throw runtime_error(os.str());
778 else
779 out0 << os.str();
780 }
781}
Index TotalNumberOfElements(const Array< Array< base > > &aa)
Determine total number of elements in an ArrayOfArray.
Definition: array.h:343
The global header file for ARTS.
void * data
int arts_omp_get_thread_num()
Wrapper for omp_get_thread_num.
Definition: arts_omp.cc:76
Header file for helper functions for OpenMP.
void forloop_agendaExecute(Workspace &ws, const Index forloop_index, const Agenda &input_agenda)
Definition: auto_md.cc:23145
void dobatch_calc_agendaExecute(Workspace &ws, Tensor7 &spectral_radiance_field, Tensor5 &radiance_field, Tensor4 &irradiance_field, Tensor5 &spectral_irradiance_field, const Index ybatch_index, const Agenda &input_agenda)
Definition: auto_md.cc:22947
void ybatch_calc_agendaExecute(Workspace &ws, Vector &y, ArrayOfVector &y_aux, Matrix &jacobian, const Index ybatch_index, const Agenda &input_agenda)
Definition: auto_md.cc:24655
void met_profile_calc_agendaExecute(Workspace &ws, Vector &y, const GriddedField3 &t_field_raw, const ArrayOfGriddedField3 &vmr_field_raw, const GriddedField3 &z_field_raw, const ArrayOfGriddedField3 &pnd_field_raw, const Vector &p_grid, const Matrix &sensor_los, const Index cloudbox_on, const ArrayOfIndex &cloudbox_limits, const Matrix &z_surface, const Agenda &input_agenda)
Definition: auto_md.cc:23962
void VectorNLogSpace(Vector &out, const Index &nelem, const Numeric &start, const Numeric &stop, const Verbosity &verbosity)
WORKSPACE METHOD: VectorNLogSpace.
void cloudboxSetManually(Index &cloudbox_on, ArrayOfIndex &cloudbox_limits, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Numeric &p1, const Numeric &p2, const Numeric &lat1, const Numeric &lat2, const Numeric &lon1, const Numeric &lon2, const Verbosity &verbosity)
WORKSPACE METHOD: cloudboxSetManually.
Definition: m_cloudbox.cc:332
The Agenda class.
Definition: agenda_class.h:44
Index nelem() const ARTS_NOEXCEPT
Number of elements.
Definition: array.h:195
Index nrows() const ARTS_NOEXCEPT
Returns the number of rows.
Definition: matpackI.cc:433
Index ncols() const ARTS_NOEXCEPT
Returns the number of columns.
Definition: matpackI.cc:436
A constant view of a Vector.
Definition: matpackI.h:489
Index nelem() const ARTS_NOEXCEPT
Returns the number of elements.
Definition: matpackI.cc:48
const Vector & get_numeric_grid(Index i) const
Get a numeric grid.
The Matrix class.
Definition: matpackI.h:1225
void resize(Index r, Index c)
Resize function.
Definition: matpackI.cc:1056
The range class.
Definition: matpackI.h:165
The Tensor4 class.
Definition: matpackIV.h:421
The Tensor5 class.
Definition: matpackV.h:506
The Tensor7 class.
Definition: matpackVII.h:2382
The Vector class.
Definition: matpackI.h:876
void resize(Index n)
Resize function.
Definition: matpackI.cc:408
Workspace class.
Definition: workspace_ng.h:40
#define abs(x)
const Numeric RAD2DEG
void ForLoop(Workspace &ws, const Agenda &forloop_agenda, const Index &start, const Index &stop, const Index &step, const Verbosity &verbosity)
WORKSPACE METHOD: ForLoop.
Definition: m_batch.cc:62
const Numeric PI
const Index GFIELD3_P_GRID
void ybatchMetProfilesClear(Workspace &ws, ArrayOfVector &ybatch, const ArrayOfArrayOfSpeciesTag &abs_species, const Agenda &met_profile_calc_agenda, const Vector &f_grid, const Matrix &met_amsu_data, const Matrix &sensor_pos, const Vector &refellipsoid, const Index &nelem_p_grid, const String &met_profile_path, const Verbosity &verbosity)
WORKSPACE METHOD: ybatchMetProfilesClear.
Definition: m_batch.cc:476
const Numeric DEG2RAD
void ybatchCalc(Workspace &ws, ArrayOfVector &ybatch, ArrayOfArrayOfVector &ybatch_aux, ArrayOfMatrix &ybatch_jacobians, const Index &ybatch_start, const Index &ybatch_n, const Agenda &ybatch_calc_agenda, const Index &robust, const Verbosity &verbosity)
WORKSPACE METHOD: ybatchCalc.
Definition: m_batch.cc:79
void ybatchMetProfiles(Workspace &ws, ArrayOfVector &ybatch, const ArrayOfArrayOfSpeciesTag &abs_species, const Agenda &met_profile_calc_agenda, const Vector &f_grid, const Matrix &met_amsu_data, const Matrix &sensor_pos, const Vector &refellipsoid, const Vector &lat_grid, const Vector &lon_grid, const Index &atmosphere_dim, const ArrayOfArrayOfSingleScatteringData &scat_data, const Index &nelem_p_grid, const String &met_profile_path, const String &met_profile_pnd_path, const Verbosity &verbosity)
WORKSPACE METHOD: ybatchMetProfiles.
Definition: m_batch.cc:238
void DOBatchCalc(Workspace &ws, ArrayOfTensor7 &dobatch_cloudbox_field, ArrayOfTensor5 &dobatch_radiance_field, ArrayOfTensor4 &dobatch_irradiance_field, ArrayOfTensor5 &dobatch_spectral_irradiance_field, const Index &ybatch_start, const Index &ybatch_n, const Agenda &dobatch_calc_agenda, const Index &robust, const Verbosity &verbosity)
WORKSPACE METHOD: DOBatchCalc.
Definition: m_batch.cc:652
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
const Joker joker
#define CREATE_OUT1
Definition: messages.h:205
#define CREATE_OUTS
Definition: messages.h:209
VectorView std(VectorView std, const Vector &y, const ArrayOfVector &ys, const Index start, const Index end_tmp)
Compute the standard deviation of the ranged ys.
Definition: raw.cc:205
constexpr Rational start(Rational Ju, Rational Jl, Polarization type) noexcept
Gives the lowest M for a polarization type of this transition.
Definition: zeemandata.h:78
This file contains declerations of functions of physical character.
Declaration of functions in rte.cc.
This file contains basic functions to handle XML data files.
void xml_write_to_file(const String &filename, const T &type, const FileType ftype, const Index no_clobber, const Verbosity &verbosity)
Write data to XML file.
Definition: xml_io.h:181
void xml_read_from_file(const String &filename, T &type, const Verbosity &verbosity)
Reads data from XML file.
Definition: xml_io.h:160
@ FILE_TYPE_ASCII
Definition: xml_io_base.h:41