diff --git a/examples/growth.ipynb b/examples/growth.ipynb
index a7a6e65..10de08c 100644
--- a/examples/growth.ipynb
+++ b/examples/growth.ipynb
@@ -697,7 +697,7 @@
     "output_dir = \"$(pwd())/output\"\n",
     "num_saves = 100\n",
     "output_times = t_init:(t_final-t_init)/num_saves:t_final\n",
-    "output_config = OutputConfig(output_dir, output_times, summary_statistics=UltraDark.Output.SummaryStatisticsMeanMaxRms)\n",
+    "output_config = OutputConfig(output_dir, output_times, summary_statistics=(Summary.SimulationTime,Summary.MeanDensity,Summary.MaxDensity,Summary.RmsDensityContrast,Summary.ScaleFactor,Summary.TimeStep))\n",
     "\n",
     "options = Config.SimulationConfig(a=a)"
    ]
diff --git a/src/UltraDark.jl b/src/UltraDark.jl
index 59438f9..cf93037 100644
--- a/src/UltraDark.jl
+++ b/src/UltraDark.jl
@@ -12,19 +12,19 @@ export Grids, PencilGrids
 export dV
 export Config, SimulationConfig, constant_scale_factor, TimeStepOptions
 export OutputConfig
-export SummaryStatistics, SummaryStatisticsMeanMaxRms
+export Summary
 
 include("grids.jl")
 include("pencil_grids.jl")
 include("phase_diff.jl")
 include("time_step.jl")
+include("summary.jl")
 include("output.jl")
 include("config.jl")
 
 import .Output: OutputConfig
 import .Output: output_grids, output_xyz
 import .Output: output_summary_row, output_summary_header
-import .Output.Summary: SummaryStatistics, SummaryStatisticsMeanMaxRms
 import .Config: SimulationConfig, constant_scale_factor, TimeStepOptions
 
 """
diff --git a/src/output.jl b/src/output.jl
index 0fd3237..c4f00a3 100644
--- a/src/output.jl
+++ b/src/output.jl
@@ -5,9 +5,7 @@ import PencilFFTs, MPI
 import NPZ
 using Statistics
 
-include("summary.jl")
-
-import .Summary: output_summary_header, output_summary_row
+import ..Summary: output_summary_header, output_summary_row
 
 export Summary
 export OutputConfig
@@ -38,13 +36,13 @@ struct OutputConfig
     rho::Bool
 
     "Type of summary statistics to collect"
-    summary_statistics::DataType
+    summary_statistics::Tuple
 end
 
 function OutputConfig(
     directory, output_times;
     box=true, slice=false, psi=true, rho=true,
-    summary_statistics=SummaryStatistics,
+    summary_statistics=(Summary.WallTime, Summary.SimulationTime,)
 )
 
     OutputConfig(directory, output_times, box, slice, psi, rho, summary_statistics)
diff --git a/src/summary.jl b/src/summary.jl
index a785d65..5ed86c1 100644
--- a/src/summary.jl
+++ b/src/summary.jl
@@ -5,7 +5,11 @@ import Dates
 import MPI
 using Statistics
 
-export SummaryStatistics, SummaryStatisticsMeanMaxRms
+export WallTime
+export SimulationTime
+export ScaleFactor
+export TimeStep
+export MeanDensity, MaxDensity, RmsDensityContrast
 
 """
     output_summary_header(output_config)
@@ -21,147 +25,319 @@ function output_summary_header(output_config)
     end
 end
 
-function generate_summary_row(summary)::String
-    line = "$(getfield(summary, 1))"
-    for i in 2:nfields(summary)
-        line = line * ",$(getfield(summary, i))"
-    end
-    line = line * "\n"
-end
-
 """
     output_summary_row(grids, output_config, t, a, Δt)
 
 Write a new row to the summary file
 """
 function output_summary_row(grids, output_config, t, a, Δt, constants)
-    summary = output_config.summary_statistics(Dates.now(), t, a, Δt, grids, constants)
-    line = generate_summary_row(summary)
+    summaries = map_summary_statistics(output_config.summary_statistics, t, a, Δt, grids, constants)
+    line = generate_summary_row(summaries)
     open(joinpath(output_config.directory, "summary.csv"), "a") do file
         write(file, line)
     end
 end
 
 """
-    output_summary_row(grids::PencilGrids, output_config, t, a, Δt)
+    map_summary_statistics(summary_statistics, sim_time, a, Δt, grids, constants)
 
-Write a new row to the summary file
+Calculate each summary statistic in summary_statistics
 
 If the grids is a PencilGrids, this uses `MPI.Reduce` to compute partial
 summaries in each task and combine them.
 """
-function output_summary_row(grids::PencilGrids, output_config, t, a, Δt, constants)
+
+function map_summary_statistics(summary_statistics, sim_time, a, Δt, grids, constants)
+    summaries = map(x->x(sim_time, a, Δt, grids, constants), summary_statistics)
+end
+
+function map_summary_statistics(summary_statistics, sim_time, a, Δt, grids::PencilGrids, constants)
     root = 0
-    summary = MPI.Reduce(
-                   output_config.summary_statistics(Dates.now(), t, a, Δt, grids, constants),
-                   pool_summarystat,
-                   root,
-                   grids.MPI_COMM
-                  )
-
-    line = generate_summary_row(summary)
-    if MPI.Comm_rank(grids.MPI_COMM) == 0
-        open(joinpath(output_config.directory, "summary.csv"), "a") do file
-            write(file, line)
-        end
-    end
+
+    local_summaries = map(x -> x(sim_time, a, Δt, grids, constants), summary_statistics)
+
+    global_summaries = map(
+                           x -> MPI.Reduce(x, pool_summarystat, root, grids.MPI_COMM,),
+                           local_summaries
+                          )
+
 end
 
 """
-    SummaryStatistics
+    column_titles(summaries)
 
-Summary statistics including the time, scale factor and current time step.
+Generate column titles from an iterable of summaries
 """
-struct SummaryStatistics
-    "wall time"
-    date::Dates.DateTime
-    "time"
-    t::Float64
-    "scale factor"
-    a::Float64
-    "time step"
-    Δt::Float64
+function column_titles(summaries)
+    mapreduce(column_title, (s1, s2)->s1 * "," * s2, summaries) * "\n"
+end
+
+"""
+    column_title(summary_struct)
+
+Get column title
+
+Assumes that the desired column title is the name of the first field.  This can
+be refined for other structs by defining a specialization for the given
+datatype.
+"""
+function column_title(summary_struct)
+    "$(fieldname(summary_struct, 1))"
 end
 
-function SummaryStatistics(wall_time, sim_time, a, Δt, grids, constants)
-    SummaryStatistics(wall_time, sim_time, a, Δt)
+"""
+    generate_summary_row(summaries)
+"""
+function generate_summary_row(summaries)::String
+    mapreduce(data, (s1, s2)-> "$s1,$s2", summaries) * "\n"
+end
+
+
+"""
+    data(summary_struct)
+
+Get column entry
+
+Assumes that the desired data is in the first field of the struct.  This can
+be refined for other structs by defining a specialization for the given
+datatype.
+"""
+function data(summary_struct)
+    getfield(summary_struct, 1)
 end
 
 """
-    SummaryStatisticsMeanMaxRms
+    WallTime
+
+The current time in the real world.
+
+
+# Examples
+
+`WallTime` created after another contains a later time.
+
+```jldoctest
+julia> using UltraDark
+
+julia> t1 = Summary.WallTime();
 
-Summary statistics including the mean density and RMS density contrast.
+julia> t2 = Summary.WallTime(0., 1., 1e-1, Grids(1.0, 16), nothing);
+
+julia> t1.date < t2.date
+true
+```
 """
-struct SummaryStatisticsMeanMaxRms
+struct WallTime
     "wall time"
     date::Dates.DateTime
+end
+
+function WallTime()
+    WallTime(Dates.now())
+end
+
+function WallTime(sim_time, a, Δt, grids, constants)
+    WallTime()
+end
+
+"""
+    pool_summarystat(S1::WallTime, S2::WallTime)
+
+MPI reduction operator for wall time
+
+Return the time of the first argument.
+
+"""
+function pool_summarystat(S1::WallTime, S2::WallTime)::WallTime
+
+    S1
+
+end
+
+struct SimulationTime
     "time"
     t::Float64
+end
+
+function SimulationTime(sim_time, a, Δt, grids, constants)
+    SimulationTime(sim_time)
+end
+
+"""
+    pool_summarystat(S1::SimulationTime, S2::SimulationTime)
+
+MPI reduction operator for simulation time
+
+Check that times are equal and return
+"""
+function pool_summarystat(S1::SimulationTime, S2::SimulationTime)::SimulationTime
+
+    if (S1.t != S2.t)
+        @error "Summaries incompatible across nodes" S1 S2
+    end
+
+    S1
+end
+
+struct ScaleFactor
     "scale factor"
     a::Float64
+end
+
+function ScaleFactor(sim_time, a, Δt, grids, constants)
+    ScaleFactor(a)
+end
+
+"""
+    pool_summarystat(S1::ScaleFactor, S2::ScaleFactor)
+
+MPI reduction operator for scale factor
+
+Check that scale factors are equal and return
+"""
+function pool_summarystat(S1::ScaleFactor, S2::ScaleFactor)::ScaleFactor
+
+    if (S1.a != S2.a)
+        @error "Summaries incompatible across nodes" S1 S2
+    end
+
+    S1
+end
+
+struct TimeStep
     "time step"
     Δt::Float64
+end
+
+function TimeStep(sim_time, a, Δt, grids, constants)
+    TimeStep(Δt)
+end
+
+"""
+    pool_summarystat(S1::TimeStep, S2::TimeStep)
+
+MPI reduction operator for scale factor
+
+Check that scale factors are equal and return
+"""
+function pool_summarystat(S1::TimeStep, S2::TimeStep)::TimeStep
+
+    if (S1.Δt != S2.Δt)
+        @error "Summaries incompatible across nodes" S1 S2
+    end
+
+    S1
+end
+
+"""
+    MeanDensity
+
+Summary statistic containing the mean density and the number of cells over
+which it was calculated.
+
+# Examples
+
+```jldoctest
+julia> using UltraDark
+
+julia> g = Grids(1.0, 16);
+
+julia> Summary.MeanDensity(g)
+UltraDark.Summary.MeanDensity(0.0, 4096)
+
+```
+"""
+struct MeanDensity
     "mean of density"
     ρx_mean::Float64
-    "max of density"
-    ρx_max::Float64
-    "RMS of density contrast"
-    δx_rms::Float64
     "number of grid points summarized"
     n::Int64
 end
 
-function SummaryStatisticsMeanMaxRms(wall_time, sim_time, a, Δt, grids, constants)
+function MeanDensity(grids)
+    ρx_mean = mean(grids.ρx)
+    n = prod(size(grids.ρx))
+
+    MeanDensity(ρx_mean, n)
+end
+
+function MeanDensity(sim_time, a, Δt, grids, constants)
     ρx_mean = mean(grids.ρx)
-    ρx_max = maximum(grids.ρx)
-    δx_rms = mean(((grids.ρx .- ρx_mean).^2))^0.5
     n = prod(size(grids.ρx))
 
-    SummaryStatisticsMeanMaxRms(wall_time, sim_time, a, Δt, ρx_mean, ρx_max, δx_rms, n)
+    MeanDensity(ρx_mean, n)
 end
 
 """
-    column_titles
+    pool_summarystat(S1::MeanDensity, S2::MeanDensity)
 
-Generate column titles for a CSV file from fields names of a struct
+MPI reduction operator for mean density
 """
-function column_titles(stat_struct)
-    reduce(replace, [":"=>"", "("=>"", ")"=>"", " "=>""], init="$(fieldnames(stat_struct))") * "\n"
+function pool_summarystat(S1::MeanDensity, S2::MeanDensity)::MeanDensity
+
+    n = S1.n + S2.n
+    ρx_mean = (S1.ρx_mean * S1.n + S2.ρx_mean * S2.n) / n
+    ρx_max = maximum(S1.ρx_max, S2.ρx_max)
+
+    MeanDensity(ρx_mean, n)
 end
 
 """
-    pool_summarystat(S1::SummaryStatistics, S2::SummaryStatistics)
+    MaxDensity
+"""
+struct MaxDensity
+    "max of density"
+    ρx_max::Float64
+end
 
-MPI reduction operator for summary statistics.
+function MaxDensity(sim_time, a, Δt, grids, constants)
+    ρx_max = maximum(grids.ρx)
+
+    MaxDensity(ρx_max)
+end
 
-Check that t, a, Δt are equal and return them.
 """
-function pool_summarystat(S1::SummaryStatistics, S2::SummaryStatistics)::SummaryStatistics
+    pool_summarystat(S1::MaxDensity, S2::MaxDensity)
 
-    if (S1.t != S2.t) || (S1.a != S2.a) || (S1.Δt != S2.Δt)
-        @error "Summaries incompatible across nodes" S1 S2
-    end
+MPI reduction operator for max density
+"""
+function pool_summarystat(S1::MaxDensity, S2::MaxDensity)::MaxDensity
 
-    S1
+    n = S1.n + S2.n
+    ρx_max = maximum(S1.ρx_max, S2.ρx_max)
+
+    MaxDensity(ρx_max, n)
 end
 
 """
-    pool_summarystat(S1::SummaryStatisticsMeanMaxRms, S2::SummaryStatisticsMeanMaxRms)
+    RmsDensityContrast
+"""
+struct RmsDensityContrast
+    "RMS of density contrast"
+    δx_rms::Float64
+    "number of grid points summarized"
+    n::Int64
+end
+
+function RmsDensityContrast(sim_time, a, Δt, grids, constants)
+    ρx_mean = mean(grids.ρx)
+    δx_rms = mean(((grids.ρx .- ρx_mean).^2))^0.5
+    n = prod(size(grids.ρx))
+
+    RmsDensityContrast(δx_rms, n)
+end
 
-MPI reduction operator for summary statistics.
 """
-function pool_summarystat(S1::SummaryStatisticsMeanMaxRms, S2::SummaryStatisticsMeanMaxRms)::SummaryStatisticsMeanMaxRms
+    pool_summarystat(S1::RmsDensityContrast, S2::RmsDensityContrast)
 
-    if (S1.t != S2.t) || (S1.a != S2.a) || (S1.Δt != S2.Δt)
-        @error "Summaries incompatible across nodes" S1 S2
-    end
+MPI reduction operator for summary statistics.
+"""
+function pool_summarystat(S1::RmsDensityContrast, S2::RmsDensityContrast)::RmsDensityContrast
 
     n = S1.n + S2.n
-    ρx_mean = (S1.ρx_mean * S1.n + S2.ρx_mean * S2.n) / n
-    ρx_max = maximum(S1.ρx_max, S2.ρx_max)
     δx_rms = ((S1.n * S1.δx_rms^2 + S2.n * S2.δx_rms^2) / n)^0.5
 
-    SummaryStatisticsMeanMaxRms(S1.date, S1.t, S1.a, S1.Δt, ρx_mean, ρx_max, δx_rms, n)
+    RmsDensityContrast(δx_rms, n)
 end
 
 end # module
diff --git a/test/runtests.jl b/test/runtests.jl
index 02b6622..3ef3775 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -13,10 +13,6 @@ end
     include("output.jl")
 end
 
-@safetestset "summary ouputs" begin
-    include("output.jl")
-end
-
 @testset "Phase grad Grids" begin
     include("phase_gradient.jl")
 end
diff --git a/test/summary.jl b/test/summary.jl
deleted file mode 100644
index 7ae1d17..0000000
--- a/test/summary.jl
+++ /dev/null
@@ -1,8 +0,0 @@
-using Test
-using UltraDark: Grids, PencilGrids, output_grids, OutputConfig
-
-for grid_type in [Grids, PencilGrids]
-    for stat in [SummaryStatistics, SummaryStatisticsMeanMaxRms]
-        @test stat(grid_type())
-    end
-end