Reject HEAT_SOURCE without a radiation model and fix its 3D ellipsoid test

.github/workflows/regression.yml

@@ -474,6 +474,15 @@ jobs:
         with:
           entrypoint: python3
           args: src/SU2/UnitTests/SU2_CFD/output/check_nemo_convergence_migration.py /github/workspace/install/bin/test_driver
+      - name: Check CConfig option errors
+        if: matrix.testdriver == 'test_driver'
+        uses: docker://ghcr.io/su2code/su2/test-su2:261011-1107
+        env:
+          OMPI_MCA_osc: pt2pt
+          LD_LIBRARY_PATH: /github/workspace/install/lib
+        with:
+          entrypoint: python3
+          args: src/SU2/UnitTests/Common/check_config_errors.py /github/workspace/install/bin/test_driver
       - name: Post Cleanup
         uses: docker://ghcr.io/su2code/su2/test-su2:261011-1107
         with:

Common/include/toolboxes/geometry_toolbox.hpp

@@ -250,5 +250,27 @@ inline void TangentProjection(Int nDim, const Mat& tensor, const Scalar* vector,
 
   for (Int iDim = 0; iDim < nDim; iDim++) proj[iDim] -= normalProj * vector[iDim];
 }
+/*!
+ * \brief Check if a point lies inside an ellipsoid (ellipse in 2D) rotated about the z axis.
+ * \note Intermediate values are stored with an explicit type, because with AD types "auto" would hold
+ *       expression templates.
+ * \param[in] nDim - Number of dimensions, only the first nDim coordinates and semi-axes are used.
+ * \param[in] coord - Coordinates of the point.
+ * \param[in] center - Center of the ellipsoid.
+ * \param[in] axes - Semi-axes of the ellipsoid.
+ * \param[in] rotZ - Rotation of the ellipsoid about the z axis, in radians.
+ * \return True if the point is inside the ellipsoid or on its surface.
+ */
+template <class T, typename Int>
+inline bool PointInEllipsoid(Int nDim, const T* coord, const T* center, const T* axes, const T& rotZ) {
+  const T dx = coord[0] - center[0];
+  const T dy = coord[1] - center[1];
+  const T x = dx * cos(rotZ) + dy * sin(rotZ);
+  const T y = -dx * sin(rotZ) + dy * cos(rotZ);
+  T check = pow(x / axes[0], 2) + pow(y / axes[1], 2);
+  if (nDim == 3) check += pow((coord[2] - center[2]) / axes[2], 2);
+  return check <= 1;
+}
+
 /// @}
 }  // namespace GeometryToolbox

Common/src/CConfig.cpp

@@ -4231,6 +4231,12 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i
 
   Radiation = (Kind_Radiation != RADIATION_MODEL::NONE);
 
+  /*--- The ellipsoidal HEAT_SOURCE is part of the radiation model. Without one, no solver applies it. ---*/
+  if (HeatSource && !Radiation) {
+    SU2_MPI::Error("HEAT_SOURCE= YES requires a RADIATION_MODEL. For other volumetric sources, use PYTHON_CUSTOM_SOURCE.",
+                   CURRENT_FUNCTION);
+  }
+
   /*--- Check for unsupported features. ---*/
 
   if ((Kind_Solver != MAIN_SOLVER::EULER && Kind_Solver != MAIN_SOLVER::NAVIER_STOKES && Kind_Solver != MAIN_SOLVER::RANS) && (TimeMarching == TIME_MARCHING::HARMONIC_BALANCE)){

SU2_CFD/src/solvers/CRadSolver.cpp

@@ -46,30 +46,16 @@ CRadSolver::CRadSolver(CGeometry* geometry, CConfig *config) : CSolver() {
 void CRadSolver::SetVolumetricHeatSource(CGeometry *geometry, CConfig *config) {
   SU2_ZONE_SCOPED
 
-  unsigned long iPoint;
-  unsigned short iDim;
-
-  su2double CP[3]={0.0,0.0,0.0};
-  su2double alpha = config->GetHeatSource_Rot_Z() * PI_NUMBER/180.0;
-  su2double OP_rot[3]={0.0,0.0,0.0};
-  const su2double *OP;
+  const su2double alpha = config->GetHeatSource_Rot_Z() * PI_NUMBER/180.0;
   const su2double *OC = config->GetHeatSource_Center();
   const su2double *Axes = config->GetHeatSource_Axes();
-  su2double check;
+
   // Reset the boolean for all points
   nodes->ResetVol_HeatSource();
   // Loop over all points and determine whether they are inside
-  for (iPoint = 0; iPoint < nPoint; iPoint ++) {
-    check = 0;
-    OP = geometry->nodes->GetCoord(iPoint);
-    // Reference point with respect to center of the ellipse
-    for (iDim = 0; iDim < nDim; iDim++) CP[iDim] = OP[iDim]-OC[iDim];
-    // Rotate point with respect to Z axis
-    OP_rot[0] = OC[0] + CP[0]*cos(alpha) + CP[1]*sin(alpha);
-    OP_rot[1] = OC[1] - CP[0]*sin(alpha) + CP[1]*cos(alpha);
-    // Check if rotated point is inside the ellipse
-    for (iDim = 0; iDim < nDim; iDim++) check += pow(OP_rot[iDim]-OC[iDim],2.0)/pow(Axes[iDim], 2.0);
-    if (check <=1) nodes->SetVol_HeatSource(iPoint);
+  for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++) {
+    const auto *OP = geometry->nodes->GetCoord(iPoint);
+    if (GeometryToolbox::PointInEllipsoid(nDim, OP, OC, Axes, alpha)) nodes->SetVol_HeatSource(iPoint);
   }
 
 }

TestCases/parallel_regression.py

@@ -403,12 +403,12 @@ def main():
     test_list.append(cylinder_lowmach)
 
     # Laminar NACA0012 with linelet preconditioner (see #2877), 4 ranks so that some partitions have no linelets
-    lam_naca0012_linelet           = TestCase('lam_naca0012_linelet')
-    lam_naca0012_linelet.cfg_dir   = "navierstokes/naca0012"
-    lam_naca0012_linelet.cfg_file  = "lam_NACA0012_LINELET.cfg"
+    lam_naca0012_linelet = TestCase('lam_naca0012_linelet')
+    lam_naca0012_linelet.cfg_dir = "navierstokes/naca0012"
+    lam_naca0012_linelet.cfg_file = "lam_NACA0012_LINELET.cfg"
     lam_naca0012_linelet.test_iter = 20
     lam_naca0012_linelet.test_vals = [-5.335771, -2.936187, -2.943926, 0.102449]
-    lam_naca0012_linelet.command   = TestCase.Command("mpirun --oversubscribe -n 4", "SU2_CFD")
+    lam_naca0012_linelet.command = TestCase.Command("mpirun --oversubscribe -n 4", "SU2_CFD")
     test_list.append(lam_naca0012_linelet)
 
     # 2D Poiseuille flow (body force driven with periodic inlet / outlet)
@@ -1352,12 +1352,12 @@ def main():
     test_list.append(transonic_stator_restart)
 
     # Multiple turbomachinery interface restart
-    multi_interface                    = TestCase('multi_interface')
-    multi_interface.cfg_dir            = "turbomachinery/multi_interface"
-    multi_interface.cfg_file           = "multi_interface_rst.cfg"
-    multi_interface.test_iter          = 5
-    multi_interface.test_vals          = [-8.634558, -8.895554, -9.348754]
-    multi_interface.test_vals_aarch64  = [-8.632227, -8.894736, -9.348706]
+    multi_interface = TestCase('multi_interface')
+    multi_interface.cfg_dir = "turbomachinery/multi_interface"
+    multi_interface.cfg_file = "multi_interface_rst.cfg"
+    multi_interface.test_iter = 5
+    multi_interface.test_vals = [-8.634558, -8.895554, -9.348754]
+    multi_interface.test_vals_aarch64 = [-8.632227, -8.894736, -9.348706]
     test_list.append(multi_interface)
 
     ######################################
@@ -1756,6 +1756,16 @@ def main():
     pywrapper_Unst_Heat_Source.command = TestCase.Command("mpirun -n 2", "python", "run.py")
     test_list.append(pywrapper_Unst_Heat_Source)
 
+    # Heat solver steady with a volumetric source in an ellipsoid, analytic solution
+    pywrapper_Steady_Heat_Source = TestCase('pywrapper_Steady_Heat_Source')
+    pywrapper_Steady_Heat_Source.cfg_dir = "py_wrapper/custom_heat_source_steady"
+    pywrapper_Steady_Heat_Source.cfg_file = "run.py"
+    pywrapper_Steady_Heat_Source.test_iter = 15
+    # Only the total heat flux is compared, the residual history of this small case depends on the partitioning.
+    pywrapper_Steady_Heat_Source.test_vals = [-64.000000]
+    pywrapper_Steady_Heat_Source.command = TestCase.Command("mpirun -n 2", "python", "run.py")
+    test_list.append(pywrapper_Steady_Heat_Source)
+
     ##############################################
     ### Method of Manufactured Solutions (MMS) ###
     ##############################################

TestCases/py_wrapper/custom_heat_source_steady/run.py

@@ -0,0 +1,157 @@
+#!/usr/bin/env python
+
+## \file run.py
+#  \brief Steady solid heat conduction with a volumetric heat source in an ellipsoid.
+#  \version 8.5.0 "Harrier"
+#
+# SU2 Project Website: https://su2code.github.io
+#
+# The SU2 Project is maintained by the SU2 Foundation
+# (http://su2foundation.org)
+#
+# Copyright 2012-2026, SU2 Contributors (cf. AUTHORS.md)
+#
+# SU2 is free software; you can redistribute it and/or
+# modify it under the terms of the GNU Lesser General Public
+# License as published by the Free Software Foundation; either
+# version 2.1 of the License, or (at your option) any later version.
+#
+# SU2 is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+# Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public
+# License along with SU2. If not, see <http://www.gnu.org/licenses/>.
+
+from mpi4py import MPI
+import math
+import sys
+import pysu2  # imports the SU2 wrapped module
+
+# Slab of width L with both x faces at T_WALL and insulated y faces.
+L = 0.016
+H = 0.004
+N_X = 17
+T_WALL = 300.0
+DENSITY = 5000.0
+CP = 500.0
+CONDUCTIVITY = 15.0
+
+# Volumetric heat source (W/m^3) inside an ellipsoid that can be rotated about the z axis.
+# The ellipsoid encloses the whole slab here, so the source is uniform.
+SOURCE_VAL = 1e6
+SOURCE_CENTER = (0.5 * L, 0.5 * H, 0.0)
+SOURCE_AXES = (0.02, 0.02, 1.0)
+SOURCE_ROTATION_Z = 0.0  # degrees
+
+settings = f"""
+SOLVER= HEAT_EQUATION
+INC_NONDIM= DIMENSIONAL
+
+FREESTREAM_TEMPERATURE= {T_WALL}
+MATERIAL_DENSITY= {DENSITY}
+SPECIFIC_HEAT_CP= {CP}
+THERMAL_CONDUCTIVITY_CONSTANT= {CONDUCTIVITY}
+
+MARKER_ISOTHERMAL= ( x_minus, {T_WALL}, x_plus, {T_WALL} )
+MARKER_HEATFLUX= ( y_minus, 0, y_plus, 0 )
+MARKER_MONITORING= ( x_minus, x_plus )
+
+PYTHON_CUSTOM_SOURCE= YES
+
+NUM_METHOD_GRAD= GREEN_GAUSS
+CFL_NUMBER= 1e8
+TIME_DISCRE_HEAT= EULER_IMPLICIT
+
+LINEAR_SOLVER= FGMRES
+LINEAR_SOLVER_PREC= ILU
+LINEAR_SOLVER_ERROR= 1E-15
+LINEAR_SOLVER_ITER= 10
+
+MESH_FORMAT= RECTANGLE
+MESH_BOX_SIZE= ( {N_X}, 5, 0 )
+MESH_BOX_LENGTH= ( {L}, {H}, 0 )
+
+ITER= 20
+CONV_FIELD= RMS_TEMPERATURE
+CONV_RESIDUAL_MINVAL= -20
+
+SCREEN_OUTPUT= INNER_ITER, RMS_TEMPERATURE, LINSOL_ITER, TOTAL_HEATFLUX
+HISTORY_OUTPUT= ITER, RMS_RES, HEAT, LINSOL
+OUTPUT_FILES= RESTART, PARAVIEW
+"""
+
+
+def HeatSource(coord):
+    """Returns the source (W/m^3) at a point, SOURCE_VAL inside the ellipsoid and zero outside."""
+    alpha = math.radians(SOURCE_ROTATION_Z)
+    dx = coord[0] - SOURCE_CENTER[0]
+    dy = coord[1] - SOURCE_CENTER[1]
+    x = dx * math.cos(alpha) + dy * math.sin(alpha)
+    y = -dx * math.sin(alpha) + dy * math.cos(alpha)
+    check = (x / SOURCE_AXES[0])**2 + (y / SOURCE_AXES[1])**2
+    if len(coord) == 3:
+        check += ((coord[2] - SOURCE_CENTER[2]) / SOURCE_AXES[2])**2
+    return SOURCE_VAL if check <= 1 else 0.0
+
+
+def main():
+    comm = MPI.COMM_WORLD
+    rank = comm.Get_rank()
+
+    if rank == 0:
+        with open('config.cfg', 'w') as f:
+            f.write(settings)
+    comm.Barrier()
+
+    # Initialize the corresponding driver of SU2, this includes solver preprocessing.
+    try:
+        driver = pysu2.CSinglezoneDriver('config.cfg', 1, comm)
+    except TypeError as exception:
+        print('A TypeError occured in pysu2.CDriver : ', exception)
+        raise
+
+    iHEATSOLVER = driver.GetSolverIndices()['HEAT']
+    nNode = driver.GetNumberNodes() - driver.GetNumberHaloNodes()
+    Coords = driver.Coordinates()
+
+    # The heat solver works with the temperature equation divided by rho*cp,
+    # so a source in W/m^3 has to be divided by rho*cp as well.
+    Source = driver.UserDefinedSource(iHEATSOLVER)
+    for iPoint in range(nNode):
+        Source.Set(iPoint, 0, HeatSource(Coords.Get(iPoint)) / (DENSITY * CP))
+
+    if rank == 0:
+        print("\n------------------------------ Begin Solver -----------------------------\n")
+    sys.stdout.flush()
+
+    driver.StartSolver()
+
+    # Compare with the exact solution T(x) = T_WALL + q*x*(L-x)/(2*k).
+    Solution = driver.Solution(iHEATSOLVER)
+    max_error = 0.0
+    for iPoint in range(nNode):
+        x = Coords.Get(iPoint)[0]
+        exact = T_WALL + SOURCE_VAL * x * (L - x) / (2 * CONDUCTIVITY)
+        max_error = max(max_error, abs(Solution.Get(iPoint)[0] - exact))
+    max_error = comm.allreduce(max_error, op=MPI.MAX)
+
+    # All the heat added by the source has to leave through the isothermal walls.
+    heat_flux = driver.GetOutputValue('TOTAL_HEATFLUX')
+
+    driver.Finalize()
+
+    # The isothermal wall is imposed weakly, which shifts the solution by q*h^2/(2*k).
+    h = L / (N_X - 1)
+    tolerance = 1.01 * SOURCE_VAL * h**2 / (2 * CONDUCTIVITY)
+    if rank == 0:
+        print(f"Maximum temperature error: {max_error:.6f} K (tolerance {tolerance:.6f} K)")
+        print(f"Total heat flux: {heat_flux:.6f} W (exact {-SOURCE_VAL * L * H:.6f} W)")
+
+    assert max_error < tolerance, f"Test FAILED, maximum temperature error = {max_error}"
+    assert abs(heat_flux / (-SOURCE_VAL * L * H) - 1) < 1e-6, f"Test FAILED, total heat flux = {heat_flux}"
+
+
+if __name__ == '__main__':
+    main()

UnitTests/Common/CConfig_tests.cpp

@@ -72,8 +72,56 @@ const std::string mixture_options =
     "PRANDTL_TURB= 0.9, 0.9\n"
     "DIFFUSIVITY_CONSTANT= 1.0e-05\n";
 
+/*--- Solid heat conduction on a rectangle mesh, used by the HEAT_SOURCE checks. ---*/
+const std::string heat_equation_options =
+    "SOLVER= HEAT_EQUATION\n"
+    "INC_NONDIM= DIMENSIONAL\n"
+    "MESH_FORMAT= RECTANGLE\n"
+    "MESH_BOX_SIZE= 5,5,0\n"
+    "MESH_BOX_LENGTH= 1,1,0\n"
+    "MARKER_ISOTHERMAL= (x_minus, 300.0, x_plus, 300.0)\n"
+    "MARKER_HEATFLUX= (y_minus, 0.0, y_plus, 0.0)\n";
+
+const std::string heat_source_options =
+    "HEAT_SOURCE= YES\n"
+    "HEAT_SOURCE_VAL= 1e6\n";
+
+/*--- Builds a config quietly, so that the test output only shows the result. ---*/
+bool GetHeatSource(const std::string& options) {
+  std::stringstream ss(options);
+  auto* orig_buf = std::cout.rdbuf(nullptr);
+  CConfig config(ss, SU2_COMPONENT::SU2_CFD, false);
+  std::cout.rdbuf(orig_buf);
+  return config.GetHeatSource();
+}
+
+/*--- Builds a config with the output left on, so that the message of SU2_MPI::Error reaches the runner. ---*/
+void BuildConfig(const std::string& options) {
+  std::stringstream ss(options);
+  CConfig config(ss, SU2_COMPONENT::SU2_CFD, false);
+}
+
 }  // namespace
 
+TEST_CASE("HEAT_SOURCE is accepted with a radiation model", "[Config]") {
+  CHECK(GetHeatSource(base_options + ideal_gas_options + heat_source_options + "RADIATION_MODEL= P1\n"));
+  CHECK_FALSE(GetHeatSource(heat_equation_options));
+}
+
+/*--- Hidden tests, SU2_MPI::Error aborts the process. check_config_errors.py runs each one on its own and checks
+ * that it fails with the expected message. ---*/
+
+TEST_CASE("HEAT_SOURCE is rejected for the heat equation", "[.HeatSourceHeatEquation]") {
+  BuildConfig(heat_equation_options + heat_source_options);
+  /*--- Returning is a successful exit, which check_config_errors.py treats as a failure. ---*/
+  SUCCEED("HEAT_SOURCE without a radiation model was not rejected");
+}
+
+TEST_CASE("HEAT_SOURCE is rejected for the incompressible solver without radiation", "[.HeatSourceIncompressible]") {
+  BuildConfig(base_options + ideal_gas_options + heat_source_options);
+  SUCCEED("HEAT_SOURCE without a radiation model was not rejected");
+}
+
 TEST_CASE("INIT_OPTION_INC defaults", "[Config]") {
   /*--- Without the option each fluid model keeps the density it used before INIT_OPTION_INC existed. ---*/