bug in curvature.h?

[Petr Karnakov]

Hi all,

It's not quite clear if following code is correct.

It seems that in some cases the value is multiplied by sigma twice.

// http://basilisk.fr/src/curvature.h?raw  
  foreach (reduction(+:sa) reduction(+:sc)) {
    
    /**
    We then construct the final curvature field using either the
    computed temporary curvature... */

    double kf;
    if (k[] < nodata)
      kf = k[];
    else if (interfacial (point, c)) {

      /**
      ...or the average of the curvatures in the $3^{d}$ neighborhood
      of interfacial cells. */
      
      double sk = 0., a = 0.;
      foreach_neighbor(1)
        if (k[] < nodata)
          sk += k[], a++;
      if (a > 0.)
        // FIRST MULTIPLICATION BY SIGMA:
        kf = sigma*sk/a, sa++;
      else

        /**
        Empty neighborhood: we try centroids as a last resort. */

        // FIRST MULTIPLICATION BY SIGMA:        kf = sigma*centroids_curvature_fit (point, c), sc++; 

    }
    else
      kf = nodata;

    /**
    We add or set *kappa*. */
    
    if (kf == nodata)
      kappa[] = nodata;
    else if (p.add)
      // SECOND MULTIPLICATION BY SIGMA:      kappa[] += sigma*kf;
    else
      // SECOND MULTIPLICATION BY SIGMA:      kappa[] = sigma*kf;      
  }

Best,

Petr

[Stephane Popinet]

Hi Petr,

You are perfectly correct. This is a serious bug in the curvature
calculation. It was introduced by this patch:

Wed Dec 13 10:21:38 CET 2017 Stephane Popinet
* Generic interfacial force implementation

and will affect all calculations done since then, involving surface
tension and under-resolved interfaces. This could in particular cause
large surface tension noise and trigger extra numerical
breakup/fragmentation.

I have just released this patch which fixes the problem:

Mon Aug 19 09:23:56 CEST 2019 Stephane Popinet <pop...@basilisk.fr>
* A serious bug fix for low-resolution interface curvature calculations

http://basilisk.fr/src/?changes=20190819072356

Although several test cases checked the curvature / surface tension
calculation, none of them used weighing by a sigma different from unity,
which explains why the bug was not picked up... I have now modified the
curvature test case so that this does not happen again.

Thanks again for finding this!

Stéphane

[Wouter Mostert]

Dear Stéphane, Petr, all,

To check my understanding of when this bug shows up and what really is meant by an under-resolved surface:
Curvature is calculated based on:
(a) A parabola/paraboloid directly via height-function calculation of the surrounding 3x3, 3x3x3 stencil
(b) A parabola/oid on a least-squares fit of height-function info if the height function orientation is not sufficiently consistent in the stencil for (a) to work

(c) A parabola/oid based on the centroids of the nearby interface fragments in the stencil if the height-function calculations yield sufficiently many points within a half-cell distance of other points. That is, if fewer than 3 points (2D) or 9 (3D) are mutually separated by more than half a cell within the stencil, then (c) 'happens'.

As I understand it, the bug corresponds with (c) being used to compute a curvature, since otherwise either (a) or (b) is used to compute the "temporary curvature" and the double-multiplication is avoided. Is this right?

If so, then an example of an under-resolved interface may be a 2D bubble or droplet fewer than 2 cells in diameter, yes? More generally, and vaguely, anywhere where the interface gets really funky on a length scale smaller than the 3x3(x3) stencil would be under-resolved.

Cheers
Wouter

[Stephane Popinet]

Hi Wouter,

> Curvature is calculated based on:
> (a) A parabola/paraboloid directly via height-function calculation of
> the surrounding 3x3, 3x3x3 stencil

The stencil is not 3x3 or 3x3x3 (in 3D), it is 3xn or 3x3xn with n
variable and arbitrarily large (larger than 3 but usually smaller than 9
though).

> (b) A parabola/oid on a least-squares fit of height-function info if the
> height function orientation is not sufficiently consistent in the
> stencil for (a) to work

Yes.

> (c) A parabola/oid based on the centroids of the nearby interface
> fragments in the stencil if the height-function calculations yield
> sufficiently many points within a half-cell distance of other points.
> That is, if fewer than 3 points (2D) or 9 (3D) are mutually separated by
> more than half a cell within the stencil, then (c) 'happens'.
>
> As I understand it, the bug corresponds with (c) being used to compute a
> curvature, since otherwise either (a) or (b) is used to compute the
> "temporary curvature" and the double-multiplication is avoided. Is this
> right?

More or less. You missed "the average of the curvatures in the $3^{d}$
neighborhood of interfacial cells" part of the algorithm.

> If so, then an example of an under-resolved interface may be a 2D bubble
> or droplet fewer than 2 cells in diameter, yes? More generally, and
> vaguely, anywhere where the interface gets really funky on a length
> scale smaller than the 3x3(x3) stencil would be under-resolved.

Yes. Unfortunately this can happen even for interfacial shapes which are
not obviously under-resolved, for example a large bubble/drop
sufficiently deformed by turbulence. In this case, the bug may cause
this bubble to shed some fragments (in particular when the surface
tension coefficient is smaller than one) which can in turn further
destabilize the bubble etc.

If you want to convince yourself of the seriousness of this bug, re-run
the atomisation example:

http://basilisk.fr/src/examples/atomisation.c

with the buggy version and you will see that the number of small
droplets dramatically increases.

cheers,

Stephane

[Wouter Mostert]

Hi Stephane,

Thanks for the explanation. I see that really all breakup dynamics will be affected (at least, if \sigma != 1).

Cheers
Wouter