[FLASH-USERS] Some problemes with FLASH3.1.1 Jeans (was Re: Negative/zero densities ...)

[Klaus Weide]

Flash Users,

This message is to explain, and help circumvent, some of the problems 
that came up earlier this month with the Jeans problem ins FLASH 3.1.1.
Some of this will probably apply to other problems that were brought up in 
the same context. In particular, private copies of Grid_markRefineDerefine 
should be checked for errors similar to those described below (near the 
end).

Seyit,

After applying all three runtime parameter changes, I was indeed also able 
to reproduce the warning messages you had reported getting from Jeans.

As you already seem to have found out yourself, the 
Grid_markRefineDerefine.F90 provided with the Jeans problem is buggy.  A 
patch is attached.  The original version resulted (with your changed 
runtime parameters) in a pattern where a region of the domain was left
at a coarser refinement than the remainder, without a physical reason.
This showed up some further problems that become obscured after 
Grid_markRefineDerefine is fixed. Let me therefore first discuss
those problems.

On Thu, 9 Apr 2009, Seyit Hocuk wrote:
> I also set convertToConsvdForMeshCalls to .True. in runtime parameters 
> instead of using convertToConsvdInMeshInterp, which should have worked better 
> because it is supposed to avoid spurius things in paramesh3+. However 
> convertToConsvdInMeshInterp calls at some point gr_sanitizeDataAfterInterp, 
> which creates all the weirdness. It is my humble opinion that something goes 
> wrong with conservation there.

It is not (runtime parameter) convertToConsvdInMeshInterp or (subroutine) 
gr_sanitizeDataAfterInterp that created any problems. 
gr_sanitizeDataAfterInterp (despite its name) only checks and reports some 
problems with data, it does not modify any data; and it happens to be 
called only when convertToConsvdInMeshInterp is TRUE. Now 
gr_sanitizeDataAfterInterp was checking maybe a bit too aggressively. That 
is, as Anshu has already written, it was checking guard cells in PARENT 
blocks for plausible values even though the values in those cells cannot 
impact the propagation of solutions *on LEAF blocks* at all.  And, as has 
just been also pointed out on FLASH-USERS, FLASH (as provided) only cares 
about evolving solutions on leaf blocks.

At the point where the checking occurs, these guard cells have indeed
just been updated by PARAMESH, by copying values from a neighboring
block; so they should have valid data IF the neighbor did.  But if
this neighbor of the PARENT (nodetype==2) block is an ANCESTOR
(nodetype==3) block, its data may not be valid.  (It can be shown that
in THIS situation at a refinement boundary, the resulting PARENT cell
values do not impact solutions on the child (nodetype==1 aka LEAF)
blocks, even though in other cases PARENT cell values can and do
impact child solutions via interpolation/prolongation.)

Normally during most simulations, ANCESTOR blocks always contain 
reasonably-looking values for all variables in their interior cells, even 
though FLASH doesn't explicitly evolve ANCESTOR (or PARENT) blocks in its 
solvers.  These values either come from the initial values established 
when Simulation_initBlock was called during initialization, or from 
higher-resolution children via LEAF->PARENT->ANCESTOR restriction done 
before checkpoints or plot files are written. They may therefore be out of 
date with respect to the current simulation time, but will not trigger 
gr_sanitizeDataAfterInterp warnings.

However, this does not apply to the Jeans simulation in question, as
far as density DENS_VAR is concerned, because of a side effect of the
Multigrid solver.  For certain kinds of boundary conditions, this
solver temporarily modifies the density values in blocks by
subtracting an offset, and later undoes this change by adding the same
value.  There is, however, an inconsistency in the sets of blocks to
which these changes are applied: the subtraction is applied to all
blocks, the addition only to leaf blocks.  I would like to emphasize
that this inconsistency causes no differences in results as long as
one only expects meaningful solution data in leaf blocks; and the only 
negative effects are unnecessary WARNINGs from gr_sanitizeDataAfterInterp
as described above.

The following one-line change will make the addition apply to the same
set of blocks as the subtraction, and will thus avoid those warnings:

--- source/Grid/GridSolvers/Multigrid/gr_hgSolve.F90    (revision 10484)
+++ source/Grid/GridSolvers/Multigrid/gr_hgSolve.F90    (working copy)
@@ -185,7 +185,7 @@
    ! when using periodic/Neumann boundary conditions).

    do m = 1, gr_hgMeshRefineMax
-     call gr_hgLevelAddScalar(m, gr_iSource, gr_hgAvgSource, MG_NODES_LEAF_ONLY) !oK
+     call gr_hgLevelAddScalar(m, gr_iSource, gr_hgAvgSource, MG_NODES_ALL_NODES) !oK
    enddo

    ! Leave boundary zones properly updated.


(Ultimately gr_sanitizeDataAfterInterp checking should be changed so that it
only checks blocks (and cells) where valid data can always be expected.)

  - * - * -

Finally, attached is a patch that fixes a few things specific to the Jeans
simulation:

-  Initialize internal energy properly in Simulation_initBlock.
-  Fix inconsistent use of lb vs gr_blkList(lb), a common mistake.
-  Changed LEAF -> ACTIVE_BLKS, the logic involving delta_max_par requires it.
-  Make sure delta_max_par, delta_max array elements are initialized.

Apply in the top directory (one level above source) with

   patch -p0 < forSeyit.diff



Klaus
-------------- next part --------------
Index: source/Simulation/SimulationMain/Jeans/Simulation_initBlock.F90
===================================================================
--- source/Simulation/SimulationMain/Jeans/Simulation_initBlock.F90	(revision 10484)
+++ source/Simulation/SimulationMain/Jeans/Simulation_initBlock.F90	(working copy)
@@ -110,8 +110,8 @@
 
         gam = sim_gamma
         ek = 0.5 * (vx*vx+vy*vy+vz*vz)
-        e  = p/(rho*(sim_gamma-1.))
-        e  = e + ek
+        ei  = p/(rho*(sim_gamma-1.))
+        e  = ei + ek
         e  = max(e, sim_smallE)
 
         startingPos(1) = 1
Index: source/Simulation/SimulationMain/Jeans/Grid_markRefineDerefine.F90
===================================================================
--- source/Simulation/SimulationMain/Jeans/Grid_markRefineDerefine.F90	(revision 10484)
+++ source/Simulation/SimulationMain/Jeans/Grid_markRefineDerefine.F90	(working copy)
@@ -41,7 +41,7 @@
   use Driver_interface, ONLY : Driver_getSimTime
   use Grid_interface, ONLY : Grid_fillGuardCells, &
     Grid_getListOfBlocks, Grid_getBlkIndexLimits, Grid_getBlkPtr, &
-    Grid_releaseBlkPtr, Grid_getLocalNumBlks
+    Grid_releaseBlkPtr
   use Grid_data, ONLY : gr_refine_cutoff,gr_derefine_cutoff, &
        gr_refine_filter, gr_numRefineVars, gr_refine_var
 
@@ -87,7 +87,7 @@
   real, save    :: nchild = 2.**NDIM
 
   call Grid_fillGuardCells(MyPE, CENTER, idir)
-  call Grid_getListOfBlocks(LEAF,gr_blkList,blockCount)
+  call Grid_getListOfBlocks(ACTIVE_BLKS,gr_blkList,blockCount)
 
   do i = 1, blockCount
      call Grid_getBlkIndexLimits(gr_blkList(i),blkLimits,blkLimitsGC)
@@ -107,13 +107,14 @@
      call gr_markRefineDerefine(MyPE,iref,ref_cut,deref_cut,ref_filter)
   end do
 
-  do lb = 1, blockCount
+  delta_max(1:lnblocks) = -huge(1.0)
 
-     call Grid_getBlkPtr(gr_blkList(lb),solnData)
+  do l = 1, blockCount
 
-     call Grid_getBlkIndexLimits(gr_blkList(lb),blkLimits,blkLimitsGC)
+     lb = gr_blkList(l)
 
-     delta_max(lb) = -huge(1.0)
+     call Grid_getBlkPtr(lb,solnData)
+     call Grid_getBlkIndexLimits(lb,blkLimits,blkLimitsGC)
 
      do k = blkLimitsGC(LOW, KAXIS), blkLimitsGC(HIGH, KAXIS)
         do j = blkLimitsGC(LOW, JAXIS), blkLimitsGC(HIGH, JAXIS)
@@ -124,7 +125,7 @@
         enddo
      enddo
 
-     call Grid_releaseBlkPtr(gr_blkList(lb),solnData)
+     call Grid_releaseBlkPtr(lb,solnData)
 
   end do
 
@@ -210,7 +211,6 @@
   ! Make sure that the refinement level of the blocks is between lrefine_min
   ! and lrefine_max
 
-  call Grid_getLocalNumBlks(lnblocks)
   do i = 1, lnblocks
 
      if ((lrefine(i) == lrefine_min) .and. (nodetype(i) == 1)) &