last version of the code breaks the wind setup when isink_radiation=1

[lsiess]

After the addition of the 4th order scheme (April 22, hash 282e51), the wind setup does not work anymore when isink_radiation = 1

[Yrisch]

Hi,
For now, you can force the code to work with leapfrog (set use_fourthorder to false). It should resolve the issue temporarily. I checked the routine called for sink radiation in the new substep force routine and it needs an update to use the extrapolation method (mandatory for 4th order). The position of gas and sink particles should be slightly shifted following Chin 2007. I will try to do that in the next few days. Which test I can try to validate the new implementation ?

[lsiess]

Hi,
An easy test is to run a wind simulation with an old version of Phantom that worked (e.g. 63d8aef) and compare the results with those obtained with the 4th order scheme.
To generate the code : writemake.sh wind > Makefile
and attached are the .in and .setup files
wind_setup.txt
wind_in.txt

Thanks!

[Yrisch]

Thanks ! I'll try to fix that soon !

[Yrisch]

Can you precise a bit more what is wrong now in this setup ? I tried both methods but I can't see a major difference...

[lsiess]

Hi Yann

If you plot the radial velocity vr = f(r)
Attached is what I obtain : in red old version in black new one

[danieljprice]

I had a look in the code and could not spot any obvious errors here, except maybe that tau is not updated? Is tau used for this test?

[lsiess]

No, isink_radiation should only add the radiation pressure force : kappa L/(4 pi r^2 c)

[danieljprice]

yes but there is the e^(-tau) term used in some cases

[lsiess]

tau should only be needed in Mat's ray tracing, i.e. if iray_resolution /= -1.
That we have not checked either!

[Yrisch]

the issue comes from the radiation pressure force that is not accounted during the middle force (using extrapolation) with FSI... I forgot to make a patch for this part of the code. I'm working on it.

[Yrisch]

Ok I made a patch, but I had to merge get_rad_accel_radiation routine with my main gas force loop. was there a specific reason that the two forces were calculated separately in the previous code ?

[Yrisch]

Should be fixed with #539

[danieljprice]

see comment on #539, it would be much nicer for this to stay inside the get_force routine if possible. Do we understand why there is a difference between calling it there and separately in the step routine?

[Yrisch]

I probably explained my patch badly... The radiation pressure forces are still computed inside the get_force routine. Before that, it was calculated outside the gas force loop (but still inside get_force). Now these radiation forces are computed inside gas loop of the get_force routine. I didn't understand why it wasn't the case in the previous version of the code (before FSI)...