dg1d: add lazy integrator (!198)

I am surprised that I did not need this before.

This is intentionally a hacky implementation as I just need that for the thesis right now.

---

This also implements
- the total rest-mass computation for the GRHD project as a first application,
- adds helper functions for the `Output` struct to read all data at once.

src/GRHD/GRHD.jl

@@ -292,6 +292,18 @@ dg1d.@parameters GRHD begin
   """
   freeze_vars = String[]
 
+  """
+  List of 0d variables to analyze during simulations.
+  Variable names enequeued here become available for `Output.variables0d`.
+
+  Available variables:
+  - `Mtot`: Total restmass density, Mtot = ∫ √(γ) D dΩ.
+  """
+  variables0d_analyze = String[]
+  @check variables0d_analyze isa Vector{String}
+  @check allunique(variables0d_analyze)
+  @check variables0d_analyze ⊂ ["Mtot"]
+
   """
   If `true` then the doublecartoon derivative terms are not accounted for in the source terms.
   This appears to be necessary in order to perform TOV star simulations with the doublecartoon

src/GRHD/callbacks.jl

@@ -21,6 +21,12 @@ function make_callback(env, P)
                                                  CallbackTiming(every_iteration=1,every_dt=0))
     push!(callbackset, cb)
   end
+  if !isempty(P.prmsdb["GRHD"]["variables0d_analyze"])
+    cbfn_analyzed0d(u, t)     = callback_analyze0d(env, P, env.mesh)
+    cb                        = FunctionCallback(cbfn_analyzed0d,
+                                                 CallbackTiming(every_iteration=1,every_dt=0))
+    push!(callbackset, cb)
+  end
   return callbackset
 end
 
@@ -475,3 +481,22 @@ function callback_analyze1d_error(env, P, mesh::Mesh1d)
     @. v_err = v .- v_ref
   end
 end
+
+
+callback_analyze0d(env, P::Project, mesh) = TODO()
+function callback_analyze0d(env, P::Project{:spherical1d}, mesh::Mesh1d)
+  vars = P.prmsdb["GRHD"]["variables0d_analyze"]
+  if "Mtot" in vars
+    @unpack D = get_dynamic_variables(mesh)
+    @unpack γrr, γθθ, γϕϕ, r = get_static_variables(mesh)
+    @unpack Mtot = get_global_variables(mesh)
+    Mtot .= dg1d.integrate(mesh) do i
+      D[i] * sqrt(γrr[i]*γθθ[i]*γϕϕ[i]) * 4*π
+    end
+    if dg1d.is_domain_symmetric(mesh)
+      # in this case we integrated the sphere twice, so halving the result
+      # should correspond to an averaging of the integral from either side
+      Mtot ./= 2
+    end
+  end
+end

src/GRHD/initialdata.jl

@@ -841,6 +841,16 @@ function initialdata_equation!(::Val{:tov}, env, P::Project{:spherical1d}, prms)
   @. init_ρ     = ρ
   @. init_ϵ     = ϵ
 
+  Mtot = dg1d.integrate(mesh) do i
+    v = D[i] * sqrt(γrr[i]*γθθ[i]*γϕϕ[i]) * 4*π
+    v
+  end
+  if dg1d.is_domain_symmetric(mesh)
+    Mtot /= 2
+  end
+  Madm = data.M
+  @show Mtot, Madm
+
 end
 
 

src/GRHD/setup.jl

@@ -481,6 +481,9 @@ function register_analysis!(mesh::Mesh1d, P)
       static_variablenames = (:c2p_reset_ϵ, :c2p_reset_atm, :c2p_limit_vr,
                               :c2p_freeze_atm, :c2p_init_admissible, :v)
   )
+  if "Mtot" in P.prmsdb["GRHD"]["variables0d_analyze"]
+    register_variables!(mesh, global_variablenames = (:Mtot,))
+  end
 end
 function register_analysis!(mesh::Mesh2d, P)
   register_variables!(mesh,
@@ -488,6 +491,10 @@ function register_analysis!(mesh::Mesh2d, P)
                               :c2p_freeze_atm, :c2p_init_admissible, :v,
                               :buf_c2p_reset_atm,)
   )
+  if "Mtot" in P.prmsdb["GRHD"]["variables0d_analyze"]
+    TODO()
+    register_variables!(mesh, global_variablenames = (:Mtot,))
+  end
 end
 
 

src/glgl.jl

@@ -62,6 +62,20 @@ function integrate(w, v, D, f)
 
   return sum
 end
+function integrate(w, v, D, fn::Function, idxs::AbstractVector{<:Integer})
+
+  N = length(w)
+  @toggled_assert (N,N) == size(D)
+  @toggled_assert N == length(idxs)
+
+  sum = 0.0
+  for i in 1:N
+    fi = fn(idxs[i])
+    sum += (w[i] + v * (D[1,i] - D[N,i])) * fi
+  end
+
+  return sum
+end
 
 
 function integrate(w, v, D, f1, f2)

src/lgl.jl

@@ -30,6 +30,14 @@ function integrate(w, f)
   @toggled_assert length(w) == length(f)
   return dot(w, f)
 end
+function integrate(w, f::Function, idxs::AbstractVector{<:Integer})
+  @toggled_assert length(w) == length(idxs)
+  sum = 0.0
+  for i in 1:length(w)
+    sum += w[i] * f(idxs[i])
+  end
+  return sum
+end
 
 
 function integrate(w, f, g)

src/mesh.jl

@@ -9,6 +9,7 @@ abstract type AbstractMesh end
 struct Mesh{Element,N_Dim,N_Sides} <: AbstractMesh
   tree::Tree{N_Dim,N_Sides} # abstract representation of mesh
   boxes::Vector{Box{N_Dim}} # physical extends of each cell
+  # TODO Fix spelling extends -> extents
   extends::NTuple{N_Dim,Tuple{Float64,Float64}} # total extend of mesh
   element::Element
   cache::Cache
@@ -286,6 +287,11 @@ function isperiodic(m::Mesh)
   return all(m.tree.periodic)
 end
 
+is_domain_symmetric(m::Mesh) = TODO()
+function is_domain_symmetric(m::Mesh1d)
+  (Lmin, Lmax), = m.extends
+  return Lmin ≈ -Lmax
+end
 
 Base.broadcastable(mesh::Mesh) = Ref(mesh)
 
@@ -389,6 +395,31 @@ function integrate_cell(u, k, mesh::Mesh2d)
 end
 
 
+function integrate(fn::Function, mesh::Mesh1d{SpectralElement})
+  @unpack invdetJ = get_static_variables(mesh)
+  Npts, K = layout(mesh)
+  mat_invdetJ = vreshape(invdetJ, (Npts,K))
+  integral = 0.0
+  for k in 1:K
+    idxs = cellindices(mesh, k)
+    integral += integrate(fn, idxs, mesh.element) / mat_invdetJ[first(idxs)]
+  end
+  return integral
+end
+function integrate(fn::Function, mesh::Mesh1d{FVElement})
+  L, = widths(mesh)
+  K, = mesh.tree.dims
+  @assert K >= 2
+  dL = L/K
+  integral = 0.0
+  fk, fkp1 = fn(1), fn(2)
+  for k in 1:K-1
+    integral += (fk+fkp1)*dL/2
+    fk, fkp1 = fkp1, fn(k+1)
+  end
+  return integral
+end
+
 function integrate_fv_central(integrand::Function, δu, uavg, mesh::Mesh{<:FVElement})
   # assume vanilla FV method so that solution in jth cell [x_(j-1/2),x_(j+1/2)] can
   # be approximated by

src/output.jl

@@ -61,7 +61,7 @@ Base.length(o::Output) = length(o.times)
 has_variable(o::Output, variable) = variable in o.variables
 
 
-function _smart_read_variable!(storage, h5_data)
+function _smart_read_variable!(storage::AbstractArray{<:Number}, h5_data)
   # figure out compatibility between storage and data
   size_data      = size(h5_data)
   ndims_data     = ndims(h5_data)
@@ -131,7 +131,7 @@ Similar to `read_variable`, but instead read data into `storage`.
 This read is smart in the sense that it recognizes the shape of the `storage` and reads
 the data accordingly.
 """
-function read_variable!(storage, o::Output, variable::AbstractString, t::Real)
+function read_variable!(storage::AbstractArray{<:Number}, o::Output, variable::AbstractString, t::Real)
   if !has_variable(o, variable)
     error("Output does not contain data for 'variable = $variable'")
   end
@@ -153,7 +153,7 @@ Similar to `read_variable`, but instead read data into `storage`.
 This read is smart in the sense that it recognizes the shape of the `storage` and reads
 the data accordingly.
 """
-function read_variable!(storage, o::Output, variable::AbstractString, index_t::Int)
+function read_variable!(storage::AbstractArray{<:Number}, o::Output, variable::AbstractString, index_t::Int)
   if !has_variable(o, variable)
     error("Output does not contain data for 'variable = $variable'")
   end
@@ -165,3 +165,28 @@ function read_variable!(storage, o::Output, variable::AbstractString, index_t::I
   _smart_read_variable!(storage, h5_variable)
   return
 end
+
+
+function read_variable(o::Output, variable::AbstractString)
+  var1 = read_variable(o, variable, 1)
+  N = length(o.times)
+  allvars = Vector{typeof(var1)}(undef, N)
+  allvars[1] = var1
+  for i in 2:N
+    allvars[i] = read_variable(o, variable, i)
+  end
+  return allvars
+end
+
+
+function read_variable!(storage::AbstractVector{<:AbstractArray{<:Number}}, o::Output,
+                        variable::AbstractString)
+  N = length(o.times)
+  if length(storage) != N
+    resize!(storage, N)
+  end
+  for i in 1:N
+    read_variable!(storage[i], o, variable, i)
+  end
+  return
+end

src/spectralelement.jl

@@ -111,6 +111,17 @@ end
 end
 
 
+@inline function integrate(fn::Function, idxs::AbstractVector{<:Integer}, el::SpectralElement)
+  @unpack quadr = el
+  @toggled_assert quadr in (:LGL, :LGL_lumped, :GLGL)
+  if quadr === :LGL || quadr === :LGL_lumped
+    return LGL.integrate(el.w, fn, idxs)
+  else # quadr === :GLGL
+    return GLGL.integrate(el.w, el.v, el.D, fn, idxs)
+  end
+end
+
+
 @inline function differentiate(u, el::SpectralElement)
   return el.D * u
 end

test/IntegrationTests/refs/grhd_bondi_accretion/grhd_bondi_accretion.toml

 [EquationOfState]
-eos = "polytrope"
-polytrope_k = 0.02143872868864732
 polytrope_gamma = "$(4/3)"
+polytrope_k = 0.02143872868864732
+eos = "polytrope"
+
+[Evolution]
+cfl = 0.2
+tend = 5.0
+
+[Output]
+variables0d = ["D_err_norm", "Mtot"]
+aligned_ts = "$(collect(range(0.5,5.0,step=0.5)))"
+variables1d = ["D_err", "D", "Sr", "τ", "p", "ρ", "ϵ", "vr"]
 
 [GRHD]
-bc = "from_id"
-id = "bondi_accretion"
-id_filename = "$(joinpath(ROOTDIR,\"initialdata\",\"bondi_accretion.h5\"))"
-c2p_enforce_causal_atm = false
 analyze_error_reference_solution = "bondi_accretion"
+id_filename = "$(joinpath(ROOTDIR,\"initialdata\",\"bondi_accretion.h5\"))"
+variables0d_analyze = ["Mtot"]
+id = "bondi_accretion"
 variables0d_analyze_error = ["D"]
 variables1d_analyze_error = ["D"]
+c2p_enforce_causal_atm = false
+bc = "from_id"
 
 [Mesh]
-range = [ 1.8, 10.0 ]
-n = 4
+periodic = false
+range = [1.8, 10.0]
 k = 35
 basis = "lgl"
-periodic = false
-
-[Output]
-variables0d = [ "D_err_norm" ]
-variables1d = [ "D_err", "D", "Sr", "τ", "p", "ρ", "ϵ", "vr" ]
-aligned_ts  = "$(collect(range(0.5,5.0,step=0.5)))"
-
-[Evolution]
-cfl = 0.2
-tend = 5.0
+n = 4

test/IntegrationTests/refs/testconfig.toml

@@ -123,7 +123,7 @@ variables1d = [ "D", "S", "tau" ]
 reltol1d = 1e-6
 
 [grhd_bondi_accretion]
-variables0d  = [ "D_err_norm" ]
+variables0d  = [ "D_err_norm", "Mtot" ]
 variables1d  = [ "D_err", "D", "Sr", "τ", "p", "ρ", "ϵ", "vr" ]
 
 [grhd_bondi_infall_avmda]