MAC velocity not getting filled correctly at the low boundary about amr-wind HOT 18 CLOSED

Passing in IndexType would require changes in amrex because that's where the function signature/call for amrex::PhysBCFunct<amrex::GpuBndryFuncFab<Functor>> physbc is happening...

from amr-wind.

Okay, that makes complete sense to me, thanks for explaining and for working on this. :)

I think the MAC velocities for valid cells should be computed fine with all the other special handling that I have put in, unless I'm missing something. Once we have this issue sorted, I will have to figure out how to apply my custom Neumann fills using the MassInflowOutflowBC operator to the MAC velocity ghost cells.

from amr-wind.

You called it, Marc. I confirmed the obvious, that indices with i = 0 are getting passed to f_user(). However, f_user is not what goes to the Rankine.H operation directly. First, it has to go through the Functor that initialized the GpuBndryFuncFab, which, for this case, is "DirichletOp" in FieldBCOps.H in AMR-Wind. This operator is what checks whether the field at the index should be filled or not, using a bool "is_bndry". This check assumes a cell-centered approach, which is why i = 0 is ignored but i = 40 isn't. After this check, either inflow_op (which would be the UDF in this case) or wall_op is called.

from amr-wind.

yes! I just found the same. The offending line is https://github.com/Exawind/amr-wind/blob/main/amr-wind/core/FieldBCOps.H#L172. This is fine for cell-centered, not fine for face-centered. Phew. Now to think of the right way to fix.

from amr-wind.

Addressed in #1093

from amr-wind.

I change the rankine.inp file to prescribe mass_inflow at xhi too, just for checking how it acts:

xlo.type = "mass_inflow"
xlo.density = 1.0
xlo.velocity.inflow_type = Rankine
xlo.temperature = 300.0
#xhi.type = "pressure_outflow"
xhi.type = "mass_inflow"
xhi.density = 1.0
xhi.velocity.inflow_type = Rankine
xhi.temperature = 300.0

I insert explicit checks within the Rankine.H UDF operator():

            amrex::IntVect ivlo(AMREX_D_DECL(0, 12, 1));
            amrex::IntVect ivhi(AMREX_D_DECL(40, 12, 1));
            if (iv == ivlo) amrex::Print() << "filling " << vel[orig_comp + comp] << " at i=0" << std::endl;
            if (iv == ivhi) amrex::Print() << "filling " << vel[orig_comp + comp] << " at i=40" << std::endl;

The output shows filling the i=40 values but not i=0. The domain has 40 cells in the x-direction and so both of them should be treated similarly as they lie on the boundary:

  MAC_projection                 9         0.01766980429       1.787159233e-14
Fill mac velocities using velocity BCs
filling 10.14358888 at i=40
filling 0.5357916479 at i=40
filling 0 at i=40
  temperature_solve              2       1.221647966e-05        1.98951966e-12

from amr-wind.

I can reproduce this. Thank you for posting the detailed information. Did you get far enough to determine a potential fix or reason for this? I will be looking at this later this week.

from amr-wind.

One thing that bugs me about this UDF is that it is assuming that it is filling a cell-centered field (e.g. the way it computes x: https://github.com/Exawind/amr-wind/blob/main/amr-wind/physics/udfs/Rankine.H#L36). So there will be a difference when you look at the value it puts there. It looks like it is filling -1 and 40:

filling 10.73181747 at i=-1
filling 1.016180833 at i=-1
filling 0 at i=-1
filling 10.14358888 at i=40
filling 0.5357916479 at i=40
filling 0 at i=40

It's not yet clear to me why it does that. Ultimately this is being called from amrex::FillPatchSingleLevel with nghost = 1 1 1, a face centered field (the mac vel), and a physbc we give it. AMReX handles the functions calls between that when we drop into physbc (our UDF). Maybe we are using this call wrong? Maybe this all makes sense to someone else but I am still digging around.

Here's some useful context for all this (when I first had a go at this madness):

• PR: #730
• PR: #743
• Issue: #735
• Issue: #174

Given the number of PRs/commits/issues etc that I went through when dealing with this stuff, I clearly lost my mind (at the time for sure, and probably still now).

from amr-wind.

One thing that bugs me about this UDF is that it is assuming that it is filling a cell-centered field

Yes, that was going to be my next concern after making sure that the correct indices are being filled. I was wondering if there is a way to pass in the IndexType of the velocity being filled in the custom velocity UDF so that the coordinates can be calculated correctly depending upon whether it is a cell-centered or face-centered velocity.

It looks like it is filling -1 and 40

Right, it fills only -1 on the lower boundary and fills both 40 and 41 on the higher boundary. But my understanding is that it should also fill 0 because that lies at the wall like 40.

Thank you for providing the additional context in terms of past work and how we seem to doing the correct thing, I have not dug around further today, I will try to make sense of it tomorrow, though it may take more time for me as I'm less familiar with this part of AW/amrex.

from amr-wind.

Okay, I think it is better to consider that feature later since it is not immediately relevant, and to focus on why i=0 is not being filled.

from amr-wind.

On that topic, I think that the unexpected behavior is not that i=0 is not being filled; rather, it's that i=40 is being filled. The boundary functions are meant to fill ghost cells, and i=0 is not a ghost cell; it's a valid cell. For a quantity on an x-face, i=40 is a valid cell, yet it's being filled. This is checked fairly easily by doing setBndry() for umac -- that operation does not touch either i=0 or i=40 values, but it does modify i=-1 and i=41.

I can see that you might want to change these valid boundary cells for the sake of your new BC, which is fine. But my understanding of the AMR-Wind implementation is to fill the ghost cells of the cell-centered velocity with correct values, and then the MAC velocities at the boundaries are calculated from those cell-centered velocities. fillpatch_sibling_fields is meant to populate the ghost cells of the MAC velocities, but not the ones at the valid boundaries.

Marc and I are getting close to understanding what needs to be done to get the correct behavior.

from amr-wind.

100% with Michael here. The curiosity is why i=40 is being filled.

  MAC_projection                 9         0.01766980429       1.791708441e-14
ghost cells here: (1,1,1) for idim: 0
--------------------------- filfcface call starts
Going to look at this tmp box: ((-1,0,0) (0,59,3) (1,0,0)) with b: ((-43,0,0) (0,59,3) (1,0,0)) and bx: ((-1,-1,-1) (41,60,4) (1,0,0))
filling 10.73181747 at i=-1
Going to look at this tmp box: ((0,-1,0) (40,-1,3) (1,0,0)) with b: ((0,-62,0) (40,-1,3) (1,0,0)) and bx: ((-1,-1,-1) (41,60,4) (1,0,0))
Going to look at this tmp box: ((0,0,-1) (40,59,-1) (1,0,0)) with b: ((0,0,-6) (40,59,-1) (1,0,0)) and bx: ((-1,-1,-1) (41,60,4) (1,0,0))
Going to look at this tmp box: ((40,0,0) (41,59,3) (1,0,0)) with b: ((40,0,0) (83,59,3) (1,0,0)) and bx: ((-1,-1,-1) (41,60,4) (1,0,0))
filling 10.14358888 at i=40
filling 10.13957005 at i=41
Going to look at this tmp box: ((0,60,0) (40,60,3) (1,0,0)) with b: ((0,60,0) (40,121,3) (1,0,0)) and bx: ((-1,-1,-1) (41,60,4) (1,0,0))
Going to look at this tmp box: ((0,0,4) (40,59,4) (1,0,0)) with b: ((0,0,4) (40,59,9) (1,0,0)) and bx: ((-1,-1,-1) (41,60,4) (1,0,0))
--------------------------- filfcface call ends

the Box tmp = b & bx; in AMReX_PhysBCFunct.H is the curious operation. For a face centered field with 1 ghost cell, we would expect i=-1 and i=41 to be filled. But the & is also including i=40. Really we would want extra checks in the UDF function: if bx.ixType() is a face and i > ihi (+1?) then fill it. You can see the fancy footwork (million if conditions) that are being done in FilfcFace in AMReX_FilFC_3D_C.H. We would want to mimic those I think. But we don't have access to bx or the ixType so we can't even begin to use those... not sure the solution to that. In any event, it doesn't much matter because the mac vel ghost cells are not needed (other than for BDS and multilevel Godunov).

All that being said. If you want i=0 or i=40 for the mac vel filled, you are going to need a different strategy... Right now I am not sure what it looks like. Other than an explicit new function call that loops on those faces you want to fill. Maybe right after that fillpatch call?

from amr-wind.

Just a comment re what mkuhn said above -- even when we fill ghost cells for cell-centered velocity there is an "understanding" in amrex that the value in that ghost cell actually represents the value on the face (even though we store it as a ghost "cell"). I wasn't sure if what you wrote above agreed or disagreed with that so just thought I'd put that here to be sure.

from amr-wind.

If you want i=0 or i=40 for the mac vel filled, you are going to need a different strategy... Right now I am not sure what it looks like. Other than an explicit new function call that loops on those faces you want to fill.

At this point, I just want symmetric behavior between the low and high boundaries. So if i=40 not being filled is the correct behavior numerically then I'd be good to go.

from amr-wind.

@asalmgren yup that's understood. Though I bet amr-wind does this inconsistently (this Rankine UDF is one of those that doesn't do that bit right I think). Going through all the UDFs that have been added in the past and making sure that they fit this assumption needs to be done. But I think it is somewhat tangential to this issue right now.

I pinged Weiqun to see if he has ideas for the i=40 index. These indices and the function calls come to us from amrex so we are a bit boxed in to what we can do/have access to.

from amr-wind.

I've gone back through this thread more carefully -- I think the missing piece here (and I'm 90% sure but not 100% sure that everything I say below is correct) is that the physbc is defined by the user and has to be appropriate for the quantity being filled. So it's not amrex's job to move data to the right place, it's the user's job when providing a physbc to make sure that physbc matches the data type. In other words, you shouldn't be using the same physbc to do a fillpatch on the cell-centered velocity as on the mac velocity

from amr-wind.

Yes I think that's about right and the same conclusion I reached last night. It would be good to achieve this without a ton of duplicate code in amr-wind. Going to have to think on this one...

My immediate concern though is with this:

Going to look at this tmp box: ((-1,0,0) (0,59,3) (1,0,0)) with b: ((-43,0,0) (0,59,3) (1,0,0)) and bx: ((-1,-1,-1) (41,60,4) (1,0,0))
filling 10.73181747 at i=-1

This box ((-1,0,0) (0,59,3) (1,0,0)) contains 0 so why does it say that it is only filling -1? I think amr-wind must be intercepting something?

from amr-wind.

Right now my thinking is this:
In FieldFillPatchOps.H:

Functor bc_functor() { return m_op(); }
Functor face_bc_functor() { return m_op.face(); }.  ----> add this and then use `face_bc_functor` in fillpatch_sibling_fields

then in FieldBCOps:

inline FaceFunctorType face() const
    {
        return FunctorType{
            m_field, m_inflow_op.device_instance_face(),
            m_wall_op.device_instance()};
    }

Though we might want using FaceFunctorType = DirichletOpFace<InflowOpType, WallOpType>; to get the right () operator that does the right is_bndry.

then in all the UDFs:

struct DeviceOpFace 2 refs|1 var
{
... something that looks like DeviceOp. Maybe it inherits it so that it is similar? Would need to make sure the members are initialized properly too and are the same as DeviceOp.
}

DeviceType device_instance_face() const { return m_op_face; }

I need to sit on this and think this through before coding anything. I am sharing this in case others have thoughts/see ways to simplify.

from amr-wind.