Merge changes from main branch

docs/cluster.md

@@ -183,14 +183,16 @@ Pertubations are controlled by the following parameters:
 <problem/cluster/init_perturb>
 # for the velocity field
 sigma_v = 0.0      # volume weighted RMS |v| in code velocity; default: 0.0 (disabled)
-k_peak_v = ???     # wavenumber in normalized units where the velocity spectrum peaks. No default value.
+l_peak_v = ???     # lengthscale (in code length units) where the velocity spectrum peaks. No default value.
+k_peak_v = ???     # (exclusive alternative to l_peak_v): wavenumber in normalized units where the velocity spectrum peaks. No default value.
 num_modes_v = 40   # (optional) number of wavemodes in spectral space; default: 40
 sol_weight_v = 1.0 # (optional) power in solenoidal (rotational) modes of the perturbation. Range between 0 (fully compressive) and 1.0 (default, fully solenoidal).
 rseed_v = 1        # (optional) integer seed for RNG for wavenumbers and amplitudes
 
 # for the magnetic field
 sigma_b = 0.0      # volume weighted RMS |B| in code magnetic; default: 0.0 (disabled)
-k_peak_b = ???     # wavenumber in normalized units where the magnetic field spectrum peaks. No default value.
+l_peak_b = ???     # lengthscale (in code length units) where the magnetic field spectrum peaks. No default value.
+k_peak_b = ???     # (exclusive alternative to l_peak_b): wavenumber in normalized units where the magnetic field spectrum peaks. No default value.
 num_modes_b = 40   # (optional) number of wavemodes in spectral space; default: 40
 rseed_b = 2        # (optional) integer seed for RNG for wavenumbers and amplitudes
 ```
@@ -306,6 +308,8 @@ kinetic_fraction = 0.3333
 ```
 These values are automatically normalized to sum to 1.0 at run time.
 
+### Thermal feedback
+
 Thermal feedback is deposited at a flat power density within a sphere of defined radius
 ```
 <problem/cluster/agn_feedback>
@@ -318,6 +322,8 @@ feedback, e.g.
 thermal_injected_mass = mdot * (1 - efficiency) * normalized_thermal_fraction;
 ```
 
+### Kinetic feedback
+
 Kinetic feedback is deposited into two disks along the axis of the jet within a
 defined radius, thickness of each disk, and an offset above and below the plane
 of the AGN disk where each disk begins.
@@ -381,7 +387,44 @@ energy density; changes in kinetic energy density will depend on the velocity
 of the ambient gas. Temperature will also change but should always increase
 with kinetic jet feedback.
 
-Magnetic feedback is injected following  ([Hui 2006](doi.org/10.1086/501499))
+#### Tracing jet material
+
+Material launched by the kinetic jet component can be traced by setting
+```
+<problem/cluster/agn_feedback>
+enable_tracer = true   # disabled by default
+```
+
+At the moment, this also requires to enable a single passive scalar by setting
+```
+<hydro>
+nscalars = 1
+```
+This may change in the future as the current implementation restricts the
+passive scalar use to tracing jet material.
+This is a design decision motivated by simplicity and not for technical
+reasons (KISS!).
+
+Note that all material launched from within the jet region is traced, i.e.,
+passive scalar concentration does not differentiate between original cell
+material and mass added through the kinetic jet feedback mechanism.
+
+
+### Magnetic feedback
+
+Multiple options exists to inject magnetic fields:
+- a simple field loop (donut) configuration (better suited for kinetically dominated jets at larger scales)
+- a more complex pinching tower model (particularly suited for magnetically dominated jets at small scales)
+
+The option is controlled by the following parameter
+```
+<problem/cluster/magnetic_tower>
+potential_type = li # alternative: "donut"
+```
+
+#### Pinching magnetic tower
+
+Magnetic feedback is injected following  ([Li 2006](doi.org/10.1086/501499))
 where the injected magnetic field follows 
 
 $$
@@ -405,7 +448,8 @@ $$
 The parameters $\alpha$ and $\ell$ can be changed with
 ```
 <problem/cluster/magnetic_tower>
-alpha = 20
+potential_type = li
+li_alpha = 20
 l_scale = 0.001
 ```
 When injected as a fraction of 
@@ -429,9 +473,50 @@ so that the total mass injected matches the accreted mass propotioned to magneti
 l_mass_scale = 0.001
 ```
 
-A magnetic tower can also be inserted at runtime and injected at a fixed
+Mass injection by the tower is enabled by default.
+It can be disabled by setting
+```
+<problem/cluster/agn_feedback>
+enable_magnetic_tower_mass_injection = false
+```
+In this case, the injected mass through kinetic and thermal feedback
+according to their ratio.
+
+#### Simple field loop (donut) feedback
+
+Magnetic energy is injected according to the following simple potential
+$$
+A_h(r, \theta, h) = B_0 L \exp^\left ( -r^2/L^2 \right)$ for $h_\mathrm{offset} \leq |h| \leq h_\mathrm{offset} + h_\mathrm{thickness}
+$$
+resultig in a magnetic field configuration of
+$$
+B_\theta(r, \theta, h) = 2 B_0 r /L \exp^\left ( -r^2/L^2 \right)$ for $h_\mathrm{offset} \leq |h| \leq h_\mathrm{offset} + h_\mathrm{thickness}
+$$
+with all other components being zero.
+
+
+```
+<problem/cluster/magnetic_tower>
+potential_type = donut
+l_scale = 0.0005         # in code length
+donut_offset = 0.001     # in code length
+donut_thickness = 0.0001 # in code length
+```
+
+It is recommended to match the donut thickness to the thickness of the kinetic jet launching region
+and to choose a lengthscale that is half the jet launching radius.
+This results in almost all injected fields being confined to the launching region (and thus being
+carried along with a dominant jet).
+
+Mass feedback for the donut model is currently identical to the Li tower above
+(and, thus, the parameters above pertaining to the mass injection equally apply here).
+
+#### Fixed magnetic feedback
+
+Magnetic feedback can also be inserted at runtime and injected at a fixed
 increase in magnetic field, and additional mass can be injected at a fixed
 rate.
+This works for all magnetic vector potentials.
 ```
 <problem/cluster/magnetic_tower>
 initial_field = 0.12431560000204142
@@ -454,3 +539,27 @@ where `power_per_bcg_mass` and `mass_rate_per_bcg_mass` is the power and mass
 per time respectively injected per BCG mass at a given radius. This SNIA
 feedback is otherwise fixed in time, spherically symmetric, and dependant on
 the BCG specified in `<problem/cluster/gravity>`.
+
+## Stellar feedback
+
+Cold, dense, potentially magnetically supported gas may accumulate in a small
+disk-like structure around the AGN depending on the specific setup.
+In reality, this gas would form stars.
+
+In absence of star particles (or similar) in the current setup, we use a simple
+formulation to convert some fraction of this gas to thermal energy.
+More specifically, gas within a given radius, above a critical number density
+and below a temperature threshold, will be converted to thermal energy with a
+given efficiency.
+The parameters can be set as follows (numerical values here just for illustration purpose -- by default they are all 0, i.e., stellar feedback is disabled).
+
+```
+<problem/cluster/stellar_feedback>
+stellar_radius = 0.025 # in code length
+efficiency = 5e-6
+number_density_threshold = 2.93799894e+74 # in code_length**-3
+temperature_threshold = 2e4 # in K
+```
+
+Note that all parameters need to be specified explicitly for the feedback to work
+(i.e., no hidden default values).

docs/pgen.md

@@ -1,11 +1,73 @@
-# Problem specific callbacks
+# Problem generators
 
-## Source functions
+Different simulation setups in AthenaPK are controlled via the so-called problem generators.
+New problem generators can easily be added and we are happy to accept and merge contibuted problem
+generators by any user via pull requests.
+
+## Addding a new problem generator
+
+In general, four small steps are required:
+
+### 1. Add a new source file to the `src/pgen` folder
+
+The file shoud includes at least the `void ProblemGenerator(Mesh *pmesh, ParameterInput *pin, MeshData<Real> *md)`
+function, which is used to initialize the data at the beginning of the simulation.
+Alternatively, the `MeshBlock` version (`void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin)`) can be used that
+operates on a single block at a time rather than on a collection of blocks -- though this is not recommended from a performance point of view.
+
+The function (and all other functions in that file) should be encapsulated in their own namespace (that, ideally, is named
+identical to the file itself) so that the functions name are unique across different problem generators.
+
+Tip: Do not write the problem generator file from scratch but mix and match from existing ones,
+e.g., start with a simple one like [orszag_tang.cpp](../src/pgen/orszag_tang.cpp).
+
+### 2. Add new function(s) to `pgen.hpp`
+
+All callback functions, i.e., at least the `ProblemGenerator` plus additional optional ones (see step 5 below),
+need to be added to [pgen.hpp](../src/pgen/pgen.hpp).
+Again, just follow the existing pattern in the file and add the new function declarations with the appropriate namespace.
+
+### 3. Add callbacks to `main.cpp`
+
+All problem specific callback functions need to be enrolled in the [`main.cpp`](../src/main.cpp) file.
+The selection (via the input file) is controlled by the `problem_id` string in the `<job>` block.
+Again, for consistency it is recommended to pick a string that matches the namespace and problem generator file.
+
+### 4. Ensure new problem generator is compiled
+
+Add the new source file to [src/pgen/CMakeLists.txt](../src/pgen/CMakeLists.txt) so that it will be compiled
+along all other problem generators.
+
+### 5. (Optional) Add more additional callback
+
+In addition to the `ProblemGenerator` that initializes the data, other callback functions exists
+that allow to modify the behavior of AthenaPK on a problem specific basis.
+See [Callback functions](#Callback-functions)] below for available options.
+
+### 6. (Optional but most likely required) Write an input file
+
+In theory, one can hardcode all paramters in the source file (like in the
+[orszag_tang.cpp](../src/pgen/orszag_tang.cpp) problem generator) but it
+prevents modification of the problem setup once the binary is compiled.
+
+The more common usecase is to create an input file that contains a problem specific
+input block.
+The convention here is to have a block named `<problem/NAME>` where `NAME` is the name
+of the problem generator (or namespace used).
+For example, the Sod shock tube problem generator processes the input file with lines like
+```
+  Real rho_l = pin->GetOrAddReal("problem/sod", "rho_l", 1.0);
+```
+to set the density on the left hand side of the domain.
+
+## Callback functions
+
+### Source functions
 
 Additional (physical) source terms (e.g., the ones typically on the right hand side of
 system of equations) can be added in various ways from an algorithmic point of view.
 
-### Unsplit
+#### Unsplit
 
 Unsplit sources are added at each stage of the (multi-stage) integration after the
 (conserved) fluxes are calculated.
@@ -23,7 +85,7 @@ by implementing a function with that signature and assigning it to the
 
 Note, there is no requirement to call the function `MyUnsplitSource`.
 
-### Split first order (generally not recommended)
+#### Split first order (generally not recommended)
 
 If for special circumstances sources need to be added in a fully operator split way,
 i.e., the source function is only called *after* the full hydro (or MHD) integration,
@@ -38,12 +100,27 @@ Note, as the name suggests, this method is only first order accurate (while ther
 is no requirement to call the function `MyFirstOrderSource`).
 
 
-### Split second order (Strang splitting)
+#### Split second order (Strang splitting)
 
-Not implemented yet.
+Strang splitting achieves second order by interleaving the main hydro/MHD update
+with a source update.
+In practice, AthenaPK calls Strang split source with `0.5*dt` before the first stage of
+each integrator and with `0.5*dt` after the last stage of the integrator.
+Note that for consistency, a Strang split source terms should update both conserved
+and primitive variables in all zones (i.e., also in the ghost zones as those
+are currently not updated after calling a source function).
+
+Users can enroll a custom source function
+```c++
+void MyStrangSplitSource(MeshData<Real> *md, const Real beta_dt);
+```
+by implementing a function with that signature and assigning it to the
+`Hydro::ProblemSourceStrangSplit` callback (currently in `main.cpp`).
+
+Note, there is no requirement to call the function `MyStrangSplitSource`.
 
 
-## Timestep restrictions
+### Timestep restrictions
 
 If additional problem specific physics are implemented (e.g., through the source
 functions above) that require a custom timestep, it can be added via the
@@ -55,4 +132,60 @@ Real ProblemEstimateTimestep(MeshData<Real> *md);
 The return value is expected to be the minimum value over all blocks in the
 contained in the `MeshData` container, cf., the hydro/mhd `EstimateTimestep` function.
 Note, that the hyperbolic CFL value is currently applied after the function call, i.e.,
-it is also applied to the problem specific function.
+it is also applied to the problem specific function.
+
+### Additional initialization on startup (adding variables/fields, custom history output, ...)
+
+Sometimes a problem generator requires a more general processing of the input file
+and/or needs to make certain global adjustments (e.g., adding a custom callback function
+for the history output or adding a new field).
+This can be achieved via modifying the AthenaPK "package" (currently called
+`parthenon::StateDescriptor`) through the following function.
+```c++
+void ProblemInitPackageData(ParameterInput *pin, parthenon::StateDescriptor *pkg)
+```
+
+For example, the [src/pgen/turbulence.cpp](../[src/pgen/turbulence.cpp]) problem generator
+add an additional field (here for an acceleration field) by calling
+```c++
+  Metadata m({Metadata::Cell, Metadata::Derived, Metadata::OneCopy},
+             std::vector<int>({3}));
+  pkg->AddField("acc", m);
+```
+in the `ProblemInitPackageData`.
+
+### Additional initialization on mesh creation/remeshing/load balancing
+
+For some problem generators it is required to initialize data on "new" blocks.
+These new blocks can, for example, be created during mesh refinement
+(or derefinement) and this data is typically not active data (like
+conserved or primitive variables as those are handled automatically)
+but more general data.
+Once example is the phase information in the turbulence driver that
+does not vary over time but spatially (and therefore at a per block level).
+
+The appropriate callback to enroll is
+```c++
+void InitMeshBlockUserData(MeshBlock *pmb, ParameterInput *pin)
+```
+
+### UserWorkAfterLoop
+
+If additional work is required once the main loop finishes (i.e., once the
+simulation reaches its final time) computations can be done inside the
+```c++
+void UserWorkAfterLoop(Mesh *mesh, ParameterInput *pin, parthenon::SimTime &tm)
+```
+callback function.
+This is, for example, done in the linear wave problem generator to calculate the
+error norms for convergence tests.
+
+### MeshBlockUserWorkBeforeOutput
+
+Sometimes it is desirable to further process data before an output file is written.
+For this the
+```c++
+void UserWorkBeforeOutput(MeshBlock *pmb, ParameterInput *pin)
+```
+callback is available, that is called once every time right before a data output
+(hdf5 or restart) is being written.

inputs/cluster/cluster.in

@@ -181,8 +181,8 @@ kinetic_jet_thickness  = 0.05
 kinetic_jet_offset = 0.05
 
 <problem/cluster/magnetic_tower>
-
-alpha = 20
+potential_type = li
+li_alpha = 20
 l_scale = 0.001
 initial_field = 0.12431560000204142
 #NOTE: Change to this line to disable initial magnetic tower

inputs/cluster/magnetic_tower.in

@@ -109,7 +109,8 @@ kinetic_fraction = 0
 thermal_fraction = 0
 
 <problem/cluster/magnetic_tower>
-alpha = 20
+potential_type = li
+li_alpha = 20
 l_scale = 0.01
 initial_field = 0.12431560000204142
 fixed_field_rate = 12.431560000204144

src/CMakeLists.txt

@@ -12,6 +12,7 @@ add_executable(
         hydro/hydro_driver.cpp
         hydro/hydro.cpp
         hydro/glmmhd/dedner_source.cpp
+        hydro/prolongation/custom_ops.hpp
         hydro/srcterms/gravitational_field.hpp
         hydro/srcterms/tabular_cooling.hpp
         hydro/srcterms/tabular_cooling.cpp