Add passive scalars

[felker]

Description

WIP. Closes #152.

Despite the discussions in that Issue, @jmstone, @tomidakn, and I decided to postpone acting on any algorithmic decisions for passive scalars in this pull request. Coming out of the UNLV Developer's Meeting, we have been uncertain about a couple of aspects:

• Do we actually need to define AthenaArray<Real> s(NSCALARS, k, j, i) as a relative concentration strictly in [0,1]? Or could we use an alternative formulation?
• Should we use / do we need a Plewa and Mueller (1999) renormalization step?

Therefore, I have implemented only the simplest interpretation of passive scalars, whereby their fluxes are computed in scalars/calculate_scalar_fluxes.cpp by:

• reconstructing sl_, sr_ face states along each direction
• upwinding these face states by the mass flux, e.g. the sign of Hydro::flux[X1DIR] is used to select s_upwind= either sl_ or sr_ and then it weights the flux by:

s_flux[X1DIR] = s_upwind*Hydro::flux[X1DIR]

Currently, there are no variable floors nor renormalization steps for the passive scalar integration. I was hoping that @jdolence and @pdmullen could weigh in with some of their hard-earned insights on these algorithmic choices. Patrick assures me that a new multimaterials implementation could easily be built on top of this branch by 1) adding user-defined source terms for volume fractions and 2) extending the EquationOfState class to allow the EOS to be dependent on volume fractions.

In order to get my passive scalar slotted cylinder test to work, I needed to partially address #247, the ability to turn off fluid evolution. I took the approach I outlined in #247 (comment)
and added

hydro/active=background

to my 2x hydro/athinput.slotted_cylinder2d* input files. This switch is stored as a public constant i the Mesh class, and its default value is true. The other future option would be false or disabled (for an empty/dummy/nonexistent Hydro class). Since we still need the Hydro fluxes for the passive scalar fluxes, the only things that background disables in time_integrator.cpp are the following function calls (inside the task functions, not in the task list construction):

• AddFluxDivergence()
• AddCoordTermsDivergence()
• AddSourceTermsHydro()
• DiffuseHydro()

@tomidakn, I could use your input on this feature. I haven't touched STS, etc. other task lists.

Testing and validation

Passive scalar advection of a slotted cylinder profile with AMR (RK3+PPM):

To-do

• Decide if/how to add an implementation of Gauss-Legendre quadrature to the Athena++ source code. It is currently required for slotted_cylinder.cpp. My original forked version used a copy of the 2x files from this open-source library: http://www.holoborodko.com/pavel/numerical-methods/numerical-integration/ .
  • I will post the author's license in a separate comment.
  • The author's code is also in the GNU Scientific Library (GSL) in a modified form in gsl/gintegration/glfixed.c under GPLv3.
  • Locally, I have placed gauss_legendre.cpp, gauss_legendre.hpp in src/utils/. The code provides:
    • precomputed high-precision abscissae and weights for some values of n
    • root-finding algorithm for finding values for other n in gauss_legendre_tbl()
    • 2x fns for integration: gauss_legendre() for 1D and gauss_legendre_2D_cube,
  • @tomidakn has mentioned that he is also considering adding tools to utils/ for helping problem generators, like B-field initialization via a user-provided vector potential (including correction for SMR/AMR, etc.).
• Improve slotted_cylinder.cpp. Implement Mesh::UserWorkAfterLoop to calculate errors like linear_wave.cpp (accounting for pmb->precon->correct_err switch).
  • Decided to forgo correct_err support and repeat calculation of GL quadrature-based initial condition. Max error increases for as the mesh is resolved given the default, unrefined input parameters.
• Extend calculate_scalar_fluxes.cpp to time/xorder=4
  • The straightforward approach is to permanently store all of the face-centered Hydro fluxes (currently computed for each direction and immediately discarded once the 4th order face-averaged fluxes are computed from them), but this would increase the memory footprint significantly for that algorithm.
  • We could conditionally store only the face-centered mass flux in each direction.
  • The best way would be to interleave the x1 Hydro and PassiveScalar flux calculations, but this would make the code unacceptably more complicated and mess up the task list dependencies.
  • 3x deep copies in hydro/calculate_fluxes.cpp could probably be improved.
• Return to Extend Kelvin-Helmholtz problem generator  #187 and implement dye concentration initial condition with this passive scalars implementation.
• Consider creating a new inputs/ subfolder (not hydro/) for zero-fluid or background-fluid input files?
• Add some flooring routines to PassiveScalars and call them in:
  • calculate_scalar_fluxes.cpp for the reconstructed states (definitely PPM, but also PLM?)
  • eos/ for the cell-centered scalars
• Add tst/regression/scripts/tests/scalars/ regression tests:
  • Mignone (2014) curvilinear limiter correction terms for PLM and PLM
  • Slotted cylinder with and without mesh refinement. Not sure what the analyze() conditions should be, but we need to make sure things at least compile and run without problems.
• Combine new reconstruct/*_simple.cpp with original *.cpp reconstruction functions in a same way that is extensible to other uses.
• Add documentation to Wiki:
  • New Outputs options: output[N]/variable = s[M], r[M]
  • New Input file options: hydro/active = true, background, disabled, problem/nu_scalar_iso, hydro/sfloor (haven't added this)
• Clean up end of TimeIntegratorTaskList ctor with conditionals that prevent deadlocking with MPI. See discussion: 742ec3f
• Consider condensing Hydro and PassiveScalars classes into an inheritance hierarchy (that would accommodate future derived classes like ChemicalSpecies, e.g.). Hydro is very nearly a proper superset of PassiveScalars, including add_flux_divergence.cpp and calculate_fluxes.cpp, once the RiemannSolver() and flooring routines are abstracted.
  • Revisit naming decisions from 742ec3f : PassiveScalars renamed to Scalars (@jmstone: mathematical connotation), Abundances, Species?
  • Revisit naming decisions from e3638f1 : 2x AddFluxDivergence() to AddHydroFluxDivergence() and AddScalarFluxDivergence()?
• Make it possible to compile without a fluid or a fixed background fluid (with or without a dummy Hydro class?). See Enumerate use-cases for zero-fluid configurations and design matching configure.py options #247.
• Future PR and project: make the following compile-time choices/configure options/macros runtime choices/input file options/variables:
  • NSCALARS
  • NGHOST

More to-do:
(added on 4/30/19 when the correct dual formalism of conserved and primitive passive scalar quantities was properly incorporated):

• Decide on naming conventions for PassiveScalars::s, r, etc. (like Athena 4.2?) Other options: rho_c and c, C and c (for color)
• Adjust timestep restriction and calculation when hydro/active=background or disabled.
• I am not incredibly confident that S/AMR and 4th order passive scalar advection are both still working perfectly after this big secondary set of changes. We need stringent convergence tests based on the passive scalars to make sure that they are fully correct.
  • In particular, src/bvals/bvals_refine.cpp involves a lot of delicate switching between primitive and conserved variables, and much of it is only activated when physical boundaries (non-periodic) are encountered on a refined mesh (which never happens in the regression tests).
• Add DIFFUSE_SCLR task to task_list.hpp and associated function. Add to STS task list too?
• Add passive scalars to the argument list for a user-defined BValFunc? This would break backwards-compatibility, so we seek a more general approach.
• Add MPI + .rst regression test. It was silently broken for a few weeks.

[felker]

retest this please

[tomidakn]

If STS is not enabled with passive scalar diffusion, we need a trap there. Also, as diffusion in hydrodynamics and in passive scalar can be different, there is a possible case where we may want to use STS for some physics and normal integrator for other physics.

As the passive scalar flux is always associated with the mass flux in hydrodynamics, technically it is possible to merge the computation and the boundary communication to hydrodynamics. For the computation (flux and divergence), I think it is better to keep them separate, mostly for the general execution model. But for the communication, you may save the MPI messages (and the tasks) if we combine the hydro and scalar communication. On the other hand, if we do so, we may lose the opportunity of computation-communication overlapping. I think the current implementation is fine, but just want to leave my comments here.

Somehow the multigrid boundary failed on Travis (Mac+Clang+OpenMPI) . Is this reproducible? And is it related to the change in the boundary?

Last but not least, Documentation must be updated.

[felker]

• For now, I have completely disabled STS + passive scalar configurations, since none of the *_SCLR tasks have been incorporated via AddTask() operations in the SuperTimeStepTaskList ctor. It should be straightforward, and @pdmullen has agreed to revisit it soon.
• Wiki (configure options, inputs, outputs, diffusion processes) has been updated, although the Passive Scalars page itself is pretty sparse.
• You are completely correct about the tight coupling of the Hydro and PassiveScalars classes and tasks. I came to the same conclusion as you did, but I think this will need to be reevaluated in greater detail when chemistry and general EOS extensions to passive scalars are considered. You could also save a ton of memory and deep copying operations by making them more tightly coupled.
• The errors on Travis CI that randomly occur on macOS VMs on the MPI grav/ tests (FFT too, I believe) have been popping up with greater frequency within the past month or so. I am not too concerned about them, since these particular MPI tests have given me trouble on macOS Travis CI since last fall, and because restarting the tests usually results in passes.

[msbc]

👍 🎉