Mean opacity for variable kappa based on piecewise-power-law approxim…

…ation (#626)

### Description

This PR implements piecewise-power-law approximation for multigroup
radiation transfer. We calculate Planck-mean, energy-mean, and flux-mean
opacity given an arbitrary function kappa(nu, T, rho) based on
piecewise-power-law fitting to kappa(nu), E(nu), and F(nu). The user can
enable this routine by setting `OpacityModel` to
`OpacityModel::piecewisePowerLaw`. See examples in
test_radhydro_shock_multigroup.cpp and test_radhydro_pulse_MG_int.hpp.
The original approach that the opacities are defined via
`ComputePlanckOpacity`, `ComputeFluxMeanOpacity`, etc are retained by
setting `OpacityModel` to `OpacityModel::user`, which is the default.

We calculate the power-law slope of a radiation quantity via the
following relations:

$$
\alpha_{i} = {\rm minmod}(s_{i-1/2}, s_{i+1/2})
$$

where

$$
s_{i+1/2} = \frac{\ln \frac{X_{i+1}}{\nu_{i+3/2} - \nu_{i+1/2}} - \ln
\frac{X_i}{\nu_{i+1/2} - \nu_{i-1/2}}}{\ln \sqrt{\nu_{i+1/2}
\nu_{i+3/2}} - \ln \sqrt{\nu_{i-1/2} \nu_{i+1/2}}}
$$

I will try to construct a better estimate of the power-law slope in a
separate PR. The user can choose to overwrite this by defining your own
`RadSystem<problem_t>::ComputeRadQuantityExponents`; see example in
test_radhydro_pulse_MG_int.cpp.

The power-law fitting of $\kappa(\nu)$ can be specified in two ways. The
first method is to specify its power-law slope and lower-bound value for
every photon bin, given $T$ and $\rho$. The second method is to specify
a precise definition of $\kappa(\nu, T, \rho)$. Then, a power-law
fitting to $\kappa(\nu)$ is done on the fly in the following way:

$$
\begin{gather}
\kappa_{L}^i = \kappa(\nu_{i-1/2}, T, \rho) \\
\alpha^i_{\kappa} = \frac{\ln(\kappa(\nu_{i+1/2}, T, \rho) /
\kappa(\nu_{i-1/2},T,\rho))}{\ln(\nu_{i+1/2} / \nu_{i-1/2})} \\
\kappa_{{\rm fit}}^i(\nu) = \kappa_{L}^i \left( \frac{\nu}{\nu_{i-1/2}}
\right)^{\alpha_{\kappa}^i}
\end{gather}
$$

This is definitely not the most accurate way to do this as it is not
volume conservative, but it's simple and this is important since it's
done for every cell. An alternative is to construct $\kappa_{L,i}$ and
$\alpha_{\kappa,i}$ in the beginning if $\kappa$ does not depend on $T$
or $\rho$.

### Related issues

N/A

### Checklist
_Before this pull request can be reviewed, all of these tasks should be
completed. Denote completed tasks with an `x` inside the square brackets
`[ ]` in the Markdown source below:_
- [x] I have added a description (see above).
- [x] I have added a link to any related issues see (see above).
- [x] I have read the [Contributing
Guide](https://github.com/quokka-astro/quokka/blob/development/CONTRIBUTING.md).
- [x] I have added tests for any new physics that this PR adds to the
code.
- [x] I have tested this PR on my local computer and all tests pass.
- [x] I have manually triggered the GPU tests with the magic comment
`/azp run`.
- [x] I have requested a reviewer for this PR.

---------

Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>

src/CMakeLists.txt

@@ -177,6 +177,7 @@ add_subdirectory(RadhydroPulse)
 add_subdirectory(RadhydroPulseDyn)
 add_subdirectory(RadhydroPulseGrey)
 add_subdirectory(RadhydroPulseMG)
+add_subdirectory(RadhydroPulseMGint)
 
 add_subdirectory(BinaryOrbitCIC)
 add_subdirectory(Cooling)

src/RadTube/test_radiation_tube.cpp

@@ -65,6 +65,7 @@ template <> struct RadSystem_Traits<TubeProblem> {
 	static constexpr double energy_unit = C::k_B;
 	static constexpr amrex::GpuArray<double, Physics_Traits<TubeProblem>::nGroups + 1> radBoundaries{0., 3.3 * T0, inf}; // Kelvin
 	static constexpr int beta_order = 1;
+	static constexpr OpacityModel opacity_model = OpacityModel::user;
 };
 
 template <>

src/RadhydroPulseMGint/CMakeLists.txt

@@ -0,0 +1,15 @@
+if (AMReX_SPACEDIM EQUAL 1)
+  add_executable(test_radhydro_pulse_MG_int test_radhydro_pulse_MG_int.cpp ../fextract.cpp ${QuokkaObjSources})
+
+  if(AMReX_GPU_BACKEND MATCHES "CUDA")
+      setup_target_for_cuda_compilation(test_radhydro_pulse_MG_int)
+  endif(AMReX_GPU_BACKEND MATCHES "CUDA")
+
+	# define variable WORK_PATH = {WORKSPACE}/sims
+	# set(WORK_PATH "${CMAKE_SOURCE_DIR}/tests/20240408-MG-int-v8.2")
+
+	# mkdir -p WORK_PATH
+	# file(MAKE_DIRECTORY ${WORK_PATH})
+
+  add_test(NAME RadhydroPulseMGint COMMAND test_radhydro_pulse_MG_int RadhydroPulse.in WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}/tests)
+endif()