Add debugging section

.. _Chap:Debugging:

Debugging

=========

Debugging is an art.  Everyone has their own favorite method.  Here we

offer a few tips we have found to be useful.

Compiling in debug mode (e.g., :cpp:`make DEBUG=TRUE`) and running with

:cpp:`amrex.fpe_trap_invalid=1` in the inputs file can be helpful.

In debug mode, many compiler debugging flags are turned on and all

:cpp:`MultiFab`s are initialized to signaling NaNs.  The

:cpp:`amrex.fpe_trap_invalid` parameter will result in backtrace files

when a floating point exception occurs.  One can then examine those

files to track down the origin of the issue.

Writing a :cpp:`MultiFab` to disk with

.. highlight:: c++

::

    VisMF::Write(const FabArray<FArrayBox>& mf, const std::string& name);

and examining it with ``Amrvis`` (section :ref:`sec:amrvis` in the AMReX documeintation) 

can be helpful as well.  

You can also use the :cpp:`print_state` routine: 

.. highlight:: c++

::

    void print_state(const MultiFab& mf, const IntVect& cell, const int n=-1);

which outputs the data for a single cell.

Valgrind is another useful debugging tool.  Note that for runs using

more than one MPI process, one can tell valgrind to output to different 

files for different processes.  For example,

.. highlight:: console

::

    mpiexec -n 4 valgrind --leak-check=yes --track-origins=yes --log-file=vallog.%p ./mfix.exe ...

   would be the case for all other geometries). But out of an intersection with

   all planar surfaces. This has the advantage of correctly describing corners.

Fluid Reconstruction

--------------------

The reconstruction algorithm is called whenever a cell in the particle's

neighbor stencil is covered. For no-slip walls, the reconstructed velocity in

that cell is linearly extrapolated from the nearest "valid" fluid cell and 0 at

the wall. This way the fluid velocity is consistent with the no-slip boundary

condition along the normal to the EB. For a planar EB wall, the following would

be enough (the level-set is called :fortran:`phi` here):

.. highlight:: fortran

::

   if ( is_covered_cell(flags(i,j,k))                         .and. &

   &   minval(abs(phi(i:i+1,j:j+1,k:k+1))) <= phi_threshold ) then

       ! Coordinates of cell center

       x_cc = ( real([i,j,k],rt) + half ) * dx

       ! Get phi at cell center

       call amrex_eb_interp_levelset(x_cc, x0, n_refine, phi, phlo, phhi, dx, phi_cc)

       ! Get normal at cell center

       call amrex_eb_normal_levelset(x_cc, x0, n_refine, phi, phlo, phhi, dx, norm_cc)

       ! Initial guess of interpolation point:

       x_i  = x_cc + two * abs(phi_cc) * norm_cc

       ! Get phi at interpolation point

       call amrex_eb_interp_levelset(x_i, x0, n_refine, phi, phlo, phhi, dx, phi_i)

       ! Compute interpolated velocity at x_i

       vel_i = trilinear_interp(vel_in, vilo, vihi, 3, x_i, x0, dx)

       ! Since interpolation point is only slightly shifted with respect to

       ! the mirror point, we approximate vel at mirror point with vel_i and

       ! then use linear interpolation between x_m and x_c

       vel_out(i,j,k,:) = vel_i * phi_cc / phi_i

If the EB represents a curved wall, the initial normal is not a good estimate

for the normal at the closest point on the wall. Therefore we start a the

covered cell, and "walking" along the EB normal on cell at a time (computed from

the level-set function), until the neighbor stencil does not include covered

cells:

.. highlight:: fortran

::

   ! Find location of interpolation point by iteration if necessary

   iter = 0

   find_xi: do

       if ( interp_stencil_is_valid(x_i, x0, dx, flags, flo, fhi) ) exit find_xi

       ! Get normal at interpolation point

       call amrex_eb_normal_levelset(x_i, x0, n_refine, phi, phlo, phhi, dx, norm_i)

       x_i = x_i + maxval(dx) * norm_i

       iter = iter + 1

       if ( iter > max_iter ) &

           call amrex_abort("reconstruct_velocity(): cannot find interpolation point")

   end do find_xi

Note that the level-set here needs to be at the same resolution as the fluid.

This is the reason why we need to keep the coarse level :cpp:`level_sets[0]`

even when running in single-level mode.

.. _AMReX EB documentation: https://amrex-codes.github.io/amrex/docs_html/EB_Chapter.html

.. _AMReX Level-Set documentation: https://amrex-codes.github.io/amrex/docs_html/EB.html#level-sets

.. _AMReX geometry documentation: https://amrex-codes.github.io/amrex/docs_html/EB.html#initializing-the-geometric-database 

   Fluids_Chapter

   Particles_Chapter

   EB

   Debugging

Notice

------