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