[FLASH-USERS] Computing dt for diffusion

[Aaron Froese]

Hi all,

I need to write a computeDt routine that can accommodate the Braginskii diffusion rates.  I have found the Diffuse_computeDt subroutine to be inadequate, leaving me manually dropping the CFL to maintain stability in the case of strong diffusion constants or large gradients.  The Diffuse_computeDt subroutine appears to determine the allowable time step using only the maximum diffusivity without checking the gradients of the properties being diffused.  Is it not possible, for example, that in the case of strong temperature gradients the entirety of a cell's thermal energy could be transported to neighbouring cells by even a moderate diffusivity?

It seems inefficient to call all the diffusion routines twice: once to determine the time step and a second time to determine the fluxes.  Would it be reasonable to base the allowable time step on the maximum flux inside the entire domain during the previous step?  Below is some code based on the updateConservedVariable subroutine that I think should give a reasonable time step limit for diffusion.  Any suggestions on how to improve it would be most appreciated.

-----------
dt_check = INF
do k = blkLimits(LOW,KAXIS), blkLimits(HIGH,KAXIS)
   do j = blkLimits(LOW,JAXIS), blkLimits(HIGH,JAXIS)
      do i = blkLimits(LOW,IAXIS), blkLimits(HIGH,IAXIS)
         do var = 1, HY_VARINUM-1
            dt_check = min( dt_check, dx*abs(U(HY_DENS+var)) / abs(xFluxR(HY_DENS_FLUX+var) - xFluxL(HY_DENS_FLUX+var)) )
#if NDIM >= 2
            dt_check = min( dt_check, dy*abs(U(HY_DENS+var)) / abs(yFluxR(HY_DENS_FLUX+var) - yFluxL(HY_DENS_FLUX+var)) )
#if NDIM == 3
            dt_check = min( dt_check, dz*abs(U(HY_DENS+var)) / abs(zFluxR(HY_DENS_FLUX+var) - zFluxL(HY_DENS_FLUX+var)) )
#endif
#endif
         enddo
      enddo
   enddo
enddo
dt_check = dt_check*safety_factor
-----------

Thanks,
Aaron