Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

periodicity of face centred data #327

Closed
darylbond opened this issue Oct 2, 2018 · 54 comments
Closed

periodicity of face centred data #327

darylbond opened this issue Oct 2, 2018 · 54 comments

Comments

@darylbond
Copy link
Contributor

I am implementing a face-centred solver and am having some issues with setting periodic domains. Essentially the problem is repeated data at the periodic boundary. An example is given below for a domain with 8 cells and 5 ghost cells initialised with the x-index of each face.

Initial face data:
[0 0 0 0 0 0 1 2 3 4 5 6 7 8 0 0 0 0 0]

Periodic face data:
[4 5 6 7 8 8 1 2 3 4 5 6 7 8 8 1 2 3 4]

What I actually want:
[3 4 5 6 7 0 1 2 3 4 5 6 7 8 1 2 3 4 5]

This is all through the FillPatch routine working on face-centred state data. For reference I am using the EB/CNS tutorial as the basis for my implementation.

Any suggestions?
@WeiqunZhang
Copy link
Member

FillPatch does not currently support face centered data. Nevertheless this can be done. You might want to take a look at https://github.com/AMReX-Astro/MAESTROeX/blob/master/Source/MaestroFillData.cpp#L233

@darylbond
Copy link
Contributor Author

@WeiqunZhang That would explain it! Many thanks.

@jared321
Copy link
Contributor

Are there plans to update FillPatch to work with face-centered variables?

@asalmgren
Copy link
Member

@jared321 -- we have been talking about this. Do you know what interpolation scheme you would want to use to fill fine face-centered data from coarse face-centered data?

@chaw0023
Copy link
Contributor

chaw0023 commented Sep 5, 2019

@jared321 -- we have been talking about this. Do you know what interpolation scheme you would want to use to fill fine face-centered data from coarse face-centered data?

Great to see activity on this front. The only requirement we have is second order accuracy be maintained during interpolation. As long as the interpolation scheme satisfies this requirement, it should be fine.

@asalmgren
Copy link
Member

Just to clarify -- no attempt with the interpolation to maintain any divergence constraint?

@kweide
Copy link
Contributor

kweide commented Sep 5, 2019

Not sure whether divergence constraint will be required. My understanding is that for FLASH5, face variables will not be used for representing magnetic field components in an MHD solver.

@asalmgren
Copy link
Member

I'm going to close this issue for now. Please re-open it (or create a new issue) when you know what interpolation scheme you would like to have as the default (we can always add more later).

@chaw0023
Copy link
Contributor

chaw0023 commented Sep 6, 2019

Just to clarify -- no attempt with the interpolation to maintain any divergence constraint?

I'm sorry I failed to mention this earlier. The divergence free constraint is critical for incompressible flow calculations so yes it is required.

@asalmgren
Copy link
Member

asalmgren commented Sep 6, 2019 via email

@asalmgren asalmgren reopened this Sep 6, 2019
@chaw0023
Copy link
Contributor

chaw0023 commented Sep 6, 2019

The existing divergence preserving quadratic interpolation method in FLASH4 as well as the stencil used is described in detail in this paper. Vanella et al., J. Comp. Phy. 229 (18) 6427-6449, 2010.

Fig. 3 (a) and 3(b) show stencil for restriction (interpolation from level l+1 to l) and prolongation (interpolation from level l to l+1) for 2D case, respectively; discussion in Sec. 3.2. 3D prolongation stencil is described in Appendix A and Fig. 21. u, v and w are the relevant face-centered variables throughout the discussion.

Figures look complex but description in Sec. 3.2 and Appendix A is well written. I am willing to take a call if oral clarification is needed.

image
(Typo in Fig.3a caption: *prolongation -> restriction)

image

@asalmgren asalmgren mentioned this issue Sep 9, 2019
Closed
@kweide
Copy link
Contributor

kweide commented May 5, 2020

Now that both this Issue and the other one that referred to it (#424) have been closed, I would just like to note (for FLASH) that we may still be interested in support for face-centered data by AMReX.

@asalmgren
Copy link
Member

asalmgren commented May 5, 2020 via email

@kweide
Copy link
Contributor

kweide commented May 5, 2020

It seems this was primarily about interpolation. @chaw0023 gave details in September 2019, see above (in Issue #327), I did not see any followup.

I don't know details about the current state of our requirements, maybe @adubey64 can say more.

@asalmgren
Copy link
Member

asalmgren commented May 5, 2020 via email

@chaw0023
Copy link
Contributor

chaw0023 commented May 7, 2020 •
edited
Loading

(I'm chipping in to explain the reason for request at that time from FLASH.)
The FillPatch method did not support face centered data at that time and the issue was opened to request this functionality. The functionality of interpolation being second order and divergence-preserving are critical as per original requirement.

Is this issue closed because it is resolved now or put on hold currently?

@asalmgren
Copy link
Member

asalmgren commented May 7, 2020 via email

@kweide
Copy link
Contributor

kweide commented May 7, 2020

Let me explain what happens when we have "ghost faces" filled when using PARAMESH, which has shaped our expectations. @chaw0023 can jump in if his requirements are different or more specific. Or if I am missing the point entirely.

We call amr_guardcell (with appropriate arguments and flags), and it will do the following, for each block that is selected (focussing on face variables only):

  • all cell faces in guard cell regions, i.e., outside of the block's valid region, [possibly only for a subset of layers,] are updated:
    1. if filling from a same-level neighboring region, just copy.
    2. if filling from a finer neighboring region, averaging happens.
    3. if filling from a coarser neighboring region, interpolation happens. This can be made divergence-free or maybe divergence-preserving.
    4. if at a domain boundary, invoke special magic.
  • That was for faces outside of the valid region. For faces on the surface of the valid region of the block, the following applies:
    1. Be default, such cell faces are never updated by guard cell filling.
    2. there is a force_consistency option provided by PARAMESH. When that is turned on, face values on the boundary of the valid region are also updated on the coarse side, based on the face values from the fine side. (Sort of like flux correction, mathematically.)
    3. There may be some additional optional action at domain boundaries.

I don't know whether this (or how much) should be implemented in FillPatch or whether it belongs elsewhere.

@asalmgren
Copy link
Member

asalmgren commented May 7, 2020 via email

@kweide
Copy link
Contributor

kweide commented May 7, 2020

No connection specific to face variables, as far as I am aware.

As far as guard cell filling more generally, its relationship to orchestration is... poorly understood by me so far. Maybe @jared321 will care to comment.

@jared321
Copy link
Contributor

jared321 commented May 7, 2020

The orchestration system has two main sub-systems - an offline toolchain and a runtime. We previously referred to the runtime as the orchestration layer.

In terms of the runtime, we hope to have the runtime make use of the phase-asynchronous GC fill whenever possible so that it will overlap the GC communication with computational work. It will be the job of the offline toolchain to figure out when this is possible and make sure that the runtime uses it in those cases. It might also be that the use of the asynchronous GC fill could be controlled by a runtime parameter.

Somehow I think that I have given you both too much and too little information. Do you have a more specific question related to the orchestration system and GC fills?

@asalmgren
Copy link
Member

asalmgren commented May 7, 2020 via email

@kweide
Copy link
Contributor

kweide commented May 7, 2020

Ideally this would all happen inside the iterator, hidden from the iterator-using code. The iterator's next() method would ship out the next block if and when there is one available, with all necessary communication (for this block) done, guard cells filled (whether there was interpolation involved or not), etc. [I guess this implies that next() could block.]

As far as iterator-using code then: some higher-level algorithms require iterating over one refinement level at a time, so an implementation of such an algorithm would request ("build") a level-specific iterator. Some algorithms just require visiting all blocks no matter what refinement or order, so it would be more natural for the implementation to request an all-level iterator. (The latter might under the hood still operate level by level; maybe already returning blocks for one level while other levels are not yet ready.)

Did I get this right, @jared321 ?

@asalmgren
Copy link
Member

asalmgren commented May 7, 2020 via email

@jared321
Copy link
Contributor

jared321 commented May 7, 2020

I have not looked specifically at how to get the phase asynchronous iterator into our code, but believe that it should be possible. My understanding of how FLASH5 would interface to this iterator matches what you write @kweide.

@asalmgren As part of the collaboration to refine the octree mode of AMReX do you think it is likely that we would evolve the current octree iterator or develop a new one?

@asalmgren
Copy link
Member

asalmgren commented May 7, 2020 via email

@kweide
Copy link
Contributor

kweide commented May 7, 2020

I'm not so sure myself what I'm asking...

But basically with "asynchronous" I meant that client code of the iterator can already process some blocks, while the iterator itself is still doing behind-the-scenes work (like waiting for communication) to deliver more blocks later.

@atmyers
Copy link
Member

atmyers commented May 8, 2020

AMReX currently has some ability to overlap communication with computation, in that you can start a FillBoundary operation, do some other stuff, and then later call finish that on the FillBoundary. It does not have the capability you talking about, however, which would be starting the FillBoundary (for example), then immediately start iterating over the blocks as they become ready. Maybe some of them require only local copies and are therefore ready immediately, and can be worked while waiting for slower MPI communication. I agree that this would be interesting and potentially useful for many AMReX applications, but I think this is really a separate issue from the FillPatching.

MAESTROeX does fill patching on face centered data. The FillPatch functions in AMReX_FillPatchUtil.cpp don't handle this case, so MAESTROeX implements a special version that does not use those functions. One way of doing this would be to implement a similar function in FLASH. That kind of makes sense to me, given that there is a particular interpolation method that you would like to use.

Alternatively, we could implement this kind of feature in AMReX and add them to the fill patch utilities.

Sorry if this redundant information, I'm just trying to summarize the current state of things.

@atmyers atmyers reopened this May 8, 2020
@kweide
Copy link
Contributor

kweide commented May 11, 2020

That was very useful information at least for me, thank you.

When you mentioned "some ability to overlap communication with computation", was that in reference to features provided by AmrTask and/or the "phase asynchronous AMR execution"? I am currently looking at that paper and trying to understand whether it is relevant to the conversation here.

@drummerdoc
Copy link
Member

I believe that the paper you are looking at describes the extent to which AMReX supports the async features you want. The "phase asynchronous" bit is a kind of nuance that allows one to limit just how deep mods for asynchrony need to reach within AMReX - it reflects a notion that many of our algorithms have been written with a level-by-level organization structure, which inherently brings with it a price, in terms of required synchronizations. I'd say the benefit is the ability to re-use a ton of code and complex algorithmic implementation, including downstream application use, which is also often organized by level.

@atmyers
Copy link
Member

atmyers commented May 11, 2020

In addition to the stuff in that paper, we also have "nowait" and "finish" versions of FillBoundary:

https://amrex-codes.github.io/amrex/doxygen/classamrex_1_1FabArray.html#ac2d0be5fcc4f1d0ca3bc90326479e6a4

So if you know you don't need multifab A in function F, you can do

A.FillBoundary_nowait(); // start exchanging ghost cells but don't block
F(); // does not need A
A.FillBoundary_finish(); // block here

// A's ghost cells are now up to date.

In principle particle redistribution and other parallel operations could do the same thing, but we don't expose an API for those right now.

@akashdhruv
Copy link
Contributor

akashdhruv commented Jan 26, 2021 •
edited
Loading

Following up on @atmyers comment: Is the requested fill patching for face centered variables now available within AMReX, or you do you suggest writing our own method for FLASH like MAESTROeX does?

There are essentially two components to this issue that we would like to address:

  1. Interpolation at fine/coarse boundaries described in Vanella and Balsara paper for face center variables. I think this maybe already available based on the discussion in issue MHD Interpolator in FillPatchUtil.cpp/InterpCrseFineBndryEMfield( ) function #397

  2. Extending the Vanella and Balsara scheme for preserving divergence free condition between parent and leaf blocks when AMR grid changes

Note: Even if the Balsara scheme for preserving divergence free condition is available for cell centered variables, that would be a great starting point for some applications within FLASH

@WeiqunZhang
Copy link
Member

A modified version of Vanella et. al.'s scheme has been implemented. What's still missing is a function to fill the physical domain boundaries. What types of physical boundary conditions do you need?

@akashdhruv
Copy link
Contributor

akashdhruv commented Jan 26, 2021 •
edited
Loading

@WeiqunZhang, can you provide more details on what has been implemented and how one can use it within FLASH? Is the divergence free interpolation/extrapolation now available for both situations described above?

In regards to the physical boundary conditions, I have primarily used Grid_bcApplyToRegion and Grid_bcApplytoRegionSpecialized to set physical boundary conditions. I don't know how that has changed with AMReX implementation. Maybe @kweide and @jared321 can provide details on that?

Following are some details on physical boundaries that would be good to have:

  1. Dirichlet boundary conditions: For guard cells that don't coincide with the boundary, a simple interpolation should be fine, for the one face centered guard cell that will coincide with the boundary, the Dirichlet value is directly applied.

  2. Neumann boundary conditions: Similar to Dirichlet, but the guard cell that coincides with the boundary in untouched.

@kweide
Copy link
Contributor

kweide commented Jan 27, 2021

@akidhruv, regarding use of the FLASH subroutines Grid_bcApplyToRegion / Grid_bcApplytoRegionSpecialized to implement boundary conditions: implementations should work with both PARAMESH and AMReX to place the desired values in guards cells. That is they are meant to work; I don't know what actually happens if one tries to apply boundary conditions to face-centered variables in a Grid implementation with Amrex. In particular, I don't know whether the callbacks (part of the FLASH code) that should be invoked by AMReX when its FillPatch functions (or similar) are called have actually been written to correctly (indirectly?) invoke Grid_bcApplyToRegion* on the face-centered cell regions involved.

@WeiqunZhang
Copy link
Member

What we have are here. https://github.com/AMReX-Codes/amrex/blob/development/Src/AmrCore/AMReX_FillPatchUtil.H#L89

It has template parameter typename BC. We still need something like https://github.com/AMReX-Codes/amrex/blob/development/Src/Base/AMReX_PhysBCFunct.H#L115

The FillPatch functions also have pre and post-interpolation hooks that can be used to call FLASH funcitons. But I don't think those alone without something like PhysBCFunct will work.

Also, we usually regard face variables on the domain face to be valid "cell".

@akidhruv Are you face variables velocity or magnetic field or something else?

@akashdhruv
Copy link
Contributor

@WeiqunZhang The face variables are velocities. We solve the fractional method for Navier-Stokes described in Vanella.et.al

Give me a couple of days and I can give you more details regarding mechanics of physical BCs

@kweide
Copy link
Contributor

kweide commented Jan 27, 2021

Also, we usually regard face variables on the domain face to be valid "cell".

So do we for FLASH, "usually". However, there is a parameter gr_pmAlwaysFillFcGcAtDomainBC that can be used to change that behavior. When this is true, the implementation routines for boundary conditions are allowed to modify such face variables on boundaries. (This was introduces specifically for Navier-Stokes purposes.)

@akashdhruv
Copy link
Contributor

akashdhruv commented Jan 28, 2021 •
edited
Loading

@WeiqunZhang, following are the details:

  1. We use face-center variables to store velocity for fractional-step Incompressible Navier-Stokes solver. There are other face center variables, but we need divergence conservation only for velocity.

  2. For physical boundary conditions, there is a callback function within FLASH-X called gr_fillPhysicalBC that is supplied to amrex_fillpatch. Maybe you can provide a stub or a generic function PhysBCFunct for the AMReX library, but FLASH-X won’t be using it.

  3. During guard cell filling we conserve fluxes for velocity that sits directly on block boundaries and require that they remain untouched by amrex_fillpatch. Please see the attached figure for details. I will write a separate method within FLASH that uses AMReX flux registers for this purpose. A flag to enable/disable guard cell filling for velocity that lie at block boundaries will be helpful. We have this option for PARAMESH as mentioned by @kweide in the previous comment.

  4. Face center points that lie outside of block boundaries will need standard guard cell treatment (see figure).

Note that 3 and 4 only apply during guard cell filling. The interaction between leaf and parent blocks during AMR grid change will use the procedure in 1 where all points need to undergo interpolation/extrapolation.

Please let me know if you have any further questions. Also please provide details on how to use the divergence free interpolation selectively for velocity and not other face center variables. Thank you!

Stencil

@kweide
Copy link
Contributor

kweide commented Jan 28, 2021

  1. During guard cell filling we conserve fluxes for velocity that sits directly on block boundaries and require that they remain untouched by amrex_fillpatch.

This is somewhat confusing. I hope the following helps to clarify.

"F ace-centered variables" and "fluxes" are quite different kinds of things. I believe that is true for all of {PARAMESH, AMReX, FLASH}. With that I mean they are different kinds of data structures. Never minds that in some places in the code, the can be used to represent the same kind of physical quantity - for example, velocity.

As far as FLASH is concerned,

  • "flux correction" (sometimes unfortunately called "flux conservation") involved "fluxes". It should only "do something" to fluxes (in the normal direction) on block boundaries.
  • "guard cell filling" involves cell-centered variables, as well as "face-centered variables" if the latter exist in a setup. (In FLASH, "guard cell filling" is responsible for doing this at all of the following kinds of locations: between blocks at the same refinement; between blocks of different resolutions, both at the coarse side (by averaging) and at the fine side (by interpolation); at domain boundaries (by applying boundary conditions (BCs) implemented in callbacks); taking periodic boundaries properly into account.)

"Guard cell filling" is, but default, only responsible for updating values outside of the valid region of blocks. I am using the term valid here in the AMReX sense; in particular, face-centered data located on cell faces that coincide with block boundaries are considered part of the valid data.

The "guard cell filling" operation of FLASH+PARAMESH provides some (optional) features that break the general rule just stated, as far as face-centered variables are concerned:

  • Cell faces on domain boundaries: can be updated by the (generally user-defined) callback that implements the BC. This is controlled by logical parameter gr_pmAlwaysFillFcGcAtDomainBC.
  • Cell faces between blocks: can be updated in order to force values that represent the same physical value to be consistent. At fine-coarse boundaries, this can have essentially the same effect as the application of "flux correction" has to fluxes: namely, coarser-resolution values are effectively replaced by finer-resolution values, because the latter are considered "better". Controlled by FLASH runtime parameter gr_pmrpForceConsistency / PARAMESH parameter force_consistency.

It is not clear to me how important those "optional" behaviors are, for the applications that @akidhruv is working on, and which of them are desireable for FLASH+AMReX. I want to emphasize that they only involve (face) locations to which fine-to-coarse interpolation should never apply.

@akashdhruv
Copy link
Contributor

akashdhruv commented Jan 28, 2021 •
edited
Loading

@kweide @WeiqunZhang for the applications that I am working on, "fluxes" and "face-centered" variables are still treated differently. The point that I was trying to make was that cell face values for velocity between block boundaries need to satisfy a conservation equation that is computed from the "fluxes". At present this has been implemented separately within the physics unit, and will continue to be that way.

After carefully studying the algorithm, I don't think it will interfere with how guard-cell filling is done at block boundaries and there won't be necessity of an "optional" parameter.

For physical boundaries the user defined callback within FLASH should take care of the conservation equation.

@akashdhruv
Copy link
Contributor

@WeiqunZhang Any updates on this issue?

@WeiqunZhang
Copy link
Member

I have not got a chance working on this, so I have asked @kngott to help on this.

@kngott
Copy link
Contributor

kngott commented Mar 15, 2021

I believe we're all set to start testing. You should be able to find the details in the PR descriptions and diffs:

A modified Vanella interpolation has been added and a FillPatchTwoLevels created to launch it and similar face-based functions: #1483
The Fortran interface for the call has just been merged: #1793

Try it out and let me know if you have any questions.

@akashdhruv
Copy link
Contributor

@kngott Awesome! I will test and get back to you.

@kweide
Copy link
Contributor

kweide commented Jan 30, 2022

Is #2452 of any relevance here?

@WeiqunZhang
Copy link
Member

#2452 is different. It uses the coarse data to override the fine data.

@akashdhruv
Copy link
Contributor

akashdhruv commented Apr 27, 2022 •
edited
Loading

I am testing this right now, and running into some issues. Need to discuss with @kweide before I can provide more details. Please don't close this issue

@akashdhruv
Copy link
Contributor

akashdhruv commented Sep 12, 2023 •
edited
Loading

I found a bug with amrex_fillpatch_two_faces when running a test in 2D mode. The multifab array, mf for face variables in each dimension returns identical values.

@akashdhruv
Copy link
Contributor

Here is a figure for reference. It looks like the data from JAXIS multifab is copied to IAXIS multifab when remaking a level.

Presentation1

@WeiqunZhang WeiqunZhang mentioned this issue Sep 13, 2023
Merged
5 tasks
WeiqunZhang added a commit that referenced this issue Sep 14, 2023
@WeiqunZhang
Fix Fortran interface for FillPatch for face variables (#3541)
5829296 
## Summary

The stride was wrong when copying MultiFab pointers from a Fortran array
to a Vector of Array of MultiFab pointers.

This also fixes cases when the user passes t_new followed by t_old,
instead of the other order. The issue was then teps was negative. This
was a minor issue, because the FillPatch function in C++ does not care
about the order.

## Additional background

#327 (comment)

## Checklist

The proposed changes:
- [x] fix a bug or incorrect behavior in AMReX
- [ ] add new capabilities to AMReX
- [ ] changes answers in the test suite to more than roundoff level
- [ ] are likely to significantly affect the results of downstream AMReX
users
- [ ] include documentation in the code and/or rst files, if appropriate
@akashdhruv akashdhruv mentioned this issue Sep 18, 2023
Merged
@akashdhruv
Copy link
Contributor

akashdhruv commented Sep 22, 2023 •
edited
Loading

I think this issue can be closed. There are still some issues with divergence free interpolation that I am investigating here:

https://github.com/Lab-Notebooks/AMReX-Divfree

When I have updates I will re-open the issue.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees
No one assigned
Labels
None yet
Projects
None yet
Milestone
No milestone
Development

No branches or pull requests
10 participants
@darylbond @drummerdoc @kweide @atmyers @akashdhruv @WeiqunZhang @asalmgren @chaw0023 @jared321 @kngott