Athena++/Atmosphere
Planetary Atmosphere Simulator
mesh.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 Mesh class
8 
9 // C headers
10 // pre-C11: needed before including inttypes.h, else won't define int64_t for C++ code
11 // #define __STDC_FORMAT_MACROS
12 
13 // C++ headers
14 #include <algorithm>
15 #include <cinttypes> // format macro "PRId64" for fixed-width integer type std::int64_t
16 #include <cmath> // std::abs(), std::pow()
17 #include <cstdint> // std::int64_t fixed-wdith integer type alias
18 #include <cstdlib>
19 #include <cstring> // std::memcpy()
20 #include <iomanip> // std::setprecision()
21 #include <iostream>
22 #include <limits>
23 #include <sstream>
24 #include <stdexcept> // runtime_error
25 #include <string> // c_str()
26 #include <vector>
27 
28 // Athena++ headers
29 #include "../athena.hpp"
30 #include "../athena_arrays.hpp"
31 #include "../bvals/bvals.hpp"
32 #include "../coordinates/coordinates.hpp"
33 #include "../eos/eos.hpp"
34 #include "../fft/athena_fft.hpp"
35 #include "../fft/turbulence.hpp"
36 #include "../field/field.hpp"
37 #include "../field/field_diffusion/field_diffusion.hpp"
38 #include "../globals.hpp"
39 #include "../gravity/fft_gravity.hpp"
40 #include "../gravity/gravity.hpp"
41 #include "../gravity/mg_gravity.hpp"
42 #include "../hydro/hydro.hpp"
43 #include "../hydro/hydro_diffusion/hydro_diffusion.hpp"
44 #include "../multigrid/multigrid.hpp"
45 #include "../outputs/io_wrapper.hpp"
46 #include "../parameter_input.hpp"
47 #include "../reconstruct/reconstruction.hpp"
48 #include "../scalars/scalars.hpp"
49 #include "../utils/buffer_utils.hpp"
50 #include "../physics/physics.hpp"
51 #include "mesh.hpp"
52 #include "mesh_refinement.hpp"
53 #include "meshblock_tree.hpp"
54 
55 // MPI/OpenMP header
56 #ifdef MPI_PARALLEL
57 #include <mpi.h>
58 #endif
59 
60 //----------------------------------------------------------------------------------------
61 // Mesh constructor, builds mesh at start of calculation using parameters in input file
62 
63 Mesh::Mesh(ParameterInput *pin, int mesh_test) :
64  // public members:
65  // aggregate initialization of RegionSize struct:
66  mesh_size{pin->GetReal("mesh", "x1min"), pin->GetReal("mesh", "x2min"),
67  pin->GetReal("mesh", "x3min"), pin->GetReal("mesh", "x1max"),
68  pin->GetReal("mesh", "x2max"), pin->GetReal("mesh", "x3max"),
69  pin->GetOrAddReal("mesh", "x1rat", RAT1),
70  pin->GetOrAddReal("mesh", "x2rat", 1.0),
71  pin->GetOrAddReal("mesh", "x3rat", 1.0),
72  pin->GetInteger("mesh", "nx1"), pin->GetInteger("mesh", "nx2"),
73  pin->GetInteger("mesh", "nx3") },
74  mesh_bcs{GetBoundaryFlag(pin->GetOrAddString("mesh", "ix1_bc", "none")),
75  GetBoundaryFlag(pin->GetOrAddString("mesh", "ox1_bc", "none")),
76  GetBoundaryFlag(pin->GetOrAddString("mesh", "ix2_bc", "none")),
77  GetBoundaryFlag(pin->GetOrAddString("mesh", "ox2_bc", "none")),
78  GetBoundaryFlag(pin->GetOrAddString("mesh", "ix3_bc", "none")),
79  GetBoundaryFlag(pin->GetOrAddString("mesh", "ox3_bc", "none"))},
80  f2(mesh_size.nx2 > 1 ? true : false), f3(mesh_size.nx3 > 1 ? true : false),
81  ndim(f3 ? 3 : (f2 ? 2 : 1)),
82  adaptive(pin->GetOrAddString("mesh", "refinement", "none") == "adaptive"
83  ? true : false),
84  multilevel((adaptive || pin->GetOrAddString("mesh", "refinement", "none") == "static")
85  ? true : false),
86  fluid_setup(GetFluidFormulation(pin->GetOrAddString("hydro", "active", "true"))),
87  start_time(pin->GetOrAddReal("time", "start_time", 0.0)), time(start_time),
88  tlim(pin->GetReal("time", "tlim")), dt(std::numeric_limits<Real>::max()),
89  dt_hyperbolic(dt), dt_parabolic(dt), dt_user(dt),
90  cfl_number(pin->GetReal("time", "cfl_number")),
91  nlim(pin->GetOrAddInteger("time", "nlim", -1)), ncycle(),
92  ncycle_out(pin->GetOrAddInteger("time", "ncycle_out", 1)),
93  dt_diagnostics(pin->GetOrAddInteger("time", "dt_diagnostics", -1)),
94  muj(), nuj(), muj_tilde(),
95  nbnew(), nbdel(),
96  step_since_lb(), gflag(), turb_flag(),
97  // private members:
98  next_phys_id_(), num_mesh_threads_(pin->GetOrAddInteger("mesh", "num_threads", 1)),
99  tree(this),
100  use_uniform_meshgen_fn_{true, true, true},
101  nreal_user_mesh_data_(), nint_user_mesh_data_(), nuser_history_output_(),
102  four_pi_G_(), grav_eps_(-1.0), grav_mean_rho_(-1.0),
103  lb_flag_(true), lb_automatic_(), lb_manual_(),
106  BoundaryFunction_{nullptr, nullptr, nullptr, nullptr, nullptr, nullptr},
107  AMRFlag_{}, UserSourceTerm_{}, UserTimeStep_{}, ViscosityCoeff_{},
108  ConductionCoeff_{}, FieldDiffusivity_{},
109  MGGravityBoundaryFunction_{MGPeriodicInnerX1, MGPeriodicOuterX1, MGPeriodicInnerX2,
110  MGPeriodicOuterX2, MGPeriodicInnerX3, MGPeriodicOuterX3} {
111  std::stringstream msg;
112  RegionSize block_size;
113  MeshBlock *pfirst{};
114  BoundaryFlag block_bcs[6];
115  std::int64_t nbmax;
116 
117  // mesh test
118  if (mesh_test > 0) Globals::nranks = mesh_test;
119 
120 #ifdef MPI_PARALLEL
121  // reserve phys=0 for former TAG_AMR=8; now hard-coded in Mesh::CreateAMRMPITag()
122  next_phys_id_ = 1;
124 #endif
125 
126  // check number of OpenMP threads for mesh
127  if (num_mesh_threads_ < 1) {
128  msg << "### FATAL ERROR in Mesh constructor" << std::endl
129  << "Number of OpenMP threads must be >= 1, but num_threads="
130  << num_mesh_threads_ << std::endl;
131  ATHENA_ERROR(msg);
132  }
133 
134  // check number of grid cells in root level of mesh from input file.
135  if (mesh_size.nx1 < 4) {
136  msg << "### FATAL ERROR in Mesh constructor" << std::endl
137  << "In mesh block in input file nx1 must be >= 4, but nx1="
138  << mesh_size.nx1 << std::endl;
139  ATHENA_ERROR(msg);
140  }
141  if (mesh_size.nx2 < 1) {
142  msg << "### FATAL ERROR in Mesh constructor" << std::endl
143  << "In mesh block in input file nx2 must be >= 1, but nx2="
144  << mesh_size.nx2 << std::endl;
145  ATHENA_ERROR(msg);
146  }
147  if (mesh_size.nx3 < 1) {
148  msg << "### FATAL ERROR in Mesh constructor" << std::endl
149  << "In mesh block in input file nx3 must be >= 1, but nx3="
150  << mesh_size.nx3 << std::endl;
151  ATHENA_ERROR(msg);
152  }
153  if (mesh_size.nx2 == 1 && mesh_size.nx3 > 1) {
154  msg << "### FATAL ERROR in Mesh constructor" << std::endl
155  << "In mesh block in input file: nx2=1, nx3=" << mesh_size.nx3
156  << ", 2D problems in x1-x3 plane not supported" << std::endl;
157  ATHENA_ERROR(msg);
158  }
159 
160  // check physical size of mesh (root level) from input file.
161  if (mesh_size.x1max <= mesh_size.x1min) {
162  msg << "### FATAL ERROR in Mesh constructor" << std::endl
163  << "Input x1max must be larger than x1min: x1min=" << mesh_size.x1min
164  << " x1max=" << mesh_size.x1max << std::endl;
165  ATHENA_ERROR(msg);
166  }
167  if (mesh_size.x2max <= mesh_size.x2min) {
168  msg << "### FATAL ERROR in Mesh constructor" << std::endl
169  << "Input x2max must be larger than x2min: x2min=" << mesh_size.x2min
170  << " x2max=" << mesh_size.x2max << std::endl;
171  ATHENA_ERROR(msg);
172  }
173  if (mesh_size.x3max <= mesh_size.x3min) {
174  msg << "### FATAL ERROR in Mesh constructor" << std::endl
175  << "Input x3max must be larger than x3min: x3min=" << mesh_size.x3min
176  << " x3max=" << mesh_size.x3max << std::endl;
177  ATHENA_ERROR(msg);
178  }
179 
180  // check the consistency of the periodic boundaries
181  if ( ((mesh_bcs[BoundaryFace::inner_x1] == BoundaryFlag::periodic
182  && mesh_bcs[BoundaryFace::outer_x1] != BoundaryFlag::periodic)
183  || (mesh_bcs[BoundaryFace::inner_x1] != BoundaryFlag::periodic
184  && mesh_bcs[BoundaryFace::outer_x1] == BoundaryFlag::periodic))
185  || (mesh_size.nx2 > 1
186  && ((mesh_bcs[BoundaryFace::inner_x2] == BoundaryFlag::periodic
187  && mesh_bcs[BoundaryFace::outer_x2] != BoundaryFlag::periodic)
188  || (mesh_bcs[BoundaryFace::inner_x2] != BoundaryFlag::periodic
189  && mesh_bcs[BoundaryFace::outer_x2] == BoundaryFlag::periodic)))
190  || (mesh_size.nx3 > 1
191  && ((mesh_bcs[BoundaryFace::inner_x3] == BoundaryFlag::periodic
192  && mesh_bcs[BoundaryFace::outer_x3] != BoundaryFlag::periodic)
193  || (mesh_bcs[BoundaryFace::inner_x3] != BoundaryFlag::periodic
194  && mesh_bcs[BoundaryFace::outer_x3] == BoundaryFlag::periodic)))) {
195  msg << "### FATAL ERROR in Mesh constructor" << std::endl
196  << "When periodic boundaries are in use, both sides must be periodic."
197  << std::endl;
198  ATHENA_ERROR(msg);
199  }
200  if ( ((mesh_bcs[BoundaryFace::inner_x1] == BoundaryFlag::shear_periodic
201  && mesh_bcs[BoundaryFace::outer_x1] != BoundaryFlag::shear_periodic)
202  || (mesh_bcs[BoundaryFace::inner_x1] != BoundaryFlag::shear_periodic
203  && mesh_bcs[BoundaryFace::outer_x1] == BoundaryFlag::shear_periodic))) {
204  msg << "### FATAL ERROR in Mesh constructor" << std::endl
205  << "When shear_periodic boundaries are in use, "
206  << "both sides must be shear_periodic." << std::endl;
207  ATHENA_ERROR(msg);
208  }
209 
210  // read and set MeshBlock parameters
211  block_size.x1rat = mesh_size.x1rat;
212  block_size.x2rat = mesh_size.x2rat;
213  block_size.x3rat = mesh_size.x3rat;
214  block_size.nx1 = pin->GetOrAddInteger("meshblock", "nx1", mesh_size.nx1);
215  if (f2)
216  block_size.nx2 = pin->GetOrAddInteger("meshblock", "nx2", mesh_size.nx2);
217  else
218  block_size.nx2 = mesh_size.nx2;
219  if (f3)
220  block_size.nx3 = pin->GetOrAddInteger("meshblock", "nx3", mesh_size.nx3);
221  else
222  block_size.nx3 = mesh_size.nx3;
223 
224  // check consistency of the block and mesh
225  if (mesh_size.nx1 % block_size.nx1 != 0
226  || mesh_size.nx2 % block_size.nx2 != 0
227  || mesh_size.nx3 % block_size.nx3 != 0) {
228  msg << "### FATAL ERROR in Mesh constructor" << std::endl
229  << "the Mesh must be evenly divisible by the MeshBlock" << std::endl;
230  ATHENA_ERROR(msg);
231  }
232 #if MAIN_TASKLIST != InversionTaskList
233  if (block_size.nx1 < 4 || (block_size.nx2 < 4 && f2)
234  || (block_size.nx3 < 4 && f3)) {
235  msg << "### FATAL ERROR in Mesh constructor" << std::endl
236  << "block_size must be larger than or equal to 4 cells." << std::endl;
237  ATHENA_ERROR(msg);
238  }
239 #endif
240 
241  // calculate the number of the blocks
242  nrbx1 = mesh_size.nx1/block_size.nx1;
243  nrbx2 = mesh_size.nx2/block_size.nx2;
244  nrbx3 = mesh_size.nx3/block_size.nx3;
245  nbmax = (nrbx1 > nrbx2) ? nrbx1:nrbx2;
246  nbmax = (nbmax > nrbx3) ? nbmax:nrbx3;
247 
248  // initialize user-enrollable functions
249  if (mesh_size.x1rat != 1.0) {
252  }
253  if (mesh_size.x2rat != 1.0) {
256  }
257  if (mesh_size.x3rat != 1.0) {
260  }
261 
262  // calculate the logical root level and maximum level
263  for (root_level=0; (1<<root_level) < nbmax; root_level++) {}
265 
266  tree.CreateRootGrid();
267 
268  // Load balancing flag and parameters
269 #ifdef MPI_PARALLEL
270  if (pin->GetOrAddString("loadbalancing","balancer","default") == "automatic")
271  lb_automatic_ = true;
272  else if (pin->GetOrAddString("loadbalancing","balancer","default") == "manual")
273  lb_manual_ = true;
274  lb_tolerance_ = pin->GetOrAddReal("loadbalancing","tolerance",0.5);
275  lb_interval_ = pin->GetOrAddReal("loadbalancing","interval",10);
276 #endif
277 
278  // SMR / AMR:
279  if (adaptive) {
280  max_level = pin->GetOrAddInteger("mesh", "numlevel", 1) + root_level - 1;
281  if (max_level > 63) {
282  msg << "### FATAL ERROR in Mesh constructor" << std::endl
283  << "The number of the refinement level must be smaller than "
284  << 63 - root_level + 1 << "." << std::endl;
285  ATHENA_ERROR(msg);
286  }
287  } else {
288  max_level = 63;
289  }
290 
291  if (EOS_TABLE_ENABLED) peos_table = new EosTable(pin);
292  InitUserMeshData(pin);
293 
294  if (multilevel) {
295  if (block_size.nx1 % 2 == 1 || (block_size.nx2 % 2 == 1 && f2)
296  || (block_size.nx3 % 2 == 1 && f3)) {
297  msg << "### FATAL ERROR in Mesh constructor" << std::endl
298  << "The size of MeshBlock must be divisible by 2 in order to use SMR or AMR."
299  << std::endl;
300  ATHENA_ERROR(msg);
301  }
302 
303  InputBlock *pib = pin->pfirst_block;
304  while (pib != nullptr) {
305  if (pib->block_name.compare(0, 10, "refinement") == 0) {
306  RegionSize ref_size;
307  ref_size.x1min = pin->GetReal(pib->block_name, "x1min");
308  ref_size.x1max = pin->GetReal(pib->block_name, "x1max");
309  if (f2) {
310  ref_size.x2min = pin->GetReal(pib->block_name, "x2min");
311  ref_size.x2max = pin->GetReal(pib->block_name, "x2max");
312  } else {
313  ref_size.x2min = mesh_size.x2min;
314  ref_size.x2max = mesh_size.x2max;
315  }
316  if (ndim == 3) {
317  ref_size.x3min = pin->GetReal(pib->block_name, "x3min");
318  ref_size.x3max = pin->GetReal(pib->block_name, "x3max");
319  } else {
320  ref_size.x3min = mesh_size.x3min;
321  ref_size.x3max = mesh_size.x3max;
322  }
323  int ref_lev = pin->GetInteger(pib->block_name, "level");
324  int lrlev = ref_lev + root_level;
325  if (lrlev > current_level) current_level = lrlev;
326  // range check
327  if (ref_lev < 1) {
328  msg << "### FATAL ERROR in Mesh constructor" << std::endl
329  << "Refinement level must be larger than 0 (root level = 0)" << std::endl;
330  ATHENA_ERROR(msg);
331  }
332  if (lrlev > max_level) {
333  msg << "### FATAL ERROR in Mesh constructor" << std::endl
334  << "Refinement level exceeds the maximum level (specify "
335  << "'maxlevel' parameter in <mesh> input block if adaptive)."
336  << std::endl;
337  ATHENA_ERROR(msg);
338  }
339  if (ref_size.x1min > ref_size.x1max || ref_size.x2min > ref_size.x2max
340  || ref_size.x3min > ref_size.x3max) {
341  msg << "### FATAL ERROR in Mesh constructor" << std::endl
342  << "Invalid refinement region is specified."<< std::endl;
343  ATHENA_ERROR(msg);
344  }
345  if (ref_size.x1min < mesh_size.x1min || ref_size.x1max > mesh_size.x1max
346  || ref_size.x2min < mesh_size.x2min || ref_size.x2max > mesh_size.x2max
347  || ref_size.x3min < mesh_size.x3min || ref_size.x3max > mesh_size.x3max) {
348  msg << "### FATAL ERROR in Mesh constructor" << std::endl
349  << "Refinement region must be smaller than the whole mesh." << std::endl;
350  ATHENA_ERROR(msg);
351  }
352  // find the logical range in the ref_level
353  // note: if this is too slow, this should be replaced with bi-section search.
354  std::int64_t lx1min = 0, lx1max = 0, lx2min = 0, lx2max = 0,
355  lx3min = 0, lx3max = 0;
356  std::int64_t lxmax = nrbx1*(1LL<<ref_lev);
357  for (lx1min=0; lx1min<lxmax; lx1min++) {
358  Real rx = ComputeMeshGeneratorX(lx1min+1, lxmax,
360  if (MeshGenerator_[X1DIR](rx, mesh_size) > ref_size.x1min)
361  break;
362  }
363  for (lx1max=lx1min; lx1max<lxmax; lx1max++) {
364  Real rx = ComputeMeshGeneratorX(lx1max+1, lxmax,
366  if (MeshGenerator_[X1DIR](rx, mesh_size) >= ref_size.x1max)
367  break;
368  }
369  if (lx1min % 2 == 1) lx1min--;
370  if (lx1max % 2 == 0) lx1max++;
371  if (f2) { // 2D or 3D
372  lxmax = nrbx2*(1LL << ref_lev);
373  for (lx2min=0; lx2min<lxmax; lx2min++) {
374  Real rx = ComputeMeshGeneratorX(lx2min+1, lxmax,
376  if (MeshGenerator_[X2DIR](rx, mesh_size) > ref_size.x2min)
377  break;
378  }
379  for (lx2max=lx2min; lx2max<lxmax; lx2max++) {
380  Real rx = ComputeMeshGeneratorX(lx2max+1, lxmax,
382  if (MeshGenerator_[X2DIR](rx, mesh_size) >= ref_size.x2max)
383  break;
384  }
385  if (lx2min % 2 == 1) lx2min--;
386  if (lx2max % 2 == 0) lx2max++;
387  }
388  if (ndim == 3) { // 3D
389  lxmax = nrbx3*(1LL<<ref_lev);
390  for (lx3min=0; lx3min<lxmax; lx3min++) {
391  Real rx = ComputeMeshGeneratorX(lx3min+1, lxmax,
393  if (MeshGenerator_[X3DIR](rx, mesh_size) > ref_size.x3min)
394  break;
395  }
396  for (lx3max=lx3min; lx3max<lxmax; lx3max++) {
397  Real rx = ComputeMeshGeneratorX(lx3max+1, lxmax,
399  if (MeshGenerator_[X3DIR](rx, mesh_size) >= ref_size.x3max)
400  break;
401  }
402  if (lx3min % 2 == 1) lx3min--;
403  if (lx3max % 2 == 0) lx3max++;
404  }
405  // create the finest level
406  if (ndim == 1) {
407  for (std::int64_t i=lx1min; i<lx1max; i+=2) {
408  LogicalLocation nloc;
409  nloc.level=lrlev, nloc.lx1=i, nloc.lx2=0, nloc.lx3=0;
410  int nnew;
411  tree.AddMeshBlock(nloc, nnew);
412  }
413  }
414  if (ndim == 2) {
415  for (std::int64_t j=lx2min; j<lx2max; j+=2) {
416  for (std::int64_t i=lx1min; i<lx1max; i+=2) {
417  LogicalLocation nloc;
418  nloc.level=lrlev, nloc.lx1=i, nloc.lx2=j, nloc.lx3=0;
419  int nnew;
420  tree.AddMeshBlock(nloc, nnew);
421  }
422  }
423  }
424  if (ndim == 3) {
425  for (std::int64_t k=lx3min; k<lx3max; k+=2) {
426  for (std::int64_t j=lx2min; j<lx2max; j+=2) {
427  for (std::int64_t i=lx1min; i<lx1max; i+=2) {
428  LogicalLocation nloc;
429  nloc.level = lrlev, nloc.lx1 = i, nloc.lx2 = j, nloc.lx3 = k;
430  int nnew;
431  tree.AddMeshBlock(nloc, nnew);
432  }
433  }
434  }
435  }
436  }
437  pib = pib->pnext;
438  }
439  }
440 
441  // initial mesh hierarchy construction is completed here
442  tree.CountMeshBlock(nbtotal);
444  tree.GetMeshBlockList(loclist, nullptr, nbtotal);
445 
446 #ifdef MPI_PARALLEL
447  // check if there are sufficient blocks
448  if (nbtotal < Globals::nranks) {
449  if (mesh_test == 0) {
450  msg << "### FATAL ERROR in Mesh constructor" << std::endl
451  << "Too few mesh blocks: nbtotal ("<< nbtotal <<") < nranks ("
452  << Globals::nranks << ")" << std::endl;
453  ATHENA_ERROR(msg);
454  } else { // test
455  std::cout << "### Warning in Mesh constructor" << std::endl
456  << "Too few mesh blocks: nbtotal ("<< nbtotal <<") < nranks ("
457  << Globals::nranks << ")" << std::endl;
458  }
459  }
460 #endif
461 
462  ranklist = new int[nbtotal];
463  nslist = new int[Globals::nranks];
464  nblist = new int[Globals::nranks];
465  costlist = new double[nbtotal];
466  if (adaptive) { // allocate arrays for AMR
467  nref = new int[Globals::nranks];
468  nderef = new int[Globals::nranks];
469  rdisp = new int[Globals::nranks];
470  ddisp = new int[Globals::nranks];
471  bnref = new int[Globals::nranks];
472  bnderef = new int[Globals::nranks];
473  brdisp = new int[Globals::nranks];
474  bddisp = new int[Globals::nranks];
475  }
476 
477  // initialize cost array with the simplest estimate; all the blocks are equal
478  for (int i=0; i<nbtotal; i++) costlist[i] = 1.0;
479 
481 
482  // Output some diagnostic information to terminal
483 
484  // Output MeshBlock list and quit (mesh test only); do not create meshes
485  if (mesh_test > 0) {
487  return;
488  }
489 
490 
491  if (SELF_GRAVITY_ENABLED == 1) {
492  gflag = 1; // set gravity flag
493  pfgrd = new FFTGravityDriver(this, pin);
494  } else if (SELF_GRAVITY_ENABLED == 2) {
495  // MGDriver must be initialzied before MeshBlocks
496  pmgrd = new MGGravityDriver(this, pin);
497  }
498  // if (SELF_GRAVITY_ENABLED == 2 && ...) // independent allocation
499  // gflag = 2;
500 
501  // create MeshBlock list for this process
502  int nbs = nslist[Globals::my_rank];
503  int nbe = nbs + nblist[Globals::my_rank] - 1;
504  // create MeshBlock list for this process
505  for (int i=nbs; i<=nbe; i++) {
506  SetBlockSizeAndBoundaries(loclist[i], block_size, block_bcs);
507  // create a block and add into the link list
508  if (i == nbs) {
509  pblock = new MeshBlock(i, i-nbs, loclist[i], block_size, block_bcs, this,
510  pin, gflag);
511  pfirst = pblock;
512  } else {
513  pblock->next = new MeshBlock(i, i-nbs, loclist[i], block_size, block_bcs,
514  this, pin, gflag);
515  pblock->next->prev = pblock;
516  pblock = pblock->next;
517  }
518  pblock->pbval->SearchAndSetNeighbors(tree, ranklist, nslist);
519  }
520  pblock = pfirst;
521 
523 
524  if (turb_flag > 0) // TurbulenceDriver depends on the MeshBlock ctor
525  ptrbd = new TurbulenceDriver(this, pin);
526 }
527 
528 //----------------------------------------------------------------------------------------
529 // Mesh constructor for restarts. Load the restart file
530 
531 Mesh::Mesh(ParameterInput *pin, IOWrapper& resfile, int mesh_test) :
532  // public members:
533  // aggregate initialization of RegionSize struct:
534  // (will be overwritten by memcpy from restart file, in this case)
535  mesh_size{pin->GetReal("mesh", "x1min"), pin->GetReal("mesh", "x2min"),
536  pin->GetReal("mesh", "x3min"), pin->GetReal("mesh", "x1max"),
537  pin->GetReal("mesh", "x2max"), pin->GetReal("mesh", "x3max"),
538  pin->GetOrAddReal("mesh", "x1rat", 1.0),
539  pin->GetOrAddReal("mesh", "x2rat", 1.0),
540  pin->GetOrAddReal("mesh", "x3rat", 1.0),
541  pin->GetInteger("mesh", "nx1"), pin->GetInteger("mesh", "nx2"),
542  pin->GetInteger("mesh", "nx3") },
543  mesh_bcs{GetBoundaryFlag(pin->GetOrAddString("mesh", "ix1_bc", "none")),
544  GetBoundaryFlag(pin->GetOrAddString("mesh", "ox1_bc", "none")),
545  GetBoundaryFlag(pin->GetOrAddString("mesh", "ix2_bc", "none")),
546  GetBoundaryFlag(pin->GetOrAddString("mesh", "ox2_bc", "none")),
547  GetBoundaryFlag(pin->GetOrAddString("mesh", "ix3_bc", "none")),
548  GetBoundaryFlag(pin->GetOrAddString("mesh", "ox3_bc", "none"))},
549  f2(mesh_size.nx2 > 1 ? true : false), f3(mesh_size.nx3 > 1 ? true : false),
550  ndim(f3 ? 3 : (f2 ? 2 : 1)),
551  adaptive(pin->GetOrAddString("mesh", "refinement", "none") == "adaptive"
552  ? true : false),
553  multilevel((adaptive || pin->GetOrAddString("mesh", "refinement", "none") == "static")
554  ? true : false),
555  fluid_setup(GetFluidFormulation(pin->GetOrAddString("hydro", "active", "true"))),
556  start_time(pin->GetOrAddReal("time", "start_time", 0.0)), time(start_time),
557  tlim(pin->GetReal("time", "tlim")), dt(std::numeric_limits<Real>::max()),
558  dt_hyperbolic(dt), dt_parabolic(dt), dt_user(dt),
559  cfl_number(pin->GetReal("time", "cfl_number")),
560  nlim(pin->GetOrAddInteger("time", "nlim", -1)), ncycle(),
561  ncycle_out(pin->GetOrAddInteger("time", "ncycle_out", 1)),
562  dt_diagnostics(pin->GetOrAddInteger("time", "dt_diagnostics", -1)),
563  muj(), nuj(), muj_tilde(),
564  nbnew(), nbdel(),
565  step_since_lb(), gflag(), turb_flag(),
566  // private members:
567  next_phys_id_(), num_mesh_threads_(pin->GetOrAddInteger("mesh", "num_threads", 1)),
568  tree(this),
569  use_uniform_meshgen_fn_{true, true, true},
570  nreal_user_mesh_data_(), nint_user_mesh_data_(), nuser_history_output_(),
571  four_pi_G_(), grav_eps_(-1.0), grav_mean_rho_(-1.0),
572  lb_flag_(true), lb_automatic_(), lb_manual_(),
575  BoundaryFunction_{nullptr, nullptr, nullptr, nullptr, nullptr, nullptr},
576  AMRFlag_{}, UserSourceTerm_{}, UserTimeStep_{}, ViscosityCoeff_{},
577  ConductionCoeff_{}, FieldDiffusivity_{},
578  MGGravityBoundaryFunction_{MGPeriodicInnerX1, MGPeriodicOuterX1, MGPeriodicInnerX2,
579  MGPeriodicOuterX2, MGPeriodicInnerX3, MGPeriodicOuterX3} {
580  std::stringstream msg;
581  RegionSize block_size;
582  BoundaryFlag block_bcs[6];
583  MeshBlock *pfirst{};
584  IOWrapperSizeT *offset{};
585  IOWrapperSizeT datasize, listsize, headeroffset;
586 
587  // mesh test
588  if (mesh_test > 0) Globals::nranks = mesh_test;
589 
590 #ifdef MPI_PARALLEL
591  // reserve phys=0 for former TAG_AMR=8; now hard-coded in Mesh::CreateAMRMPITag()
592  next_phys_id_ = 1;
594 #endif
595 
596  // check the number of OpenMP threads for mesh
597  if (num_mesh_threads_ < 1) {
598  msg << "### FATAL ERROR in Mesh constructor" << std::endl
599  << "Number of OpenMP threads must be >= 1, but num_threads="
600  << num_mesh_threads_ << std::endl;
601  ATHENA_ERROR(msg);
602  }
603 
604  // get the end of the header
605  headeroffset = resfile.GetPosition();
606  // read the restart file
607  // the file is already open and the pointer is set to after <par_end>
608  IOWrapperSizeT headersize = sizeof(int)*3+sizeof(Real)*2
609  + sizeof(RegionSize)+sizeof(IOWrapperSizeT);
610  char *headerdata = new char[headersize];
611  if (Globals::my_rank == 0) { // the master process reads the header data
612  if (resfile.Read(headerdata, 1, headersize) != headersize) {
613  msg << "### FATAL ERROR in Mesh constructor" << std::endl
614  << "The restart file is broken." << std::endl;
615  ATHENA_ERROR(msg);
616  }
617  }
618 #ifdef MPI_PARALLEL
619  // then broadcast the header data
620  MPI_Bcast(headerdata, headersize, MPI_BYTE, 0, MPI_COMM_WORLD);
621 #endif
622  IOWrapperSizeT hdos = 0;
623  std::memcpy(&nbtotal, &(headerdata[hdos]), sizeof(int));
624  hdos += sizeof(int);
625  std::memcpy(&root_level, &(headerdata[hdos]), sizeof(int));
626  hdos += sizeof(int);
628  std::memcpy(&mesh_size, &(headerdata[hdos]), sizeof(RegionSize));
629  hdos += sizeof(RegionSize);
630  std::memcpy(&time, &(headerdata[hdos]), sizeof(Real));
631  hdos += sizeof(Real);
632  std::memcpy(&dt, &(headerdata[hdos]), sizeof(Real));
633  hdos += sizeof(Real);
634  std::memcpy(&ncycle, &(headerdata[hdos]), sizeof(int));
635  hdos += sizeof(int);
636  std::memcpy(&datasize, &(headerdata[hdos]), sizeof(IOWrapperSizeT));
637  hdos += sizeof(IOWrapperSizeT); // (this updated value is never used)
638 
639  delete [] headerdata;
640 
641  // initialize
643  offset = new IOWrapperSizeT[nbtotal];
644  costlist = new double[nbtotal];
645  ranklist = new int[nbtotal];
646  nslist = new int[Globals::nranks];
647  nblist = new int[Globals::nranks];
648 
649  block_size.nx1 = pin->GetOrAddInteger("meshblock", "nx1", mesh_size.nx1);
650  block_size.nx2 = pin->GetOrAddInteger("meshblock", "nx2", mesh_size.nx2);
651  block_size.nx3 = pin->GetOrAddInteger("meshblock", "nx3", mesh_size.nx3);
652 
653  // calculate the number of the blocks
654  nrbx1 = mesh_size.nx1/block_size.nx1;
655  nrbx2 = mesh_size.nx2/block_size.nx2;
656  nrbx3 = mesh_size.nx3/block_size.nx3;
657 
658  // initialize user-enrollable functions
659  if (mesh_size.x1rat != 1.0) {
662  }
663  if (mesh_size.x2rat != 1.0) {
666  }
667  if (mesh_size.x3rat != 1.0) {
670  }
671 
672  // Load balancing flag and parameters
673 #ifdef MPI_PARALLEL
674  if (pin->GetOrAddString("loadbalancing", "balancer", "default") == "automatic")
675  lb_automatic_ = true;
676  else if (pin->GetOrAddString("loadbalancing", "balancer", "default") == "manual")
677  lb_manual_ = true;
678  lb_tolerance_ = pin->GetOrAddReal("loadbalancing", "tolerance", 0.5);
679  lb_interval_ = pin->GetOrAddReal("loadbalancing", "interval", 10);
680 #endif
681 
682  // SMR / AMR
683  if (adaptive) {
684  max_level = pin->GetOrAddInteger("mesh", "numlevel", 1) + root_level - 1;
685  if (max_level > 63) {
686  msg << "### FATAL ERROR in Mesh constructor" << std::endl
687  << "The number of the refinement level must be smaller than "
688  << 63 - root_level + 1 << "." << std::endl;
689  ATHENA_ERROR(msg);
690  }
691  } else {
692  max_level = 63;
693  }
694 
695  if (EOS_TABLE_ENABLED) peos_table = new EosTable(pin);
696  InitUserMeshData(pin);
697 
698  // read user Mesh data
699  IOWrapperSizeT udsize = 0;
700  for (int n=0; n<nint_user_mesh_data_; n++)
701  udsize += iuser_mesh_data[n].GetSizeInBytes();
702  for (int n=0; n<nreal_user_mesh_data_; n++)
703  udsize += ruser_mesh_data[n].GetSizeInBytes();
704  if (udsize != 0) {
705  char *userdata = new char[udsize];
706  if (Globals::my_rank == 0) { // only the master process reads the ID list
707  if (resfile.Read(userdata, 1, udsize) != udsize) {
708  msg << "### FATAL ERROR in Mesh constructor" << std::endl
709  << "The restart file is broken." << std::endl;
710  ATHENA_ERROR(msg);
711  }
712  }
713 #ifdef MPI_PARALLEL
714  // then broadcast the ID list
715  MPI_Bcast(userdata, udsize, MPI_BYTE, 0, MPI_COMM_WORLD);
716 #endif
717 
718  IOWrapperSizeT udoffset=0;
719  for (int n=0; n<nint_user_mesh_data_; n++) {
720  std::memcpy(iuser_mesh_data[n].data(), &(userdata[udoffset]),
721  iuser_mesh_data[n].GetSizeInBytes());
722  udoffset += iuser_mesh_data[n].GetSizeInBytes();
723  }
724  for (int n=0; n<nreal_user_mesh_data_; n++) {
725  std::memcpy(ruser_mesh_data[n].data(), &(userdata[udoffset]),
726  ruser_mesh_data[n].GetSizeInBytes());
727  udoffset += ruser_mesh_data[n].GetSizeInBytes();
728  }
729  delete [] userdata;
730  }
731 
732  // read the ID list
733  listsize = sizeof(LogicalLocation)+sizeof(Real);
734  //allocate the idlist buffer
735  char *idlist = new char[listsize*nbtotal];
736  if (Globals::my_rank == 0) { // only the master process reads the ID list
737  if (resfile.Read(idlist, listsize, nbtotal) != static_cast<unsigned int>(nbtotal)) {
738  msg << "### FATAL ERROR in Mesh constructor" << std::endl
739  << "The restart file is broken." << std::endl;
740  ATHENA_ERROR(msg);
741  }
742  }
743 #ifdef MPI_PARALLEL
744  // then broadcast the ID list
745  MPI_Bcast(idlist, listsize*nbtotal, MPI_BYTE, 0, MPI_COMM_WORLD);
746 #endif
747 
748  int os = 0;
749  for (int i=0; i<nbtotal; i++) {
750  std::memcpy(&(loclist[i]), &(idlist[os]), sizeof(LogicalLocation));
751  os += sizeof(LogicalLocation);
752  std::memcpy(&(costlist[i]), &(idlist[os]), sizeof(double));
753  os += sizeof(double);
754  if (loclist[i].level > current_level) current_level = loclist[i].level;
755  }
756  delete [] idlist;
757 
758  // calculate the header offset and seek
759  headeroffset += headersize + udsize + listsize*nbtotal;
760  if (Globals::my_rank != 0)
761  resfile.Seek(headeroffset);
762 
763  // rebuild the Block Tree
764  tree.CreateRootGrid();
765  for (int i=0; i<nbtotal; i++)
766  tree.AddMeshBlockWithoutRefine(loclist[i]);
767  int nnb;
768  // check the tree structure, and assign GID
769  tree.GetMeshBlockList(loclist, nullptr, nnb);
770  if (nnb != nbtotal) {
771  msg << "### FATAL ERROR in Mesh constructor" << std::endl
772  << "Tree reconstruction failed. The total numbers of the blocks do not match. ("
773  << nbtotal << " != " << nnb << ")" << std::endl;
774  ATHENA_ERROR(msg);
775  }
776 
777 #ifdef MPI_PARALLEL
778  if (nbtotal < Globals::nranks) {
779  if (mesh_test == 0) {
780  msg << "### FATAL ERROR in Mesh constructor" << std::endl
781  << "Too few mesh blocks: nbtotal ("<< nbtotal <<") < nranks ("
782  << Globals::nranks << ")" << std::endl;
783  ATHENA_ERROR(msg);
784  } else { // test
785  std::cout << "### Warning in Mesh constructor" << std::endl
786  << "Too few mesh blocks: nbtotal ("<< nbtotal <<") < nranks ("
787  << Globals::nranks << ")" << std::endl;
788  delete [] offset;
789  return;
790  }
791  }
792 #endif
793 
794  if (adaptive) { // allocate arrays for AMR
795  nref = new int[Globals::nranks];
796  nderef = new int[Globals::nranks];
797  rdisp = new int[Globals::nranks];
798  ddisp = new int[Globals::nranks];
799  bnref = new int[Globals::nranks];
800  bnderef = new int[Globals::nranks];
801  brdisp = new int[Globals::nranks];
802  bddisp = new int[Globals::nranks];
803  }
804 
806 
807  // Output MeshBlock list and quit (mesh test only); do not create meshes
808  if (mesh_test > 0) {
810  delete [] offset;
811  return;
812  }
813 
814  if (SELF_GRAVITY_ENABLED == 1) {
815  gflag = 1; // set gravity flag
816  pfgrd = new FFTGravityDriver(this, pin);
817  } else if (SELF_GRAVITY_ENABLED == 2) {
818  // MGDriver must be initialzied before MeshBlocks
819  pmgrd = new MGGravityDriver(this, pin);
820  }
821  // if (SELF_GRAVITY_ENABLED == 2 && ...) // independent allocation
822  // gflag=2;
823 
824  // allocate data buffer
825  int nb = nblist[Globals::my_rank];
826  int nbs = nslist[Globals::my_rank];
827  int nbe = nbs + nb - 1;
828  char *mbdata = new char[datasize*nb];
829  // load MeshBlocks (parallel)
830  if (resfile.Read_at_all(mbdata, datasize, nb, headeroffset+nbs*datasize) !=
831  static_cast<unsigned int>(nb)) {
832  msg << "### FATAL ERROR in Mesh constructor" << std::endl
833  << "The restart file is broken or input parameters are inconsistent."
834  << std::endl;
835  ATHENA_ERROR(msg);
836  }
837  for (int i=nbs; i<=nbe; i++) {
838  // Match fixed-width integer precision of IOWrapperSizeT datasize
839  std::uint64_t buff_os = datasize * (i-nbs);
840  SetBlockSizeAndBoundaries(loclist[i], block_size, block_bcs);
841  // create a block and add into the link list
842  if (i == nbs) {
843  pblock = new MeshBlock(i, i-nbs, this, pin, loclist[i], block_size,
844  block_bcs, costlist[i], mbdata+buff_os, gflag);
845  pfirst = pblock;
846  } else {
847  pblock->next = new MeshBlock(i, i-nbs, this, pin, loclist[i], block_size,
848  block_bcs, costlist[i], mbdata+buff_os, gflag);
849  pblock->next->prev = pblock;
850  pblock = pblock->next;
851  }
852  pblock->pbval->SearchAndSetNeighbors(tree, ranklist, nslist);
853  }
854  pblock = pfirst;
855  delete [] mbdata;
856  // check consistency
857  if (datasize != pblock->GetBlockSizeInBytes()) {
858  msg << "### FATAL ERROR in Mesh constructor" << std::endl
859  << "The restart file is broken or input parameters are inconsistent."
860  << std::endl;
861  ATHENA_ERROR(msg);
862  }
863 
865 
866  // clean up
867  delete [] offset;
868 
869  if (turb_flag > 0) // TurbulenceDriver depends on the MeshBlock ctor
870  ptrbd = new TurbulenceDriver(this, pin);
871 }
872 
873 //----------------------------------------------------------------------------------------
874 // destructor
875 
877  while (pblock->prev != nullptr) // should not be true
878  delete pblock->prev;
879  while (pblock->next != nullptr)
880  delete pblock->next;
881  delete pblock;
882  delete [] nslist;
883  delete [] nblist;
884  delete [] ranklist;
885  delete [] costlist;
886  delete [] loclist;
887  if (SELF_GRAVITY_ENABLED == 1) delete pfgrd;
888  else if (SELF_GRAVITY_ENABLED == 2) delete pmgrd;
889  if (turb_flag > 0) delete ptrbd;
890  if (adaptive) { // deallocate arrays for AMR
891  delete [] nref;
892  delete [] nderef;
893  delete [] rdisp;
894  delete [] ddisp;
895  delete [] bnref;
896  delete [] bnderef;
897  delete [] brdisp;
898  delete [] bddisp;
899  }
900  // delete user Mesh data
901  if (nreal_user_mesh_data_>0) delete [] ruser_mesh_data;
902  if (nuser_history_output_ > 0) {
903  delete [] user_history_output_names_;
904  delete [] user_history_func_;
905  delete [] user_history_ops_;
906  }
907  if (nint_user_mesh_data_>0) delete [] iuser_mesh_data;
908  if (EOS_TABLE_ENABLED) delete peos_table;
909 }
910 
911 //----------------------------------------------------------------------------------------
913 // \brief print the mesh structure information
914 
916  RegionSize block_size;
917  BoundaryFlag block_bcs[6];
918  FILE *fp = nullptr;
919 
920  // open 'mesh_structure.dat' file
921  if (f2) {
922  if ((fp = std::fopen("mesh_structure.dat","wb")) == nullptr) {
923  std::cout << "### ERROR in function Mesh::OutputMeshStructure" << std::endl
924  << "Cannot open mesh_structure.dat" << std::endl;
925  return;
926  }
927  }
928 
929  // Write overall Mesh structure to stdout and file
930  std::cout << std::endl;
931  std::cout << "Root grid = " << nrbx1 << " x " << nrbx2 << " x " << nrbx3
932  << " MeshBlocks" << std::endl;
933  std::cout << "Total number of MeshBlocks = " << nbtotal << std::endl;
934  std::cout << "Number of physical refinement levels = "
935  << (current_level - root_level) << std::endl;
936  std::cout << "Number of logical refinement levels = " << current_level << std::endl;
937 
938  // compute/output number of blocks per level, and cost per level
939  int *nb_per_plevel = new int[max_level];
940  int *cost_per_plevel = new int[max_level];
941  for (int i=0; i<=max_level; ++i) {
942  nb_per_plevel[i] = 0;
943  cost_per_plevel[i] = 0;
944  }
945  for (int i=0; i<nbtotal; i++) {
946  nb_per_plevel[(loclist[i].level - root_level)]++;
947  cost_per_plevel[(loclist[i].level - root_level)] += costlist[i];
948  }
949  for (int i=root_level; i<=max_level; i++) {
950  if (nb_per_plevel[i-root_level] != 0) {
951  std::cout << " Physical level = " << i-root_level << " (logical level = " << i
952  << "): " << nb_per_plevel[i-root_level] << " MeshBlocks, cost = "
953  << cost_per_plevel[i-root_level] << std::endl;
954  }
955  }
956 
957  // compute/output number of blocks per rank, and cost per rank
958  std::cout << "Number of parallel ranks = " << Globals::nranks << std::endl;
959  int *nb_per_rank = new int[Globals::nranks];
960  int *cost_per_rank = new int[Globals::nranks];
961  for (int i=0; i<Globals::nranks; ++i) {
962  nb_per_rank[i] = 0;
963  cost_per_rank[i] = 0;
964  }
965  for (int i=0; i<nbtotal; i++) {
966  nb_per_rank[ranklist[i]]++;
967  cost_per_rank[ranklist[i]] += costlist[i];
968  }
969  for (int i=0; i<Globals::nranks; ++i) {
970  std::cout << " Rank = " << i << ": " << nb_per_rank[i] <<" MeshBlocks, cost = "
971  << cost_per_rank[i] << std::endl;
972  }
973 
974  // output relative size/locations of meshblock to file, for plotting
975  double real_max = std::numeric_limits<double>::max();
976  double mincost = real_max, maxcost = 0.0, totalcost = 0.0;
977  for (int i=root_level; i<=max_level; i++) {
978  for (int j=0; j<nbtotal; j++) {
979  if (loclist[j].level == i) {
980  SetBlockSizeAndBoundaries(loclist[j], block_size, block_bcs);
981  std::int64_t &lx1 = loclist[j].lx1;
982  std::int64_t &lx2 = loclist[j].lx2;
983  std::int64_t &lx3 = loclist[j].lx3;
984  int &ll = loclist[j].level;
985  mincost = std::min(mincost,costlist[i]);
986  maxcost = std::max(maxcost,costlist[i]);
987  totalcost += costlist[i];
988  std::fprintf(fp,"#MeshBlock %d on rank=%d with cost=%g\n", j, ranklist[j],
989  costlist[j]);
990  std::fprintf(
991  fp, "# Logical level %d, location = (%" PRId64 " %" PRId64 " %" PRId64")\n",
992  ll, lx1, lx2, lx3);
993  if (ndim == 2) {
994  std::fprintf(fp, "%g %g\n", block_size.x1min, block_size.x2min);
995  std::fprintf(fp, "%g %g\n", block_size.x1max, block_size.x2min);
996  std::fprintf(fp, "%g %g\n", block_size.x1max, block_size.x2max);
997  std::fprintf(fp, "%g %g\n", block_size.x1min, block_size.x2max);
998  std::fprintf(fp, "%g %g\n", block_size.x1min, block_size.x2min);
999  std::fprintf(fp, "\n\n");
1000  }
1001  if (ndim == 3) {
1002  std::fprintf(fp, "%g %g %g\n", block_size.x1min, block_size.x2min,
1003  block_size.x3min);
1004  std::fprintf(fp, "%g %g %g\n", block_size.x1max, block_size.x2min,
1005  block_size.x3min);
1006  std::fprintf(fp, "%g %g %g\n", block_size.x1max, block_size.x2max,
1007  block_size.x3min);
1008  std::fprintf(fp, "%g %g %g\n", block_size.x1min, block_size.x2max,
1009  block_size.x3min);
1010  std::fprintf(fp, "%g %g %g\n", block_size.x1min, block_size.x2min,
1011  block_size.x3min);
1012  std::fprintf(fp, "%g %g %g\n", block_size.x1min, block_size.x2min,
1013  block_size.x3max);
1014  std::fprintf(fp, "%g %g %g\n", block_size.x1max, block_size.x2min,
1015  block_size.x3max);
1016  std::fprintf(fp, "%g %g %g\n", block_size.x1max, block_size.x2min,
1017  block_size.x3min);
1018  std::fprintf(fp, "%g %g %g\n", block_size.x1max, block_size.x2min,
1019  block_size.x3max);
1020  std::fprintf(fp, "%g %g %g\n", block_size.x1max, block_size.x2max,
1021  block_size.x3max);
1022  std::fprintf(fp, "%g %g %g\n", block_size.x1max, block_size.x2max,
1023  block_size.x3min);
1024  std::fprintf(fp, "%g %g %g\n", block_size.x1max, block_size.x2max,
1025  block_size.x3max);
1026  std::fprintf(fp, "%g %g %g\n", block_size.x1min, block_size.x2max,
1027  block_size.x3max);
1028  std::fprintf(fp, "%g %g %g\n", block_size.x1min, block_size.x2max,
1029  block_size.x3min);
1030  std::fprintf(fp, "%g %g %g\n", block_size.x1min, block_size.x2max,
1031  block_size.x3max);
1032  std::fprintf(fp, "%g %g %g\n", block_size.x1min, block_size.x2min,
1033  block_size.x3max);
1034  std::fprintf(fp, "%g %g %g\n", block_size.x1min, block_size.x2min,
1035  block_size.x3min);
1036  std::fprintf(fp, "\n\n");
1037  }
1038  }
1039  }
1040  }
1041 
1042  // close file, final outputs
1043  if (f2) std::fclose(fp);
1044  std::cout << "Load Balancing:" << std::endl;
1045  std::cout << " Minimum cost = " << mincost << ", Maximum cost = " << maxcost
1046  << ", Average cost = " << totalcost/nbtotal << std::endl << std::endl;
1047  std::cout << "See the 'mesh_structure.dat' file for a complete list"
1048  << " of MeshBlocks." << std::endl;
1049  std::cout << "Use 'python ../vis/python/plot_mesh.py' or gnuplot"
1050  << " to visualize mesh structure." << std::endl << std::endl;
1051 
1052  delete [] nb_per_plevel;
1053  delete [] cost_per_plevel;
1054  delete [] nb_per_rank;
1055  delete [] cost_per_rank;
1056 
1057  return;
1058 }
1059 
1060 //----------------------------------------------------------------------------------------
1061 // \!fn void Mesh::NewTimeStep()
1062 // \brief function that loops over all MeshBlocks and find new timestep
1063 // this assumes that phydro->NewBlockTimeStep is already called
1064 
1066  MeshBlock *pmb = pblock;
1067 
1068  // prevent timestep from growing too fast in between 2x cycles (even if every MeshBlock
1069  // has new_block_dt > 2.0*dt_old)
1070  dt = static_cast<Real>(2.0)*dt;
1071  // consider first MeshBlock on this MPI rank's linked list of blocks:
1072  dt = std::min(dt, pmb->new_block_dt_);
1075  dt_user = pmb->new_block_dt_user_;
1076  pmb = pmb->next;
1077 
1078  while (pmb != nullptr) {
1079  dt = std::min(dt, pmb->new_block_dt_);
1083  pmb = pmb->next;
1084  }
1085 
1086 #ifdef MPI_PARALLEL
1087  // pack array, MPI allreduce over array, then unpack into Mesh variables
1088  Real dt_array[4] = {dt, dt_hyperbolic, dt_parabolic, dt_user};
1089  MPI_Allreduce(MPI_IN_PLACE, dt_array, 4, MPI_ATHENA_REAL, MPI_MIN, MPI_COMM_WORLD);
1090  dt = dt_array[0];
1091  dt_hyperbolic = dt_array[1];
1092  dt_parabolic = dt_array[2];
1093  dt_user = dt_array[3];
1094 #endif
1095 
1096  if (time < tlim && (tlim - time) < dt) // timestep would take us past desired endpoint
1097  dt = tlim - time;
1098 
1099  return;
1100 }
1101 
1102 //----------------------------------------------------------------------------------------
1104 // \brief Enroll a user-defined boundary function
1105 
1106 void Mesh::EnrollUserBoundaryFunction(BoundaryFace dir, BValFunc my_bc) {
1107  std::stringstream msg;
1108  if (dir < 0 || dir > 5) {
1109  msg << "### FATAL ERROR in EnrollBoundaryCondition function" << std::endl
1110  << "dirName = " << dir << " not valid" << std::endl;
1111  ATHENA_ERROR(msg);
1112  }
1113  if (mesh_bcs[dir] != BoundaryFlag::user) {
1114  msg << "### FATAL ERROR in EnrollUserBoundaryFunction" << std::endl
1115  << "The boundary condition flag must be set to the string 'user' in the "
1116  << " <mesh> block in the input file to use user-enrolled BCs" << std::endl;
1117  ATHENA_ERROR(msg);
1118  }
1119  BoundaryFunction_[static_cast<int>(dir)]=my_bc;
1120  return;
1121 }
1122 
1123 //----------------------------------------------------------------------------------------
1125 // MGBoundaryFunc my_bc)
1126 // \brief Enroll a user-defined Multigrid boundary function
1127 
1129  std::stringstream msg;
1130  if (dir < 0 || dir > 5) {
1131  msg << "### FATAL ERROR in EnrollBoundaryCondition function" << std::endl
1132  << "dirName = " << dir << " not valid" << std::endl;
1133  ATHENA_ERROR(msg);
1134  }
1135  MGGravityBoundaryFunction_[static_cast<int>(dir)] = my_bc;
1136  return;
1137 }
1138 
1139 // DEPRECATED(felker): provide trivial overloads for old-style BoundaryFace enum argument
1141  EnrollUserBoundaryFunction(static_cast<BoundaryFace>(dir), my_bc);
1142  return;
1143 }
1144 
1146  EnrollUserMGGravityBoundaryFunction(static_cast<BoundaryFace>(dir), my_bc);
1147  return;
1148 }
1149 
1150 //----------------------------------------------------------------------------------------
1152 // \brief Enroll a user-defined function for checking refinement criteria
1153 
1155  if (adaptive)
1156  AMRFlag_ = amrflag;
1157  return;
1158 }
1159 
1160 //----------------------------------------------------------------------------------------
1162 // \brief Enroll a user-defined function for Mesh generation
1163 
1165  std::stringstream msg;
1166  if (dir < 0 || dir >= 3) {
1167  msg << "### FATAL ERROR in EnrollUserMeshGenerator function" << std::endl
1168  << "dirName = " << dir << " not valid" << std::endl;
1169  ATHENA_ERROR(msg);
1170  }
1171  if (dir == X1DIR && mesh_size.x1rat > 0.0) {
1172  msg << "### FATAL ERROR in EnrollUserMeshGenerator function" << std::endl
1173  << "x1rat = " << mesh_size.x1rat <<
1174  " must be negative for user-defined mesh generator in X1DIR " << std::endl;
1175  ATHENA_ERROR(msg);
1176  }
1177  if (dir == X2DIR && mesh_size.x2rat > 0.0) {
1178  msg << "### FATAL ERROR in EnrollUserMeshGenerator function" << std::endl
1179  << "x2rat = " << mesh_size.x2rat <<
1180  " must be negative for user-defined mesh generator in X2DIR " << std::endl;
1181  ATHENA_ERROR(msg);
1182  }
1183  if (dir == X3DIR && mesh_size.x3rat > 0.0) {
1184  msg << "### FATAL ERROR in EnrollUserMeshGenerator function" << std::endl
1185  << "x3rat = " << mesh_size.x3rat <<
1186  " must be negative for user-defined mesh generator in X3DIR " << std::endl;
1187  ATHENA_ERROR(msg);
1188  }
1189  use_uniform_meshgen_fn_[dir] = false;
1190  MeshGenerator_[dir] = my_mg;
1191  return;
1192 }
1193 
1194 //----------------------------------------------------------------------------------------
1196 // \brief Enroll a user-defined source function
1197 
1199  UserSourceTerm_ = my_func;
1200  return;
1201 }
1202 
1203 //----------------------------------------------------------------------------------------
1205 // \brief Enroll a user-defined time step function
1206 
1208  UserTimeStep_ = my_func;
1209  return;
1210 }
1211 
1212 //----------------------------------------------------------------------------------------
1214 // \brief set the number of user-defined history outputs
1215 
1218  user_history_output_names_ = new std::string[n];
1221  for (int i=0; i<n; i++) user_history_func_[i] = nullptr;
1222 }
1223 
1224 //----------------------------------------------------------------------------------------
1226 // const char *name, UserHistoryOperation op)
1227 // \brief Enroll a user-defined history output function and set its name
1228 
1229 void Mesh::EnrollUserHistoryOutput(int i, HistoryOutputFunc my_func, const char *name,
1230  UserHistoryOperation op) {
1231  std::stringstream msg;
1232  if (i >= nuser_history_output_) {
1233  msg << "### FATAL ERROR in EnrollUserHistoryOutput function" << std::endl
1234  << "The number of the user-defined history output (" << i << ") "
1235  << "exceeds the declared number (" << nuser_history_output_ << ")." << std::endl;
1236  ATHENA_ERROR(msg);
1237  }
1238  user_history_output_names_[i] = name;
1239  user_history_func_[i] = my_func;
1240  user_history_ops_[i] = op;
1241 }
1242 
1243 //----------------------------------------------------------------------------------------
1245 // \brief Enroll a user-defined metric for arbitrary GR coordinates
1246 
1248  UserMetric_ = my_func;
1249  return;
1250 }
1251 
1252 //----------------------------------------------------------------------------------------
1254 // \brief Enroll a user-defined magnetic field diffusivity function
1255 
1257  ViscosityCoeff_ = my_func;
1258  return;
1259 }
1260 
1261 //----------------------------------------------------------------------------------------
1263 // \brief Enroll a user-defined thermal conduction function
1264 
1266  ConductionCoeff_ = my_func;
1267  return;
1268 }
1269 
1270 //----------------------------------------------------------------------------------------
1272 // \brief Enroll a user-defined magnetic field diffusivity function
1273 
1275  FieldDiffusivity_ = my_func;
1276  return;
1277 }
1278 //----------------------------------------------------------------------------------------
1280 // \brief Allocate Real AthenaArrays for user-defned data in Mesh
1281 
1283  if (nreal_user_mesh_data_ != 0) {
1284  std::stringstream msg;
1285  msg << "### FATAL ERROR in Mesh::AllocateRealUserMeshDataField"
1286  << std::endl << "User Mesh data arrays are already allocated" << std::endl;
1287  ATHENA_ERROR(msg);
1288  }
1291  return;
1292 }
1293 
1294 //----------------------------------------------------------------------------------------
1296 // \brief Allocate integer AthenaArrays for user-defned data in Mesh
1297 
1299  if (nint_user_mesh_data_ != 0) {
1300  std::stringstream msg;
1301  msg << "### FATAL ERROR in Mesh::AllocateIntUserMeshDataField"
1302  << std::endl << "User Mesh data arrays are already allocated" << std::endl;
1303  ATHENA_ERROR(msg);
1304  }
1307  return;
1308 }
1309 
1310 
1311 //----------------------------------------------------------------------------------------
1312 // \!fn void Mesh::ApplyUserWorkBeforeOutput(ParameterInput *pin)
1313 // \brief Apply MeshBlock::UserWorkBeforeOutput
1314 
1316  MeshBlock *pmb = pblock;
1317  while (pmb != nullptr) {
1318  pmb->UserWorkBeforeOutput(pin);
1319  pmb = pmb->next;
1320  }
1321 }
1322 
1323 //----------------------------------------------------------------------------------------
1324 // \!fn void Mesh::Initialize(int res_flag, ParameterInput *pin)
1325 // \brief initialization before the main loop
1326 
1327 void Mesh::Initialize(int res_flag, ParameterInput *pin) {
1328  bool iflag = true;
1329  int inb = nbtotal;
1330  int nthreads = GetNumMeshThreads();
1332  std::vector<MeshBlock*> pmb_array(nmb);
1333 
1334  do {
1335  // initialize a vector of MeshBlock pointers
1337  if (static_cast<unsigned int>(nmb) != pmb_array.size()) pmb_array.resize(nmb);
1338  MeshBlock *pmbl = pblock;
1339  for (int i=0; i<nmb; ++i) {
1340  pmb_array[i] = pmbl;
1341  pmbl = pmbl->next;
1342  }
1343 
1344  if (res_flag == 0) {
1345 #pragma omp parallel for num_threads(nthreads)
1346  for (int i=0; i<nmb; ++i) {
1347  MeshBlock *pmb = pmb_array[i];
1348  pmb->ProblemGenerator(pin);
1349  pmb->pbval->CheckUserBoundaries();
1350  pmb->phydro->CheckHydro();
1351  pmb->pphy->Initialize(pmb->phydro->w);
1352  }
1353  }
1354 
1355  // add initial perturbation for decaying or impulsive turbulence
1356  if (((turb_flag == 1) || (turb_flag == 2)) && (res_flag == 0))
1357  ptrbd->Driving();
1358 
1359  // Create send/recv MPI_Requests for all BoundaryData objects
1360 #pragma omp parallel for num_threads(nthreads)
1361  for (int i=0; i<nmb; ++i) {
1362  MeshBlock *pmb = pmb_array[i];
1363  // BoundaryVariable objects evolved in main TimeIntegratorTaskList:
1364  pmb->pbval->SetupPersistentMPI();
1365  // other BoundaryVariable objects:
1366  if (SELF_GRAVITY_ENABLED == 1)
1367  pmb->pgrav->gbvar.SetupPersistentMPI();
1368  }
1369 
1370  // solve gravity for the first time
1371  if (SELF_GRAVITY_ENABLED == 1)
1372  pfgrd->Solve(1, 0);
1373  else if (SELF_GRAVITY_ENABLED == 2)
1374  pmgrd->Solve(1);
1375 
1376 #pragma omp parallel num_threads(nthreads)
1377  {
1378  MeshBlock *pmb;
1379  BoundaryValues *pbval;
1380 
1381  // prepare to receive conserved variables
1382 #pragma omp for private(pmb,pbval)
1383  for (int i=0; i<nmb; ++i) {
1384  pmb = pmb_array[i]; pbval = pmb->pbval;
1385  if (SHEARING_BOX) {
1386  pbval->ComputeShear(time);
1387  }
1388  pbval->StartReceiving(BoundaryCommSubset::mesh_init);
1389  }
1390 
1391  // send conserved variables
1392 #pragma omp for private(pmb,pbval)
1393  for (int i=0; i<nmb; ++i) {
1394  pmb = pmb_array[i]; pbval = pmb->pbval;
1395  pmb->phydro->hbvar.SwapHydroQuantity(pmb->phydro->u,
1397  pmb->phydro->hbvar.SendBoundaryBuffers();
1398  if (MAGNETIC_FIELDS_ENABLED)
1399  pmb->pfield->fbvar.SendBoundaryBuffers();
1400  // and (conserved variable) passive scalar masses:
1401  if (NSCALARS > 0)
1402  pmb->pscalars->sbvar.SendBoundaryBuffers();
1403  }
1404 
1405  // wait to receive conserved variables
1406 #pragma omp for private(pmb,pbval)
1407  for (int i=0; i<nmb; ++i) {
1408  pmb = pmb_array[i]; pbval = pmb->pbval;
1409  pmb->phydro->hbvar.ReceiveAndSetBoundariesWithWait();
1410  if (MAGNETIC_FIELDS_ENABLED)
1411  pmb->pfield->fbvar.ReceiveAndSetBoundariesWithWait();
1412  if (NSCALARS > 0)
1413  pmb->pscalars->sbvar.ReceiveAndSetBoundariesWithWait();
1414  if (SHEARING_BOX) {
1416  }
1417  pbval->ClearBoundary(BoundaryCommSubset::mesh_init);
1418  }
1419 
1420  // With AMR/SMR GR send primitives to enable cons->prim before prolongation
1421  if (GENERAL_RELATIVITY && multilevel) {
1422  // prepare to receive primitives
1423 #pragma omp for private(pmb,pbval)
1424  for (int i=0; i<nmb; ++i) {
1425  pmb = pmb_array[i]; pbval = pmb->pbval;
1426  pbval->StartReceiving(BoundaryCommSubset::gr_amr);
1427  }
1428 
1429  // send primitives
1430 #pragma omp for private(pmb,pbval)
1431  for (int i=0; i<nmb; ++i) {
1432  pmb = pmb_array[i]; pbval = pmb->pbval;
1433  pmb->phydro->hbvar.SwapHydroQuantity(pmb->phydro->w,
1435  pmb->phydro->hbvar.SendBoundaryBuffers();
1436  }
1437 
1438  // wait to receive AMR/SMR GR primitives
1439 #pragma omp for private(pmb,pbval)
1440  for (int i=0; i<nmb; ++i) {
1441  pmb = pmb_array[i]; pbval = pmb->pbval;
1442  pmb->phydro->hbvar.ReceiveAndSetBoundariesWithWait();
1443  pbval->ClearBoundary(BoundaryCommSubset::gr_amr);
1444  pmb->phydro->hbvar.SwapHydroQuantity(pmb->phydro->u,
1446  }
1447  } // multilevel
1448 
1449  // perform fourth-order correction of midpoint initial condition:
1450  // (correct IC on all MeshBlocks or none; switch cannot be toggled independently)
1451  bool correct_ic = pmb_array[0]->precon->correct_ic;
1452  if (correct_ic)
1453  CorrectMidpointInitialCondition(pmb_array, nmb);
1454 
1455  // Now do prolongation, compute primitives, apply BCs
1456  Hydro *ph;
1457  Field *pf;
1458  PassiveScalars *ps;
1459 #pragma omp for private(pmb,pbval,ph,pf,ps)
1460  for (int i=0; i<nmb; ++i) {
1461  pmb = pmb_array[i];
1462  pbval = pmb->pbval, ph = pmb->phydro, pf = pmb->pfield, ps = pmb->pscalars;
1463  if (multilevel)
1464  pbval->ProlongateBoundaries(time, 0.0);
1465 
1466  int il = pmb->is, iu = pmb->ie,
1467  jl = pmb->js, ju = pmb->je,
1468  kl = pmb->ks, ku = pmb->ke;
1469  if (pbval->nblevel[1][1][0] != -1) il -= NGHOST;
1470  if (pbval->nblevel[1][1][2] != -1) iu += NGHOST;
1471  if (pmb->block_size.nx2 > 1) {
1472  if (pbval->nblevel[1][0][1] != -1) jl -= NGHOST;
1473  if (pbval->nblevel[1][2][1] != -1) ju += NGHOST;
1474  }
1475  if (pmb->block_size.nx3 > 1) {
1476  if (pbval->nblevel[0][1][1] != -1) kl -= NGHOST;
1477  if (pbval->nblevel[2][1][1] != -1) ku += NGHOST;
1478  }
1479  pmb->peos->ConservedToPrimitive(ph->u, ph->w1, pf->b,
1480  ph->w, pf->bcc, pmb->pcoord,
1481  il, iu, jl, ju, kl, ku);
1482  if (NSCALARS > 0) {
1483  // r1/r_old for GR is currently unused:
1484  pmb->peos->PassiveScalarConservedToPrimitive(ps->s, ph->w, ps->r, ps->r,
1485  pmb->pcoord,
1486  il, iu, jl, ju, kl, ku);
1487  }
1488  // --------------------------
1489  int order = pmb->precon->xorder;
1490  if (order == 4) {
1491  // fourth-order EOS:
1492  // for hydro, shrink buffer by 1 on all sides
1493  if (pbval->nblevel[1][1][0] != -1) il += 1;
1494  if (pbval->nblevel[1][1][2] != -1) iu -= 1;
1495  if (pbval->nblevel[1][0][1] != -1) jl += 1;
1496  if (pbval->nblevel[1][2][1] != -1) ju -= 1;
1497  if (pbval->nblevel[0][1][1] != -1) kl += 1;
1498  if (pbval->nblevel[2][1][1] != -1) ku -= 1;
1499  // for MHD, shrink buffer by 3
1500  // TODO(felker): add MHD loop limit calculation for 4th order W(U)
1501  pmb->peos->ConservedToPrimitiveCellAverage(ph->u, ph->w1, pf->b,
1502  ph->w, pf->bcc, pmb->pcoord,
1503  il, iu, jl, ju, kl, ku);
1504  if (NSCALARS > 0) {
1506  ps->s, ps->r, ps->r, pmb->pcoord, il, iu, jl, ju, kl, ku);
1507  }
1508  }
1509  // --------------------------
1510  // end fourth-order EOS
1511 
1512  // Swap Hydro and (possibly) passive scalar quantities in BoundaryVariable
1513  // interface from conserved to primitive formulations:
1515  if (NSCALARS > 0)
1516  ps->sbvar.var_cc = &(ps->r);
1517 
1518  pbval->ApplyPhysicalBoundaries(time, 0.0);
1519  }
1520 
1521  // Calc initial diffusion coefficients
1522 #pragma omp for private(pmb,ph,pf)
1523  for (int i=0; i<nmb; ++i) {
1524  pmb = pmb_array[i]; ph = pmb->phydro, pf = pmb->pfield;
1525  if (ph->hdif.hydro_diffusion_defined)
1526  ph->hdif.SetDiffusivity(ph->w, pf->bcc);
1527  if (MAGNETIC_FIELDS_ENABLED) {
1528  if (pf->fdif.field_diffusion_defined)
1529  pf->fdif.SetDiffusivity(ph->w, pf->bcc);
1530  }
1531  }
1532 
1533  if (!res_flag && adaptive) {
1534 #pragma omp for
1535  for (int i=0; i<nmb; ++i) {
1536  pmb_array[i]->pmr->CheckRefinementCondition();
1537  }
1538  }
1539  } // omp parallel
1540 
1541  if (!res_flag && adaptive) {
1542  iflag = false;
1543  int onb = nbtotal;
1545  if (nbtotal == onb) {
1546  iflag = true;
1547  } else if (nbtotal < onb && Globals::my_rank == 0) {
1548  std::cout << "### Warning in Mesh::Initialize" << std::endl
1549  << "The number of MeshBlocks decreased during AMR grid initialization."
1550  << std::endl
1551  << "Possibly the refinement criteria have a problem." << std::endl;
1552  }
1553  if (nbtotal > 2*inb && Globals::my_rank == 0) {
1554  std::cout
1555  << "### Warning in Mesh::Initialize" << std::endl
1556  << "The number of MeshBlocks increased more than twice during initialization."
1557  << std::endl
1558  << "More computing power than you expected may be required." << std::endl;
1559  }
1560  }
1561  } while (!iflag);
1562 
1563  // calculate the first time step
1564 #pragma omp parallel for num_threads(nthreads)
1565  for (int i=0; i<nmb; ++i) {
1566  pmb_array[i]->phydro->NewBlockTimeStep();
1567  }
1568 
1569  NewTimeStep();
1570  return;
1571 }
1572 
1573 //----------------------------------------------------------------------------------------
1575 // \brief return the MeshBlock whose gid is tgid
1576 
1578  MeshBlock *pbl = pblock;
1579  while (pbl != nullptr) {
1580  if (pbl->gid == tgid)
1581  break;
1582  pbl = pbl->next;
1583  }
1584  return pbl;
1585 }
1586 
1587 //----------------------------------------------------------------------------------------
1588 // \!fn void Mesh::SetBlockSizeAndBoundaries(LogicalLocation loc,
1589 // RegionSize &block_size, BundaryFlag *block_bcs)
1590 // \brief Set the physical part of a block_size structure and block boundary conditions
1591 
1593  BoundaryFlag *block_bcs) {
1594  std::int64_t &lx1 = loc.lx1;
1595  int &ll = loc.level;
1596  std::int64_t nrbx_ll = nrbx1 << (ll - root_level);
1597 
1598  // calculate physical block size, x1
1599  if (lx1 == 0) {
1600  block_size.x1min = mesh_size.x1min;
1601  block_bcs[BoundaryFace::inner_x1] = mesh_bcs[BoundaryFace::inner_x1];
1602  } else {
1604  block_size.x1min = MeshGenerator_[X1DIR](rx, mesh_size);
1605  block_bcs[BoundaryFace::inner_x1] = BoundaryFlag::block;
1606  }
1607  if (lx1 == nrbx_ll - 1) {
1608  block_size.x1max = mesh_size.x1max;
1609  block_bcs[BoundaryFace::outer_x1] = mesh_bcs[BoundaryFace::outer_x1];
1610  } else {
1612  block_size.x1max = MeshGenerator_[X1DIR](rx, mesh_size);
1613  block_bcs[BoundaryFace::outer_x1] = BoundaryFlag::block;
1614  }
1615 
1616  // calculate physical block size, x2
1617  if (mesh_size.nx2 == 1) {
1618  block_size.x2min = mesh_size.x2min;
1619  block_size.x2max = mesh_size.x2max;
1620  block_bcs[BoundaryFace::inner_x2] = mesh_bcs[BoundaryFace::inner_x2];
1621  block_bcs[BoundaryFace::outer_x2] = mesh_bcs[BoundaryFace::outer_x2];
1622  } else {
1623  std::int64_t &lx2 = loc.lx2;
1624  nrbx_ll = nrbx2 << (ll - root_level);
1625  if (lx2 == 0) {
1626  block_size.x2min = mesh_size.x2min;
1627  block_bcs[BoundaryFace::inner_x2] = mesh_bcs[BoundaryFace::inner_x2];
1628  } else {
1630  block_size.x2min = MeshGenerator_[X2DIR](rx, mesh_size);
1631  block_bcs[BoundaryFace::inner_x2] = BoundaryFlag::block;
1632  }
1633  if (lx2 == (nrbx_ll) - 1) {
1634  block_size.x2max = mesh_size.x2max;
1635  block_bcs[BoundaryFace::outer_x2] = mesh_bcs[BoundaryFace::outer_x2];
1636  } else {
1638  block_size.x2max = MeshGenerator_[X2DIR](rx, mesh_size);
1639  block_bcs[BoundaryFace::outer_x2] = BoundaryFlag::block;
1640  }
1641  }
1642 
1643  // calculate physical block size, x3
1644  if (mesh_size.nx3 == 1) {
1645  block_size.x3min = mesh_size.x3min;
1646  block_size.x3max = mesh_size.x3max;
1647  block_bcs[BoundaryFace::inner_x3] = mesh_bcs[BoundaryFace::inner_x3];
1648  block_bcs[BoundaryFace::outer_x3] = mesh_bcs[BoundaryFace::outer_x3];
1649  } else {
1650  std::int64_t &lx3 = loc.lx3;
1651  nrbx_ll = nrbx3 << (ll - root_level);
1652  if (lx3 == 0) {
1653  block_size.x3min = mesh_size.x3min;
1654  block_bcs[BoundaryFace::inner_x3] = mesh_bcs[BoundaryFace::inner_x3];
1655  } else {
1657  block_size.x3min = MeshGenerator_[X3DIR](rx, mesh_size);
1658  block_bcs[BoundaryFace::inner_x3] = BoundaryFlag::block;
1659  }
1660  if (lx3 == (nrbx_ll) - 1) {
1661  block_size.x3max = mesh_size.x3max;
1662  block_bcs[BoundaryFace::outer_x3] = mesh_bcs[BoundaryFace::outer_x3];
1663  } else {
1665  block_size.x3max = MeshGenerator_[X3DIR](rx, mesh_size);
1666  block_bcs[BoundaryFace::outer_x3] = BoundaryFlag::block;
1667  }
1668  }
1669 
1670  block_size.x1rat = mesh_size.x1rat;
1671  block_size.x2rat = mesh_size.x2rat;
1672  block_size.x3rat = mesh_size.x3rat;
1673 
1674  return;
1675 }
1676 
1677 
1678 void Mesh::CorrectMidpointInitialCondition(std::vector<MeshBlock*> &pmb_array, int nmb) {
1679  MeshBlock *pmb;
1680  Hydro *ph;
1681  PassiveScalars *ps;
1682 #pragma omp for private(pmb, ph, ps)
1683  for (int nb=0; nb<nmb; ++nb) {
1684  pmb = pmb_array[nb];
1685  ph = pmb->phydro;
1686  ps = pmb->pscalars;
1687 
1688  // Assume cell-centered analytic value is computed at all real cells, and ghost
1689  // cells with the cell-centered U have been exchanged
1690  int il = pmb->is, iu = pmb->ie, jl = pmb->js, ju = pmb->je,
1691  kl = pmb->ks, ku = pmb->ke;
1692 
1693  // Laplacian of cell-averaged conserved variables, scalar concentrations
1694  AthenaArray<Real> delta_cons_, delta_s_;
1695 
1696  // Allocate memory for 4D Laplacian
1697  int ncells4 = NHYDRO;
1698  int nl = 0;
1699  int nu = ncells4 - 1;
1700  delta_cons_.NewAthenaArray(ncells4, pmb->ncells3, pmb->ncells2, pmb->ncells1);
1701 
1702  // Compute and store Laplacian of cell-averaged conserved variables
1703  pmb->pcoord->Laplacian(ph->u, delta_cons_, il, iu, jl, ju, kl, ku, nl, nu);
1704 
1705  // TODO(felker): assuming uniform mesh with dx1f=dx2f=dx3f, so this factors out
1706  // TODO(felker): also, this may need to be dx1v, since Laplacian is cell-center
1707  Real h = pmb->pcoord->dx1f(il); // pco->dx1f(i); inside loop
1708  Real C = (h*h)/24.0;
1709 
1710  // Compute fourth-order approximation to cell-centered conserved variables
1711  for (int n=nl; n<=nu; ++n) {
1712  for (int k=kl; k<=ku; ++k) {
1713  for (int j=jl; j<=ju; ++j) {
1714  for (int i=il; i<=iu; ++i) {
1715  // We do not actually need to store all cell-centered cons. variables,
1716  // but the ConservedToPrimitivePointwise() implementation operates on 4D
1717  ph->u(n,k,j,i) = ph->u(n,k,j,i) + C*delta_cons_(n,k,j,i);
1718  }
1719  }
1720  }
1721  }
1722 
1723  // If NSCALARS < NHYDRO, could reuse delta_cons_ allocated memory...
1724  int ncells4_s = NSCALARS;
1725  int nl_s = 0;
1726  int nu_s = ncells4_s -1;
1727  if (NSCALARS > 0) {
1728  delta_s_.NewAthenaArray(ncells4_s, pmb->ncells3, pmb->ncells2, pmb->ncells1);
1729  pmb->pcoord->Laplacian(ps->s, delta_s_, il, iu, jl, ju, kl, ku, nl_s, nu_s);
1730  }
1731 
1732  // Compute fourth-order approximation to cell-centered conserved variables
1733  for (int n=nl_s; n<=nu_s; ++n) {
1734  for (int k=kl; k<=ku; ++k) {
1735  for (int j=jl; j<=ju; ++j) {
1736  for (int i=il; i<=iu; ++i) {
1737  ps->s(n,k,j,i) = ps->s(n,k,j,i) + C*delta_s_(n,k,j,i);
1738  }
1739  }
1740  }
1741  }
1742  } // end loop over MeshBlocks
1743 
1744  // begin second exchange of ghost cells with corrected cell-averaged <U>
1745  // ----------------- (mostly copied from above section in Mesh::Initialize())
1746  BoundaryValues *pbval;
1747  // prepare to receive conserved variables
1748 #pragma omp for private(pmb,pbval)
1749  for (int i=0; i<nmb; ++i) {
1750  pmb = pmb_array[i]; pbval = pmb->pbval;
1751  if (SHEARING_BOX) {
1752  pbval->ComputeShear(time);
1753  }
1754  // no need to re-SetupPersistentMPI() the MPI requests for boundary values
1755  pbval->StartReceiving(BoundaryCommSubset::mesh_init);
1756  }
1757 
1758 #pragma omp for private(pmb,pbval)
1759  for (int i=0; i<nmb; ++i) {
1760  pmb = pmb_array[i];
1761  pmb->phydro->hbvar.SwapHydroQuantity(pmb->phydro->u,
1763  pmb->phydro->hbvar.SendBoundaryBuffers();
1764  if (MAGNETIC_FIELDS_ENABLED)
1765  pmb->pfield->fbvar.SendBoundaryBuffers();
1766  // and (conserved variable) passive scalar masses:
1767  if (NSCALARS > 0)
1768  pmb->pscalars->sbvar.SendBoundaryBuffers();
1769  }
1770 
1771  // wait to receive conserved variables
1772 #pragma omp for private(pmb,pbval)
1773  for (int i=0; i<nmb; ++i) {
1774  pmb = pmb_array[i]; pbval = pmb->pbval;
1775  pmb->phydro->hbvar.SwapHydroQuantity(pmb->phydro->u,
1777  pmb->phydro->hbvar.ReceiveAndSetBoundariesWithWait();
1778  if (MAGNETIC_FIELDS_ENABLED)
1779  pmb->pfield->fbvar.ReceiveAndSetBoundariesWithWait();
1780  if (NSCALARS > 0)
1781  pmb->pscalars->sbvar.ReceiveAndSetBoundariesWithWait();
1782  if (SHEARING_BOX) {
1784  }
1785  pbval->ClearBoundary(BoundaryCommSubset::mesh_init);
1786  } // end second exchange of ghost cells
1787  return;
1788 }
1789 
1790 // Public function for advancing next_phys_id_ counter
1791 // E.g. if chemistry or radiation elects to communicate additional information with MPI
1792 // outside the framework of the BoundaryVariable classes
1793 
1794 // Store signed, but positive, integer corresponding to the next unused value to be used
1795 // as unique ID for a BoundaryVariable object's single set of MPI calls (formerly "enum
1796 // AthenaTagMPI"). 5 bits of unsigned integer representation are currently reserved
1797 // for this "phys" part of the bitfield tag, making 0, ..., 31 legal values
1798 
1799 int Mesh::ReserveTagPhysIDs(int num_phys) {
1800  // TODO(felker): add safety checks? input, output are positive, obey <= 31= MAX_NUM_PHYS
1801  int start_id = next_phys_id_;
1802  next_phys_id_ += num_phys;
1803  return start_id;
1804 }
1805 
1806 // private member fn, called in Mesh() ctor
1807 
1808 // depending on compile- and runtime options, reserve the maximum number of "int physid"
1809 // that might be necessary for each MeshBlock's BoundaryValues object to perform MPI
1810 // communication for all BoundaryVariable objects
1811 
1812 // TODO(felker): deduplicate this logic, which combines conditionals in MeshBlock ctor
1813 
1815 #ifdef MPI_PARALLEL
1816  // if (FLUID_ENABLED) {
1817  // Advance Mesh's shared counter (initialized to next_phys_id=1 if MPI)
1818  // Greedy reservation of phys IDs (only 1 of 2 needed for Hydro if multilevel==false)
1819  ReserveTagPhysIDs(HydroBoundaryVariable::max_phys_id);
1820  // }
1821  if (MAGNETIC_FIELDS_ENABLED) {
1822  ReserveTagPhysIDs(FaceCenteredBoundaryVariable::max_phys_id);
1823  }
1824  if (SELF_GRAVITY_ENABLED) {
1825  ReserveTagPhysIDs(CellCenteredBoundaryVariable::max_phys_id);
1826  }
1827  if (NSCALARS > 0) {
1828  ReserveTagPhysIDs(CellCenteredBoundaryVariable::max_phys_id);
1829  }
1830 #endif
1831  return;
1832 }
1833 
1834 //----------------------------------------------------------------------------------------
1836 // \brief Parses input string to return scoped enumerator flag specifying boundary
1837 // condition. Typically called in Mesh() ctor initializer list
1838 
1839 FluidFormulation GetFluidFormulation(const std::string& input_string) {
1840  if (input_string == "true") {
1841  return FluidFormulation::evolve;
1842  } else if (input_string == "disabled") {
1844  } else if (input_string == "background") {
1846  } else {
1847  std::stringstream msg;
1848  msg << "### FATAL ERROR in GetFluidFormulation" << std::endl
1849  << "Input string=" << input_string << "\n"
1850  << "is an invalid fluid formulation" << std::endl;
1851  ATHENA_ERROR(msg);
1852  }
1853 }
1854 
1856  const int dt_precision = std::numeric_limits<Real>::max_digits10 - 1;
1857  const int ratio_precision = 3;
1858  if (ncycle_out != 0) {
1859  if (ncycle % ncycle_out == 0) {
1860  std::cout << "cycle=" << ncycle << std::scientific
1861  << std::setprecision(dt_precision)
1862  << " time=" << time << " dt=" << dt;
1863  if (dt_diagnostics != -1) {
1864  if (STS_ENABLED) {
1865  if (UserTimeStep_ == nullptr)
1866  std::cout << "=dt_hyperbolic";
1867  // remaining dt_parabolic diagnostic output handled in STS StartupTaskList
1868  } else {
1869  Real ratio = dt / dt_hyperbolic;
1870  std::cout << "\ndt_hyperbolic=" << dt_hyperbolic << " ratio="
1871  << std::setprecision(ratio_precision) << ratio
1872  << std::setprecision(dt_precision);
1873  ratio = dt / dt_parabolic;
1874  std::cout << "\ndt_parabolic=" << dt_parabolic << " ratio="
1875  << std::setprecision(ratio_precision) << ratio
1876  << std::setprecision(dt_precision);
1877  }
1878  if (UserTimeStep_ != nullptr) {
1879  Real ratio = dt / dt_user;
1880  std::cout << "\ndt_user=" << dt_user << " ratio="
1881  << std::setprecision(ratio_precision) << ratio
1882  << std::setprecision(dt_precision);
1883  }
1884  } // else (empty): dt_diagnostics = -1 -> provide no additional timestep diagnostics
1885  std::cout << std::endl;
1886  }
1887  }
1888  return;
1889 }
void(*)(FieldDiffusion *pfdif, MeshBlock *pmb, const AthenaArray< Real > &w, const AthenaArray< Real > &bmag, int is, int ie, int js, int je, int ks, int ke) FieldDiffusionCoeffFunc
Definition: athena.hpp:204
UserHistoryOperation
Definition: athena.hpp:167
void(*)(AthenaArray< Real > &dst, Real time, int nvar, int is, int ie, int js, int je, int ks, int ke, int ngh, Real x0, Real y0, Real z0, Real dx, Real dy, Real dz) MGBoundaryFunc
Definition: athena.hpp:191
double Real
Definition: athena.hpp:29
Real(*)(Real x, RegionSize rs) MeshGenFunc
Definition: athena.hpp:178
CoordinateDirection
Definition: athena.hpp:155
@ X2DIR
Definition: athena.hpp:155
@ X3DIR
Definition: athena.hpp:155
@ X1DIR
Definition: athena.hpp:155
void(*)(Real x1, Real x2, Real x3, ParameterInput *pin, AthenaArray< Real > &g, AthenaArray< Real > &g_inv, AthenaArray< Real > &dg_dx1, AthenaArray< Real > &dg_dx2, AthenaArray< Real > &dg_dx3) MetricFunc
Definition: athena.hpp:187
void(*)(HydroDiffusion *phdif, MeshBlock *pmb, const AthenaArray< Real > &w, const AthenaArray< Real > &bc, int is, int ie, int js, int je, int ks, int ke) ViscosityCoeffFunc
Definition: athena.hpp:195
FluidFormulation
Definition: athena.hpp:166
int(*)(MeshBlock *pmb) AMRFlagFunc
Definition: athena.hpp:177
Real(*)(MeshBlock *pmb) TimeStepFunc
Definition: athena.hpp:182
void(*)(MeshBlock *pmb, Coordinates *pco, AthenaArray< Real > &prim, FaceField &b, Real time, Real dt, int is, int ie, int js, int je, int ks, int ke, int ngh) BValFunc
Definition: athena.hpp:176
void(*)(MeshBlock *pmb, const Real time, const Real dt, const AthenaArray< Real > &prim, const AthenaArray< Real > &bcc, AthenaArray< Real > &cons) SrcTermFunc
Definition: athena.hpp:181
Real(*)(MeshBlock *pmb, int iout) HistoryOutputFunc
Definition: athena.hpp:183
void(*)(HydroDiffusion *phdif, MeshBlock *pmb, const AthenaArray< Real > &w, const AthenaArray< Real > &bc, int is, int ie, int js, int je, int ks, int ke) ConductionCoeffFunc
Definition: athena.hpp:199
std::size_t GetSizeInBytes() const
AthenaArray< Real > dx1f
Definition: coordinates.hpp:41
virtual void Laplacian(const AthenaArray< Real > &s, AthenaArray< Real > &delta_s, const int il, const int iu, const int jl, const int ju, const int kl, const int ku, const int nl, const int nu)
void ConservedToPrimitiveCellAverage(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)
void PassiveScalarConservedToPrimitiveCellAverage(AthenaArray< Real > &s, const AthenaArray< Real > &r_old, AthenaArray< Real > &r, Coordinates *pco, int il, int iu, int jl, int ju, int kl, int ku)
Definition: eos_scalars.cpp:68
void PassiveScalarConservedToPrimitive(AthenaArray< Real > &s, const AthenaArray< Real > &w, const AthenaArray< Real > &r_old, AthenaArray< Real > &r, Coordinates *pco, int il, int iu, int jl, int ju, int kl, int ku)
Definition: eos_scalars.cpp:32
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)
void SwapHydroQuantity(AthenaArray< Real > &var_hydro, HydroBoundaryQuantity hydro_type)
Definition: hydro.hpp:34
AthenaArray< Real > w
Definition: hydro.hpp:46
HydroDiffusion hdif
Definition: hydro.hpp:72
AthenaArray< Real > u
Definition: hydro.hpp:46
HydroBoundaryVariable hbvar
Definition: hydro.hpp:70
AthenaArray< Real > w1
Definition: hydro.hpp:47
void CheckHydro()
Definition: hydro.cpp:206
InputBlock * pnext
std::string block_name
Reconstruction * precon
Definition: mesh.hpp:124
Physics * pphy
Definition: mesh.hpp:138
BoundaryValues * pbval
Definition: mesh.hpp:123
Gravity * pgrav
Definition: mesh.hpp:130
void UserWorkBeforeOutput(ParameterInput *pin)
MeshBlock * prev
Definition: mesh.hpp:145
void ProblemGenerator(ParameterInput *pin)
Real new_block_dt_
Definition: mesh.hpp:171
RegionSize block_size
Definition: mesh.hpp:94
Hydro * phydro
Definition: mesh.hpp:128
EquationOfState * peos
Definition: mesh.hpp:133
int ie
Definition: mesh.hpp:100
int js
Definition: mesh.hpp:100
Real new_block_dt_hyperbolic_
Definition: mesh.hpp:171
int ncells1
Definition: mesh.hpp:97
int je
Definition: mesh.hpp:100
Real new_block_dt_parabolic_
Definition: mesh.hpp:171
int ke
Definition: mesh.hpp:100
PassiveScalars * pscalars
Definition: mesh.hpp:132
int ncells3
Definition: mesh.hpp:97
Real new_block_dt_user_
Definition: mesh.hpp:172
int is
Definition: mesh.hpp:100
int ks
Definition: mesh.hpp:100
Coordinates * pcoord
Definition: mesh.hpp:122
int gid
Definition: mesh.hpp:101
int ncells2
Definition: mesh.hpp:97
std::size_t GetBlockSizeInBytes()
Definition: meshblock.cpp:545
MeshBlock * next
Definition: mesh.hpp:145
Field * pfield
Definition: mesh.hpp:129
int GetNumMeshBlocksThisRank(int my_rank)
Definition: mesh.hpp:236
int * rdisp
Definition: mesh.hpp:296
void EnrollFieldDiffusivity(FieldDiffusionCoeffFunc my_func)
Definition: mesh.cpp:1274
void OutputCycleDiagnostics()
Definition: mesh.cpp:1855
int * nslist
Definition: mesh.hpp:292
int nrbx1
Definition: mesh.hpp:305
int nint_user_mesh_data_
Definition: mesh.hpp:311
friend class BoundaryValues
Definition: mesh.hpp:209
int * nref
Definition: mesh.hpp:295
ConductionCoeffFunc ConductionCoeff_
Definition: mesh.hpp:334
void InitUserMeshData(ParameterInput *pin)
double lb_tolerance_
Definition: mesh.hpp:322
FFTGravityDriver * pfgrd
Definition: mesh.hpp:262
Real dt_parabolic
Definition: mesh.hpp:248
int nreal_user_mesh_data_
Definition: mesh.hpp:311
Real dt
Definition: mesh.hpp:248
int lb_interval_
Definition: mesh.hpp:323
const int ndim
Definition: mesh.hpp:245
UserHistoryOperation * user_history_ops_
Definition: mesh.hpp:315
BoundaryFlag mesh_bcs[6]
Definition: mesh.hpp:243
void EnrollUserMGGravityBoundaryFunction(BoundaryFace dir, MGBoundaryFunc my_bc)
Definition: mesh.cpp:1128
int num_mesh_threads_
Definition: mesh.hpp:291
void EnrollUserTimeStepFunction(TimeStepFunc my_func)
Definition: mesh.cpp:1207
void CalculateLoadBalance(double *clist, int *rlist, int *slist, int *nlist, int nb)
const bool multilevel
Definition: mesh.hpp:246
int * bnderef
Definition: mesh.hpp:297
Real dt_user
Definition: mesh.hpp:248
FieldDiffusionCoeffFunc FieldDiffusivity_
Definition: mesh.hpp:335
void EnrollUserBoundaryFunction(BoundaryFace face, BValFunc my_func)
Definition: mesh.cpp:1106
int nrbx2
Definition: mesh.hpp:305
int current_level
Definition: mesh.hpp:290
int next_phys_id_
Definition: mesh.hpp:289
int * ranklist
Definition: mesh.hpp:292
int ncycle_out
Definition: mesh.hpp:249
const bool adaptive
Definition: mesh.hpp:246
int * ddisp
Definition: mesh.hpp:296
void EnrollUserMetric(MetricFunc my_func)
Definition: mesh.cpp:1247
int nbtotal
Definition: mesh.hpp:251
~Mesh()
Definition: mesh.cpp:876
int * nblist
Definition: mesh.hpp:292
Mesh(ParameterInput *pin, int test_flag=0)
Definition: mesh.cpp:63
Real tlim
Definition: mesh.hpp:248
SrcTermFunc UserSourceTerm_
Definition: mesh.hpp:329
bool lb_manual_
Definition: mesh.hpp:321
AMRFlagFunc AMRFlag_
Definition: mesh.hpp:328
int max_level
Definition: mesh.hpp:290
int * nderef
Definition: mesh.hpp:295
void NewTimeStep()
Definition: mesh.cpp:1065
TurbulenceDriver * ptrbd
Definition: mesh.hpp:261
int ReserveTagPhysIDs(int num_phys)
Definition: mesh.cpp:1799
void ResetLoadBalanceVariables()
int * brdisp
Definition: mesh.hpp:298
void EnrollUserExplicitSourceFunction(SrcTermFunc my_func)
Definition: mesh.cpp:1198
bool lb_automatic_
Definition: mesh.hpp:321
MeshBlockTree tree
Definition: mesh.hpp:302
void Initialize(int res_flag, ParameterInput *pin)
Definition: mesh.cpp:1327
HistoryOutputFunc * user_history_func_
Definition: mesh.hpp:331
void AllocateIntUserMeshDataField(int n)
Definition: mesh.cpp:1298
const bool f3
Definition: mesh.hpp:244
int dt_diagnostics
Definition: mesh.hpp:249
MGGravityDriver * pmgrd
Definition: mesh.hpp:263
void EnrollUserHistoryOutput(int i, HistoryOutputFunc my_func, const char *name, UserHistoryOperation op=UserHistoryOperation::sum)
Definition: mesh.cpp:1229
Real time
Definition: mesh.hpp:248
void EnrollUserRefinementCondition(AMRFlagFunc amrflag)
Definition: mesh.cpp:1154
int * bnref
Definition: mesh.hpp:297
friend class MeshBlock
Definition: mesh.hpp:206
int nrbx3
Definition: mesh.hpp:305
MGBoundaryFunc MGGravityBoundaryFunction_[6]
Definition: mesh.hpp:336
BValFunc BoundaryFunction_[6]
Definition: mesh.hpp:327
int root_level
Definition: mesh.hpp:290
EosTable * peos_table
Definition: mesh.hpp:256
void SetBlockSizeAndBoundaries(LogicalLocation loc, RegionSize &block_size, BoundaryFlag *block_bcs)
Definition: mesh.cpp:1592
void ReserveMeshBlockPhysIDs()
Definition: mesh.cpp:1814
const bool f2
Definition: mesh.hpp:244
double * costlist
Definition: mesh.hpp:293
void LoadBalancingAndAdaptiveMeshRefinement(ParameterInput *pin)
void EnrollUserMeshGenerator(CoordinateDirection dir, MeshGenFunc my_mg)
Definition: mesh.cpp:1164
void OutputMeshStructure(int dim)
Definition: mesh.cpp:915
friend class TurbulenceDriver
Definition: mesh.hpp:219
int nuser_history_output_
Definition: mesh.hpp:313
std::string * user_history_output_names_
Definition: mesh.hpp:314
void EnrollViscosityCoefficient(ViscosityCoeffFunc my_func)
Definition: mesh.cpp:1256
void AllocateUserHistoryOutput(int n)
Definition: mesh.cpp:1216
friend class MGGravityDriver
Definition: mesh.hpp:221
MeshBlock * pblock
Definition: mesh.hpp:259
int gflag
Definition: mesh.hpp:254
MetricFunc UserMetric_
Definition: mesh.hpp:332
Real dt_hyperbolic
Definition: mesh.hpp:248
friend class FFTGravityDriver
Definition: mesh.hpp:218
void ApplyUserWorkBeforeOutput(ParameterInput *pin)
Definition: mesh.cpp:1315
MeshBlock * FindMeshBlock(int tgid)
Definition: mesh.cpp:1577
RegionSize mesh_size
Definition: mesh.hpp:242
int ncycle
Definition: mesh.hpp:249
AthenaArray< int > * iuser_mesh_data
Definition: mesh.hpp:266
int turb_flag
Definition: mesh.hpp:255
int GetNumMeshThreads() const
Definition: mesh.hpp:237
bool use_uniform_meshgen_fn_[3]
Definition: mesh.hpp:310
void EnrollConductionCoefficient(ConductionCoeffFunc my_func)
Definition: mesh.cpp:1265
AthenaArray< Real > * ruser_mesh_data
Definition: mesh.hpp:265
LogicalLocation * loclist
Definition: mesh.hpp:301
int * bddisp
Definition: mesh.hpp:298
TimeStepFunc UserTimeStep_
Definition: mesh.hpp:330
ViscosityCoeffFunc ViscosityCoeff_
Definition: mesh.hpp:333
MeshGenFunc MeshGenerator_[3]
Definition: mesh.hpp:326
void CorrectMidpointInitialCondition(std::vector< MeshBlock * > &pmb_array, int nmb)
Definition: mesh.cpp:1678
void AllocateRealUserMeshDataField(int n)
Definition: mesh.cpp:1282
Real GetReal(std::string block, std::string name)
int GetOrAddInteger(std::string block, std::string name, int value)
int GetInteger(std::string block, std::string name)
std::string GetOrAddString(std::string block, std::string name, std::string value)
Real GetOrAddReal(std::string block, std::string name, Real value)
InputBlock * pfirst_block
CellCenteredBoundaryVariable sbvar
Definition: scalars.hpp:51
AthenaArray< Real > r
Definition: scalars.hpp:37
AthenaArray< Real > s
Definition: scalars.hpp:35
double min(double x1, double x2, double x3)
Definition: core.h:15
double max(double x1, double x2, double x3)
Definition: core.h:18
FluidFormulation GetFluidFormulation(const std::string &input_string)
Definition: mesh.cpp:1839
Real DefaultMeshGeneratorX3(Real x, RegionSize rs)
Definition: mesh.hpp:459
Real DefaultMeshGeneratorX2(Real x, RegionSize rs)
Definition: mesh.hpp:442
Real UniformMeshGeneratorX1(Real x, RegionSize rs)
Definition: mesh.hpp:476
Real UniformMeshGeneratorX2(Real x, RegionSize rs)
Definition: mesh.hpp:485
Real DefaultMeshGeneratorX1(Real x, RegionSize rs)
Definition: mesh.hpp:424
Real UniformMeshGeneratorX3(Real x, RegionSize rs)
Definition: mesh.hpp:493
Real ComputeMeshGeneratorX(std::int64_t index, std::int64_t nrange, bool sym_interval)
Definition: mesh.hpp:402
int nranks
Definition: globals.cpp:24
int my_rank
Definition: globals.cpp:23
std::int64_t lx3
Definition: athena.hpp:64
std::int64_t lx2
Definition: athena.hpp:64
std::int64_t lx1
Definition: athena.hpp:64
Real x3rat
Definition: athena.hpp:86
int nx3
Definition: athena.hpp:88
Real x1max
Definition: athena.hpp:85
Real x1min
Definition: athena.hpp:84
Real x2rat
Definition: athena.hpp:86
Real x3min
Definition: athena.hpp:84
Real x1rat
Definition: athena.hpp:86
int nx2
Definition: athena.hpp:88
Real x3max
Definition: athena.hpp:85
Real x2min
Definition: athena.hpp:84
int nx1
Definition: athena.hpp:88
Real x2max
Definition: athena.hpp:85