Skip to content

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

Open
pragyaangaur wants to merge 10 commits into
su2code:developfrom
pragyaangaur:heat-solver-volumetric-source
Open

pragyaangaur wants to merge 10 commits into
su2code:developfrom
pragyaangaur:heat-solver-volumetric-source

Conversation

@pragyaangaur

@pragyaangaur pragyaangaur commented Sep 29, 2026 •

Copy link
Copy Markdown

Proposed Changes

HEAT_SOURCE is only applied by the P1 radiation model. Without a radiation model the option is accepted and then ignored. A solid zone run with SOLVER= HEAT_EQUATION and HEAT_SOURCE= YES converges and exits normally, but the solid is never heated. In the incompressible solver the same setting would read a radiation solver that does not exist. This PR makes that case an error and fixes the ellipsoid test that the P1 source uses in 3D.

The changes are listed below.

  • CConfig stops with an error when HEAT_SOURCE= YES is set without a RADIATION_MODEL. The message points to PYTHON_CUSTOM_SOURCE for other volumetric sources.
  • The ellipsoid test moves from CRadSolver::SetVolumetricHeatSource into GeometryToolbox::PointInEllipsoid. The old code left the rotated z coordinate at zero, so in 3D the z term used the distance from the center to the plane z = 0 and ignored the z coordinate of the point. The helper measures z from the center. 2D results are unchanged.
  • UnitTests/Common/check_config_errors.py runs two hidden Catch2 tests in CConfig_tests.cpp, one for a HEAT_EQUATION solid and one for an incompressible case without radiation, and checks that each stops with the HEAT_SOURCE error. It is registered in meson and in the unit test job of the regression workflow, in the same way as the NEMO convergence field probes. A visible test checks that HEAT_SOURCE is still accepted with RADIATION_MODEL= P1.
  • A unit test for PointInEllipsoid covers the 2D, rotated and 3D cases. The 3D case fails with the old logic.
  • A new python wrapper test case, py_wrapper/custom_heat_source_steady, applies a volumetric source in an ellipsoid to a solid slab through PYTHON_CUSTOM_SOURCE. It divides the source in W/m^3 by rhocp, because the heat solver works with the temperature equation divided by rhocp. The script checks the result against the exact parabolic solution and the heat balance, and it is added to parallel_regression.py.
  • config_template.cfg states that the option needs a radiation model and gives the units of HEAT_SOURCE_VAL.

An earlier version of this PR also applied the source in CHeatSolver. That part was removed after review, and the new python wrapper test case replaces it.

Related Work

The py_wrapper/custom_heat_source test case shows the Python route for a volumetric source in the heat solver. No test case sets HEAT_SOURCE= YES without a radiation model, so the new check does not affect any existing regression test.

PR Checklist

  • I am submitting my contribution to the develop branch.
  • My contribution generates no new compiler warnings (try with --warnlevel=3 when using meson).
  • My contribution is commented and consistent with SU2 style (https://su2code.github.io/docs_v7/Style-Guide/).
  • I used the pre-commit hook to prevent dirty commits and used pre-commit run --all to format old commits.
  • I have added a test case that demonstrates my contribution, if necessary.
  • I have updated appropriate documentation (Tutorials, Docs Page, config_template.cpp), if necessary.

The default release build gives no new warnings in the changed files. I have not yet built with --warnlevel=3.

HEAT_SOURCE was only used by the incompressible energy equation through the P1 radiation solver. A solid zone run with SOLVER= HEAT_EQUATION read the option, converged and exited normally, but the source never entered the residual, so the solid was not heated.

CHeatSolver now marks the points inside the source ellipsoid when it is built and adds the source in Source_Residual, scaled by the same reference as a wall heat flux. It is applied to solid conduction only, because in the weakly coupled fluid case the source belongs to the flow solver's energy equation.

The ellipsoid test moves into GeometryToolbox::PointInEllipsoid and CRadSolver uses it too. The old radiation code left the rotated z coordinate at zero. In 3D its z term therefore used the distance from the center to the plane z = 0 and ignored the z coordinate of the point. The shared helper measures z from the center, which fixes 3D sources in the radiation solver as well.

A new regression case, solid_heat_conduction/volumetric_source, heats a slab with a uniform source and uses the built-in rectangle mesh, so it needs no mesh file. The computed profile has the exact parabolic shape, the heat leaving through the isothermal walls equals the source (64 W per unit depth), and the peak is within 0.033 K of the analytic 302.133 K. That offset comes from the weak isothermal wall and falls by four times when the mesh spacing is halved. A unit test covers the ellipsoid helper in 2D, rotated and 3D.
@bigfooted

Copy link
Copy Markdown
Contributor

What is the use of this volumetric heat source? It looks very specific and at the same time not very useful. Why not let the python wrapper deal with this? With pywrapper you could define any kind of source term that you want. I'd rather see this hardcoded source removed and replaced by a python wrapper option.
For P1 radiation it's different because it's a well known radiation model that users might just want to have available as an immediate option.

Following review, the hard-coded source in CHeatSolver and its regression case are removed. PYTHON_CUSTOM_SOURCE already provides volumetric sources for the heat solver and for other solvers.

HEAT_SOURCE is only applied by the P1 radiation model. Without a radiation model, a heat equation run accepted the option and ignored it, and the incompressible solver would read a radiation solver that does not exist. CConfig now stops with an error that points to PYTHON_CUSTOM_SOURCE.

PointInEllipsoid now stores its intermediate values with an explicit type. With AD types, "auto" held CoDiPack expression templates and the forward and reverse builds failed to compile.
@pragyaangaur pragyaangaur changed the title Apply HEAT_SOURCE in the heat equation solver Reject HEAT_SOURCE without a radiation model and fix its 3D ellipsoid test Sep 30, 2026
@pragyaangaur

Copy link
Copy Markdown
Author

Thanks, that makes sense. PYTHON_CUSTOM_SOURCE already covers this for the heat solver, so I have removed the hard-coded source from CHeatSolver along with its regression case.

What is left fixes the problem that led me here. With SOLVER= HEAT_EQUATION, HEAT_SOURCE= YES was accepted and then ignored, so the run converged with no heating and gave no warning. CConfig now stops with an error when HEAT_SOURCE is set without a RADIATION_MODEL, and the message points to PYTHON_CUSTOM_SOURCE. The same check also covers the incompressible solver, which would otherwise read a radiation solver that does not exist.

I also kept a fix for the P1 source itself. In 3D the old ellipsoid test ignored the z coordinate of the point, because the rotated z value was left at zero. The test is now in GeometryToolbox::PointInEllipsoid with a unit test, and 2D results do not change.

@bigfooted

Copy link
Copy Markdown
Contributor

Please write the python wrapper replacement and test it.

This is the python wrapper replacement for the hard-coded source that was removed from CHeatSolver. The script py_wrapper/custom_heat_source_steady/run.py applies a source in an ellipsoid to a solid slab through PYTHON_CUSTOM_SOURCE and checks the result against the exact parabolic solution.

The heat solver works with the temperature equation divided by rho*cp, so the script divides the source in W/m^3 by rho*cp before it sets it. With that scaling the heat leaving through the isothermal walls equals the source, 64 W per unit depth, and the nodal temperatures match the exact solution to within the offset of the weak isothermal wall, q*h^2/(2*k) = 0.0333 K. The script asserts both and is added to parallel_regression.py.
@pragyaangaur

Copy link
Copy Markdown
Author

Done. I added py_wrapper/custom_heat_source_steady/run.py. It applies a volumetric source in an ellipsoid (center, axes and rotation about z) to a solid slab through PYTHON_CUSTOM_SOURCE, and it replaces the hard-coded source that I removed from CHeatSolver.

The case is a slab with both x faces at 300 K, insulated y faces and a uniform source of 1e6 W/m^3. The script divides the source by rhocp before it sets it, because the heat solver works with the temperature equation divided by rhocp. It then asserts two things. The heat leaving through the isothermal walls equals the source, which is 64 W per unit depth. The nodal temperatures match the exact solution T(x) = 300 + q x (L - x) / (2 k) to within q h^2 / (2 k) = 0.0333 K, which is the offset of the weakly imposed isothermal wall.

I ran it on 1 and 2 ranks, and both pass. The test is added to parallel_regression.py as pywrapper_Steady_Heat_Source.

The residual of this small case at a fixed iteration depends on how the mesh is partitioned, so the stored values did not match on the CI machines. The run now uses a fixed number of iterations and the regression entry compares only the total heat flux, which is the same on 1 to 4 ranks. The script still asserts the temperature field against the exact solution.
@pragyaangaur

Copy link
Copy Markdown
Author

The regression failure was in the new test . I had stored a residual value at a fixed iteration, and for this small case that value depends on the partitioning. The entry now compares only the total heat flux, and the script still puts the temperature field against the exact solution.

@rois1995

Copy link
Copy Markdown
Contributor

Looks good. Please add a regression test for the new configuration check.

SU2_MPI::Error ends the process, so the rejected configurations are hidden Catch2 tests and UnitTests/Common/check_config_errors.py runs each one on its own. The runner passes only when the test driver exits with an error and its output contains the HEAT_SOURCE message. It covers a HEAT_EQUATION solid and an incompressible case without radiation, and it is registered in meson and in the unit test job of the regression workflow, in the same way as the NEMO convergence field probes.

A visible test checks that HEAT_SOURCE is still accepted with RADIATION_MODEL= P1, and that a heat equation case without HEAT_SOURCE builds normally.
@pragyaangaur

Copy link
Copy Markdown
Author

Thanks. I added a test for the new configuration check. SU2_MPI::Error ends the process, so the regression harness cannot check it, and I followed the pattern used for the NEMO convergence field probes instead. Two hidden Catch2 tests in CConfig_tests.cpp build a HEAT_EQUATION solid with HEAT_SOURCE and an incompressible case with HEAT_SOURCE but no radiation model. UnitTests/Common/check_config_errors.py runs each one in its own process and passes only when it stops with the HEAT_SOURCE message. The runner is registered in meson and in the unit test job of the regression workflow. A visible test checks that HEAT_SOURCE is still accepted with RADIATION_MODEL= P1.

I checked that the runner fails when the check in CConfig is turned off.

The 5 new CodeFactor issues are in lines from #2899 and #2961 that came in with the merge of develop. This PR does not change those lines.

@rois1995

Copy link
Copy Markdown
Contributor

Can you also fix the CodeFactor errors? At least we won't get those in the next PRs.

CodeFactor flags the aligned assignments of lam_naca0012_linelet and multi_interface in parallel_regression.py. This changes whitespace only.
@pragyaangaur

Copy link
Copy Markdown
Author

Yes, done!

Comment on lines 46 to +58
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);

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

the pointinEllipsoid function can be removed entirely and replaced by the python wrapper.
also, this entire function can be removed.

Copy link
Copy Markdown
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks, that makes sense. Removing HEAT_SOURCE completely is a bigger change than this PR, so I would like to do it in a separate PR, which I will open in a few hours. That PR will remove the HEAT_SOURCE options, SetVolumetricHeatSource, PointInEllipsoid and the code in CIncEulerSolver that uses them, and the config check will then tell users that HEAT_SOURCE has been removed and that they should use PYTHON_CUSTOM_SOURCE. The steady python wrapper test from this PR already shows that replacement.

If you are fine with that, could you approve this PR and resolve this thread so that it can be merged? I will then open the removal PR on top of it.

This branch has not been deployed

No deployments
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Projects

None yet

Development

Successfully merging this pull request may close these issues.

3 participants