diff --git a/src/GRHD/callbacks.jl b/src/GRHD/callbacks.jl
index 13d3051153af3517c4076adc7cafef9257c7015a..840fe3550a4f6b867e14b7a42bc67fbb341b8851 100644
--- a/src/GRHD/callbacks.jl
+++ b/src/GRHD/callbacks.jl
@@ -52,8 +52,8 @@ end
 
 compute_atmosphere_mask(env, P, mesh) = nothing
 function compute_atmosphere_mask(env, P, mesh::Mesh1d{SpectralElement})
-  enforce        = P.prmsdb["GRHD"]["c2p_enforce_atm"]
-  enforce_causal = P.prmsdb["GRHD"]["c2p_enforce_causal_atm"]
+  enforce        = P.prms.c2p_enforce_atm
+  enforce_causal = P.prms.c2p_enforce_causal_atm
   (!enforce && !enforce_causal) && return
   if formulation(P) === :valencia1d
     broadcast_volume!(cons2prim_valencia1d, P.equation, env.mesh)
diff --git a/src/GRHD/rhs.jl b/src/GRHD/rhs.jl
index 71827f7c05c7497d8d3540a24f52ecdac754bf36..dd7ed4804fab56e678f0c073baff0dcd4c3db264 100644
--- a/src/GRHD/rhs.jl
+++ b/src/GRHD/rhs.jl
@@ -93,9 +93,9 @@ end
 
 
 function rhs!(mesh::Mesh1d{FVElement}, P::Project)
-  if P.prmsdb["GRHD"]["hrsc"] === "fv_central"
+  if P.prms.hrsc === :fv_central
     rhs_fv_central!(mesh, P)
-  elseif P.prmsdb["GRHD"]["hrsc"] === "slope_limiter"
+  elseif P.prms.hrsc === :slope_limiter
     rhs_fv_slope_limiter!(mesh, P)
   else
     TODO(P.prms.hrsc)
@@ -117,7 +117,7 @@ function rhs_fv_central!(mesh, P::Project{:spherical1d})
           bdry_D, bdry_Sr, bdry_Ï„ = get_bdry_variables(cache)
   @unpack dt                      = get_global_variables(cache)
 
-  if P.prmsdb["GRHD"]["c2p_enforce_causal_atm"]
+  if P.prms.c2p_enforce_causal_atm
     broadcast_volume!(cons2prim_spherical1d_freeze_flags, equation, mesh)
   else
     broadcast_volume!(cons2prim_spherical1d, equation, mesh)
@@ -130,7 +130,7 @@ function rhs_fv_central!(mesh, P::Project{:spherical1d})
   fv_update_step!(rhs_Sr, Sr, flr_Sr, src_Sr, bdry_Sr, nflr_Sr, dt[1], mesh)
   fv_update_step!(rhs_Ï„,  Ï„,  flr_Ï„,  src_Ï„,  bdry_Ï„,  nflr_Ï„,  dt[1], mesh)
 
-  if !P.prmsdb["GRHD"]["atm_evolve"]
+  if !P.prms.atm_evolve
     broadcast_volume!(determine_atmosphere, P.equation, mesh)
     broadcast_volume!(stop_atmosphere_evolution_spherical1d, P.equation, mesh)
   end
@@ -154,7 +154,7 @@ function rhs_fv_central!(mesh, P::Project{:rescaled_spherical1d})
 
   broadcast_volume!(unscale_conservatives, equation, mesh)
   impose_symmetry!(P, mesh)
-  if P.prmsdb["GRHD"]["c2p_enforce_causal_atm"]
+  if P.prms.c2p_enforce_causal_atm
     broadcast_volume!(cons2prim_rescaled_spherical1d_freeze_flags, equation, mesh)
   else
     broadcast_volume!(cons2prim_rescaled_spherical1d, equation, mesh)
@@ -168,7 +168,7 @@ function rhs_fv_central!(mesh, P::Project{:rescaled_spherical1d})
   fv_update_step!(rhs_rSr, rSr, flr_rSr, src_rSr, bdry_Sr, nflr_rSr, dt, mesh)
   fv_update_step!(rhs_rτ,  rτ,  flr_rτ,  src_rτ,  bdry_τ,  nflr_rτ,  dt, mesh)
 
-  if !P.prmsdb["GRHD"]["atm_evolve"]
+  if !P.prms.atm_evolve
     broadcast_volume!(determine_atmosphere, P.equation, mesh)
     broadcast_volume!(stop_atmosphere_evolution_rescaled_spherical1d, P.equation, mesh)
   end
@@ -191,7 +191,7 @@ function rhs_fv_central!(mesh, P::Project{:doublecartoon})
           bdry_rD, bdry_rSx, bdry_rτ = get_bdry_variables(cache)
 
   broadcast_volume!(unscale_conservatives_doublecartoon, equation, mesh)
-  if P.prmsdb["GRHD"]["c2p_enforce_causal_atm"]
+  if P.prms.c2p_enforce_causal_atm
     broadcast_volume!(cons2prim_doublecartoon_freeze_flags, equation, mesh)
   else
     broadcast_volume!(cons2prim_doublecartoon, equation, mesh)
@@ -202,17 +202,17 @@ function rhs_fv_central!(mesh, P::Project{:doublecartoon})
   broadcast_bdry!(fv_bdry_flux_doublecartoon_tov_mirrored, equation, mesh)
 
   dt = timestep(mesh, P, P.hrsc)
-  if "rD" ∉ P.prmsdb["GRHD"]["freeze_vars"]
+  if :rD ∉ P.prms.freeze_vars
     fv_update_step!(rhs_rD,  rD,  flx_rDx,  src_rD,  bdry_rD,  nflx_rD,  dt, mesh)
   end
-  if "rSx" ∉ P.prmsdb["GRHD"]["freeze_vars"]
+  if :rSx ∉ P.prms.freeze_vars
     fv_update_step!(rhs_rSx, rSx, flx_rSxx, src_rSx, bdry_rSx, nflx_rSx, dt, mesh)
   end
-  if "rτ" ∉ P.prmsdb["GRHD"]["freeze_vars"]
+  if :rτ ∉ P.prms.freeze_vars
     fv_update_step!(rhs_rτ,  rτ,  flx_rτx,  src_rτ,  bdry_rτ,  nflx_rτ,  dt, mesh)
   end
 
-  if !P.prmsdb["GRHD"]["atm_evolve"]
+  if !P.prms.atm_evolve
     broadcast_volume!(determine_atmosphere, P.equation, mesh)
     broadcast_volume!(stop_atmosphere_evolution_doublecartoon, P.equation, mesh)
   end
@@ -272,36 +272,36 @@ function rhs_fv_slope_limiter!(mesh, P::Project{:spherical1d})
 
   # SSPRK3 stepping from Hesthaven 2018 book, sec 10.2 (cf Script 10.13)
   step_fv_slope_limiter!(mesh, P, (D,Sr,Ï„))
-  if "D" ∉ P.prmsdb["GRHD"]["freeze_vars"]
+  if :D ∉ P.prms.freeze_vars
     @. D1  = D  + dt*rhs_D
   end
-  if "Sr" ∉ P.prmsdb["GRHD"]["freeze_vars"]
+  if :Sr ∉ P.prms.freeze_vars
     @. Sr1 = Sr + dt*rhs_Sr
   end
-  if "τ" ∉ P.prmsdb["GRHD"]["freeze_vars"]
+  if :τ ∉ P.prms.freeze_vars
     @. Ï„1  = Ï„  + dt*rhs_Ï„
   end
 
   step_fv_slope_limiter!(mesh, P, (D1,Sr1,Ï„1))
-  if "D" ∉ P.prmsdb["GRHD"]["freeze_vars"]
+  if :D ∉ P.prms.freeze_vars
     @. D2  = (3*D  + D1  + dt*rhs_D)/4
   end
-  if "Sr" ∉ P.prmsdb["GRHD"]["freeze_vars"]
+  if :Sr ∉ P.prms.freeze_vars
     @. Sr2 = (3*Sr + Sr1 + dt*rhs_Sr)/4
   end
-  if "τ" ∉ P.prmsdb["GRHD"]["freeze_vars"]
+  if :τ ∉ P.prms.freeze_vars
     @. Ï„2  = (3*Ï„  + Ï„1  + dt*rhs_Ï„)/4
   end
 
   step_fv_slope_limiter!(mesh, P, (D2,Sr2,Ï„2))
   # atm final result is published via rhs_ variables
-  if "D" ∉ P.prmsdb["GRHD"]["freeze_vars"]
+  if :D ∉ P.prms.freeze_vars
     @. rhs_D  = (D  + 2*D2  + 2*dt*rhs_D)/3
   end
-  if "Sr" ∉ P.prmsdb["GRHD"]["freeze_vars"]
+  if :Sr ∉ P.prms.freeze_vars
     @. rhs_Sr = (Sr + 2*Sr2 + 2*dt*rhs_Sr)/3
   end
-  if "τ" ∉ P.prmsdb["GRHD"]["freeze_vars"]
+  if :τ ∉ P.prms.freeze_vars
     @. rhs_Ï„  = (Ï„  + 2*Ï„2  + 2*dt*rhs_Ï„)/3
   end
 
@@ -355,7 +355,7 @@ function step_fv_slope_limiter!(mesh, P::Project{:spherical1d}, state)
   D .= sD
   Sr .= sSr
   τ .= sτ
-  if P.prmsdb["GRHD"]["c2p_enforce_causal_atm"]
+  if P.prms.c2p_enforce_causal_atm
     broadcast_volume!(cons2prim_spherical1d_freeze_flags, equation, mesh)
   else
     broadcast_volume!(cons2prim_spherical1d, equation, mesh)
@@ -373,8 +373,8 @@ function step_fv_slope_limiter!(mesh, P::Project{:spherical1d}, state)
   L, = dg1d.widths(mesh)
   dl = L/K
 
-  id = P.prmsdb["GRHD"]["id"]
-  bc = P.prmsdb["GRHD"]["bc"]
+  id = P.prms.id
+  bc = P.prms.bc
   if id == "bondi_accretion" || id == "bondi_accretion_infall"
     if bc == "from_id"
       @unpack init_ρ, init_vr, init_ϵ, r = get_static_variables(mesh)
@@ -475,7 +475,7 @@ function step_fv_slope_limiter!(mesh, P::Project{:spherical1d}, state)
   vars = (;bγrr, bα, bβru, bρ, bvr, bϵ, bρL, bvrL, bϵL, bpL,
            bρR, bvrR, bϵR, bpR, bdρ_m, bdρ_p, bdvr_m, bdvr_p, bdϵ_m, bdϵ_p, bde_max_v,
            nflx_D_L, nflx_Sr_L, nflx_Ï„_L, nflx_D_R, nflx_Sr_R, nflx_Ï„_R, eos)
-  numflux_method = P.prmsdb["GRHD"]["fv_numerical_flux"]
+  numflux_method = P.prms.fv_numerical_flux
   if numflux_method == "llf"
     numflux_local_lax_friedrich!(vars, K, P)
   else
@@ -484,17 +484,17 @@ function step_fv_slope_limiter!(mesh, P::Project{:spherical1d}, state)
 
   # evaluate rhs
   # we have a + in the braces, because we use the nx with fluxes above
-  if "D" ∉ P.prmsdb["GRHD"]["freeze_vars"]
+  if :D ∉ P.prms.freeze_vars
     @turbo for k in 1:K
       rhs_D[k]  = -(nflx_D_R[k]  + nflx_D_L[k] )/dl + src_D[k]
     end
   end
-  if "Sr" ∉ P.prmsdb["GRHD"]["freeze_vars"]
+  if :Sr ∉ P.prms.freeze_vars
     @turbo for k in 1:K
       rhs_Sr[k] = -(nflx_Sr_R[k] + nflx_Sr_L[k])/dl + src_Sr[k]
     end
   end
-  if "τ" ∉ P.prmsdb["GRHD"]["freeze_vars"]
+  if :τ ∉ P.prms.freeze_vars
     @turbo for k in 1:K
       rhs_Ï„[k]  = -(nflx_Ï„_R[k]  + nflx_Ï„_L[k] )/dl + src_Ï„[k]
     end
@@ -507,7 +507,7 @@ function step_fv_slope_limiter!(mesh, P::Project{:spherical1d}, state)
     rhs_D[1] = rhs_Sr[1] = rhs_Ï„[1] = rhs_D[end] = rhs_Sr[end] = rhs_Ï„[end] = 0
   end
 
-  if !P.prmsdb["GRHD"]["atm_evolve"]
+  if !P.prms.atm_evolve
     broadcast_volume!(determine_atmosphere, P.equation, mesh)
     broadcast_volume!(stop_atmosphere_evolution_spherical1d, P.equation, mesh)
   end
@@ -578,36 +578,36 @@ function rhs_fv_slope_limiter!(mesh, P::Project{:doublecartoon})
 
   # SSPRK3 stepping from Hesthaven 2018 book, sec 10.2 (cf Script 10.13)
   step_fv_slope_limiter!(mesh, P, (rD,rSx,rτ))
-  if "rD" ∉ P.prmsdb["GRHD"]["freeze_vars"]
+  if :rD ∉ P.prms.freeze_vars
     @. rD1  = rD  + dt*rhs_rD
   end
-  if "rSx" ∉ P.prmsdb["GRHD"]["freeze_vars"]
+  if :rSx ∉ P.prms.freeze_vars
     @. rSx1 = rSx + dt*rhs_rSx
   end
-  if "rτ" ∉ P.prmsdb["GRHD"]["freeze_vars"]
+  if :rτ ∉ P.prms.freeze_vars
     @. rτ1  = rτ  + dt*rhs_rτ
   end
 
   step_fv_slope_limiter!(mesh, P, (rD1,rSx1,rτ1))
-  if "rD" ∉ P.prmsdb["GRHD"]["freeze_vars"]
+  if :rD ∉ P.prms.freeze_vars
     @. rD2  = (3*rD  + rD1  + dt*rhs_rD)/4
   end
-  if "rSx" ∉ P.prmsdb["GRHD"]["freeze_vars"]
+  if :rSx ∉ P.prms.freeze_vars
     @. rSx2 = (3*rSx + rSx1 + dt*rhs_rSx)/4
   end
-  if "rτ" ∉ P.prmsdb["GRHD"]["freeze_vars"]
+  if :rτ ∉ P.prms.freeze_vars
     @. rτ2  = (3*rτ  + rτ1  + dt*rhs_rτ)/4
   end
 
   step_fv_slope_limiter!(mesh, P, (rD2,rSx2,rτ2))
   # atm final result is published via rhs_ variables
-  if "rD" ∉ P.prmsdb["GRHD"]["freeze_vars"]
+  if :rD ∉ P.prms.freeze_vars
     @. rhs_rD  = (rD  + 2*rD2  + 2*dt*rhs_rD)/3
   end
-  if "rSx" ∉ P.prmsdb["GRHD"]["freeze_vars"]
+  if :rSx ∉ P.prms.freeze_vars
     @. rhs_rSx = (rSx + 2*rSx2 + 2*dt*rhs_rSx)/3
   end
-  if "rτ" ∉ P.prmsdb["GRHD"]["freeze_vars"]
+  if :rτ ∉ P.prms.freeze_vars
     @. rhs_rτ  = (rτ  + 2*rτ2  + 2*dt*rhs_rτ)/3
   end
 
@@ -648,7 +648,7 @@ function step_fv_slope_limiter!(mesh, P::Project{:doublecartoon}, state)
   rSx .= sSx
   rτ .= sτ
   broadcast_volume!(unscale_conservatives_doublecartoon, equation, mesh)
-  if P.prmsdb["GRHD"]["c2p_enforce_causal_atm"]
+  if P.prms.c2p_enforce_causal_atm
     broadcast_volume!(cons2prim_doublecartoon_freeze_flags, equation, mesh)
   else
     broadcast_volume!(cons2prim_doublecartoon, equation, mesh)
@@ -667,7 +667,7 @@ function step_fv_slope_limiter!(mesh, P::Project{:doublecartoon}, state)
   L, = dg1d.widths(mesh)
   dl = L/K
 
-  bc = P.prmsdb["GRHD"]["bc"]
+  bc = P.prms.bc
   if bc == "tov_symmetric_domain"
     bdry_γxx_l, bdry_γxx_ll           = γxx[1],        γxx[2]
     bdry_α_l, bdry_α_ll               = α[1],          α[2]
@@ -748,7 +748,7 @@ function step_fv_slope_limiter!(mesh, P::Project{:doublecartoon}, state)
   vars = (;bγxx, bα, bβux, bsqrtdetγ, bρ, bvx, bϵ, bρL, bvxL, bϵL, bpL,
            bρR, bvxR, bϵR, bpR, bdρ_m, bdρ_p, bdvx_m, bdvx_p, bdϵ_m, bdϵ_p, bde_max_v,
            nflx_rD_L, nflx_rSx_L, nflx_rτ_L, nflx_rD_R, nflx_rSx_R, nflx_rτ_R, eos)
-  numflux_method = P.prmsdb["GRHD"]["fv_numerical_flux"]
+  numflux_method = P.prms.fv_numerical_flux
   if numflux_method == "llf"
     numflux_local_lax_friedrich!(vars, K, P)
   elseif numflux_method == "roe"
@@ -761,17 +761,17 @@ function step_fv_slope_limiter!(mesh, P::Project{:doublecartoon}, state)
 
   # evaluate rhs
   # we have a + in the braces, because we use the nx with fluxes above
-  if "rD" ∉ P.prmsdb["GRHD"]["freeze_vars"]
+  if :rD ∉ P.prms.freeze_vars
     @turbo for k in 1:K
       rhs_rD[k]  = -(nflx_rD_R[k]  + nflx_rD_L[k] )/dl + src_rD[k]
     end
   end
-  if "rSr" ∉ P.prmsdb["GRHD"]["freeze_vars"]
+  if :rSr ∉ P.prms.freeze_vars
     @turbo for k in 1:K
       rhs_rSx[k] = -(nflx_rSx_R[k] + nflx_rSx_L[k])/dl + src_rSx[k]
     end
   end
-  if "rτ" ∉ P.prmsdb["GRHD"]["freeze_vars"]
+  if :rτ ∉ P.prms.freeze_vars
     @turbo for k in 1:K
       rhs_rτ[k]  = -(nflx_rτ_R[k]  + nflx_rτ_L[k] )/dl + src_rτ[k]
     end
@@ -784,7 +784,7 @@ function step_fv_slope_limiter!(mesh, P::Project{:doublecartoon}, state)
     rhs_rD[1] = rhs_rSx[1] = rhs_rτ[1] = rhs_rD[end] = rhs_rSx[end] = rhs_rτ[end] = 0
   end
 
-  if !P.prmsdb["GRHD"]["atm_evolve"]
+  if !P.prms.atm_evolve
     broadcast_volume!(determine_atmosphere, P.equation, mesh)
     broadcast_volume!(stop_atmosphere_evolution_doublecartoon, P.equation, mesh)
   end
@@ -1186,7 +1186,7 @@ function rhs!(mesh::Mesh1d, P::Project{:valencia1d}, hrsc::Nothing)
 
   broadcast_volume!(unscale_conservatives_valencia1d, equation, mesh)
   impose_symmetry!(P, mesh)
-  if P.prmsdb["GRHD"]["c2p_enforce_causal_atm"]
+  if P.prms.c2p_enforce_causal_atm
     broadcast_volume!(cons2prim_valencia1d_freeze_flags, equation, mesh)
   else
     broadcast_volume!(cons2prim_valencia1d, equation, mesh)
@@ -1220,7 +1220,7 @@ function rhs!(mesh::Mesh1d, P::Project{:valencia1d}, hrsc::Nothing)
   compute_rhs_weak_form!(rhs_rSr, flr_rSr, src_rSr, nflr_rSr, mesh)
   compute_rhs_weak_form!(rhs_rτ,  flr_rτ,  src_rτ,  nflr_rτ,  mesh)
 
-  if !P.prmsdb["GRHD"]["atm_evolve"]
+  if !P.prms.atm_evolve
     broadcast_volume!(determine_atmosphere, P.equation, mesh)
     broadcast_volume!(stop_atmosphere_evolution_valencia1d, P.equation, mesh)
   end
@@ -1243,7 +1243,7 @@ function rhs!(mesh::Mesh1d, P::Project{:spherical1d}, hrsc::Nothing)
           bdry_D, bdry_Sr, bdry_Ï„,
           bdry_max_v, bdry_vr, bdry_p = get_bdry_variables(cache)
 
-  if P.prmsdb["GRHD"]["c2p_enforce_causal_atm"]
+  if P.prms.c2p_enforce_causal_atm
     broadcast_volume!(cons2prim_spherical1d_freeze_flags, equation, mesh)
   else
     broadcast_volume!(cons2prim_spherical1d, equation, mesh)
@@ -1277,7 +1277,7 @@ function rhs!(mesh::Mesh1d, P::Project{:spherical1d}, hrsc::Nothing)
   compute_rhs_weak_form!(rhs_Sr, flr_Sr, src_Sr, nflr_Sr, mesh)
   compute_rhs_weak_form!(rhs_Ï„,  flr_Ï„,  src_Ï„,  nflr_Ï„,  mesh)
 
-  if !P.prmsdb["GRHD"]["atm_evolve"]
+  if !P.prms.atm_evolve
     broadcast_volume!(determine_atmosphere, P.equation, mesh)
     broadcast_volume!(stop_atmosphere_evolution_spherical1d, P.equation, mesh)
   end
@@ -1302,7 +1302,7 @@ function rhs!(mesh::Mesh1d, P::Project{:rescaled_spherical1d}, hrsc::Nothing)
 
   broadcast_volume!(unscale_conservatives, equation, mesh)
   impose_symmetry!(P, mesh)
-  if P.prmsdb["GRHD"]["c2p_enforce_causal_atm"]
+  if P.prms.c2p_enforce_causal_atm
     broadcast_volume!(cons2prim_rescaled_spherical1d_freeze_flags, equation, mesh)
   else
     broadcast_volume!(cons2prim_rescaled_spherical1d, equation, mesh)
@@ -1330,7 +1330,7 @@ function rhs!(mesh::Mesh1d, P::Project{:rescaled_spherical1d}, hrsc::Nothing)
   compute_rhs_weak_form!(rhs_rSr, flr_rSr, src_rSr, nflr_rSr, mesh)
   compute_rhs_weak_form!(rhs_rτ,  flr_rτ,  src_rτ,  nflr_rτ,  mesh)
 
-  if !P.prmsdb["GRHD"]["atm_evolve"]
+  if !P.prms.atm_evolve
     broadcast_volume!(determine_atmosphere, P.equation, mesh)
     broadcast_volume!(stop_atmosphere_evolution_rescaled_spherical1d, P.equation, mesh)
   end
@@ -1354,7 +1354,7 @@ function rhs!(mesh::Mesh1d, P::Project{:doublecartoon}, hrsc::Nothing)
           bdry_max_v, bdry_vx, bdry_p = get_bdry_variables(cache)
 
   broadcast_volume!(unscale_conservatives_doublecartoon, equation, mesh)
-  if P.prmsdb["GRHD"]["c2p_enforce_causal_atm"]
+  if P.prms.c2p_enforce_causal_atm
     broadcast_volume!(cons2prim_doublecartoon_freeze_flags, equation, mesh)
   else
     broadcast_volume!(cons2prim_doublecartoon, equation, mesh)
@@ -1392,17 +1392,17 @@ function rhs!(mesh::Mesh1d, P::Project{:doublecartoon}, hrsc::Nothing)
     end
   end
 
-  if "rD" ∉ P.prmsdb["GRHD"]["freeze_vars"]
+  if :rD ∉ P.prms.freeze_vars
     compute_rhs_weak_form!(rhs_rD,  flx_rDx,  src_rD,  nflx_rD,  mesh)
   end
-  if "rSx" ∉ P.prmsdb["GRHD"]["freeze_vars"]
+  if :rSx ∉ P.prms.freeze_vars
     compute_rhs_weak_form!(rhs_rSx, flx_rSxx, src_rSx, nflx_rSx, mesh)
   end
-  if "rτ" ∉ P.prmsdb["GRHD"]["freeze_vars"]
+  if :rτ ∉ P.prms.freeze_vars
     compute_rhs_weak_form!(rhs_rτ,  flx_rτx,  src_rτ,  nflx_rτ,  mesh)
   end
 
-  if !P.prmsdb["GRHD"]["atm_evolve"]
+  if !P.prms.atm_evolve
     broadcast_volume!(determine_atmosphere, P.equation, mesh)
     broadcast_volume!(stop_atmosphere_evolution_doublecartoon, P.equation, mesh)
   end
@@ -1431,7 +1431,7 @@ function rhs!(mesh, P::Project{:spherical1d}, hrsc::HRSC.AbstractReconstruction)
   HRSC.reconstruct!(Sr, flag, hrsc, isperiodic=isperiodic(mesh), limit_with_boundaries=false)
   HRSC.reconstruct!(Ï„,  flag, hrsc, isperiodic=isperiodic(mesh), limit_with_boundaries=false)
 
-  if P.prmsdb["GRHD"]["c2p_enforce_causal_atm"]
+  if P.prms.c2p_enforce_causal_atm
     broadcast_volume!(cons2prim_freeze_flags, equation, mesh)
   else
     broadcast_volume!(cons2prim, equation, mesh)
@@ -1460,7 +1460,7 @@ function rhs!(mesh, P::Project{:spherical1d}, hrsc::HRSC.AbstractReconstruction)
   compute_rhs_weak_form!(rhs_Sr, flr_Sr, src_Sr, nflr_Sr, mesh)
   compute_rhs_weak_form!(rhs_Ï„,  flr_Ï„,  src_Ï„,  nflr_Ï„,  mesh)
 
-  if !P.prmsdb["GRHD"]["atm_evolve"]
+  if !P.prms.atm_evolve
     broadcast_volume!(determine_atmosphere, P.equation, mesh)
     broadcast_volume!(stop_atmosphere_evolution_spherical1d, P.equation, mesh)
   end
@@ -1490,7 +1490,7 @@ function rhs!(mesh, P::Project{:doublecartoon}, hrsc::HRSC.AbstractReconstructio
   # HRSC.reconstruct!(rτ,  flag, hrsc, isperiodic=isperiodic(mesh), limit_with_boundaries=false)
 
   broadcast_volume!(unscale_conservatives_doublecartoon, equation, mesh)
-  if P.prmsdb["GRHD"]["c2p_enforce_causal_atm"]
+  if P.prms.c2p_enforce_causal_atm
     broadcast_volume!(cons2prim_doublecartoon_freeze_flags, equation, mesh)
   else
     broadcast_volume!(cons2prim_doublecartoon, equation, mesh)
@@ -1528,17 +1528,17 @@ function rhs!(mesh, P::Project{:doublecartoon}, hrsc::HRSC.AbstractReconstructio
     end
   end
 
-  if "rD" ∉ P.prmsdb["GRHD"]["freeze_vars"]
+  if :rD ∉ P.prms.freeze_vars
     compute_rhs_weak_form!(rhs_rD,  flx_rDx,  src_rD,  nflx_rD,  mesh)
   end
-  if "rSx" ∉ P.prmsdb["GRHD"]["freeze_vars"]
+  if :rSx ∉ P.prms.freeze_vars
     compute_rhs_weak_form!(rhs_rSx, flx_rSxx, src_rSx, nflx_rSx, mesh)
   end
-  if "rτ" ∉ P.prmsdb["GRHD"]["freeze_vars"]
+  if :rτ ∉ P.prms.freeze_vars
     compute_rhs_weak_form!(rhs_rτ,  flx_rτx,  src_rτ,  nflx_rτ,  mesh)
   end
 
-  if !P.prmsdb["GRHD"]["atm_evolve"]
+  if !P.prms.atm_evolve
     broadcast_volume!(determine_atmosphere, P.equation, mesh)
     broadcast_volume!(stop_atmosphere_evolution_doublecartoon, P.equation, mesh)
   end
@@ -1568,7 +1568,7 @@ function rhs!(mesh, P::Project{:valencia1d}, hrsc::HRSC.AbstractArtificialViscos
 
   broadcast_volume!(unscale_conservatives_valencia1d, equation, mesh)
   impose_symmetry!(P, mesh)
-  if P.prmsdb["GRHD"]["c2p_enforce_causal_atm"]
+  if P.prms.c2p_enforce_causal_atm
     broadcast_volume!(cons2prim_valencia1d_freeze_flags, equation, mesh)
   else
     broadcast_volume!(cons2prim_valencia1d, equation, mesh)
@@ -1632,7 +1632,7 @@ function rhs!(mesh, P::Project{:valencia1d}, hrsc::HRSC.AbstractArtificialViscos
   compute_rhs_weak_form!(rhs_rSr, flr_rSr, src_rSr, nflr_rSr, mesh)
   compute_rhs_weak_form!(rhs_rτ,  flr_rτ,  src_rτ,  nflr_rτ,  mesh)
 
-  if !P.prmsdb["GRHD"]["atm_evolve"]
+  if !P.prms.atm_evolve
     broadcast_volume!(determine_atmosphere, P.equation, mesh)
     broadcast_volume!(stop_atmosphere_evolution_valencia1d, P.equation, mesh)
   end
@@ -1659,7 +1659,7 @@ function rhs!(mesh, P::Project{:spherical1d}, hrsc::HRSC.AbstractArtificialVisco
           bdry_rhs_D, bdry_rhs_Sr, bdry_rhs_Ï„,
           bdry_smoothed_mu                    = get_bdry_variables(cache)
 
-  if P.prmsdb["GRHD"]["c2p_enforce_causal_atm"]
+  if P.prms.c2p_enforce_causal_atm
     broadcast_volume!(cons2prim_spherical1d_freeze_flags, equation, mesh)
   else
     broadcast_volume!(cons2prim_spherical1d, equation, mesh)
@@ -1676,8 +1676,8 @@ function rhs!(mesh, P::Project{:spherical1d}, hrsc::HRSC.AbstractArtificialVisco
   broadcast_volume!(flux_source_spherical1d, equation, mesh)
   impose_symmetry_sources!(P, mesh)
   broadcast_faces!(llf_spherical1d, equation, mesh)
-  id = P.prmsdb["GRHD"]["id"]
-  bc = P.prmsdb["GRHD"]["bc"]
+  id = P.prms.id
+  bc = P.prms.bc
   if id == "bondi_accretion" || id == "bondi_accretion_infall"
     if bc == "from_id"
       # already accounts for inflow boundary
@@ -1697,13 +1697,13 @@ function rhs!(mesh, P::Project{:spherical1d}, hrsc::HRSC.AbstractArtificialVisco
     TODO("unknown boundary conditions for initial data $id: $bc")
   end
 
-  if "D" ∉ P.prmsdb["GRHD"]["freeze_vars"]
+  if :D ∉ P.prms.freeze_vars
     compute_rhs_weak_form!(rhs_D,  flr_D,  src_D,  nflr_D,  mesh)
   end
-  if "Sr" ∉ P.prmsdb["GRHD"]["freeze_vars"]
+  if :Sr ∉ P.prms.freeze_vars
     compute_rhs_weak_form!(rhs_Sr, flr_Sr, src_Sr, nflr_Sr, mesh)
   end
-  if "τ" ∉ P.prmsdb["GRHD"]["freeze_vars"]
+  if :τ ∉ P.prms.freeze_vars
     compute_rhs_weak_form!(rhs_Ï„,  flr_Ï„,  src_Ï„,  nflr_Ï„,  mesh)
   end
 
@@ -1748,17 +1748,17 @@ function rhs!(mesh, P::Project{:spherical1d}, hrsc::HRSC.AbstractArtificialVisco
 
   broadcast_faces!(av_nflux_covariant_spherical1d, equation, mesh)
 
-  if "D" ∉ P.prmsdb["GRHD"]["freeze_vars"]
+  if :D ∉ P.prms.freeze_vars
     compute_rhs_weak_form!(rhs_D,  flr_D,  src_D,  nflr_D,  mesh)
   end
-  if "Sr" ∉ P.prmsdb["GRHD"]["freeze_vars"]
+  if :Sr ∉ P.prms.freeze_vars
     compute_rhs_weak_form!(rhs_Sr, flr_Sr, src_Sr, nflr_Sr, mesh)
   end
-  if "τ" ∉ P.prmsdb["GRHD"]["freeze_vars"]
+  if :τ ∉ P.prms.freeze_vars
     compute_rhs_weak_form!(rhs_Ï„,  flr_Ï„,  src_Ï„,  nflr_Ï„,  mesh)
   end
 
-  if !P.prmsdb["GRHD"]["atm_evolve"]
+  if !P.prms.atm_evolve
     broadcast_volume!(determine_atmosphere, P.equation, mesh)
     broadcast_volume!(stop_atmosphere_evolution_spherical1d, P.equation, mesh)
   end
@@ -1804,7 +1804,7 @@ function rhs!(mesh::Mesh1d, P::Project{:doublecartoon}, hrsc::AbstractArtificial
           bdry_rhs_rD, bdry_rhs_rSx, bdry_rhs_rτ = get_bdry_variables(cache)
 
   broadcast_volume!(unscale_conservatives_doublecartoon, equation, mesh)
-  if P.prmsdb["GRHD"]["c2p_enforce_causal_atm"]
+  if P.prms.c2p_enforce_causal_atm
     broadcast_volume!(cons2prim_doublecartoon_freeze_flags, equation, mesh)
   else
     broadcast_volume!(cons2prim_doublecartoon, equation, mesh)
@@ -1842,13 +1842,13 @@ function rhs!(mesh::Mesh1d, P::Project{:doublecartoon}, hrsc::AbstractArtificial
     end
   end
 
-  if "rD" ∉ P.prmsdb["GRHD"]["freeze_vars"]
+  if :rD ∉ P.prms.freeze_vars
     compute_rhs_weak_form!(rhs_rD,  flx_rDx,  src_rD,  nflx_rD,  mesh)
   end
-  if "rSx" ∉ P.prmsdb["GRHD"]["freeze_vars"]
+  if :rSx ∉ P.prms.freeze_vars
     compute_rhs_weak_form!(rhs_rSx, flx_rSxx, src_rSx, nflx_rSx, mesh)
   end
-  if "rτ" ∉ P.prmsdb["GRHD"]["freeze_vars"]
+  if :rτ ∉ P.prms.freeze_vars
     compute_rhs_weak_form!(rhs_rτ,  flx_rτx,  src_rτ,  nflx_rτ,  mesh)
   end
 
@@ -1896,17 +1896,17 @@ function rhs!(mesh::Mesh1d, P::Project{:doublecartoon}, hrsc::AbstractArtificial
 
   broadcast_faces!(av_nflux_covariant_doublecartoon, equation, mesh)
 
-  if "rD" ∉ P.prmsdb["GRHD"]["freeze_vars"]
+  if :rD ∉ P.prms.freeze_vars
     compute_rhs_weak_form!(rhs_rD,  flx_rDx,  src_rD,  nflx_rD,  mesh)
   end
-  if "rSx" ∉ P.prmsdb["GRHD"]["freeze_vars"]
+  if :rSx ∉ P.prms.freeze_vars
     compute_rhs_weak_form!(rhs_rSx, flx_rSxx, src_rSx, nflx_rSx, mesh)
   end
-  if "rτ" ∉ P.prmsdb["GRHD"]["freeze_vars"]
+  if :rτ ∉ P.prms.freeze_vars
     compute_rhs_weak_form!(rhs_rτ,  flx_rτx,  src_rτ,  nflx_rτ,  mesh)
   end
 
-  if !P.prmsdb["GRHD"]["atm_evolve"]
+  if !P.prms.atm_evolve
     broadcast_volume!(determine_atmosphere, P.equation, mesh)
     broadcast_volume!(stop_atmosphere_evolution_doublecartoon, P.equation, mesh)
   end
@@ -1938,7 +1938,7 @@ function rhs!(mesh::Mesh2d, P::Project{:cartoon})
           bdry_max_v,  bdry_vx, bdry_vz, bdry_p = get_bdry_variables(cache)
 
   broadcast_volume!(unscale_conservatives_cartoon, equation, mesh)
-  if P.prmsdb["GRHD"]["c2p_enforce_causal_atm"]
+  if P.prms.c2p_enforce_causal_atm
     broadcast_volume!(cons2prim_cartoon_freeze_flags, equation, mesh)
   else
     broadcast_volume!(cons2prim_cartoon, equation, mesh)
@@ -1966,7 +1966,7 @@ function rhs!(mesh::Mesh2d, P::Project{:cartoon})
   compute_rhs_weak_form!(rhs_rSz, flx_rSzx, flx_rSzz, src_rSz, nflx_rSz, mesh)
   compute_rhs_weak_form!(rhs_rτ,  flx_rτx,  flx_rτz,  src_rτ,  nflx_rτ,  mesh)
 
-  if !P.prmsdb["GRHD"]["atm_evolve"]
+  if !P.prms.atm_evolve
     broadcast_volume!(determine_atmosphere, P.equation, mesh)
     broadcast_volume!(stop_atmosphere_evolution_cartoon, P.equation, mesh)
   end
diff --git a/src/GRHD/setup.jl b/src/GRHD/setup.jl
index 2c9f178dbf44b4041b42554e6cad44a984fa8238..2dc1b712a426b325bde49fdfedec0bd46e3fcc1f 100644
--- a/src/GRHD/setup.jl
+++ b/src/GRHD/setup.jl
@@ -34,11 +34,20 @@ function Project(env::Environment, prms)
   cold_K, cold_Γ = prms["GRHD"]["c2p_cold_eos_parameters"]
   cold_eos = Polytrope(cold_K, cold_Γ)
   c2p_set_atmosphere_on_failure = prms["GRHD"]["c2p_set_atmosphere_on_failure"]
+  c2p_enforce_atm = Bool(prms["GRHD"]["c2p_enforce_atm"])
+  c2p_enforce_causal_atm = Bool(prms["GRHD"]["c2p_enforce_causal_atm"])
+  atm_evolve = Bool(prms["GRHD"]["atm_evolve"])
+  hrsc = Symbol(prms["GRHD"]["hrsc"])
+  freeze_vars = Symbol.(prms["GRHD"]["freeze_vars"])
+  fv_numerical_flux = prms["GRHD"]["fv_numerical_flux"]
+  id = prms["GRHD"]["id"]
+  bc = prms["GRHD"]["bc"]
   fixedprms = (; av_regularization=:covariant, id_smooth=true,
                  bernstein, slope_limiter_method, slope_limiter_tvb_M,
-                 c2p_dynamic_atm,
-                 c2p_set_atmosphere_on_failure,
-                 av_drag, problem)
+                 c2p_dynamic_atm, atm_evolve,
+                 c2p_set_atmosphere_on_failure, c2p_enforce_causal_atm, c2p_enforce_atm,
+                 av_drag, problem, freeze_vars,
+                 id, bc, hrsc, fv_numerical_flux)
 
   # construct Project
   eos = EquationOfState.make_EquationOfState(prms["EquationOfState"])