monitor functions

[dmingram66]

First of all can I apologise for accidentally posting this as an issue #5505,

I want to write a function which can interrogate the AMR data structure to evaluate one of the variables at a point closest to a specified co-ordinate. I'm new to AMREX and am working with IAMReX. The function needs to be called at every coarse time step but find the variable on the finest grid. I'm having problems with the approach of getting the data using

MultiFab& phi = get_new_data(4);
Real y0 = find_level(0.0, phi);  // scan the x=0 column to find where phi(x)=0

causing a segmentation fault in find_level when it starts to loop over the patches in the MultiFab

find_level(const Real x_ordinate, MultiFab& phi)
{
  // loop over all the patches
#ifdef AMREX_USE_OMP
#pragma omp parallel if (Gpu::notInLaunchRegion())
#endif
    for (MFIter mfi(phi,TilingIfNotGPU()); mfi.isValid(); ++mfi)
    {
     ...

[WeiqunZhang]

https://amrex-codes.github.io/amrex/doxygen/namespaceamrex.html#a77c1d037f86dc7d1fcfe255428503c99

You can try the get_cell_data function.

[WeiqunZhang]

I asked codex to write a function using amrex::get_cell_data and amrex::KeyValue<int,Real>. Here is what it wrote. I think it's correct.

Real get_data_at (RealVect const& x, Vector<MultiFab const*> const& phi,
                  Vector<Geometry> const& geom)
{
    AMREX_ALWAYS_ASSERT(phi.size() == geom.size());

    // KeyValuePair<int,Real> compares on the first field, so store level there.                                                                                                   
    KeyValuePair<int,Real> finest{-1, Real(0.0)};

    for (int lev = 0; lev < static_cast<int>(phi.size()); ++lev) {
        if (phi[lev] == nullptr) {
            continue;
        }

        if (!geom[lev].ProbDomain().contains(x)) {
            continue;
        }

        auto const cell = geom[lev].CellIndex(x.dataPtr());
        if (!geom[lev].Domain().contains(cell)) {
            continue;
        }

        auto const data = amrex::get_cell_data(*phi[lev], cell);
        if (!data.empty()) {
            finest = {lev, data[0]};
        }
    }

    ParallelAllReduce::Max(finest, ParallelContext::CommunicatorSub());

    AMREX_ALWAYS_ASSERT_WITH_MESSAGE(finest.first() >= 0,
                                     "get_data_at: no data found at x");
    return finest.second();
}

[WeiqunZhang]

A more efficient version is

Real get_data_at (RealVect const& x, Vector<MultiFab const*> const& phi,
                  Vector<Geometry> const& geom)
{
    AMREX_ALWAYS_ASSERT(phi.size() == geom.size());

    // KeyValuePair<int,Real> compares on the first field, so store level there. 
    KeyValuePair<int,Real> finest{-1, Real(0.0)};

    for (int lev = static_cast<int>(phi.size()) - 1; lev >= 0; --lev) {
        if (phi[lev] == nullptr) {
            continue;
        }

        if (!geom[lev].ProbDomain().contains(x)) {
            continue;
        }

        auto const cell = geom[lev].CellIndex(x.dataPtr());
        if (!geom[lev].Domain().contains(cell)) {
            continue;
        }

        auto const data = amrex::get_cell_data(*phi[lev], cell);
        if (!data.empty()) {
            finest = {lev, data[0]};
            break;
        }
    }

    ParallelAllReduce::Max(finest, ParallelContext::CommunicatorSub());

    AMREX_ALWAYS_ASSERT_WITH_MESSAGE(finest.first() >= 0,
                                     "get_data_at: no data found at x");
    return finest.second();
}

[dmingram66]

I now need to find the real world (x,y,z) coordinates of a point in a box on the finest level of refinement. Trying to use

geom.CellCenter()

results in a real coordinate which is not inside the computational domain (I assume it's using the wrong value of dx).

Google suggests that

geom[level].CellCenter()

will work but I just get errors saying geom doesn't support [level].

Help would be appreciated.

I have the following code working...

void
NavierStokes::findFreeSurfaceLevel()
{
  // Determine the finest level and the current time.
  const int finest_level = parent->finestLevel();
  const Real time = state[State_Type].curTime();

  // Now we need to get all the boxes which make up the computational
  // domain.  We are looking for the finest box which contains the
  // left hand boundary of the domain.  This box will have a
  // smalEnd(0) ordinate i=0. We start by getting the finest
  // refinement level from the MultiFab.
  NavierStokes& ns_level = getLevel(finest_level);
  MultiFab &mf = ns_level.get_new_data(State_Type);

  // get the box array containing the boxes in the MultiFAB
  const amrex::BoxArray& ba = mf.boxArray();
  const auto prob_lo = geom.ProbLo(1);
  const auto dx = geom.CellSize(finest_level); 

  // Now loop over them, we have to do this since we're not allowed nested MFItter calls
  // which are also used inside the min and max functions.
  for (int i = 0; i < ba.size(); i++) {
    // This is the valid Box of the current FArrayBox. By "valid",
    // we mean the original ungrown Box in BoxArray.
    const Box& box = ba[i];
    
    // is this one of the boxes we are looking for?
    Dim3 lo = lbound(box);
    Dim3 hi = ubound(box);
    if (lo.x == 0) {
	  // get the minimum and maximum values of phi, if the box
	  // contains the free surface then the product of these will
	  // be negative.
    	  Real phi0 = mf.min(box, phicomp);
	  Real phi1 = mf.max(box, phicomp);

	  if (phi0*phi1<0) {
	    // Now we have a box which contains the free surface. We
	    // now have to find the y-ordinate associated with the
	    // location of the free surface at i=0.
	    FArrayBox& fab = mf[i];
	    auto const& a = fab.const_array();

	    int jFace = -1;
	    int i = lo.x; 
	    for (int k = lo.z; k <= hi.z; ++k) {
	      for (int j = lo.y; j <= hi.y; ++j) {
		if ((a(i,j,k,phicomp) * a(i,j+1,k,phicomp)) < 0) {
		  jFace = j;
		}
	      }
	    }

	    // compute the factional position of the zero
	    Real y_zero = (0.0 - a(i,jFace+1, lo.z,phicomp))
	      / ( a(i,jFace, lo.z,phicomp) - a(i,jFace+1, lo.z,phicomp));

	    // Now need to find the ordinate of the bottom of the cell which is located at
            //  (I,jFace,lo.z,phicomp) and dx[finest_level] so that we can compute it.

            // debugging info...
	    Print() << "Box found " <<" [" << box.smallEnd() << ", " << box.bigEnd() << "].\n";
	    Print() << "Phi:  "<< jFace << " " << y_zero << ".\n";
	  }
    }
  }
}

[dmingram66]

It seems that

  const Real dy = parent->Geom(finest_level).CellSize()[1];    
  const Real y0 = geom.ProbLo(1);

allowing

	    // compute the factional position of the zero within the cell
	    Real y_zero = (0.0 - a(i,jFace+1, lo.z,phicomp))
	      / ( a(i,jFace, lo.z,phicomp) - a(i,jFace+1, lo.z,phicomp))
	      * dy;

	    // add the y ordinate of the cell center.  Please note
	    // that if EB are being used this will be in accurate as
	    // we are not calculating the cell Centroid.
	    y_zero =+ y0 + (jFace+0.5)*dy;

works when run on a single MPI process, but the whole routine causes amrex to crash (segfault) when called in parallel.

[WeiqunZhang]

I have trouble understanding what you are trying to do, so I chatted with codex.

Me: do you understand what the user is asking in AMReX Discussion 5508?

Codex: • Yes.

They are trying to write a monitor/free-surface locator in IAMReX/NavierStokes. The practical request is:

They want to find where phi crosses zero along the left boundary, on the finest AMR level, then convert that AMR cell/index location into physical coordinates.

The latest issue is not really get_cell_data anymore. It is two things:

1. They need the finest-level geometry:
   geom in their member function appears to be the current level’s Geometry, not an array. So geom[level] fails. geom.CellSize(finest_level) is also wrong because the argument is direction, not
   AMR level. They found the right idea with parent->Geom(finest_level).CellSize()[1].

2. Their parallel segfault is probably because they loop over mf.boxArray() globally and then do mf[i]. In parallel, not every rank owns every box. That works with one MPI rank, but can crash
   when i refers to a box not local to that process. They should iterate local FABs with MFIter, or check the DistributionMapping owner before touching mf[i].

So the answer should probably explain: use parent->Geom(finest_level) for finest-level physical coordinates, and do not index MultiFab by global BoxArray index unless the FAB is local.

Me: › What does phi crosses zero mean?

Codex: • “phi crosses zero” means the scalar field phi changes sign between two neighboring sample points.

For example, if along a vertical line:

phi(i,j,k) = -0.2
phi(i,j+1,k) = 0.3

then the product is negative:

phi(i,j,k) * phi(i,j+1,k) < 0

so somewhere between those two points, phi = 0.

In their code, that zero location is being treated as the free surface/interface location. They then linearly interpolate between the two neighboring values to estimate the physical y
coordinate where phi equals zero.

Me: › What if there are k's that satisfies the condition?

Codex: • Then the current code is ambiguous and likely wrong.

Right now it does:

int jFace = -1;

for (int k = lo.z; k <= hi.z; ++k) {
for (int j = lo.y; j <= hi.y; ++j) {
if ((a(i,j,k,phicomp) * a(i,j+1,k,phicomp)) < 0) {
jFace = j;
}
}
}

If multiple k values satisfy the condition, jFace just gets overwritten. The final value is the last crossing found in loop order, and the code later uses lo.z instead of the k where the
crossing was found. So in 3D, that is not well-defined.

They need to decide what quantity they actually want:

• One free-surface height for each k: store (jFace, kFace) or compute y_zero[k].
• A single monitor value: define a reduction, e.g. min height, max height, average height, or height at a specific z.
• If physically there should be only one crossing, assert/error when more than one k or more than one j is found.

At minimum, the search should track both indices:

int jFace = -1;
int kFace = -1;

if ((a(i,j,k,phicomp) * a(i,j+1,k,phicomp)) < 0) {
jFace = j;
kFace = k;
}

and then interpolate using kFace, not lo.z.

[WeiqunZhang]

So maybe instead of trying to fix the code (your assumed solution), you can explain what you try to do (the underlying goal)?
https://en.wikipedia.org/wiki/XY_problem

[dmingram66]

I'm impressed that codex does such a good job.

In a heterogeneous multi-fluid simulations it is common to have a marker which can be used to identify the location of the free surface. This could be the volume fraction [in Volume of Fluid methods] or an advected signed distance function [in Level-set methods]. I'm working with IAMReX which implements the level-set method so have an advected scalar phi, which is zero on the free surface. i.e,

$$\frac{\partial \phi}{\partial t}+u\nabla \phi=0,$$ where the set of points on the free surface is, $$\Gamma = \left\{ (x,y) | \phi(x,y)=0 \right\}.$$

The test problem is 2D, though real problems could be 3D. The problem (in 2D) is:

$$\text{given } x=x_s \text{ find } y | \phi(x_s,y)=0.$$

For the test problem $$x_s=0$$, though in a more complex problem we could have a list of locations $$x_s=x_0, x_1, x_2, \ldots $$. Essential I'm trying to create the numerical analog of a fluid level sensor located at $$x=x_s$$. Clearly the 3D problem involves finding y given x and z.

Having found the required value of i,
$$I=floor\left(\frac{x_s-x_0}{\Delta x}\right)$$
we need to find j.

Because phi(i,j,k) contains the centroid values of phi, so we can identify the location at which phi=0, by finding the cells where phi changes sign [i.e. checking for phi(i,j,k)*phi(i,j+1,k)<0 ]. Linear interpolation can then be used to solve for phi=0. This problem is only solved on the finest level box containing the desired coordinate. With the resulting y value and the time being output to a file.

A patch only needs to be searched if it contains the desired coordinate and it contains both fluids, so ubound(box,phi)*lbound(box,phi)<0

[WeiqunZhang]

To fix the segfault, I think you should use MFIter instead of looping of BoxArray.

Real phi0 = mf.min(box, phicomp);

Here the reduce should be on fab not multifab. fab.min(amrex::makeSlab(box,0,lo.x), phicomp)