Force vector as a function of scalar gradients?
Closed this issue · 11 comments
Hi,
I am trying to add a body force to the momentum equation, which is proportional to the gradients of an advected scalar (in this case Tracer). What I have noticed is that the getForce() function is called several times in each time step, across several files, depending on the procedure (scalar advection, velocity advection, initial time step estimate etc), with each function call using different ngrow, scomp, scalScomp etc.
I assume that I only need to add this body force in the functions where getForce() is called relating to the velocity updates. Hence at these points, I have tried to make a new MultiFab for the scalars that has more ghost cells than the force vector, e.g. in NavierStokesBase::initial_velocity_diffusion_update
MultiFab& U_old = get_old_data(State_Type);
MultiFab UOldBorder(grids, dmap, NUM_STATE, 4);
MultiFab::Copy(UOldBorder, U_old, 0, 0, NUM_STATE, U_old.nGrow());
UOldBorder.FillBoundary(geom.periodicity());
Then MFIter-ating over UOldBorder as the Scal argument for getForce() . Then in FORT_MAKEFORCE:
do k = f_lo(3), f_hi(3)
do j = f_lo(2), f_hi(2)
do i = f_lo(1), f_hi(1)
dTdx = (scal(i+1,j,k,nTracScal)-scal(i-1,j,k,nTracScal))/(2*hx)
dTdy = (scal(i,j+1,k,nTracScal)-scal(i,j-1,k,nTracScal))/(2*hy)
force(i,j,k,nXvel) = someFunctionOf(dTdx)
force(i,j,k,nYvel) = someFunctionOf(dTdy)
enddo
enddo
enddo
When I try this I get strange behaviour, e.g. the pressure field is divided into square patches with the surrounding cells of the squares being wrong values. This may indicate that the way I am filling ghost cells or indexing the arrays is wrong, but I use a similar procedure for tagging based on density gradients which works okay (this may just be something to debug on my end).
I am somewhat confused by the various different places the function is called and the different values of scomp, ncomp etc., so any suggestions would be much appreciated. Is there something I need to do in addition to what I have described above in order to obtain the correct forcing, say changing some other arguments of getForce? Is there an example of a similar method in any other branches/AMReX code? Thanks
Quick question: have you confirmed that the grow cells of scale have been properly filled to allow you gradient calc?
Thanks for your response.
I have tried to confirm this quickly, and found that the ghost cells have not been filled correctly (with the density being negative even in places), seeing the below output from getForceVerbose:
A - Predict velocity:
Calling getForce...
NavierStokesBase::getForce(): Entered...
time = 0
scomp = 0
ncomp = 2
ngrow = 1
scalScomp = 0
Doing all components
NavierStokesBase::getForce(): Force Domain:
(-1,-1) - (128,128)
NavierStokesBase::getForce(): Vel Domain:
(-3,-3) - (130,130)
NavierStokesBase::getForce(): Scal Domain:
(-4,-4) - (131,131)
Vel 0 min/max 0 / 0
Vel 1 min/max 0 / 0
Scal 0 min/max 0 / 0.9999999986
Scal 1 min/max 0 / 0.9999999971
...
B - velocity advection:
Calling getForce...
NavierStokesBase::getForce(): Entered...
time = 0
scomp = 0
ncomp = 2
ngrow = 1
scalScomp = 0
Doing all components
NavierStokesBase::getForce(): Force Domain:
(-1,-1) - (128,128)
NavierStokesBase::getForce(): Vel Domain:
(-3,-3) - (130,130)
NavierStokesBase::getForce(): Scal Domain:
(-4,-4) - (131,131)
Vel 0 min/max 0 / 0
Vel 1 min/max 0 / 0
Scal 0 min/max -0.04675671551 / 0.04582436846
Scal 1 min/max -0.03736710405 / 0.0439111718
Is what I have done above to fill out the ghost cells incorrect? I used mf.FillBoundary(geom.periodicity()) as in the documentation.
Two thoughts, one is pure speculation...
- For the Godunov integrator, getForce is needed both to predict the time-centered state to face centers, and to evaluate the (conservative) update. For method-of-lines (used with EB), it is not needed to get face values because no extrapolation in time is done.
- (speculation) Internal to the code, the advective updates are done over tiles of the grid, so the box outlines are likely due to these tiles. You need that no matter how the work is divided that you get the same answer, meaning in your case that you need the gradient correctly evaluated in the grow cells needed for the forcing, and since you are using a second-order centered gradient, you will need yet one more valid cell for the scalar, which you may not have.
If 2 is correct, you need to find a way to get a bigger patch of data into the routine so that you can properly evaluate the gradient, but as you point out, there are so many calls to getForce and you don't have (easy) control over how the data is processed before calling it. Here's what I would do (warning, it's a bit of a hack)...
I would make a new MultiFab owned by the NavierStokes class to hold the gradient, and I'd fill it properly with as many grow cells as required (see below). Then, inside of getForce, if I was filling forces for velocity, I'd magically get data from that structure to construct the force. If that worked, we can talk about "proper" ways to do this.
To make the gradient vector "properly", things can get a bit messy because you need to account for physical boundary conditions, as well as fine-fine and coarse-fine. The "easiest" way to avoid most of the problems is to use FillPatch to create a version of the scalar out to nGrow+1 cells and then apply the gradient to that, but modify your centered stencil to account for the fact the ghost cell values at EXT_DIR boundaries really live on faces. Short of that detail, the FillPatch operation will internally take care of the FillBoundary(geom.Periodicity()) bit, but also the interpolation at coarse-fine that you will need done.
Also... if you try my approach, you have to remember to fill the gradient multifab prior to calling routines that will call getForce ... that's the hacky part.
If this fixes your issue, "we" (where I'm not sure who is included in the "we") can think about whether it makes sense to modify the routines to take a MultiFab for the force. It would change a lot of software, but it seems like the way its currently written really places quite a burden on the "user" of the code that wants to add their own custom forcing.
The alternative here is to derive a new class from NavierStokes to do "your own physics" where you can over-write all the functions from NavierStokes that you don't like. But that sure seems like a really big hammer for small nail.
Thanks for this.
Just to confirm a couple of things:
In order to get the correct forcing for the velocity, I would have to compute the scalar gradients in each of these functions:
NavierStokesBase::estTimeStepNavierStokesBase::velocity_advectionNavierStokesBase::velocity_advection_updateNavierStokesBase::initial_velocity_diffusion_updateNavierStokes::predict_velocityMacProj::mac_sync_compute
As per your suggestion, what I am doing is creating a MultiFab, scalGrad, with the same number of ghost cells as the state vector. Then I create a state vector with more ghost cells like so:
MultiFab& U_new = get_new_data(State_Type);
MultiFab SBorder(grids, dmap, NUM_STATE,U_new.nGrow()+2);
MultiFab scalGrad(grids, dmap, 2, U_new.nGrow());
MultiFab::Copy(SBorder, U_new, 0, 0, NUM_STATE, U_new.nGrow());
FillPatch(*this, SBorder, U_new.nGrow()+2, state[State_Type].curTime(), State_Type, 0, NUM_STATE);
Then I have a function, computeScalarGradient that takes in the bx of the state vector and fills scalGrad based on central differences. scalGrad[mfi] is then passed as an additional argument to getForce, and the gradients are accessed at the corresponding index of the force vector.
If this is what you had in mind (feel free to iterate/correct) then I will focus on the development (I anticipate that this will not be trivial to debug).
Not quite. U_new in your example likely has 1 grow cell, but the advection probably needs 4 (check me on that...). Typically then, FillPatch would be called to fill SBorder to 4 grow cells. The force might be needed out 4 grow cells as well (again, need to check the requirements of the function). If so, you will then need to allocate scalGrad with 4 grow cells. However, to fill scalGrad including its 4 grow cells you will need the scalar for it to have 5 grow cells.
The main error above is that you base the number of grow cells for scalGrad on the number of grow cells that have been allocated in the state MultiFab, but that is not directly associated with the number required for the operator. FillPatchIterator is used to load up valid+grow cells, and works even when the number required is larger than what is in the state because it allocates its own temporary to hold the result. FillPatch copies from that temporary into the container you pass through the argument list (again, not the state).
Even if you find that you only need, say, 2 grow cells for the force, this 2 is 2 beyond the size of the final result from the function (e.g., hyperbolic flux divergence). It's NOT the number of grow cells in the state + 1 (well, it might be, but that's a coincidence and the code would be misleading). For example, if you wanted to extend your code to incorporate embedded boundaries, you may need to add another 2 grow cells....independent of how many are in the state.
I hope this helps
Okay, thanks for pointing this out.
So to try this again, the code would be something like:
MultiFab& U_new = get_new_data(State_Type);
MultiFab SBorder(grids, dmap, NUM_STATE,5);
MultiFab scalGrad(grids, dmap, 2, 4);
MultiFab::Copy(SBorder, U_new, 0, 0, NUM_STATE, U_new.nGrow());
FillPatch(*this, SBorder, 5, state[State_Type].curTime(), State_Type, 0, NUM_STATE);
Then, before getForce is called, there would have to be another tiling loop to compute the scalar gradients based on the tilebox of the MFIters of scalGrad, to fill the MultiFab with data based on SBorder?
Sounds good.
Sorry to bother you again, I just need to clarify a couple of things once more in light of some issues I am having.
This relates again to the different procedures in which the forcing is called.
Take for example F - velocity advection update (half time). The scalar Fab used in getForce is the average of the old and new time states. I have assumed that I also have to find scalar gradients at the old and new time states and average them, like below:
for(MFIter SGmfi(scalOldGradFab,true); SGmfi.isValid(); ++SGmfi)
{
const Box& bx = SGmfi.tilebox();
computeScalarGradient(SOldBorder[SGmfi], scalOldGradFab[SGmfi], bx, dx);
}
for(MFIter SGmfi(scalNewGradFab,true); SGmfi.isValid(); ++SGmfi)
{
const Box& bx = SGmfi.tilebox();
computeScalarGradient(SNewBorder[SGmfi], scalNewGradFab[SGmfi], bx, dx);
}
MultiFab scalGradHalf(grids, dmap, 3,4);
scalGradHalf.setVal(0);
MultiFab::LinComb(scalGradHalf, 0.5, scalNewGradFab, 0, 0.5, scalOldGradFab, 0, 0, 3, 4);
And how about other instances, such as in A - predict velocity where you have the state MultiFabs used in getForce derived from FillPatchIterators:
FillPatchIterator U_fpi(*this,visc_terms,Godunov::hypgrow(),prev_time,State_Type,Xvel,BL_SPACEDIM);
MultiFab& Umf=U_fpi.get_mf();
...
FillPatchIterator S_fpi(*this,visc_terms,1,prev_time,State_Type,Density,NUM_SCALARS);
MultiFab& Smf=S_fpi.get_mf();
Do I compute gradients based on the data from Umf, Smf, or get_old_data(State_Type)?
I ask because I seem to get differing answers for values within my scalar gradient MultiFab depending on the particular stage of the force calculation (which is one of the hazards of computing the same variable multiple times I guess) - in particular the above steps (A and F) are when things seem to go awry. Thanks
For your "A", you would use Smf, except that you have only asked for 1 grow cell to be filled for some reason. It was my understanding that you needed enough grow cells in the scalar to be able to compute Godunov::hypgrow() grow cells of it's gradient, meaning that you need to create Smf with Godunov::hypgrow()+1 grow cells, then evaluate the scalGrad with that.
The other question is a bit subtle. You need scalNewGrad over the entire box, which means you need scalNew out to one grow cell. This leads to requirements that need to be pulled through everything. So, if you needed 4 grow cells to get the hyperbolic terms filled with no grow cells, but are now added a forcing that needs the final result out 1 grow cell, you then need 6 grow cells - one extra to end up with 1 instead of 0 grow cells in the result, and another because you need to compute the gradient of the scalar in order to evaluate the force in the predictor.
If you incorporating EB, you need yet another 2 grow cells.