How to get a multifab to interpolate to cell faces at domain boundaries

[tpg2114]

I'm struggling to get my boundary conditions implemented on domain boundaries and am looking for some advice on how to get it to work. Our code is cell centered and when I create a new refinement level at the domain boundaries, I want to set the ghost cell center of the new level such that the cell face on the domain boundary has a value interpolated from the coarser level. I am using the Fortran interface and have tried something like this in my function to remake a level:

! All my lo_bc and hi_bc are set to amrex_bc_ext_dir
! Build my cell centered new level -- the ghost cells on domain boundaries will be unset because my physical BC routine doesn't do anything at this stage.
CALL amrex_fillpatch(new_vars, <long list of arguments using the coarser level and current level>)

! Create a multifab for the I-face-centered values, ultimately I only want the I-face coincident with the domain boundary
CALL amrex_multifab_build(new_vars_Ifaces, boxarray, distromap, nvars, nghosts, (/ .TRUE., .FALSE., .FALSE. /))

CALL amrex_fillpatch(new_vars_Ifaces, <long list of arguments using the coarser level and current level cell-centered data>)

If I then look at the solution value in new_vars_Ifaces on the face coincident with the domain boundary, it looks like it just took the value from the first interior cell center and copied it to the face value. What I wanted is that face value to include both sides of the face, something like 0.5*(ghost_cell + first interior cell). However, it looks like the amrex_fillpatch routines ignore ghost cells when setting the face values.

Is there some way to get this to work the way I expected to, so the amrex_fillpatch will do the interpolation to the cell faces for all cells, regardless of where the domain boundary is?

Alternatively, I'm open to frame-challenges -- am I going about this the wrong way? My boundary condition ghost cell needs to be set based on the face value and the first interior cell, but I can't figure out how to get the face value in a straightforward way.

Thanks,

Tim

[asalmgren]

The ext_dir type of boundary condition treats the Dirichlet values as living on faces.

…

On Tue, Nov 9, 2021 at 9:32 AM tpg2114 ***@***.***> wrote: I'm struggling to get my boundary conditions implemented on domain boundaries and am looking for some advice on how to get it to work. Our code is cell centered and when I create a new refinement level at the domain boundaries, I want to set the ghost cell center of the new level such that the cell face on the domain boundary has a value interpolated from the coarser level. I am using the Fortran interface and have tried something like this in my function to remake a level: ` ! All my lo_bc and hi_bc are set to amrex_bc_ext_dir ! Build my cell centered new level -- the ghost cells on domain boundaries will be unset because my physical BC routine doesn't do anything at this stage. CALL amrex_fillpatch(new_vars, ) ! Create a multifab for the I-face-centered values, ultimately I only want the I-face coincident with the domain boundary CALL amrex_multifab_build(new_vars_Ifaces, boxarray, distromap, nvars, nghosts, (/ .TRUE., .FALSE., .FALSE. /)) CALL amrex_fillpatch(new_vars_Ifaces, ) ` If I then look at the solution value in new_vars_Ifaces on the face coincident with the domain boundary, it looks like it just took the value from the first interior cell center and copied it to the face value. What I wanted is that face value to include both sides of the face, something like 0.5*(ghost_cell + first interior cell). However, it looks like the amrex_fillpatch routines ignore ghost cells when setting the face values. Is there some way to get this to work the way I expected to, so the amrex_fillpatch will do the interpolation to the cell faces for *all* cells, regardless of where the domain boundary is? Alternatively, I'm open to frame-challenges -- am I going about this the wrong way? My boundary condition ghost cell needs to be set based on the face value and the first interior cell, but I can't figure out how to get the face value in a straightforward way. Thanks, Tim — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub <#2468>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ACRE6YR46WFVOGKBMS263O3ULFLLZANCNFSM5HV3S3FQ> . Triage notifications on the go with GitHub Mobile for iOS <https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675> or Android <https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub>.

-- Ann Almgren Senior Scientist; CCSE Group Lead Pronouns: she/her/hers

[tpg2114]

So if I want amrex_fillpatch to compute I-face-centered values (e.g. nodal=(/ T, F, F /)) from a cell-centered multifab at all I-faces including those that happen at domain boundaries and/or between all types of ghost cells, what type of boundary condition do I need to tell AMReX to apply in that case?

I was under the impression ext_dir was just a no-op BC, but all of the calls to amrex_fillpatch refuse to do anything in my domain boundary ghost cells.

If we forget for the moment that this is to apply a boundary condition specifically, is there any way to get AMReX to reconstruct to faces for all cells in a multifab?

[drummerdoc]

It is a little confusing. So maybe a second pass with different phrasing might help. In AMReX applications, ext_dir, indicates that the user will provide a boundary value that is to be applied on the boundary face. The data structure that is used to communicate this information to the code is actually cell-centered - that's the confusing part - but there is a 1:1 correspondence between the data on each face and the container that holds that data. In the guts of the advection and diffusion operators, for example, the stencils are indeed modified to account for the centering shift of this point. Note that the other boundary condition types (e.g., reflect_even, foextrap) all fill data at the center of the grow cell outside the domain, and the stencils for the AMReX operators know that as well. So, after a fill-patch operation, the way that the data that lives in grow cells is used depends still on the BC type.

I hope that helps.

[tpg2114]

Okay, that helps slightly but I think my use-case here might be atypical... so maybe I can fill in some more context for why I'm trying to do what I'm stuck on.

I have a legacy, structured multiblock code that is not based on AMR or AMReX. I have bolted in AMReX, but all of our numerics and all of our physics (including BCs) happen with our legacy code. I am basically using AMReX to handle the gridding, refinement, and communication -- I copy data from my legacy structure into AMReX, do the communication, copy back to my legacy structure, and once in awhile call amrex_regrid. If AMReX does regrid, I tear down all my out-of-date multiblock structure (the data still lives since it was copied to AMReX), then I rebuild my new multiblock structure and copy data back over from the AMReX and then advance. Each box in AMReX is one block in my legacy structure.

So ultimately, the only boundary condition I need AMReX to treat on its own is periodic. Everything else, I just want AMReX to not mess up the ghost cells because those ghost cells have data copied into them from our legacy structures and those have the BC values I want. So amrex_bc_ext_dir worked, since it didn't touch my ghost cell values.

For most boundary conditions, this is no problem because our legacy code sets the ghost value based on the interior (walls, traditional Dirichlet or Neumann). So zeros or junk values in ghost cells will just be overwritten by proper values.

However, our implementation of the Navier-Stokes Characteristic BCs is based on applying a source term to the ghost cell based on the type/reflectivity of the boundary condition. It's essentially like a control problem rather than a typical Dirichlet or Neumann BC. What this means is that the value in the ghost cell has a time-history and it isn't as simple as just setting it based on the current state of the interior -- it's value is the time-integrated history of the initial value plus all the source terms added so far.

I'm getting in trouble now because I can't get the ghost cells set correctly when I am creating a new level or remaking a level when the new boxes contain the domain ghost cells of a subsonic boundary condition. Let's forget about grid refinement for a minute because I run into problems just on remaking a level where processors swap existing boxes.

What I need to do is find what the face value is on the domain boundary (if I were writing this by hand, I would do U_face = 0.5*(U_ghost + U_interior) on the existing multifab, and then I can compute the ghost cell value on the new multifab as U_ghost = 2*U_face - U_interior. I could hand-roll all the math and communication needed to do that I suppose, by first computing that U_face value in all existing boxes that hold a domain boundary, then by finding out which existing boxes need to talk to which new boxes, and then figure out all the processors involved and build up a message buffer.

But hand-rolling that seems overly complicated and that sounds like something the multifab should support. What I really want is to take an existing multifab and have it compute the I-face-centered values, then transfer that to the new multifab. It doesn't need to know or care about boundary conditions, or where those faces are physically -- I want it to do that cell-center-to-cell-face operation at every cell regardless of those things. Just a purely mathematical operation on the cells that are there.

Is this something that can be done or do I need to hand-roll things? I was looking again at the documentation for the BC types and maybe I can (ab)use the amrex_bc_int_dir setting? I thought that was only for periodic boundary conditions since that's how all the examples use it, but the documentation does say "including periodic boundary conditions" so perhaps it would work for what I want so long as periodicity is 0?

[drummerdoc]

Hmmm, interesting on several counts.

I wonder though if instead of bastardizing the BC stuff for your purposes, if you instead created the BoxArray that defines the valid region to include the bits grown over the physical boundary. Then you could simply treat the AMReX data structure as a data structure and not as something more complicated using fill-patch stuff where you have to side-step the BCs. You could copy into it and from it and not worry about masking out or cleverly undoing what was done. Does this make sense?

As for the NSCBC, a separate and uncalled for question that comes from the fact that I also need to tackle these kinds of boundaries...do you have the problem that accumulating stuff since the IC can accumulate too many errors? I'm curious to learn more about that part, if you have a paper or something...again, unrelated to your actual question here, but I'm an opportunist :)

[tpg2114]

Okay, so I was poking around in the Fortran interface files to see if I could find anything useful to help me and I found amrex_average_cellcenters_to_face in amrex_multifabutil_mod... I had no idea that subroutine existed, but it appears to do what I was trying to do -- average the left and right cells centers of every cell, no matter where it appears. It takes geom as an argument but ignores it in 3D.

I did have to re-wrap it to add an ncomp argument to it, because the default value is 1 and the Fortran interface doesn't expose that option. So I can now get my domain-face values on my existing levels/boxarrays. Now I need to figure out how to project that onto the new boxarrays being created for the new level (or remade existing level).

If I used fillpatch I will be right back at the physical boundary condition issue I had before. I'll dig through what fillpatch does internally and see if it's something I can reproduce in a fillpatch_ignore_bcs or something that just takes the two levels of multifabs and fills in the values of a new multifab for every cell including all ghosts...

[drummerdoc]

Sure! Thanks. Marc.Day@nrel.gov

[tpg2114]

Alright, I appreciate the tips and suggestions @drummerdoc and @asalmgren -- I managed to get it to work the way I wanted! On the off chance anybody else has a similar desire, here's a summary for what I worked out:

• The following works as expected using both amrex_bc_int_dir and amrex_bc_ext_dir, so I don't quite appreciate what those really do I suppose. But it probably works because my fillphys_bc is basically a no-op and I set my domain ghost cells other ways.
• I found amrex_average_cellcenter_to_face and traced it through to the implementation where it does the $U_F = 0.5(U_L + U_R)$ that I was looking for. I saw it also ignores the geometry passed in, so I suspected it wouldn't care about domain boundaries or ghost cells or the like -- turns out that hunch was correct.
• The function to do the averaging defaults to ncomp=1 but the Fortran interface does not expose the ncomp argument. So I re-wrapped it in the sources in my code to add the ncomp argument. Here is a quick patch I knocked up that adds it as an optional argument to the AMReX Fortran interface (so it should be backwards compatible) if it is of interest.
• With that, a very rough pseudo-code of my remake grid looks like (and new grid from scratch is similar -- both were originally identical to the tutorials):

if remaking level > 0:
    fill patch of vars on remade level using vars on old level-1 and old level

    if not all periodic:
        build multifabs for the I, J, K face-centered vars on old level-1
        build multifabs for the I, J, K face-centered vars on old level
        build multifabs for the I, J, K face-centered vars on remade level

        average_cellcenter_to_face on old level-1
        average_cellcenter_to_face on old level

        fill patch of I, J, K face-centered vars on remade level using face-centered vars on old level-1 and old level

        allow multiple mfiters
        loop over the remade cell-centered mfiter and remade face-centered mfiter
             if on non-periodic domain boundary, reconstruct ghost cells using remade cell-centered and remade face-centered values

        clean up
     replace old level with remade level

If there are other/better ways to accomplish the same, I'm all ears!

---

Of course, proof is in the pudding -- the test case is a left-ward traveling acoustic pulse with a quiescent background with a perfectly-reflective, Navier-Stokes Characteristic Boundary Condition inflow. The image is the flow rate at multiple snapshots in time, with the AMReX-powered AMR result and our traditional, structured multiblock result plotted together. The AMR has 3 levels and dynamically tracks the pulse, with the finest level resolution equal to the structured multiblock resolution (which is uniform).

It may not look like it, but each of the pulses is actually 2 lines -- one for AMR, one for structured multiblock, but it's impossible to tell them apart. I consider that a win!

---

Lastly, since it is tangential but was asked -- the NSCBC are very very sensitive to error accumulation. In fact, I discovered that if the initial velocity in the inflow ghost cell is -1e-11 instead of 0.0, the green lines in the image (flow rate when the pulse has bounced off the inflow) peaks at -0.05 instead of -0.01. In other words, there is a 500% error at that point in time due to a 1e-11 difference in one cell at initialization.

We also discovered about 5 years ago that floating point truncation errors can make or break the BC. If we do not compile the boundary condition routines with strict IEEE 754 compliance with double precision (-fp-model strict with the Intel compiler), the boundary conditions will drift and frequently drift far enough to diverge the simulation entirely. I didn't track it down precisely, but I suspect what is happening is that individual terms have very large magnitudes but their addition or difference are not that big. However, fast-math approximations or loose compliance with underflow and truncation rules lead those errors to build to the point that the answer is either completely wrong or the simulation is unrealizable (negative temperature or pressure). This happens for simple flows -- vortex or turbulence convection for example.

If it is as I suspect, there are probably ways to re-write the code such that the individual terms are not computed directly but their sum or difference is -- however, given that developer and SME time is far far more expensive than CPU time, we opted for code clarity/fidelity to the equations as written and eat the performance penalty of strict IEEE 754 compliance.