search

danieljprice / phantom

Phantom Smoothed Particle Hydrodynamics and Magnetohydrodynamics code
https://phantomsph.github.io
Other
103 stars 223 forks source link

Sinks cause fatal errors in cloud collisions #584

Open teogeo1996 opened 1 month ago

teogeo1996 commented 1 month ago

There is an issue with sink formation and evolution where, there needs to be a very specific combination of sink parameters for it not to throw a fatal error for the decrease of bin size beyond the allowed level. This affects mostly simulations with magnetic fields but it can also be present in simulation without MHD.

Yrisch commented 4 weeks ago

Hello, Could you please give more details on the setup used, the set of parameters that seems to work and the error logs when it crashes. Thank you.

teogeo1996 commented 4 weeks ago

Hi,

Regarding the issue with the sinks, I am simulating cloud-cloud collisions. I am using 1 million particles in the simulations each with 0.001 solar mass. The clouds are uniform density sphere and there is also a hot background with 1/100 of the sphere density and 100*T0. I run two different sets of simulations, one with and one without magnetic fields. My initial set of sink parameters is:

 isink_potential =           0    ! sink potential (0=1/r,1=surf)

   icreate_sinks =           1    ! allow automatic sink particle creation

    rho_crit_cgs =   1.000E-14    ! density above which sink particles are created (g/cm^3)

          r_crit =       0.001    ! critical radius for point mass creation (no new sinks < r_crit from existing sink)

           h_acc =       0.0003    ! accretion radius for new sink particles

 f_crit_override =        100.    ! unconditional sink formation if rho > f_crit_override*rho_crit

  h_soft_sinkgas =       0.002    ! softening length for new sink particles

 h_soft_sinksink =       0.002    ! softening length between sink particles

           f_acc =       0.800    ! particles < f_acc*h_acc accreted without checks

  r_merge_uncond =       0.000    ! sinks will unconditionally merge within this separation

    r_merge_cond =       0.000    ! sinks will merge if bound within this radius

However this combination throws a Fatal error of bin size being too small, with the magnetic field case throwing the error much faster than the non-MHD case. I was not using the f_crit_override variable initially however including it still got me a fatal bin error.

The only set of simulations that runs is:

 isink_potential =           0    ! sink potential (0=1/r,1=surf)

   icreate_sinks =           1    ! allow automatic sink particle creation

    rho_crit_cgs =   1.000E-15    ! density above which sink particles are created (g/cm^3)

          r_crit =       0.000    ! critical radius for point mass creation (no new sinks < r_crit from existing sink)

           h_acc =       0.001    ! accretion radius for new sink particles

 f_crit_override =        100.    ! unconditional sink formation if rho > f_crit_override*rho_crit

  h_soft_sinkgas =       0.004    ! softening length for new sink particles

 h_soft_sinksink =       0.004    ! softening length between sink particles

           f_acc =       1.000    ! particles < f_acc*h_acc accreted without checks

  r_merge_uncond =       0.003    ! sinks will unconditionally merge within this separation

    r_merge_cond =       0.000    ! sinks will merge if bound within this radius

However most of the sink particles are forced once the density exceeds f_crit_override*rho_crit_cgs rather than passing the tests.

I started by using eos=8, the barotropic EOS. With this EOS I was getting to densities up to 10(-8) without particles passing test, and thus cause a fatal error. After using f_crit_override, I still got fatal bin errors though. I changed to a dual isothermal EOS that set the temperature to 10K for the cloud and 1000K for the background. This formed some stars but still most of the stars were formed through f_crit_override and eventually gave me a fatal bin error. I need to say that I have experienced this in phantom again with bin size increasing to exceed allowed limit in other cases so it is not new. It bothered me that sink formation required that much tuning to work, which is why I decided to open an issue as from what I get in the code, given that I choose an h_acc larger than the jean’s length for rho_crit_cgs, at some point I should have sinks forming if the density becomes too large (and surely not as large as 10(-8) g cm**-3. There might be an issue in the Etherm or Eton check ? I was not able to identify one however.

Cheers Theo

On 13 Aug 2024, at 09:10, Yann bernard @.***> wrote:

External email to Cardiff University - Take care when replying/opening attachments or links. Nid ebost mewnol o Brifysgol Caerdydd yw hwn - Cymerwch ofal wrth ateb/agor atodiadau neu ddolenni.

Hello, Could you please give more details on the setup used, the set of parameters that seems to work and the error logs when it crashes. Thank you.

— Reply to this email directly, view it on GitHubhttps://github.com/danieljprice/phantom/issues/584#issuecomment-2285626659, or unsubscribehttps://github.com/notifications/unsubscribe-auth/A6TZEC4TMUNQZZDZYGR6K43ZRG5OJAVCNFSM6AAAAABML2XIQOVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDEOBVGYZDMNRVHE. You are receiving this because you authored the thread.Message ID: @.***>

Yrisch commented 4 weeks ago

I think the issue comes from r_crit. When you're setting it to 0.001, it means that no sinks can form into a sphere of that radius around each already formed sinks. However density in accretion flows near sinks can rise quite a lot which then causes the reduction of the step size and even a crash if the density is higher than the maximum bin. I had same behaviors in my runs. It also explains why the second setup works. With no constraint, sinks can form wherever the conditions are fulfilled and no particles reach super high density that will need very small time steps.

So the real question is why density in accretion flow can rise to so high values (much higher than rho_crit_cgs) near sinks vicinity ?

I fixed an issue in my last pull request where particles could only accrete where the particle with index 1 was awake which is really strange... If accretion is slowed down by that, it can artificially increase the density in region near sink boundaries causing the crash. Have you tried your simulations with the very last version of phantom ?

teogeo1996 commented 4 weeks ago

I understand this issue which is why I opted to change the r_crit. However, if the density rises to such high levels close to the sink, shouldn’t f_crit_override create a sink no matter what the r_crit value resctrits? As from my understanding of the pimass.F90, it overrides all other checks. I think other than changing r_crit, I needed to change f_acc to unity to stop this from happening.

Still the issue remains that if you do not use f_crit_override you essentially cannot pass the checks and form a sink for a parcel of gas that is collapsing to densities up to 10-8 g cm-3 when the r_crit_cgs is 10**-14 in cgs units. For a parcel of gas to collapse to such high densities it should be able to pass the tests set.

Cheers Theo

On 13 Aug 2024, at 14:14, Yann bernard @.***> wrote:

External email to Cardiff University - Take care when replying/opening attachments or links. Nid ebost mewnol o Brifysgol Caerdydd yw hwn - Cymerwch ofal wrth ateb/agor atodiadau neu ddolenni.

I think the issue comes from r_crit. When you're setting it to 0.001, it means that no sinks can form into a sphere of that radius around each already formed sinks. However density in accretion flows near sinks can rise quite a lot which then causes the reduction of the step size and even a crash if the density is higher than the maximum bin. I had same behaviors in my runs. It also explains why the second setup works. With no constraint, sinks can form wherever the conditions are fulfilled and no particles reach super high density that will need very small time steps.

So the real question is why density in accretion flow can rise to so high values (much higher than rho_crit_cgs) near sinks vicinity ?

I fixed an issue in my last pull request where particles could only accrete where the particle with index 1 was awake which is really strange... If accretion is slowed down by that, it can artificially increase the density in region near sink boundaries causing the crash. Have you tried your simulations with the very last version of phantom ?

— Reply to this email directly, view it on GitHubhttps://github.com/danieljprice/phantom/issues/584#issuecomment-2286226454, or unsubscribehttps://github.com/notifications/unsubscribe-auth/A6TZEC3GXQII3CHAQ2KJJC3ZRIBFBAVCNFSM6AAAAABML2XIQOVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDEOBWGIZDMNBVGQ. You are receiving this because you authored the thread.Message ID: @.***>

Yrisch commented 4 weeks ago

No, r_crit is used in the force routine when phantom try to find the densest particle to transform in sink (in force.F90). Overriding methods happen after this choice in ptmass_create. So it never tries to create new sink in other sink vicinity even with overriding.

Yrisch commented 4 weeks ago

What are the most common reasons why sinks don't want to form in your simulation ? In my case, I often had 'not all active particle' fail when using IND_TIMESTEPS. Is it also your case ?

teogeo1996 commented 4 weeks ago

For me it was the same. However I assumed that once the particles were synchronised that wouldn’t be an issue. The second most common issue was that the coronal parameter was larger than 1 and that etot>0

From: Yann bernard @.> Sent: Tuesday, August 13, 2024 3:15 PM To: danieljprice/phantom @.> Cc: Theotokis Georgatos @.>; Author @.> Subject: Re: [danieljprice/phantom] Sinks cause fatal errors (Issue #584)

External email to Cardiff University - Take care when replying/opening attachments or links. Nid ebost mewnol o Brifysgol Caerdydd yw hwn - Cymerwch ofal wrth ateb/agor atodiadau neu ddolenni.

What are the most common reasons why sinks don't want to form in your simulation ? In my case, I often had 'not all active particle' fail when using IND_TIMESTEPS. Is it also your case ?

— Reply to this email directly, view it on GitHubhttps://github.com/danieljprice/phantom/issues/584#issuecomment-2286369458, or unsubscribehttps://github.com/notifications/unsubscribe-auth/A6TZECZ3VDC7MJIP7BFUWGLZRIIHRAVCNFSM6AAAAABML2XIQOVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDEOBWGM3DSNBVHA. You are receiving this because you authored the thread.Message ID: @.***>

Yrisch commented 4 weeks ago

I don't think that these tests can hide issues. They have been extensively tested since the beginning of the code. I re-checked everything and it looks OK to me at least...

I advise you to retry your computations with my last patch if it wasn't the case yet. If it still crashes, you may try to recheck your initial conditions. Finally, it also could be a resolution issue. You might need more particles in your simulation even if you are resolving minimum Jeans mass with your accretion radius...

teogeo1996 commented 4 weeks ago

Ok sure I will. Thank you. I am assuming that the patch mentioned has been integrated into the main branch correct?

On 14 Aug 2024, at 09:34, Yann bernard @.***> wrote:

External email to Cardiff University - Take care when replying/opening attachments or links. Nid ebost mewnol o Brifysgol Caerdydd yw hwn - Cymerwch ofal wrth ateb/agor atodiadau neu ddolenni.

I don't think that these tests can hide issues. They have been extensively tested since the beginning of the code. I re-checked everything and it looks OK to me at least...

I advise you to retry your computations with my last patch if it wasn't the case yet. If it still crashes, you may try to recheck your initial conditions. Finally, it also could be a resolution issue. You might need more particles in your simulation even if you are resolving minimum Jeans mass with your accretion radius...

— Reply to this email directly, view it on GitHubhttps://github.com/danieljprice/phantom/issues/584#issuecomment-2288168432, or unsubscribehttps://github.com/notifications/unsubscribe-auth/A6TZEC6SP7DL6LAUFX2LABTZRMJB7AVCNFSM6AAAAABML2XIQOVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDEOBYGE3DQNBTGI. You are receiving this because you authored the thread.Message ID: @.***>

teogeo1996 commented 2 weeks ago

Hello again,

I am still facing this issue even with the new patch. I have increased resolution and still I cannot get passed the fatal timstep error. No matter what I do, the code just crashes. I do not understand why this is the case. I am using facc=1.0 so everything that goes into the accretion radius is accreted. I have also introduced a merging radius about 4 times the accretion radius to ensure I don’t get sinks too close together. I have set rcrit =0 to avoid having over densities close to the sink that are overlooked. I have increased the accretion radius to a crazy number (0.1 pc) and even reduced the sink formation threshold from 10-14 cos to 10-17 cgs. No matter what I do, the code just keeps on crashing. Can I get some support on this as this may be an issue that the tests have not caught. I get some simulations saying that the Courant tilmestep is too small, and some that the dt_clean tilmestep is too small, and some the dt_force tilmestep. After getting a log for all time steps I can see that the code is forcing a control time step on specific particles and that there are about 5-6 particles that go on to get very small time steps.

Cheers Theo

On 14 Aug 2024, at 09:34, Yann bernard @.***> wrote:

External email to Cardiff University - Take care when replying/opening attachments or links. Nid ebost mewnol o Brifysgol Caerdydd yw hwn - Cymerwch ofal wrth ateb/agor atodiadau neu ddolenni.

I don't think that these tests can hide issues. They have been extensively tested since the beginning of the code. I re-checked everything and it looks OK to me at least...

I advise you to retry your computations with my last patch if it wasn't the case yet. If it still crashes, you may try to recheck your initial conditions. Finally, it also could be a resolution issue. You might need more particles in your simulation even if you are resolving minimum Jeans mass with your accretion radius...

— Reply to this email directly, view it on GitHubhttps://github.com/danieljprice/phantom/issues/584#issuecomment-2288168432, or unsubscribehttps://github.com/notifications/unsubscribe-auth/A6TZEC6SP7DL6LAUFX2LABTZRMJB7AVCNFSM6AAAAABML2XIQOVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDEOBYGE3DQNBTGI. You are receiving this because you authored the thread.Message ID: @.***>

teogeo1996 commented 2 weeks ago

To further shed light on this, I am simulating cloud collisions. At the centre of the collision front I have fragmentation into sinks. Gas keeps flowing on the collision front and at some point I get this error. After further study, I find that some particles inside the collision front acquire extremely large velocities in a small period of time (they go from 1-2 km/s to about 200-300 km/s). This is what causes the tilmestep to decrease into very small values. If I opt to use a very small dtmax it catches this and overcomes the issue but I can’t help but feel like this is not natural behaviour and that something is wrong. The only thing that could cause this is sink-gas interactions but still, it doesn’t make sense. Could there be an issue with gas interacting with multiple sinks? I can provide the files that generate the error if nescecery.

Cheers Theo

On 25 Aug 2024, at 18:55, Theotokis Georgatos @.***> wrote:

Hello again,

I am still facing this issue even with the new patch. I have increased resolution and still I cannot get passed the fatal timstep error. No matter what I do, the code just crashes. I do not understand why this is the case. I am using facc=1.0 so everything that goes into the accretion radius is accreted. I have also introduced a merging radius about 4 times the accretion radius to ensure I don’t get sinks too close together. I have set rcrit =0 to avoid having over densities close to the sink that are overlooked. I have increased the accretion radius to a crazy number (0.1 pc) and even reduced the sink formation threshold from 10-14 cos to 10-17 cgs. No matter what I do, the code just keeps on crashing. Can I get some support on this as this may be an issue that the tests have not caught. I get some simulations saying that the Courant tilmestep is too small, and some that the dt_clean tilmestep is too small, and some the dt_force tilmestep. After getting a log for all time steps I can see that the code is forcing a control time step on specific particles and that there are about 5-6 particles that go on to get very small time steps.

Cheers Theo

On 14 Aug 2024, at 09:34, Yann bernard @.***> wrote:

External email to Cardiff University - Take care when replying/opening attachments or links. Nid ebost mewnol o Brifysgol Caerdydd yw hwn - Cymerwch ofal wrth ateb/agor atodiadau neu ddolenni.

I don't think that these tests can hide issues. They have been extensively tested since the beginning of the code. I re-checked everything and it looks OK to me at least...

I advise you to retry your computations with my last patch if it wasn't the case yet. If it still crashes, you may try to recheck your initial conditions. Finally, it also could be a resolution issue. You might need more particles in your simulation even if you are resolving minimum Jeans mass with your accretion radius...

— Reply to this email directly, view it on GitHubhttps://github.com/danieljprice/phantom/issues/584#issuecomment-2288168432, or unsubscribehttps://github.com/notifications/unsubscribe-auth/A6TZEC6SP7DL6LAUFX2LABTZRMJB7AVCNFSM6AAAAABML2XIQOVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDEOBYGE3DQNBTGI. You are receiving this because you authored the thread.Message ID: @.***>

Yrisch commented 2 weeks ago

Hi, It looks strange indeed... Yes, the setup and in files could be helpful to reproduce this issue. If the computation is long before crashing could you also send the latest full dump to relaunch the simulation just before the crash please. I'll test that on my side. ..

teogeo1996 commented 2 weeks ago

Attachment available until 26 Sep 2024 Hello again,

This is the closest I can get to the crash. I also noticed that there might be an issue when you choose to not save every dtmax but instead every 2 (e.g) dtmax. The code throws an error and doesn’t write a dump file. I have zipped the .in and the data file in this email.

Cheers Theo

Click to Downloadhttps://www.icloud.com/attachment/?u=https%3A%2F%2Fcvws.icloud-content.com%2FB%2FAbXqzdLtKkKqKyvZq85VYsXqOsVkAXNwIK6kPgV-DJVXBxFHl2_-uEpZ%2F%24%7Bf%7D%3Fo%3DAm2ZCOvwIkeA7dg_fZIE4NxKFT3ROlGPeNNlQTgcD8Fk%26v%3D1%26x%3D3%26a%3DCAog1k2oSnnbN5HeM8KXmDwpdQDuWgJzomtuRtQ6b4FA818SehCliZOpmTIYpZmO_aIyIgEAKgkC6AMA_yBYJgVSBOo6xWRaBP64SllqJ02Kja6uH_2pt4RcClGILRhec1ZUWkbACuUjdcGlBKlDa4Si-lFXRHInKilJ0Kmi_jq74EvoFFAdk-poBvrkqxS6Iwqhle-U8KZnRaSSY7S6%26e%3D1727376100%26fl%3D%26r%3D4C4180E5-CDF9-4E49-8D06-8FF074EFF5A3-1%26k%3D%24%7Buk%7D%26ckc%3Dcom.apple.largeattachment%26ckz%3D7BCC543A-EAFF-49AA-B046-209F01F3B2D6%26p%3D136%26s%3DanxJVElGzHntcoRY-uN8XSy1nmE&uk=JOlX2v_jyiHlYpl7OaGpew&f=Archive.zip&sz=110376677 Archive.zip 110.4 MB

On 26 Aug 2024, at 12:00, Yann bernard @.***> wrote:

External email to Cardiff University - Take care when replying/opening attachments or links. Nid ebost mewnol o Brifysgol Caerdydd yw hwn - Cymerwch ofal wrth ateb/agor atodiadau neu ddolenni.

Hi, It looks strange indeed... Yes, the setup and in files could be helpful to reproduce this issue. If the computation is long before crashing could you also send the latest full dump to relaunch the simulation just before the crash please. I'll test that on my side. ..

— Reply to this email directly, view it on GitHubhttps://github.com/danieljprice/phantom/issues/584#issuecomment-2309932772, or unsubscribehttps://github.com/notifications/unsubscribe-auth/A6TZEC5DMBJKWX3IHCPJD43ZTMDDXAVCNFSM6AAAAABML2XIQOVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDGMBZHEZTENZXGI. You are receiving this because you authored the thread.Message ID: @.***>

teogeo1996 commented 2 weeks ago

Hello again,

I also noticed that there might be an issue when you choose to not save every dtmax but instead every 2 (e.g) dtmax. The code throws an error and doesn’t write a dump file. I can provide you with the files that throw the error, but I can't upload them here since there is a file size limit even if they are zipped. If you provide me with an email I can send you a link to an one drive folder with them.

Cheers Theo