code_saturne User's Forum

Bonjour
Nous pensons avoir identifie un bug assez technique dans le traitement du radiatif en géométrie 2d axi dans Syrthes 5.0 qui conduit a faire diverger des cas par accumulations d erreurs surestimées qui ne divergeraient pas autrement. Nous avons patche notre version.
Nous ne trouvons ps de canal pour faire remonter l info aux développeurs.
D'ou ce message.
Y a t il un process pour ce genre de remontée d informations (comme sur github)?
Merci

JFL

Hello,

There is currently no external GitHub of Gitlab mirror for Syrthes, and the support address may not be checked as regularly as that for code_saturne (development resources on Syrthes are currently quite limited, but some members of the code_saturne dev team have write access to the internal Syrthes git repository).

So the best solution (at least for Syrthes) is to post your patch here, or send it to saturne-support@edf.fr, so we can integrate it (or at least try).

It is preferable to use English on this forum (so a broader community can benefit), but you can use French on the support address.

Thanks for your contribution !

Yvan

Ok Please see below
Just to be clear these bugs were found by an AI Agent that we tasked to audit the reason some of our cases were crashing without any apparent reasons.

After we applied the patch all our concerned cases converged nicely and provided coherent results where they crashed before.
Any feedback will be appreciated
JFL

# SYRTHES 5.0 — Three Bug Reports
version: SYRTHES 5.0 (kernel source)

Configurations affected: 2D axisymmetric confined radiation (`PRISE EN COMPTE DU RAYONNEMENT CONFINE= OUI`, `iaxisy=1`), single-band gray model, sequential (non-MPI-parallel radiation) path with external view factor file (`lec_fdf=1`).

Patches for all three files are attached. A pre-built patched library (`libsyrthes_cfd.patched.a`) can be also provided for immediate use without recompilation.

---

## Bug 1 — `syr_calfdf.c`: wrong `alter_axi1` direction flag after view factor computation

### File and location

`syr_calfdf.c`, function `calfdf()`, the line immediately after the call to `facforme_2a`:

```c
/* BUGGY (line ~157) */
if (maillnodray.iaxisy==1) alter_axi1(1, maillnodray, defsymper);
```
### Root cause

`alter_axi1(flag, ...)` applies a coordinate permutation on the 2D axisymmetric radiation mesh:
- `flag=1` — forward transform: `(x, r) → (r, -x)`
- `flag=0` — inverse transform: `(r, -x) → (x, r)` (restore)

The forward transform (`flag=1`) is applied before `facforme_2a` to align the mesh for the integration kernel. After integration the mesh must be restored to its original `(x, r)` layout. The line above fires after `facforme_2a` returns and should therefore use `flag=0`.

With `flag=1` at this location the transform is applied a second time in the same direction, leaving all nodes at `(-x, -r)` — a reflection through the origin that has no physical meaning.

### Observable symptom

The solid/radiation mesh matching step (`corresp_solray` or equivalent) computes correspondence distances on the order of the domain diagonal (~1 m for a typical metre-scale case) instead of zero. SYRTHES prints something like:

```
* corresp_solray : distance max entre maillages = 9.87e-01
```

All coupling exchanges between the solid solver and the radiation mesh are then meaningless.

### Fix

```c
/* FIXED */
if (maillnodray.iaxisy==1) alter_axi1(0, maillnodray, defsymper);
```

### Minimal verification

Set up any 2D axisymmetric radiation case (even a single annular ring with 4 elements). With the bug, `syrthes.log` will report a non-zero correspondence distance. With the fix, the distance should be exactly zero:

```
* corresp_solray : distance max entre maillages = 0.00e+00
distance moy entre maillages = 0.00e+00
```
---

## Bug 2 — `syr_resray.c`: `pow(rr, 0.25)` returns NaN when irradiation sum is negative

### File and location

`syr_resray.c`, function `fi2teq()`, single-band (`phyray.bandespec.nb==1`) branch, equivalent temperature computation loop:

```c
/* BUGGY (line ~683) */
rr=0;
for (nn=0; nn<listefdf->nb; nn++)
rr += listefdf->fdf[nn] * radios[0][listefdf->numenface[nn]];
rr += diagfdf * radios[0];
if (propinf.actif)
rr += fdfnp1*sigma*propinf.EO[0] * pow(propinf.TO+tkel, 4);

trayeq = pow(rr / (sufray*sigma), 0.25); /* BUG: rr may be < 0 */
```

### Root cause

`rr` is the total incoming radiative power on face `i`, accumulated from view-factor-weighted radiosity values. The Gaussian quadrature in `facforme_2a` (2D axisymmetric kernel) can produce small negative view factor entries due to near-singularity integration error (this is expected and acknowledged in the `smoogc` regularisation philosophy). When an element has mostly near-zero neighbors, the sum of negative view factor contributions can make `rr` slightly negative (e.g., `rr ≈ -1e-3 W/m²`).

In C, `pow(x, 0.25)` with `x < 0` is a domain error and returns `NaN` (per IEEE 754 and C99 §7.12.7.4). This is not a platform-specific behaviour: it is guaranteed on all conforming implementations.

The NaN then propagates:
1. `trayeq[i] = NaN`
2. → `transRS` copies it to `scoupr.t` (the surface temperature used in conjugate exchange)
3. → the linearised radiation exchange coefficient `hray` computed from `scoupr.t` becomes NaN
4. → NaN enters the diagonal of the coupled system matrix
5. → GRCONJ diverges on the first iteration (residual = NaN from step 1)

The solver prints:

```
GRCONJ : iter= 1 res= NaN
GRCONJ : iter= 2 res= NaN
...
```

and the simulation produces all-NaN temperature fields.

### Fix

Clamp the unphysical negative irradiation to zero before the `pow` call:

```c
/* FIXED */
if (rr < 0.0) rr = 0.0;
trayeq[i] = pow(rr / (sufray[i]*sigma), 0.25);
```

This is physically correct: negative incoming radiation is non-physical and is the result of numerical noise. Clamping to zero is equivalent to assuming zero incoming flux, which is conservative and safe.

### Note on scope

The parallel counterpart `fi2teq_parall()` computes the irradiation sum via `fdf_Y_parall` and stores it in `trayeq[i]` before the `pow` call (line ~922). It is also vulnerable if the accumulated sum is negative, and should receive the same guard. The multiband branches in both functions accumulate `rr` over spectral bands; the risk there is smaller (multiple positive contributions) but not zero.

### Minimal verification

Take any 2D axisymmetric single-band radiation case with a coarse radiation mesh (few elements, high L/gap aspect ratio). Add a `printf` before the `pow` call:

```c
printf("face %d: rr=%e\n", i, rr);
```

At least one face will print a small negative value. With the fix, GRCONJ converges normally:

```
GRCONJ : 1 res= 3.2e-02
GRCONJ : 2 res= 4.7e-04
...
GRCONJ : converge
```

---

## Bug 3 — `syr_init_radiation.c`: `iniparallray` not called in sequential path, causing NULL-pointer crash in `lire_fdf`

### File and location

`syr_init_radiation.c`, function `init_radiation_seq()`, before the first mesh-level operation (`lire_divers_ray`):

```c
/* BUGGY — iniparallray missing entirely */
void init_radiation_seq(...)
{
...
lire_divers_ray(*maillnodray, raysymper, volconnexe, &(phyray->bandespec), propinf, mask);
...
/* later, if lec_fdf==1: */
lire_fdf(listefdf, diagfdf, propinf, sdparall_ray); /* CRASH: sdparall_ray->nelrayloc is NULL */
}
```

### Root cause

`iniparallray(maillnodray, sdparall_ray)` allocates and fills two arrays inside `sdparall_ray`:
- `sdparall_ray->nelrayloc[npartsray]` — number of radiation elements local to each partition
- `sdparall_ray->ieledeb[npartsray+1]` — element index offsets

`lire_fdf()` (called when `lec_fdf=1`, i.e., `LECTURE DES FACTEURS DE FORME SUR FICHIER= OUI`) dereferences both arrays unconditionally:

```c
/* inside lire_fdf, ~line 513 */
nelemloc = sdparall_ray->nelrayloc[sdparall_ray->rangray];
```

The parallel path `init_radiation_parall()` correctly calls `iniparallray` at its start (line ~318). The sequential path `init_radiation_seq()` does not, leaving `nelrayloc` and `ieledeb` as NULL pointers. The code path through `calfdf` (view-factor computation) happens to call `iniparallray` internally as well, so the bug is invisible when view factors are computed at runtime (`lec_fdf=0`). It only triggers when an external view factor file is loaded (`lec_fdf=1`), because in that case `calfdf` is skipped entirely and `lire_fdf` is called directly.

### Observable symptom

Segmentation fault (SIGSEGV) during radiation initialisation, immediately after printing:

```
* calfdf : lecture des facteurs de forme sur fichier...
```

The backtrace points into `lire_fdf` at the `sdparall_ray->nelrayloc[...]` dereference.

### Fix

Add `iniparallray` at the top of `init_radiation_seq()`, mirroring what `init_radiation_parall()` already does:

```c
/* FIXED */
void init_radiation_seq(...)
{
...
/* iniparallray allocates nelrayloc and ieledeb required by lire_fdf
when lec_fdf=1. The parallel path already does this; add it here
for the sequential path. */
iniparallray(*maillnodray, sdparall_ray);

lire_divers_ray(*maillnodray, raysymper, volconnexe, &(phyray->bandespec), propinf, mask);
...
}
```

### Minimal verification

Enable external view factor loading in the `.syd` file:

```
LECTURE DES FACTEURS DE FORME SUR FICHIER= OUI
```

Provide a valid `.fdf` file. Without the fix the run crashes immediately with SIGSEGV. With the fix, SYRTHES reads the file and continues normally:

```
* lire_fdf : lecture des facteurs de forme...
* lire_fdf : 58 facteurs de forme lus
```

---

## Summary table

| # | File | Function | Trigger condition | Symptom |
|---|------|----------|-------------------|---------|
| 1 | `syr_calfdf.c` | `calfdf` | Any 2D axisym. radiation case (`iaxisy=1`) | Mesh correspondence distance ~1 m; broken coupling |
| 2 | `syr_resray.c` | `fi2teq` | Single-band gray, slight negative VF entries | GRCONJ diverges with NaN residuals |
| 3 | `syr_init_radiation.c` | `init_radiation_seq` | External VF file (`lec_fdf=1`) in sequential mode | SIGSEGV in `lire_fdf` |

Bugs 1 and 2 are independent and can each occur alone. Bug 3 only manifests when `lec_fdf=1`; it is masked by bug 1 if both are present because the crash happens before the coupling distances are evaluated.
: