[FLASH-USERS] cell density goes to zero when refining several times

[Andrea Gatto]

Dear Christoph, Klaus, et al,

I try to be more specific.

I am using a (split) Bouchut 5 waves solver, since I'd like to include
magnetic fields in future. I'm using MODE_DENS_EI.

I have a simulation of gas spanning over a wide ranges of density and
temperature, with a box size of the order of 10^20-10^21 cm.
The gas is shaped by stellar feedback.
This setup causes the medium to be extremely turbulent and intermittent.

Specifically, I'm interested in capturing the correct physics processes at
small scales. However, since the gas is extremely turbulent, increasing
lrefine_max doesn't help me since the 2nd derivative-based refinement
criterion causes the whole (or majority of the) domain to be refined to
the this maximum level, which makes the simulation too slow and impossible
to run.

Moreover, I'm not really interested in resolving shocks/discontinuity et al.
Even if I have a big dense cloud with uniform thermodynamical properties,
I'd like to have this refined to high resolution.

For this reason, I've modified the Grid unit allowing for different
refinement criteria, with each one of them having its own maximum
refinement level.
These criteria are: refinement on variables threshold and jeans lenght
refinement (based on the default routine present in the sink particle
unit).

What I do in practice is to rearrange the refine, derefine, and stay flags
of the blocks in the following way in Grid_markRefineDerefine.

1) I call gr_markRefineDerefine and I let the code mark the blocks
following the default 2nd derivative-based criterion. I rearrange the
flags so that this method refines only up to a level lrefine_2ndD. Each
block above this refinement level is marked automatically for
derefinement.

2) I call the new criteria, e.g. gr_markVarThreshold, which reorganize the
refine, derefine and stay flags according to what I want to a higher
maximum refinement level.
For instance, I could say I want a block refined if the density of one of
its cell is above 1.e-22 gr/cm^3, and I want it to derefine if the
densities are below 1.e-24 gr/cm^-3.
This routines are able to refine up to a higher maximum refinement level
with respect to lrefine_2ndD.

This method is extremely powerful, as I'm able to reach high resolutions
exactly where I want without having too many blocks.
For instance, in my tests I see a number of blocks of 40 000 - 60 000
against the 2 - 3 millions I receive if I allow the 2nd derivative-based
criterion to refine to the highest refinement level.

Now, I know this is a bad idea, but it's the only way to reach high
resolutions where I want.
In principle I could play a bit with the 2nd derivative-based criterion
parameters, but still this would not allow me to reach high resolutions in
the gas parcels I'm interested in.

I made a number of tests with this new method, from a single supernova
explosion in a uniform medium to ram-pressure stripping of a dense blob
moving at high velocities, and the blocks are refined and derefined
correctly following my criteria without any problem.

The problem arises when I have a turbulent medium. The method refines and
derefines correctly, but at a certain point the solver crashes, as I
described in the previous email.
I've found that a single (and the neighbor cells) has a density which goes
to zero, and internal and kinetic energies that explode.
Every time I've checked I've found that this cell is a guard cell between
two different blocks at different high refinement levels.

This problem arises regardless the method I choose (jeans, threshold or a
combination of them), with or without sink particles, with or without
gravity, using a periodic box or an elongated box along z with periodic
boundary conditions in x and y and diode in z.
As long as I have multiple maximum refinement levels in a turbulent
environment, the code crashes. I also tried different combinations of
eintSwitch and density/energy floors, but the problem still remains.

Given the "universality" of this problem, I suppose this has to do with
the amr routines which handle the density/energy/velocity filling from
parent to child blocks and vice versa.

I'll try to use the routine Klaus has suggested to have a deeper look on
the problem, however I'd really appreciate any help in this regard.

Sorry for the lenghty email, but I think now you understand what I'm
trying to do.

Thank you.

Best,

Andrea


> On Wed, 28 May 2014, Andrea Gatto wrote:
>> Dear all,
>> I'm running a simulation with FLASH 4.0.1 of a turbulent interstellar
medium.
>> For different reasons, some regions of my domain are refined several
times
>> during the evolution (nrefs=2), so that every two timesteps blocks in
the
>> same region are refined, derefined, refined again etc.
>> After few hundreds timesteps, however, the code crashes due to a steep
increase in internal and kinetic energy.
>> Every time I checked, I found that a single/few cells belonging to one
of
>> this highly time-varying blocks show a density that rapidly goes to
zero
>> and a subsequent explosion in temperature and velocity.
>> I tried several solutions, such as the inclusion of guard cells filling
routines in different parts of the code, but the problem remains. This
doesn't happen when the blocks are not refined so often.
>> Do you have any idea on how to solve this problem?
> Andrea,
> First of all, I think a situation like this should be avoided if at all
possible. I.e., don't refine-then-derefine the same region so many times.
> Second, the behavior will likely depend on the Hydro version you are
using
> and its runtime parameters, like "order", "RiemannSolver", etc. if you
are
> using unsplit Hydro. Others will be better able to give guidance here.
Third, the following MAY be helpful for circumventing or analyzing the
problematic behavior:
> In FLASH4.2, we have added two runtime parameter to make it easier to
deal with situations where some invalid values start appearing in
density
> (also, total and internal energy). Note that the mode that "fixes" data
by flooring, gr_sanitizeDataMode=3, has not been tested much in
practice.
> Excerpt from source/Grid/GridMain/paramesh/paramesh4/Config :
> ========================================================================
D gr_sanitizeDataMode What to do when gr_sanitizeDataAfterInterp is called
> D & to check for acceptable values in the dens, ener, and eint
> cell-centered variables
> D & after a Grid operation may have resulted in grid interpolation. D &
0: Do nothing.
> D & 1: Check (if variable is not masked out) and report (see
> sanitizeVerbosity).
> D & 2: Check (ignoring variable mask) and report (see
sanitizeVerbosity).
> D & 3: Check (if variable is not masked out) and fix (apply floor
value).
> D & 4: Check (if variable is not masked out) and abort if cell is found
below floor value.
> PARAMETER gr_sanitizeDataMode   INTEGER 1       [0,1,2,3,4]
> D gr_sanitizeVerbosity How to write information about unacceptable D &
values in the dens, ener, and eint cell-centered variables if
gr_sanitizeDataAfterInterp
> D & finds value that are below the acceptable floor.
> D & This reporting is in addition to other actions selected with
gr_sanitizeDataMode=3 or 4.
> D & 0: Be quiet.
> D & 1: Only write a log file message per block if unacceptable value
found
> on MASTER_PE.
> D & 4: As 1, and each proc writes a line to standard output for each
block
> with bad values.
> D & 5: As 4, and each proc writes lines showing the values in all cells
of
> the block (in 1D/2D)
> D & 5: or a 2D slice (in 3D).
> PARAMETER gr_sanitizeVerbosity  INTEGER 5       [0,1,4,5]
> ========================================================================
HTH,
> Klaus