comphy-lab / comphy-lab/MultiRheoFlow

log-conform: omega = (...)/(Lambda.y - Lambda.x) guarded only at 1e-20; reaches 2e5 in isotropic far-field cells

Open
#8 0 comments 0 reactions 0 assignees View on GitHub
Dominant language
C
Stars
5
Forks
2
Avg merge
1h 16m
Merged PRs (30d)
1

Description

## Summary

In `log-conform-viscoelastic-scalar-2D.h` (and the same construction in the 3D header), the rotation rate

```c
double omega = (Lambda.y*M.x.y + Lambda.x*M.y.x)/(Lambda.y - Lambda.x);
OM = (R.x.x*R.y.y - R.x.y*R.y.x)*omega;
```

is guarded only by `fabs(Lambda.x - Lambda.y) <= 1e-20`. In a bubble-bursting campaign at `Oh = 0.024`, `De = 0.02–0.2`, `Ec = 0.009`, level 12, we measured `|OM|` up to `2.4e5` in far-field cells where `A = I` to `5e-9` — i.e. the eigenvalue gap is small but far above the `1e-20` threshold, so the guard never fires and the division amplifies round-off.

## Why it is harmless today

Those cells take the `sq(A12) < 1e-15` branch, which sets `R = I` and hence `Psi12 = 0` exactly. The enormous `OM` then multiplies `(Psi22 - Psi11)` and `s * OM` terms that are identically zero or cancel, so no damage propagates. The header is correct *as long as* that branch and the `OM` construction stay coupled. The cancellation is incidental rather than designed, and any future change to the `A12` threshold, the `R = I` shortcut, or the `Psi` update ordering would expose the division.

## Suggested hardening

Either of these keeps behaviour identical in the well-conditioned case and removes the latent hazard:

1. Skip the `omega` branch when the conformation is isotropic to working precision, e.g. `fabs(Lambda.x - Lambda.y) <= 1e-12 * max(fabs(Lambda.x), fabs(Lambda.y))` (relative, not absolute), and set `OM = 0` there — the rotation of an isotropic tensor is undefined and physically irrelevant.
2. Or compute `omega` from the eigenvector-projected velocity gradient directly, which avoids dividing by the eigenvalue gap at all (Fattal & Kupferman's original construction admits this form).

## Context

Found while diagnosing a divergence in `critical-bursting-bubble-ve` that turned out to be a different unbounded term (the axisymmetric hoop source `Psi_qq += dt*2*u_r/r`, for which we have added a `CFL_conform` stability bound in our vendored copy; will propose separately once its A/B test finishes). This one is reported because it was measured along the way and is a maintainability hazard rather than a live bug.

Header provenance in our copy: vendored from MultiRheoFlow, reconciled with `7d9c3df` (2026-08-30).

Contributor guide

Open the contributing guide

Assessment

This issue has not been assessed yet.

Get new issues in your inbox

A short digest of beginner-friendly GitHub issues.