Adding a sink particle to the grid

The easiest solution to get what I am looking for (a sink at the origin that accretes material in a spherically symmetric manner after the center most cell satisfies the Fedderath conditions) is to hack astrobear rather than apply fixes and checks to the code for particles near reflecting boundary conditions. This is something I would be interested in pursuing, but due to time wanted to just get the ball rolling on getting a sink at the origin to check flow properties post formation..

Given the sparse documentation on the wiki to do this, I thought I'd better outline the procedure I learned about here, so it can be ported over to an actual wiki page later..

The hack I did has the following attributes:

  1. Initialize a fixed sink and associated point gravity object
  2. Kill any new checks for sinks so that as simulation progresses, no others will form (for the BE collapse case no other sinks SHOULD form anyway)
  3. Modify accretion so that new values of the mass is increased by factor of 8 to account for accretion from the reflecting ghost zones along x,y,z.
  4. Kill all motion from sink by removing particle momentum

As I outlined in previous blog posts on the wiki, amr_control.f90 calls various routines for checking for new sinks, initializing and advancing them, etc.

In particular the first pertinent routine it calls is particle_preupdate.. In this routine, there are loops that 1) check for new particles and 2) collects new particles. I commented both of these out and took some code from collects new particles and placed in my problem.f90.

Inside collects new particles, there are the following lines:

NULLIFY(NewParticle) CALL CreateParticle(NewParticle) NewParticle%xloc(1:nDim)=GlobalParticles(1:nDim,i) NewParticle%moments=GlobalParticles(nDim+1:,i) NewParticle%iAccrete=FEDERRATH_ACCRETION CALL AddSinkParticle(NewParticle) CALL CreatePointGravityObject(NewParticle%PointGravityObj) NewParticle%PointGravityObj%x0(1:nDim)=GlobalParticles(1:nDim,i) NewParticle%PointGravityObj%soft_length=1d0*sink_dx NewParticle%PointGravityObj%soft_function=SPLINESOFT

These lines initialize a new particle and associate with it a point gravity object. Now, there is a distinction between sink particles and point gravity objects. To my understanding, sinks can accrete mass and move, but point gravity objects can only exert a gravitational force. Thus when one wishes to accrete mass onto a gravitational point, one should initialize a sink + point gravity object. (Whereas if you only cared about exerting a force, but no need for accretion onto the particle — accumulation of matter in the surrounding cells would be fine — you could use just a point grav obj).

So, I added these lines to my problem.f90 in the routine problem module init, below the ambient and clump objects. To use these routines, I also added a USE Particle Declarations at the top of the file.

IF (.NOT.LRESTART) THEN

NULLIFY(Particle)

CALL CreateParticle(Particle)

Particle%lfixed=.true. !keeps particle fixed at xloc

Particle%xloc(1:ndim)=0d0

Particle%iAccrete=1 !Fedderath accretion, as given in particle def

Particle%moments(1:3)=0

!check print statement for moments

CALL AddSinkParticle(Particle) !accesses particle_

CALL CreatePointGravityObject(Particle%PointGravityObj)

!accesses pointgravitysrc

Particle%PointGravityObj%x0(1:ndim)=0

Particle%PointGravityObj%soft_length=1d0*sink_dx

!not sure about this

Particle%PointGravityObj%soft_function=SPLINESOFT

!default mass = 0 as defined in the type

! starts gaining mass when surrounding gas is jeans unstable for fedd accretion (should coincide with same time a sink formed in previous simulations..)

END IF

The grid is then initialized with this point particle at the origin. It lives there, massless and with 0 momentum throughout the simulation until it begins to accrete. I found that I did not have to explicitly 0 out the momentum though, as checks were already in place that did this for a fixed particle.

To modify the accretion, in finalize accretion routine, I multiplied dq*8:

Particle%Q(1:NrHydroVars)=Particle%Q(1:NrHydroVars)+Particle%dQ(1:NrHydroVars)*8d0

Now the sink will be accreting 8 octants worth of material :)

Now that the edits have been made, the overall structure of the code goes like this —

1) sink + point gravity obj is initialized on the grid

2) code goes through pre-update that no longer checks for new particles, but for the one particle in the particle list (added to it in the problem module) it does accretion, synchs accretion, and finalizes accretion. In the do accretion routine, the code checks whether the jeans criteria has been violated around the sink particle. If it has, then the accretion takes place. Fedderath accretion removes just the amount of gas from the cells to make the cells no longer jeans unstable, and adds that gas to the particle. This is the type of accretion I am using. Away from the sink, one expects a bondi flow to be reached..

3) Code Finalizes accretions where it zeroes out particle's momentum, and adds new values to the q array for accreted mass, energy (multiplied by the factor of 8).

And this all repeats for the next time step.

Comments

No comments.