Re: failed to converge in neighbour iteration in density()

From: Volker Springel <volker_at_MPA-Garching.MPG.DE>
Date: Mon, 04 May 2009 14:23:20 +0200

Sander Valcke wrote:
> Dear gadget list,
>
> I recently had a problem with the smoothing length update for SPH
> particles: for a certain (SPH-only) setup there are some particles for
> which the smoothing length updater does not converge. Some relevant
> output (the lines with (0) were added to the code manually, showing some
> variables of an erroneous particle):
>
> (0) ngb: 61.5331 h: 0.149935 left: 0 right: 0
> ngb iteration 1: need to repeat for 0008918068 particles.
> (0) ngb: 47.4712 h: 0.118996 left: 0 right: 0.149935
> ngb iteration 2: need to repeat for 0002706616 particles.
> (0) ngb: 39.3642 h: 0.0731442 left: 0 right: 0.118996
> ngb iteration 3: need to repeat for 0000786578 particles.
> (0) ngb: 0 h: 0 left: 0 right: 0.0731442
> ngb iteration 4: need to repeat for 0000246135 particles.
> (0) ngb: 0 h: 0 left: 0 right: 0.0731442
> ngb iteration 5: need to repeat for 0000005784 particles.
> (0) ngb: 0 h: 0 left: 0 right: 0.0731442
> ngb iteration 6: need to repeat for 0000001390 particles.
> (0) ngb: 0 h: 0 left: 0 right: 0.0731442
> ngb iteration 7: need to repeat for 0000000290 particles.
> (0) ngb: 0 h: 0 left: 0 right: 0.0731442
> ngb iteration 8: need to repeat for 0000000032 particles.
> (0) ngb: 0 h: 0 left: 0 right: 0.0731442
>
> You can see that this particle is assigned an h of 0 which, together
> with left being zero and right having the value it has, results in
> Gadget using "SphP[i].Hsml /= 1.26;" to update h, remaining at 0.
>
> Why did h become zero? Because of a large value of DhsmlDensityFactor in
>
> SphP[i].Hsml *=
> 1 - (SphP[i].NumNgb -
> All.DesNumNgb) / (NUMDIMS * SphP[i].NumNgb) *
> SphP[i].DhsmlDensityFactor;
>
> h actually became negative, it is then reset to 0 using
>
> if(SphP[i].Hsml < All.MinGasHsml)
> SphP[i].Hsml = All.MinGasHsml;
>
> What works for me is to set MinGasHsmlFractional to a non-zero
> (arbitrarily small) value in the parameterfile. What I have not looked
> into is why SphP[i].DhsmlDensityFactor is larger than it is supposed to,
> it is probably due to having a/some very close neighbour(s).
>
> Just thought I'd mention this in case someone had the same problem
> occur, or in case I'm missing something here.
>

Hi Sander,

Thanks for describing this issue in a clear way. Gadget2 applies Newton-Raphson
iteration to solve for the new Hsml. In rare cases it can however happen that
the present point is outside the convergence radius of the Newton-Raphson
scheme, i.e. one is too far away from the zero of the equation. Then the delta_h
correction calculated for the step can become inaccurate and may even lead one
into invalid terrain (for example into the negative regime).

One simple fix that I have adopted in newer versions of the code is to limit the
maximum allowed Newton Raphson change per iteration. To this end one can for
example replace

  SphP[i].Hsml *= 1 - (SphP[i].NumNgb -
                      All.DesNumNgb) / (NUMDIMS * SphP[i].NumNgb) *
                SphP[i].DhsmlDensityFactor;

with something like

   fac = 1 - (SphP[i].NumNgb -
                      All.DesNumNgb) / (NUMDIMS * SphP[i].NumNgb) *
                SphP[i].DhsmlDensityFactor;

   if(fac > 1.26)
        fac = 1.26;
   if(fac < 1/1.26)
        fac = 1/1.26;

   SphP[i].Hsml *= fac;


A too large Newton-Raphson step then keeps its "direction", but it should not be
able to produce a pathological overshoot any more.

Volker
Received on 2009-05-04 14:23:20