Athena++/Atmosphere
Planetary Atmosphere Simulator
meshblock.cpp
Go to the documentation of this file.
1 //========================================================================================
2 // Athena++ astrophysical MHD code
3 // Copyright(C) 2014 James M. Stone <jmstone@princeton.edu> and other code contributors
4 // Licensed under the 3-clause BSD License, see LICENSE file for details
5 //========================================================================================
7 // \brief implementation of functions in MeshBlock class
8 
9 // C headers
10 
11 // C++ headers
12 #include <algorithm> // sort()
13 #include <cstdlib>
14 #include <cstring> // memcpy()
15 #include <ctime> // clock(), CLOCKS_PER_SEC, clock_t
16 #include <iomanip>
17 #include <iostream>
18 #include <sstream>
19 #include <stdexcept> // runtime_error
20 #include <string> // c_str()
21 
22 // Athena++ headers
23 #include "../athena.hpp"
24 #include "../athena_arrays.hpp"
25 #include "../bvals/bvals.hpp"
26 #include "../coordinates/coordinates.hpp"
27 #include "../eos/eos.hpp"
28 #include "../fft/athena_fft.hpp"
29 #include "../field/field.hpp"
30 #include "../globals.hpp"
31 #include "../gravity/gravity.hpp"
32 #include "../gravity/mg_gravity.hpp"
33 #include "../hydro/hydro.hpp"
34 #include "../parameter_input.hpp"
35 #include "../reconstruct/reconstruction.hpp"
36 #include "../scalars/scalars.hpp"
37 #include "../utils/buffer_utils.hpp"
38 #include "mesh.hpp"
39 #include "mesh_refinement.hpp"
40 #include "meshblock_tree.hpp"
41 #include "../thermodynamics/thermodynamics.hpp"
42 #include "../chemistry/chemistry.hpp"
43 #include "../radiation/radiation.hpp"
44 #include "../physics/physics.hpp"
45 #include "../diagnostics/diagnostics.hpp"
46 #include "../debugger/debugger.hpp"
47 #include "../particles/particles.hpp"
48 #include "../inversion/inversion.hpp"
49 #include "../communicator/communicator.hpp"
50 
51 //----------------------------------------------------------------------------------------
52 // MeshBlock constructor: constructs coordinate, boundary condition, hydro, field
53 // and mesh refinement objects.
54 
55 MeshBlock::MeshBlock(int igid, int ilid, LogicalLocation iloc, RegionSize input_block,
56  BoundaryFlag *input_bcs, Mesh *pm, ParameterInput *pin,
57  int igflag, bool ref_flag) :
58  pmy_mesh(pm), loc(iloc), block_size(input_block),
59  gid(igid), lid(ilid), gflag(igflag), nuser_out_var(), prev(nullptr), next(nullptr),
60  new_block_dt_{}, new_block_dt_hyperbolic_{}, new_block_dt_parabolic_{},
61  new_block_dt_user_{},
62  nreal_user_meshblock_data_(), nint_user_meshblock_data_(), cost_(1.0) {
63  // initialize grid indices
64  is = NGHOST;
65  ie = is + block_size.nx1 - 1;
66 
67  ncells1 = block_size.nx1 + 2*NGHOST;
68  ncc1 = block_size.nx1/2 + 2*NGHOST;
69  if (pmy_mesh->f2) {
70  js = NGHOST;
71  je = js + block_size.nx2 - 1;
72  ncells2 = block_size.nx2 + 2*NGHOST;
73  ncc2 = block_size.nx2/2 + 2*NGHOST;
74  } else {
75  js = je = 0;
76  ncells2 = 1;
77  ncc2 = 1;
78  }
79 
80  if (pmy_mesh->f3) {
81  ks = NGHOST;
82  ke = ks + block_size.nx3 - 1;
83  ncells3 = block_size.nx3 + 2*NGHOST;
84  ncc3 = block_size.nx3/2 + 2*NGHOST;
85  } else {
86  ks = ke = 0;
87  ncells3 = 1;
88  ncc3 = 1;
89  }
90 
91  if (pm->multilevel) {
92  cnghost = (NGHOST + 1)/2 + 1;
93  cis = NGHOST; cie = cis + block_size.nx1/2 - 1;
94  cjs = cje = cks = cke = 0;
95  if (pmy_mesh->f2) // 2D or 3D
96  cjs = NGHOST, cje = cjs + block_size.nx2/2 - 1;
97  if (pmy_mesh->f3) // 3D
98  cks = NGHOST, cke = cks + block_size.nx3/2 - 1;
99  }
100 
101  // (probably don't need to preallocate space for references in these vectors)
102  vars_cc_.reserve(3);
103  vars_fc_.reserve(3);
104 
105  // construct objects stored in MeshBlock class. Note in particular that the initial
106  // conditions for the simulation are set in problem generator called from main, not
107  // in the Hydro constructor
108 
109  // mesh-related objects
110  // debug
111  pdebug = new Debugger(this);
112 
113  // Boundary
114  pbval = new BoundaryValues(this, input_bcs, pin);
115 
116  // Coordinates
117  if (std::strcmp(COORDINATE_SYSTEM, "cartesian") == 0) {
118  pcoord = new Cartesian(this, pin, false);
119  } else if (std::strcmp(COORDINATE_SYSTEM, "cylindrical") == 0) {
120  pcoord = new Cylindrical(this, pin, false);
121  } else if (std::strcmp(COORDINATE_SYSTEM, "spherical_polar") == 0) {
122  pcoord = new SphericalPolar(this, pin, false);
123  } else if (std::strcmp(COORDINATE_SYSTEM, "minkowski") == 0) {
124  pcoord = new Minkowski(this, pin, false);
125  } else if (std::strcmp(COORDINATE_SYSTEM, "schwarzschild") == 0) {
126  pcoord = new Schwarzschild(this, pin, false);
127  } else if (std::strcmp(COORDINATE_SYSTEM, "kerr-schild") == 0) {
128  pcoord = new KerrSchild(this, pin, false);
129  } else if (std::strcmp(COORDINATE_SYSTEM, "gr_user") == 0) {
130  pcoord = new GRUser(this, pin, false);
131  }
132 
133  // Reconstruction: constructor may implicitly depend on Coordinates, and PPM variable
134  // floors depend on EOS, but EOS isn't needed in Reconstruction constructor-> this is ok
135  precon = new Reconstruction(this, pin);
136 
137  if (pm->multilevel) pmr = new MeshRefinement(this, pin);
138 
139  // physics-related, per-MeshBlock objects: may depend on Coordinates for diffusion
140  // terms, and may enroll quantities in AMR and BoundaryVariable objs. in BoundaryValues
141 
142  // TODO(felker): prepare this section of the MeshBlock ctor to become more complicated
143  // for several extensions:
144  // 1) allow solver to compile without a Hydro class (or with a Hydro class for the
145  // background fluid that is not dynamically evolved)
146  // 2) MPI ranks containing MeshBlocks that solve a subset of the physics, e.g. Gravity
147  // but not Hydro.
148  // 3) MAGNETIC_FIELDS_ENABLED, SELF_GRAVITY_ENABLED, NSCALARS, (future) FLUID_ENABLED,
149  // etc. become runtime switches
150 
151  // if (FLUID_ENABLED) {
152  // if (this->hydro_block)
153  phydro = new Hydro(this, pin);
154  // } else
155  // }
156  // Regardless, advance MeshBlock's local counter (initialized to bvars_next_phys_id=1)
157  // Greedy reservation of phys IDs (only 1 of 2 needed for Hydro if multilevel==false)
158  pbval->AdvanceCounterPhysID(HydroBoundaryVariable::max_phys_id);
159  // }
160  if (MAGNETIC_FIELDS_ENABLED) {
161  // if (this->field_block)
162  pfield = new Field(this, pin);
163  pbval->AdvanceCounterPhysID(FaceCenteredBoundaryVariable::max_phys_id);
164  }
165  if (SELF_GRAVITY_ENABLED) {
166  // if (this->grav_block)
167  pgrav = new Gravity(this, pin);
168  pbval->AdvanceCounterPhysID(CellCenteredBoundaryVariable::max_phys_id);
169  if (SELF_GRAVITY_ENABLED == 2)
170  pmg = new MGGravity(pmy_mesh->pmgrd, this);
171  }
172  if (NSCALARS > 0) {
173  // if (this->scalars_block)
174  pscalars = new PassiveScalars(this, pin);
175  pbval->AdvanceCounterPhysID(CellCenteredBoundaryVariable::max_phys_id);
176  }
177  // KGF: suboptimal solution, since developer must copy/paste BoundaryVariable derived
178  // class type that is used in each PassiveScalars, Gravity, Field, Hydro, ... etc. class
179  // in order to correctly advance the BoundaryValues::bvars_next_phys_id_ local counter.
180 
181  // TODO(felker): check that local counter pbval->bvars_next_phys_id_ agrees with shared
182  // Mesh::next_phys_id_ counter (including non-BoundaryVariable / per-MeshBlock reserved
183  // values). Compare both private member variables via BoundaryValues::CheckCounterPhysID
184 
185  peos = new EquationOfState(this, pin);
186  pthermo = new Thermodynamics(this, pin);
187  prad = new Radiation(this, pin);
188  pphy = new Physics(this, pin);
189  // initialize particles and remove the first one (header)
190  ppart = new Particles(this, pin);
191  if (ppart->next != nullptr) {
192  ppart = ppart->next;
193  delete ppart->prev;
194  } else {
195  delete ppart;
196  ppart = nullptr;
197  }
198  pchem = new Chemistry(this, pin);
199  pdiag = new Diagnostics(this, pin);
200  pfit = nullptr;
201  pcomm = new Communicator(this);
202 
203  // Create user mesh data
205 
206  return;
207 }
208 
209 //----------------------------------------------------------------------------------------
210 // MeshBlock constructor for restarts
211 
212 MeshBlock::MeshBlock(int igid, int ilid, Mesh *pm, ParameterInput *pin,
213  LogicalLocation iloc, RegionSize input_block,
214  BoundaryFlag *input_bcs,
215  double icost, char *mbdata, int igflag) :
216  pmy_mesh(pm), loc(iloc), block_size(input_block),
217  gid(igid), lid(ilid), gflag(igflag), nuser_out_var(), prev(nullptr), next(nullptr),
218  new_block_dt_{}, new_block_dt_hyperbolic_{}, new_block_dt_parabolic_{},
219  new_block_dt_user_{},
220  nreal_user_meshblock_data_(), nint_user_meshblock_data_(), cost_(icost) {
221  // initialize grid indices
222  is = NGHOST;
223  ie = is + block_size.nx1 - 1;
224 
225  ncells1 = block_size.nx1 + 2*NGHOST;
226  ncc1 = block_size.nx1/2 + 2*NGHOST;
227  if (pmy_mesh->f2) {
228  js = NGHOST;
229  je = js + block_size.nx2 - 1;
230  ncells2 = block_size.nx2 + 2*NGHOST;
231  ncc2 = block_size.nx2/2 + 2*NGHOST;
232  } else {
233  js = je = 0;
234  ncells2 = 1;
235  ncc2 = 1;
236  }
237 
238  if (pmy_mesh->f3) {
239  ks = NGHOST;
240  ke = ks + block_size.nx3 - 1;
241  ncells3 = block_size.nx3 + 2*NGHOST;
242  ncc3 = block_size.nx3/2 + 2*NGHOST;
243  } else {
244  ks = ke = 0;
245  ncells3 = 1;
246  ncc3 = 1;
247  }
248 
249  if (pm->multilevel) {
250  cnghost = (NGHOST + 1)/2 + 1;
251  cis = NGHOST; cie = cis + block_size.nx1/2 - 1;
252  cjs = cje = cks = cke = 0;
253  if (pmy_mesh->f2) // 2D or 3D
254  cjs = NGHOST, cje = cjs + block_size.nx2/2 - 1;
255  if (pmy_mesh->f3) // 3D
256  cks = NGHOST, cke = cks + block_size.nx3/2 - 1;
257  }
258 
259  // (re-)create mesh-related objects in MeshBlock
260  // debug
261  pdebug = new Debugger(this);
262 
263  // Boundary
264  pbval = new BoundaryValues(this, input_bcs, pin);
265 
266  // Coordinates
267  if (std::strcmp(COORDINATE_SYSTEM, "cartesian") == 0) {
268  pcoord = new Cartesian(this, pin, false);
269  } else if (std::strcmp(COORDINATE_SYSTEM, "cylindrical") == 0) {
270  pcoord = new Cylindrical(this, pin, false);
271  } else if (std::strcmp(COORDINATE_SYSTEM, "spherical_polar") == 0) {
272  pcoord = new SphericalPolar(this, pin, false);
273  } else if (std::strcmp(COORDINATE_SYSTEM, "minkowski") == 0) {
274  pcoord = new Minkowski(this, pin, false);
275  } else if (std::strcmp(COORDINATE_SYSTEM, "schwarzschild") == 0) {
276  pcoord = new Schwarzschild(this, pin, false);
277  } else if (std::strcmp(COORDINATE_SYSTEM, "kerr-schild") == 0) {
278  pcoord = new KerrSchild(this, pin, false);
279  } else if (std::strcmp(COORDINATE_SYSTEM, "gr_user") == 0) {
280  pcoord = new GRUser(this, pin, false);
281  }
282 
283  // Reconstruction (constructor may implicitly depend on Coordinates)
284  precon = new Reconstruction(this, pin);
285 
286  if (pm->multilevel) pmr = new MeshRefinement(this, pin);
287 
288  // (re-)create physics-related objects in MeshBlock
289 
290  // if (FLUID_ENABLED) {
291  // if (this->hydro_block)
292  phydro = new Hydro(this, pin);
293  // } else
294  // }
295  // Regardless, advance MeshBlock's local counter (initialized to bvars_next_phys_id=1)
296  // Greedy reservation of phys IDs (only 1 of 2 needed for Hydro if multilevel==false)
297  pbval->AdvanceCounterPhysID(HydroBoundaryVariable::max_phys_id);
298  // }
299  if (MAGNETIC_FIELDS_ENABLED) {
300  // if (this->field_block)
301  pfield = new Field(this, pin);
302  pbval->AdvanceCounterPhysID(FaceCenteredBoundaryVariable::max_phys_id);
303  }
304  if (SELF_GRAVITY_ENABLED) {
305  // if (this->grav_block)
306  pgrav = new Gravity(this, pin);
307  pbval->AdvanceCounterPhysID(CellCenteredBoundaryVariable::max_phys_id);
308  if (SELF_GRAVITY_ENABLED == 2)
309  pmg = new MGGravity(pmy_mesh->pmgrd, this);
310  }
311 
312  if (NSCALARS > 0) {
313  // if (this->scalars_block)
314  pscalars = new PassiveScalars(this, pin);
315  pbval->AdvanceCounterPhysID(CellCenteredBoundaryVariable::max_phys_id);
316  }
317 
318  peos = new EquationOfState(this, pin);
319  pthermo = new Thermodynamics(this, pin);
320  prad = new Radiation(this, pin);
321  pphy = new Physics(this, pin);
322  // initialize particles and remove the first one (header)
323  ppart = new Particles(this, pin);
324  if (ppart->next != nullptr) {
325  ppart = ppart->next;
326  delete ppart->prev;
327  } else {
328  delete ppart;
329  ppart = nullptr;
330  }
331  pchem = new Chemistry(this, pin);
332  pdiag = new Diagnostics(this, pin);
333  pfit = nullptr;
334  pcomm = new Communicator(this);
335 
337 
338  std::size_t os = 0;
339  // NEW_OUTPUT_TYPES:
340 
341  // load hydro and field data
342  std::memcpy(phydro->u.data(), &(mbdata[os]), phydro->u.GetSizeInBytes());
343  // outflow boundary condition
344  if (pbval->block_bcs[inner_x1] == BoundaryFlag::outflow) {
345  int kl = block_size.nx3 == 1 ? ks : ks-NGHOST;
346  int ku = block_size.nx3 == 1 ? ke : ke+NGHOST;
347  int jl = block_size.nx2 == 1 ? js : js-NGHOST;
348  int ju = block_size.nx2 == 1 ? je : je+NGHOST;
350  phydro->w, pfield->bcc, pcoord,
351  is-NGHOST, is-1, jl, ju, kl, ku);
352  phydro->w1 = phydro->w;
353  }
354  // load it into the other memory register(s) too
355  std::memcpy(phydro->u1.data(), &(mbdata[os]), phydro->u1.GetSizeInBytes());
356  os += phydro->u.GetSizeInBytes();
357  if (GENERAL_RELATIVITY) {
358  std::memcpy(phydro->w.data(), &(mbdata[os]), phydro->w.GetSizeInBytes());
359  os += phydro->w.GetSizeInBytes();
360  std::memcpy(phydro->w1.data(), &(mbdata[os]), phydro->w1.GetSizeInBytes());
361  os += phydro->w1.GetSizeInBytes();
362  }
363  if (MAGNETIC_FIELDS_ENABLED) {
364  std::memcpy(pfield->b.x1f.data(), &(mbdata[os]), pfield->b.x1f.GetSizeInBytes());
365  std::memcpy(pfield->b1.x1f.data(), &(mbdata[os]), pfield->b1.x1f.GetSizeInBytes());
366  os += pfield->b.x1f.GetSizeInBytes();
367  std::memcpy(pfield->b.x2f.data(), &(mbdata[os]), pfield->b.x2f.GetSizeInBytes());
368  std::memcpy(pfield->b1.x2f.data(), &(mbdata[os]), pfield->b1.x2f.GetSizeInBytes());
369  os += pfield->b.x2f.GetSizeInBytes();
370  std::memcpy(pfield->b.x3f.data(), &(mbdata[os]), pfield->b.x3f.GetSizeInBytes());
371  std::memcpy(pfield->b1.x3f.data(), &(mbdata[os]), pfield->b1.x3f.GetSizeInBytes());
372  os += pfield->b.x3f.GetSizeInBytes();
373  }
374 
375  // (conserved variable) Passive scalars:
376  if (NSCALARS > 0) {
377  std::memcpy(pscalars->s.data(), &(mbdata[os]), pscalars->s.GetSizeInBytes());
378  // load it into the other memory register(s) too
379  std::memcpy(pscalars->s1.data(), &(mbdata[os]), pscalars->s1.GetSizeInBytes());
380  os += pscalars->s.GetSizeInBytes();
381  }
382 
383  // load user MeshBlock data
384  for (int n=0; n<nint_user_meshblock_data_; n++) {
385  std::memcpy(iuser_meshblock_data[n].data(), &(mbdata[os]),
386  iuser_meshblock_data[n].GetSizeInBytes());
388  }
389  for (int n=0; n<nreal_user_meshblock_data_; n++) {
390  std::memcpy(ruser_meshblock_data[n].data(), &(mbdata[os]),
391  ruser_meshblock_data[n].GetSizeInBytes());
393  }
394 
395  // load radiation data
396  os += prad->loadRestartData(&(mbdata[os]));
397 
398  // load physics data
399  os += pphy->LoadRestartData(&(mbdata[os]));
400 
401  // load particle data
402  if (ppart != nullptr)
403  os += ppart->LoadRestartData(&(mbdata[os]));
404 
405  return;
406 }
407 
408 //----------------------------------------------------------------------------------------
409 // MeshBlock destructor
410 
412  if (prev != nullptr) prev->next = next;
413  if (next != nullptr) next->prev = prev;
414 
415  delete pcoord;
416  delete precon;
417  if (pmy_mesh->multilevel) delete pmr;
418 
419  delete phydro;
420  if (MAGNETIC_FIELDS_ENABLED) delete pfield;
421  delete peos;
422  if (SELF_GRAVITY_ENABLED) delete pgrav;
423  if (NSCALARS > 0) delete pscalars;
424 
425  // BoundaryValues should be destructed AFTER all BoundaryVariable objects are destroyed
426  delete pbval;
427  // delete user output variables array
428  if (nuser_out_var > 0) {
429  delete [] user_out_var_names_;
430  delete [] user_out_var_longnames_;
431  delete [] user_out_var_units_;
432  }
433  // delete user MeshBlock data
436 
437  delete pthermo;
438  delete prad;
439  delete pphy;
440 
441  while (pdiag->prev != nullptr)
442  delete pdiag->prev;
443  while (pdiag->next != nullptr)
444  delete pdiag->next;
445  delete pdiag;
446 
447  while (pdebug->prev != nullptr)
448  delete pdebug->prev;
449  while (pdebug->next != nullptr)
450  delete pdebug->next;
451  delete pdebug;
452 
453  if (ppart != nullptr) {
454  while (ppart->prev != nullptr)
455  delete ppart->prev;
456  while (ppart->next != nullptr)
457  delete ppart->next;
458  delete ppart;
459  }
460  delete pchem;
461 
462  if (pfit != nullptr)
463  delete pfit;
464  delete pcomm;
465 }
466 
467 //----------------------------------------------------------------------------------------
469 // \brief Allocate Real AthenaArrays for user-defned data in MeshBlock
470 
472  if (nreal_user_meshblock_data_ != 0) {
473  std::stringstream msg;
474  msg << "### FATAL ERROR in MeshBlock::AllocateRealUserMeshBlockDataField"
475  << std::endl << "User MeshBlock data arrays are already allocated" << std::endl;
476  ATHENA_ERROR(msg);
477  }
480  return;
481 }
482 
483 //----------------------------------------------------------------------------------------
485 // \brief Allocate integer AthenaArrays for user-defned data in MeshBlock
486 
488  if (nint_user_meshblock_data_ != 0) {
489  std::stringstream msg;
490  msg << "### FATAL ERROR in MeshBlock::AllocateIntusermeshblockDataField"
491  << std::endl << "User MeshBlock data arrays are already allocated" << std::endl;
492  ATHENA_ERROR(msg);
493  return;
494  }
497  return;
498 }
499 
500 //----------------------------------------------------------------------------------------
502 // \brief Allocate user-defined output variables
503 
505  if (n <= 0) return;
506  if (nuser_out_var != 0) {
507  std::stringstream msg;
508  msg << "### FATAL ERROR in MeshBlock::AllocateUserOutputVariables"
509  << std::endl << "User output variables are already allocated." << std::endl;
510  ATHENA_ERROR(msg);
511  return;
512  }
513  nuser_out_var = n;
514  user_out_var.NewAthenaArray(nuser_out_var, ncells3, ncells2, ncells1);
515  user_out_var_names_ = new std::string[n];
516  user_out_var_longnames_ = new std::string[n];
517  user_out_var_units_ = new std::string[n];
518  return;
519 }
520 
521 
522 //----------------------------------------------------------------------------------------
524 // \brief set the user-defined output variable name
525 
526 void MeshBlock::SetUserOutputVariableName(int n, const char *name,
527  const char *long_name, const char *units) {
528  if (n >= nuser_out_var) {
529  std::stringstream msg;
530  msg << "### FATAL ERROR in MeshBlock::SetUserOutputVariableName"
531  << std::endl << "User output variable is not allocated." << std::endl;
532  ATHENA_ERROR(msg);
533  return;
534  }
535  user_out_var_names_[n] = name;
536  user_out_var_longnames_[n] = long_name;
537  user_out_var_units_[n] = units;
538  return;
539 }
540 
541 //----------------------------------------------------------------------------------------
543 // \brief Calculate the block data size required for restart.
544 
546  std::size_t size;
547  // NEW_OUTPUT_TYPES:
548  size = phydro->u.GetSizeInBytes();
549  if (GENERAL_RELATIVITY) {
550  size += phydro->w.GetSizeInBytes();
551  size += phydro->w1.GetSizeInBytes();
552  }
553  if (MAGNETIC_FIELDS_ENABLED)
554  size += (pfield->b.x1f.GetSizeInBytes() + pfield->b.x2f.GetSizeInBytes()
555  + pfield->b.x3f.GetSizeInBytes());
556  if (SELF_GRAVITY_ENABLED)
557  size += pgrav->phi.GetSizeInBytes();
558  if (NSCALARS > 0)
559  size += pscalars->s.GetSizeInBytes();
560 
561  // calculate user MeshBlock data size
562  for (int n=0; n<nint_user_meshblock_data_; n++)
563  size += iuser_meshblock_data[n].GetSizeInBytes();
564  for (int n=0; n<nreal_user_meshblock_data_; n++)
565  size += ruser_meshblock_data[n].GetSizeInBytes();
566 
567  // radiation data
568  size += prad->getRestartDataSizeInBytes();
569 
570  // physics data
571  size += pphy->RestartDataSizeInBytes();
572 
573  // particle data
574  if (ppart != nullptr)
575  size += ppart->RestartDataSizeInBytes();
576 
577  return size;
578 }
579 
580 //----------------------------------------------------------------------------------------
582 // \brief stop time measurement and accumulate it in the MeshBlock cost
583 
585  if (pmy_mesh->lb_manual_) {
586  cost_ = std::min(cost, TINY_NUMBER);
587  pmy_mesh->lb_flag_ = true;
588  }
589 }
590 
591 //----------------------------------------------------------------------------------------
593 // \brief reset the MeshBlock cost for automatic load balancing
594 
596  if (pmy_mesh->lb_automatic_) cost_ = TINY_NUMBER;
597 }
598 
599 //----------------------------------------------------------------------------------------
601 // \brief start time measurement for automatic load balancing
602 
604  if (pmy_mesh->lb_automatic_) {
605 #ifdef OPENMP_PARALLEL
606  lb_time_ = omp_get_wtime();
607 #else
608  lb_time_ = static_cast<double>(clock());
609 #endif
610  }
611 }
612 
613 //----------------------------------------------------------------------------------------
615 // \brief stop time measurement and accumulate it in the MeshBlock cost
616 
618  if (pmy_mesh->lb_automatic_) {
619 #ifdef OPENMP_PARALLEL
620  lb_time_ = omp_get_wtime() - lb_time_;
621 #else
622  lb_time_ = static_cast<double>(clock()) - lb_time_;
623 #endif
624  cost_ += lb_time_;
625  }
626 }
627 
628 
630  vars_cc_.push_back(pvar_cc);
631  return;
632 }
633 
634 
636  vars_fc_.push_back(pvar_fc);
637  return;
638 }
639 
640 
641 // TODO(felker): consider merging the MeshRefinement::pvars_cc/fc_ into the
642 // MeshBlock::pvars_cc/fc_. Would need to weaken the MeshBlock std::vector to use tuples
643 // of pointers instead of a std::vector of references, so that:
644 // - nullptr can be passed for the second entry if multilevel==false
645 // - we can rebind the pointers to Hydro for GR purposes in bvals_refine.cpp
646 // If GR, etc. in the future requires additional flexiblity from non-refinement load
647 // balancing, we will need to use ptrs instead of references anyways, and add:
648 
649 // void MeshBlock::SetHydroData(HydroBoundaryQuantity hydro_type)
650 // Hydro *ph = pmy_block_->phydro;
651 // // hard-coded assumption that, if multilevel, then Hydro is always present
652 // // and enrolled in mesh refinement in the first pvars_cc_ vector entry
653 // switch (hydro_type) {
654 // case (HydroBoundaryQuantity::cons): {
655 // pvars_cc_.front() = &ph->u;
656 // break;
657 // }
658 // case (HydroBoundaryQuantity::prim): {
659 // pvars_cc_.front() = &ph->w;
660 // break;
661 // }
662 // }
663 // return;
664 // }
std::size_t GetSizeInBytes() const
Diagnostics * next
Definition: diagnostics.hpp:18
Diagnostics * prev
Definition: diagnostics.hpp:18
void ConservedToPrimitive(AthenaArray< Real > &cons, const AthenaArray< Real > &prim_old, const FaceField &b, AthenaArray< Real > &prim, AthenaArray< Real > &bcc, Coordinates *pco, int il, int iu, int jl, int ju, int kl, int ku)
AthenaArray< Real > u1
Definition: hydro.hpp:47
AthenaArray< Real > w
Definition: hydro.hpp:46
AthenaArray< Real > u
Definition: hydro.hpp:46
AthenaArray< Real > w1
Definition: hydro.hpp:47
MeshRefinement * pmr
Definition: mesh.hpp:125
void AllocateUserOutputVariables(int n)
Definition: meshblock.cpp:504
std::string * user_out_var_units_
Definition: mesh.hpp:115
Reconstruction * precon
Definition: mesh.hpp:124
double cost_
Definition: mesh.hpp:193
Physics * pphy
Definition: mesh.hpp:138
void StopTimeMeasurement()
Definition: meshblock.cpp:617
BoundaryValues * pbval
Definition: mesh.hpp:123
Debugger * pdebug
Definition: mesh.hpp:140
int nint_user_meshblock_data_
Definition: mesh.hpp:176
Gravity * pgrav
Definition: mesh.hpp:130
int cie
Definition: mesh.hpp:102
int ncc2
Definition: mesh.hpp:99
friend class BoundaryValues
Definition: mesh.hpp:72
MeshBlock * prev
Definition: mesh.hpp:145
Communicator * pcomm
Definition: mesh.hpp:143
Inversion * pfit
Definition: mesh.hpp:142
void RegisterMeshBlockData(AthenaArray< Real > &pvar_cc)
Definition: meshblock.cpp:629
friend class Hydro
Definition: mesh.hpp:76
RegionSize block_size
Definition: mesh.hpp:94
Hydro * phydro
Definition: mesh.hpp:128
AthenaArray< Real > * ruser_meshblock_data
Definition: mesh.hpp:118
EquationOfState * peos
Definition: mesh.hpp:133
int cjs
Definition: mesh.hpp:102
Mesh * pmy_mesh
Definition: mesh.hpp:92
MGGravity * pmg
Definition: mesh.hpp:131
int ie
Definition: mesh.hpp:100
int cnghost
Definition: mesh.hpp:102
void AllocateRealUserMeshBlockDataField(int n)
Definition: meshblock.cpp:471
AthenaArray< Real > user_out_var
Definition: mesh.hpp:112
int js
Definition: mesh.hpp:100
std::vector< std::reference_wrapper< AthenaArray< Real > > > vars_cc_
Definition: mesh.hpp:177
int ncells1
Definition: mesh.hpp:97
Diagnostics * pdiag
Definition: mesh.hpp:139
int je
Definition: mesh.hpp:100
void AllocateIntUserMeshBlockDataField(int n)
Definition: meshblock.cpp:487
AthenaArray< int > * iuser_meshblock_data
Definition: mesh.hpp:119
std::string * user_out_var_longnames_
Definition: mesh.hpp:114
int ncc3
Definition: mesh.hpp:99
Particles * ppart
Definition: mesh.hpp:141
int ke
Definition: mesh.hpp:100
int cis
Definition: mesh.hpp:102
PassiveScalars * pscalars
Definition: mesh.hpp:132
int ncells3
Definition: mesh.hpp:97
void StartTimeMeasurement()
Definition: meshblock.cpp:603
double lb_time_
Definition: mesh.hpp:193
Radiation * prad
Definition: mesh.hpp:137
int ncc1
Definition: mesh.hpp:99
Chemistry * pchem
Definition: mesh.hpp:136
int is
Definition: mesh.hpp:100
int cks
Definition: mesh.hpp:102
int nuser_out_var
Definition: mesh.hpp:111
void ResetTimeMeasurement()
Definition: meshblock.cpp:595
int ks
Definition: mesh.hpp:100
Coordinates * pcoord
Definition: mesh.hpp:122
int cke
Definition: mesh.hpp:102
int cje
Definition: mesh.hpp:102
int nreal_user_meshblock_data_
Definition: mesh.hpp:176
void InitUserMeshBlockData(ParameterInput *pin)
void SetUserOutputVariableName(int n, const char *name, const char *long_name="", const char *units="")
Definition: meshblock.cpp:526
Thermodynamics * pthermo
Definition: mesh.hpp:135
std::vector< std::reference_wrapper< FaceField > > vars_fc_
Definition: mesh.hpp:178
int ncells2
Definition: mesh.hpp:97
std::size_t GetBlockSizeInBytes()
Definition: meshblock.cpp:545
void SetCostForLoadBalancing(double cost)
Definition: meshblock.cpp:584
MeshBlock(int igid, int ilid, LogicalLocation iloc, RegionSize input_size, BoundaryFlag *input_bcs, Mesh *pm, ParameterInput *pin, int igflag, bool ref_flag=false)
Definition: meshblock.cpp:55
MeshBlock * next
Definition: mesh.hpp:145
std::string * user_out_var_names_
Definition: mesh.hpp:113
Field * pfield
Definition: mesh.hpp:129
Definition: mesh.hpp:203
bool lb_flag_
Definition: mesh.hpp:321
const bool multilevel
Definition: mesh.hpp:246
bool lb_manual_
Definition: mesh.hpp:321
bool lb_automatic_
Definition: mesh.hpp:321
const bool f3
Definition: mesh.hpp:244
MGGravityDriver * pmgrd
Definition: mesh.hpp:263
const bool f2
Definition: mesh.hpp:244
AthenaArray< Real > s
Definition: scalars.hpp:35
AthenaArray< Real > s1
Definition: scalars.hpp:35
size_t loadRestartData(char *psrc)
Definition: radiation.cpp:277
size_t getRestartDataSizeInBytes()
Definition: radiation.cpp:249
double min(double x1, double x2, double x3)
Definition: core.h:15
int nx3
Definition: athena.hpp:88
int nx2
Definition: athena.hpp:88
int nx1
Definition: athena.hpp:88