[FLASH-USERS] Peaks in RHD solution

[Klaus Weide]

On Fri, 28 Apr 2017, ascenzi wrote:

> I am using FLASH4.4 to simulate an expanding ball of matter with a maximum
> speed of 1.e-4c.
> 
> After a preliminary (1D) HD simulation I set up the same problem in RHD
> (taking into account c=1 units), but I obtained a solution with a lot of
> spikes, in spite of the smooth solution of the former case.
> 
> When I implemented the same problem in 2D (cylindrical coordinates) I found
> the same behavior passing from the HD to RHD case (figure in attachment).
> 
> Do you know what could be the origin of these peaks?

Hi Stefano,

First please following comment by a FLASH developer in HydroMain/split/RHD/Config:

# DWL: PLUTO would be  a better choice for doing Relativistic Hydro than FLASH.

You should not expect too much from the relativistic Hydro code in FLASH;
It is not being maintained by anyone AFAIK. (If you have bee able to 
fix something, please share any patches!)


> I noticed that if I rise the maximum speed to 1.e-2 c  I obtain a smooth
> solution. Is it possible that the RHD solver does not work well because when
> v=1.e-4 c I am too far from the relativistic limit?

I don't understand enough about it to say whether that is likely.


> I wanted also to try to solve this problem with the artificial viscosity. Do
> you know if it is implemented also for RHD solver or just for the HD unsplit
> solver ?
> By the way I tried to rise and lower the parameter cvisc in the 1D HD case but
> I have not seen any difference in the solution. How can I be sure that the
> artificial viscosity is active ?

To find out whether a particular runtime parameter is know to a FLASH 
build, look for that runtime parameter by name in
  (a) setup_params   (after running ./setup) or
  (b) default.par    (after running ./setup) or
  (c) the main FLASH log file (after running flash4).

If it is not listed, it cannot possibly be recognized at runtime (and you 
should get a warning about that if it occurs in your .par file.)

In your case, I think you will find neither "cvisc" nor "use_avisc".


HOWEVER, searching for "viscosity" in source/physics/Hydro/HydroMain/split/RHD/*.F90
yields this:

! ---------------------------------------------------
!     Add artificial (Lapidus) Viscosity
! ---------------------------------------------------

  .....

  do i = sweepBegin, sweepEnd
    a = 0.1d0*max(0.e0, Vc(VELN_VAR, i-1) - Vc(VELN_VAR,i))
    flux(DENS_FLUX,i) = flux(DENS_FLUX,i) - a*(Uc(DENS_VAR,i) - Uc(DENS_VAR,i-1))
    .....
  end do


So it appears that artificial viscosity is implemented and is always 
enabled, with a fixed strength corresponding to cvisc = 0.1 .

Klaus