Improvements to TOV evolutions, pt 2

Follow up to !107 (merged).

Atm it seems that we can stably, and also fairly good, evolve a stable TOV in Cowling approximation with resolution that is comparable (in one dimension) with that of Deppe's paper, that is, n = 5, k = 16-24. What is interesting is that it seems like we don't need an AV for the stable star.

The key change seems to have been the 'freezing' of the atmosphere, combined with the velocity limiter from the Dumbser paper. Freezing atmosphere means that we restrict the points at which we are allowed to impose the atmosphere within ERK substeps. Which points are allowed to evolve during substeps is determined from the initial data of each step. The idea is to enforce a finite propagation speed for the atmosphere treatment.

Whether that works also in dynamical scenarios needs to be checked. With regards to the finite propagation speed. I think one can cook up a nice argument in relation to a CFL like condition to show that freezing the atmosphere makes sense (=^= light cones).

src/GRHD/callbacks.jl

@@ -18,8 +18,34 @@ end
 #######################################################################
 
 
+function update_atm_domain_of_dependence!(mesh::Mesh1d)
+  @unpack c2p_freeze_atm, c2p_reset_atm = get_static_variables(mesh)
+  freeze, reset = c2p_freeze_atm, c2p_reset_atm
+  for i in 2:length(reset)-1
+    rl, rc, rr = reset[i-1], reset[i], reset[i+1]
+    if rl == rc == rr
+      # atmosphere can't spread within next timestep
+      freeze[i-1] = freeze[i] = freeze[i+1] = 1
+      continue
+    else
+      freeze[i-1] = freeze[i] = freeze[i+1] = 0
+    end
+  end
+end
+
+
 function callback_equation(state_u, state_t, env, P, isperiodic, eq::Equation)
   wrap_dynamic_variables!(env.cache, state_u)
+  if P.prms.c2p_freeze_atm_reset
+    if formulation(P) === :spherical1d
+      broadcast_volume!(cons2prim, eq, env.mesh)
+    elseif formulation(P) === :rescaled_spherical1d
+      broadcast_volume!(cons2prim_rescaled_spherical1d, eq, env.mesh)
+    else
+      TODO()
+    end
+    update_atm_domain_of_dependence!(env.mesh)
+  end
   callback_equation(state_t, isperiodic, P.equation, env.mesh)
 end
 function callback_equation(state_t, isperiodic, eq::Equation, mesh::Mesh1d)
@@ -101,14 +127,20 @@ function callback_hrsc(state_u, state_t, env, P, isperiodic, hrsc::HRSC.Abstract
   end
   # display(extrema(cellmax_v))
 
+  @unpack abslog_D, abslog_τ = get_static_variables(cache)
+  @. abslog_D = log10.(abs.(D))
+  @. abslog_τ = log10.(abs.(τ))
+
   fill!(tmp_mu, 0.0)
   HRSC.compute_viscosity!(
       mu, tmp_mu,
+      (abslog_D, abslog_τ), cellmax_v,
+      # (log10.(abs.(D)),log10.(abs.(τ))), cellmax_v,
       # (log10.(abs.(D)),Sr,log10.(abs.(τ))), cellmax_v,
       # (log10.(abs.(D)),), cellmax_v,
       # (log10.(abs.(D)),log10.(abs.(Sr)),log10.(abs.(τ))), cellmax_v,
       # (D,Sr,τ), cellmax_v,
-      (D,τ), cellmax_v,
+      # (D,τ), cellmax_v,
       # (vr,), cellmax_v,
       hrsc)
   # ρmin = equation.atm_density * equation.atm_factor
@@ -132,20 +164,20 @@ function callback_hrsc(state_u, state_t, env, P, isperiodic, hrsc::HRSC.Abstract
   # end
   mu[1] = mu[2] = mu[end] = 0 # ignore boundaries
 
-  prms = parameters(:HRSC)
-  prms["method"] = "filter"
-  prms["filter_method"] = "bernstein"
-  dg1d.check_parameters(:HRSC, prms)
-  bernstein_hrsc = HRSC.make_HRSC(mesh, prms)
-  flag = mu .> 0.0
-  @unpack D, Sr, τ = get_dynamic_variables(mesh)
-  logD = log.(abs.(D))
-  logτ = log.(abs.(τ))
-  HRSC.reconstruct!(logD,  flag, bernstein_hrsc, isperiodic=false)
-  HRSC.reconstruct!(Sr, flag, bernstein_hrsc, isperiodic=false)
-  HRSC.reconstruct!(logτ,  flag, bernstein_hrsc, isperiodic=false)
-  D .= exp.(logD)
-  τ .= exp.(logτ)
+  # prms = parameters(:HRSC)
+  # prms["method"] = "filter"
+  # prms["filter_method"] = "bernstein"
+  # dg1d.check_parameters(:HRSC, prms)
+  # bernstein_hrsc = HRSC.make_HRSC(mesh, prms)
+  # flag = mu .> 0.0
+  # @unpack D, Sr, τ = get_dynamic_variables(mesh)
+  # logD = log.(abs.(D))
+  # logτ = log.(abs.(τ))
+  # HRSC.reconstruct!(logD,  flag, bernstein_hrsc, isperiodic=false)
+  # HRSC.reconstruct!(Sr, flag, bernstein_hrsc, isperiodic=false)
+  # HRSC.reconstruct!(logτ,  flag, bernstein_hrsc, isperiodic=false)
+  # D .= exp.(logD)
+  # τ .= exp.(logτ)
   # broadcast_volume!(cons2prim, P.equation, mesh)
   # broadcast_volume!(maxspeed, P.equation, mesh)
   # broadcast_volume!(flux_source_spherical1d, P.equation, mesh)
@@ -170,9 +202,18 @@ function callback_hrsc(state_u, state_t, env, P, isperiodic, hrsc::HRSC.Smoothed
   @unpack cache, mesh = env
   @unpack smoothed_mu = get_static_variables(cache)
   hrsc.smoother(smoothed_mu, mu, mesh, false)
-  max_D = maximum(D)
-  min_D = minimum(D)
-  smoothed_mu .*= (D .- min_D) / max_D
+  # min_D, max_D = extrema(D)
+  # max_D = maximum(D)
+  # min_D = minimum(D)
+  # smoothed_mu .*= (D .- min_D) / max_D
+  # @. smoothed_mu *= D * min_D / max_D
+  K, = env.mesh.tree.dims
+  mat_D = dg1d.vreshape(D, layout(env.mesh))
+  mat_smoothed_mu = dg1d.vreshape(smoothed_mu, layout(env.mesh))
+  for k in 1:K
+    min_D, max_D = extrema(@views mat_D[:,k])
+    @. @views mat_smoothed_mu[:,k] *= mat_D[:,k] * min_D / max_D
+  end
 end