Question about Grid_DepositMustRefineParticles.C

[Nathan Goldbaum]

Hi all,

I'm currently trying to get sink particles working under the active particles framework. This is the first time I've dived deep into the enzo source, so I'm hoping the combined wisdom of the enzo developers list can help set me straight on one issue I'm having.

Right now I'm trying to get sink particle accretion working. One problem I've run into is that there doesn't seem to be an obvious way of making sure that the accretion zone around a sink particle is always at the maximum refinement level.

My first inclination was to simply make the sink particles a type of must refine particle and do the refinement using Grid_DepositMustRefineParticles.C or, if that's not possible with the way active particles are set up, use a new analogous refinement criterion. My understanding was that Grid_DMRP deposits the must refine particles onto the Flagging field using a nearest grid point algorithm. Taking a closer look at Grid_DMRP, it does in fact do this. There also seems to be an option to also flag a buffer region surrounding the must refine particle. By default, the buffer size is set to 1 but for problem types 106 and 107, the buffer size is set to 16 cells. Does anyone know what would happen if I changed buffersize? I have a feeling that the answer depends on the hydro method in use. I suspect that if I try to flag cells that are outside of the grid under consideration, I will cause a crash.

I could also see a more expensive way to do the flagging:

for each sink particle
for each grid on this level
if the grid overlaps with the sink accretion zone
flag the cells that are in the accretion zone

Before I go ahead and do this search, I want to make sure that the 'easy' way is actually not viable.

Thanks for your help with this,

Cheers,

Nathan Goldbaum

[Greg Bryan]

Hi Nathan -- I think it would work to make the sink particles a type of must refine particle. I haven't thought about how this works with the new ActiveParticle framework, but I don't see any obvious reason why it wouldn't work.

The one caveat is that the must refine particle is really a local refinement criteria. I see the buffer zone stuff you refer to; I didn't know about this, but I don't think it will work to flag cells beyond the edge of the current grid. At least for standard particles (non-Active), it will not. I don't think it will crash (since the code properly handles the edge of the grid), it just won't flag cells outside of the grid with the particle. At least, that's my reading of this code.

Your more expensive way seems like it would work, and allow a general way to specify some refinement region. If the number of sink particles is small, it doesn't seem that expensive.

Cheers,
Greg

p.s. Thanks for all your work helping to implement the ActiveParticles.

> --
> You received this message because you are subscribed to the Google Groups "enzo-dev" group.
> To post to this group, send email to enzo...@googlegroups.com.
> To unsubscribe from this group, send email to enzo-dev+u...@googlegroups.com.
> For more options, visit this group at http://groups.google.com/group/enzo-dev?hl=en.
>

[Nathan Goldbaum]

Hi Greg,

Thanks for your reply. I think I will have to go with the second choice as it's very important that the accretion zone is maximally resolved at all times. Seeing as there is currently no mechanism in enzo to have a refinement criterion that may affect grids on other processors, it looks like I'll have to add that logic myself.

The active particle have the ability to control their own formation, do local feedback, i/o, and do distributed feedback and other non-local operations at the beginning and end of EvolveLevel. To do what I want, I think we need one more function that will get called immediately after DepositParticleMassFlaggingField in RebuildHeirarchy. Each active particle type will then have the ability to modify the ParticleMassFlaggingField and thus force any region in the simulation to be resolved as needed.

If anyone sees possible issues with me modifying RebuildHierarchy in such a way, please let me know.

Cheers,

-Nathan

[David Collins]

Hi, Nathan--

I had one thought on this problem, for what it's worth. I'd propose
to do this with a boundary condition call on the FlaggingField, that
also includes ghost zones on the flagging field. So the flagging
algorithm would also flag into ghost zones, then a call to
CheckForOverlap would ensure that active zones on one grid that
overlap flagged boundary zones on another are also flagged. This
would involve breaking the RebuildHierarchy loop
for( j < NumberOfGrids )
FindSubgrids( Grid[j] )
loop into two, one that flags the zones and one that creates subgrids,
with a call to CheckForOverlap(grid1, grid2, *(your_new_function) in
between. This will scale with the rest of the machinery in the code,
so if I want to do a billion star particle run I can.

I'm mostly making this up on the spot, so feel free to ignore if this
doesn't do what you need.

d.

--
Sent from my computer.

[Greg Bryan]

Dave -- Nice idea! Although, it wouldn't work if the number of
required buffer zones in the flagging field was larger than the number
of ghost zones (as would be the case in the given example, since 16 >
3).

Greg

> --
> Sent from my computer.
>

[David Collins]

> Dave -- Nice idea!  Although, it wouldn't work if the number of required
> buffer zones in the flagging field was larger than the number of ghost zones
> (as would be the case in the given example, since 16 > 3).

Ah, right-- 16 ghost zones is a bit much :)

d.

>> --
>> Sent from my computer.
>>
>> --
>> You received this message because you are subscribed to the Google Groups
>> "enzo-dev" group.
>> To post to this group, send email to enzo...@googlegroups.com.
>> To unsubscribe from this group, send email to
>> enzo-dev+u...@googlegroups.com.

>> For more options, visit this group at
>> http://groups.google.com/group/enzo-dev?hl=en.
>>
>
> --
> You received this message because you are subscribed to the Google Groups
> "enzo-dev" group.
> To post to this group, send email to enzo...@googlegroups.com.
> To unsubscribe from this group, send email to
> enzo-dev+u...@googlegroups.com.
> For more options, visit this group at
> http://groups.google.com/group/enzo-dev?hl=en.
>

--
Sent from my computer.

[Nathan Goldbaum]

Hi David,

Thanks for the idea and I think that might actually work for my use case as I only need 4 buffer zones rather than 16. I'm not sure why someone required 16 buffer zones for problem types 106 and 107 - I certainly don't need that many.

If you have a chance, it would be great if you could get on the enzo irc channel and chat with me about this - I'm still a little bit hazy on how I'd implement your suggestion.

[David Collins]

Hi, Nathan--

> Thanks for the idea and I think that might actually work for my use case as I only need 4 buffer zones rather than 16.  I'm not sure why someone > required 16 buffer zones for problem types 106 and 107 - I certainly don't need that many.
> If you have a chance, it would be great if you could get on the enzo irc channel and chat with me about this - I'm still a little bit hazy on how I'd
> implement your suggestion.

Here's a cartoon of this idea. Now that I've written it, I have
doubts about it working. I think it's a necessary but possibly not
sufficient solution. Or, it might not work at all. I'll be on IRC
probably tomorrow, but it's easier for me to put all the pieces in an
email first.

The creation of new (empty) subgrids done in two steps, one grid at a
time. It's all done in FindSubgrids.C, which takes a single
HiererarchyEnergy as input, and hangs the new grids off of it in the
NextGridNextLevel linked list. The first stage is the flagging, which
is done in Grid_SetFlaggingField. The second stage creates a
ProtoSubgrid based on FlaggingField, and then hacks it into bite sized
pieces, with a new proto-subgrid for each bit sized piece. This
iterates until the computer gets bored, then the final set of
ProtoSubgrids are converted into Subgrids, and hung onto the linked
list. (I'll skip those details for now, not relevant here).

My suggest has this broken into two distinct loops, with a sort of
"boundary set" in between. I'm doing this off the top of my head, so
don't take this syntax too literally

//Flag Zones Only.
for( i=0; i<NumberOfSubgrids; i++)
FlagZones( Grid[i] )

//Set "boundary" conditions on the flagging field. here I gloss over
the communication and some arguments to the function call
for(i=0; i<NumberOfGrids; i++)
for( j=0; j< SiblingList.NumberOfSiblings[i]; j++)
Grids[i]->GridData->CheckForOverlap( SiblingList[ i ].GridList[ j
], POINTER_TO_NEW_FUNCTION);

//make new subgrids
for( i=0; i< NubmerOfSubgrids )
MakeNewSubgrids( Grid[i] )

So basically, FlagZones is everyting in FindSubgrids up to and
including the call Grid->SetFlaggingField, and MakeNewSubgrids is
everything after. CheckForOverlap is a pretty great function that
takes a pair of grids and looks for overlap, including the periodic
shifts, and if overlap is found, apply POINTER_TO_NEW_FUNCTION to the
pair of grids. The creation of FlagZones and MakeNewSubgrids should
be as easy as copy-paste (and ensuring that FlaggingField isn't
deleted). The middle bit is a little more complex.

For NEW_FUNCTION, I'd digest what Grid_CopyZonesFromGrid does, with a
couple modifications. Grid1->CopyZonesFromGrid( Grid2 ) looks at all
the points where Grid1 (active and ghost zones) overlaps the active
zones of Grid2, and copies that data from Grid2 to Grid2. For your
thing, you'd need to examine the ghost zones of Grid2 as well, and do
a max of the values of the flagging field in the overlap region.

There is some communication necessary-- SetBoundaryConditions.C is a
good place to learn how that's done. It's basically three steps:
post receive, post send, sit around and wait (done in
CommunicationReceiveHandler.C) You'll need to add an option to the
catcher in the latter routine to call your new function.

CommuniationSendRegion (called from CopyZonesFromGrid, it does most of
the heavy lifting) will accept a single field (in your case
FlaggingField), you might need to point a new BaryonField
(BaryonField[ NumberOfBaryonFields++ ]) to FlaggingField for this part
to maximize your code reuse.

So it's basically 3 steps;
1.) hack out the first half of FindSubgrids, put it in it's own loop
2.) hack out the second half of FindSubgrids, put it in it's own loop
3.) Create Grid_CopyFlaggingField to pull the relevant flagging field
from one grid to another
a.) clone CopyZonesFromGrid to examine the "other" grids ghost
zones, FlaggingField= max(FlaggingField,OtherFlaggingField)
b.) write the three-stage set of loops to handle the communication:
two double loops over CheckForOverlap, one call to RecieveHandler
c.) install CopyFlaggingField in RecieveHandler

Things are always way less complicated in my head before i actually
write out all the steps. Also, it's not clear that this is the best
solution, since you need 4 zones and there are typically only 3 ghost
zones, so that will cost you in time and memory. However, I think
that this is a necessary step on the way to proper hierarchy aware
refinement, which I think (and this has been the topic of conversation
for a while) is necessary for ensuring proper nesting 100% of the
time.

You might be able to use part of this to call a grid-pair-wise call to
SetFlaggingField, in the event that ghost zones aren't enough.

My major misgiving is the event in which you have to communicate to a
grid's parent that it needs to refine. That's a harder problem for
Enzo, and I haven't yet convinced myself that it's never going to
happen.

More soon,
d.

[Greg Bryan]

Although note that the number of buffer zones required (4) is larger than the number of ghost zones used by default (3), so more ghost zones would be required.

Greg