- Mail actions: [ respond to this message ] [ mail a new topic ]
- Contemporary messages sorted: [ by date ] [ by thread ] [ by subject ] [ by author ] [ by messages with attachments ]

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

Date: Mon, 04 May 2009 14:23:20 +0200

Sander Valcke wrote:

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

*
This archive was generated by hypermail 2.3.0
: 2023-01-10 10:01:31 CET
*