[FLASH-BUGS] 3 flash2.1 bugs

[Erik-Jan Rijkhorst]

Dear Flash developers,

I found 3 bugs in Flash2.1 which I've listed below.

Some info on the system I'm using:
(http://www.sara.nl/userinfo/teras/description/index.html)
Machine:          SGI Origin 3800
Compiler:         MIPSpro Fortran 90 compiler
Operating system: IRIX64 p1 6.5 04101930 IP35

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

Hard coded value in source/mesh/amr/paramesh2.0/ref_marking.F90:

...
! compute the error
                 num = 0.
                 denom = 0.

                 do kk = 1, ndim2
                    num = num + delu2(kk)**2
                    denom = denom + (delu3(kk) + &
                         &  (epsil*delu4(kk)+1.e-20))**2
                 end do
...

When using very low densities (of the order of 1e-26) the code won't
refine. I solved this by putting my own version of ref_marking.F90 in the
setups/problem directory and replaced the value of 1.e-20 by 1.e-200. It
would be nice to change this to a parameter in a future Flash release?

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

Multiplication by cf_area is done twice at step 7b in
source/hydro/explicit/delta_form/kurganov/hydro_3d.F90
(you'll only notice this when running 3D simulations...)

...
!!!!!!!!! Step 7b: Apply geometry factors.
!!!!!!!!!          Update deltaU.
        do j=1,ny
        do i=1,nx
          cf_area = (            ( xr(i)-xl(i) )                )     &
                  * (            ( yr(j)-yl(j) )                )
          flux_block(:,i,j,:) = flux_block(:,i,j,:)* cf_area
          do k=1,nz
            c_volume_i = 1.0/( cf_area*( zr(k) - zl(k) ) )
            flux_block(:,i,j,k) = flux_block(:,i,j,k)* cf_area
            du_local(:,i,j,k) =  du_local(:,i,j,k  ) -                &
                c_volume_i*(   flux_block(:,i,j,k  )                  &
                             - flux_block(:,i,j,k-1) )
          enddo
        enddo
        enddo
...

I think this should be

...
!!!!!!!!! Step 7b: Apply geometry factors.
!!!!!!!!!          Update deltaU.
        do j=1,ny
        do i=1,nx
          cf_area = (            ( xr(i)-xl(i) )                )     &
                  * (            ( yr(j)-yl(j) )                )
          flux_block(:,i,j,:) = flux_block(:,i,j,:)* cf_area
          do k=1,nz
            c_volume_i = 1.0/( cf_area*( zr(k) - zl(k) ) )
            du_local(:,i,j,k) =  du_local(:,i,j,k  ) -                &
                c_volume_i*(   flux_block(:,i,j,k  )                  &
                             - flux_block(:,i,j,k-1) )
          enddo
        enddo
        enddo
...

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

Wrong index in source/formulation/state_form/du_update_zface.F90:
I used this formulation to implement an approximate Rieman solver that we
have here at Leiden.
I copied source/hydro/explicit/delta_form/kurganov to my own
source/hydro/explicit/delta_form/roe directory (see previous bug).
Note: you won't notice this bug until you start running 3D sims...

...
  do 200 j=1,ny
!   Transpose ys into eos1d format.
    do n=1,ns
      ys_line(:,n) = ys(n,:,j)
    enddo
    call eos1d(2, 1, nx,        &
                   rho(:,k),    &
                     t(:,k),    &
                     p(:,k),    &
                  eint(:,k),    &
               gamma_c(:,k),    &
             ys_line(:,:  ),    &
               nx,ns          )
200 continue
...

I changed this to:

...
  do 200 j=1,ny
!   Transpose ys into eos1d format.
    do n=1,ns
      ys_line(:,n) = ys(n,:,j)
    enddo
    call eos1d(2, 1, nx,        &
                   rho(:,j),    &
                     t(:,j),    &
                     p(:,j),    &
                  eint(:,j),    &
               gamma_c(:,j),    &
             ys_line(:,:  ),    &
               nx,ns          )
200 continue
...

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

cheers,

Erik-Jan Rijkhorst