[FLASH-USERS] Missing term in analytic Sedov (+ Nupdate?)

[Massimo Gaspari]

Hi guys,

[pre-scriptum: have you tried to compute Nupdate in just one of your simulation (as pointed out in my previous post)? I am really interested to know this number that you can compute in a couple of seconds (please write the used modules, Hydro, Gravity, etc.). Nupdate is simply the number of cells updated per second and per processor. Could you be so kind to provide this number?]

I was doing a benchmark to test FLASH (+Hydro) with the Sedov problem and I found that the computed analytic solution was incorrect.There are two corrections to do:

1) Major. In sedov_vfunc of Simulation_initBlock.F90 it is missing Vpfe in the last line (see Landau & Lifshitz, Fluid Mechanics, p.405)):

  xi = nu1*alog(5.-(3.*gamma-1.)*Vfpe) + &        nu2*alog(gamma*Vfpe-1.) - &        nu1*alog(gam7/gamp1) - nu2*alog(gamm1/gamp1) - &        2.*alog(0.5*gamp1 *Vfpe) !<--- Vfpe WAS MISSING (see Landau p. 405)

2) Minor. The solution for xi --> 0. gives estimate for P about two orders of magnitude lower (in set_analytic_sedov). Therefore, I have taken the solution given in the book of Sedov (1959), which provides the correct asymptotic value for pressure (constant).This should be used at least when:
        if (xi .le. 0.075) then     ! Use asymptotic r->0 solution. 
Note that ".le. 1.e-6" is too low.

OLD (what is the reference?):
           kappa = ( (0.5*gamp1/gamma)**2. * &                 (gamp1/gam7* &                 (5.-(3.*gamma-1.)/gamma))**nu1 )**(1./nu2) * &                 gamm1/gamp1 / gamma           
           epsilon = k**(nu5+1.) * (k*gamma*kappa)**nu3 * &                 (gamp1/gam7*(3.*gamma-1.))**nu4 * &                 ((2.*gamma+1)/gamma/(3.*gamma-1.)) * &                 (gamm1/gamma)           zeta = gamm1*gamm1/(2.*gamma*gamma*kappa)           p(i) = rho_ambient/gamma * 0.16*(R0/t)**2 * epsilon*zeta

NEW (ok):   
           kappa = (0.32/2.**(6./5.)) *&                   ((gamma+1.)**((gamma+1.)/(3.*gamma-1.))/gamma**(4.*gamma/(3.*gamma-1.))) *&                   ((2.*gamma+1.)/(7.-gamma))**((13.*gamma**2-7.*gamma+12.)/&                    (5.*(2.-gamma)*(3.*gamma-1.))) *&                   beta**2 !--> Sedov assumed beta =1, see Eq. 11.11           p(i)  = kappa*rho_ambient*(E/rho_ambient)**(2./5.) *t**(-6./5.)




Now the analytic and simulated profiles are consistent (within errors).The plots in the user's guide seem somehow correct (e.g. Fig. 25.9); did you use other routines?

I would also advise to use higher n_iter_max (>1000) and to add only V = 1./gamma + 1.e-18, instead ofV = 1./gamma + 1.e-8 (in compute_sedov_v).


I hope this can help if you are testing Sedov problem.(Let me know if I missed something)

Cheers,


   Max
 		 	   		  
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://flash.rochester.edu/pipermail/flash-users/attachments/20120615/2ec241ad/attachment.htm>