work in progress

Project.toml

@@ -15,3 +15,4 @@ Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
 Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Trixi = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb"
+TruncatedStacktraces = "781d530d-4396-4725-bb49-402e4bee1e77"

src/DiscSimulations.jl

@@ -1,29 +1,22 @@
-module DiscSimulationsjl
+module DiscSimulations
 
 using LinearAlgebra
 using Plots
 using Printf
-
-include("solution.jl")
-
-include("burgers/utils.jl")
-#include("burgers/spectral.jl")
-#include("burgers/simple.jl")
-include("simpleGravitySim/simulation.jl")
-include("younsi-2012.jl")
-include("Euler_with_source.jl")
-include("burgers/trixi.jl")
-
-struct Parameters{}
-    N
-    xmin
-    xmax
+using TruncatedStacktraces
+TruncatedStacktraces.VERBOSE[] = true
+
+#struct Parameters{T,DensityFuncType,SourceFuncType}
+struct Parameters{T}
+    N::Int
+    xmin::T
+    xmax::T
     ρ
     source
 end
 
 abstract type SimulationType end
-main(p::Parameters, ::SimulationType) = error("Unknown simulation type.")
+main(p::Parameters, ::SimulationType) = error("Unknown simulation type: $(typeof(t))")
 
 struct BurgersSimulation <: SimulationType end
 struct EulerWithSourceSimulation <: SimulationType end
@@ -33,6 +26,18 @@ main(p::Parameters, ::BurgersSimulation) = BurgerTrixi.solve_disc(p.N, p.xmax)
 main(p::Parameters, ::EulerWithSourceSimulation) = EulerSource.solve_euler(p.source)
 main(p::Parameters, ::SimpleGravitySimulation) = SimpleGravitySimulation.solve_gravity_sim(p.ρ)
 
+include("solution.jl")
+
+include("burgers/utils.jl")
+include("burgers/trixi.jl")
+#include("burgers/spectral.jl")
+#include("burgers/simple.jl")
+include("simpleGravitySim/simpleGravity.jl")
+include("simpleGravitySim/simulation.jl")
+include("younsi-2012.jl")
+include("Euler_with_source.jl")
+
+
 export Parameters, SimulationType, BurgersSimulation, EulerWithSourceSimulation, SimpleGravitySimulation
 
 end # module DiscSimulations

src/Euler_with_source.jl

@@ -20,7 +20,7 @@ using Plots
     du1 = -cos(ϕ)*ρ*(G*mₛ/(r^2))
     du2 = -sin(ϕ)*ρ*(G*mₛ/(r^2))
     du3 = -(ux*cos(ϕ) + uy*sin(ϕ))*ρ*(G*mₛ/(r^2))
-    return SVector(du1, du2, du3)
+    return SVector(0, du1, du3)
 end
 
 source_zero(u, x, t, equations) = SVector(0, 0, 0)

src/burgers/trixi.jl

@@ -5,15 +5,10 @@ using DifferentialEquations
 using Plots
 using Printf
 
-include("utils.jl")
-using .BurgerUtils
+import DiscSimulations:
+    STANDARD_BURGER_INIT, STANDARD_BURGER_TSPAN, DiscSolution, OneDimension
 
-include("../solution.jl")
-
-#import DiscSimulations:
-#    BurgerUtils.STANDARD_BURGER_INIT, BurgerUtils.STANDARD_BURGER_TSPAN, DiscSolution, OneDimension
-
-function setup(N, x_max, init = BurgerUtils.STANDARD_BURGER_INIT, t_span = BurgerUtils.STANDARD_BURGER_TSPAN)
+function setup(N, x_max, init = STANDARD_BURGER_INIT, t_span = STANDARD_BURGER_TSPAN)
     x_min = 0.0
 
     equations = InviscidBurgersEquation1D()
@@ -22,15 +17,15 @@ function setup(N, x_max, init = BurgerUtils.STANDARD_BURGER_INIT, t_span = Burge
     mesh = TreeMesh(x_min, x_max, initial_refinement_level = 3, n_cells_max = N)
     semi = SemidiscretizationHyperbolic(mesh, equations, (x, t, eqs) -> init(x), solver)
 
-    problem = semidiscretize(semi, t_span)
+    ode = semidiscretize(semi, t_span)
 
-    semi, problem
+    semi, ode
 end
 
-function solve_disc(N, x_max, init = BurgerUtils.STANDARD_BURGER_INIT, t_span = BurgerUtils.STANDARD_BURGER_TSPAN)
+function solve_disc(N, x_max, init = STANDARD_BURGER_INIT, t_span = STANDARD_BURGER_TSPAN)
     semi, problem = setup(N, x_max, init, t_span)
     sol = solve(problem, RDPK3SpFSAL49(), abstol = 1.0e-7, reltol = 1.0e-7)
-    return DiscSolution(sol, semi)
+    return DiscSolution(sol, semi, OneDimension())
 end
 
 end # module BurgerTrixi

src/burgers/utils.jl

@@ -1,5 +1,3 @@
-module BurgerUtils
-
 function central_difference(Δu, N)
     lower = -ones(Float64, N - 1)
     mid = zeros(Float64, N)
@@ -11,8 +9,4 @@ function STANDARD_BURGER_INIT(x)
     @. 1 - cos(x)
 end
 
-const STANDARD_BURGER_TSPAN = (0.0, 1.0)
-
-end # module BurgerUtils
-
-export BurgerUtils
+const STANDARD_BURGER_TSPAN = (0.0, 1.0)

src/simpleGravitySim/simulation.jl

@@ -1,7 +1,7 @@
 module simpleGravitySimulation
 
-include("./simpleGravity.jl")
-using .simpleGravity
+import DiscSimulations:
+    simpleGravity
 using Trixi
 using OrdinaryDiffEq
 using Plots

src/solution.jl

@@ -7,33 +7,40 @@ struct OneDimension <: AbstractDimension end
 struct DiscSolution{Dimension,SolutionType,DiscretiziationType}
     sol::SolutionType
     semi::DiscretiziationType
-    #to add: keyword args
 end
 
 function DiscSolution(sol, semi, dimension_type::AbstractDimension)
     DiscSolution{typeof(dimension_type),typeof(sol),typeof(semi)}(sol, semi)
 end
 
 #1D version of plot
-function dataplot(solution::DiscSolution{OneDimension})
-    return PlotData1D(solution.sol(t), solution.semi)
+function dataplot(solution::DiscSolution{<:OneDimension}, solver, t)
+    return PlotData1D(solver(t), solution.semi)
 end
 
 #2D version of plot
-function dataplot(solution::DiscSolution{TwoDimension})
-    return PlotData2D(solution.sol(t), solution.semi)
+function dataplot(solution::DiscSolution{<:TwoDimension}, solver, t)
+    return PlotData2D(solver(t), solution.semi)
 end
 
-function plotgif(solution, tmin, tmax)
+dataplot(solution::DiscSolution{<:AbstractDimension}, t) = error("Unknown simulation type: $(typeof(t))")
+
+function plotgif(solution::DiscSolution, tmin, tmax)#, kwargs...)
     ts = range(tmin, tmax, 150)
+    solver = @time solution.sol
     frames = Plots.@animate for t in ts
-        pd = dataplot(solution)
+        pd = dataplot(solution, solver, t)
         plot(
             pd,
             label = Printf.@sprintf("t = %1.2f", t),
-            ylims = (-0.1, 2.1),
-            legend = :topright,
+            #ylims = (-1, 1.1), 
+            legend = :topright
+            #kwargs
         )
     end
     gif(frames, "temp.gif", fps = 10) |> display
 end
+
+function Base.show(io::IO, sol::DiscSolution)
+    println(io, string(typeof(sol)))
+end

src/test.jl

@@ -1,9 +1,11 @@
 include("DiscSimulations.jl")
-using .DiscSimulationsjl
+using .DiscSimulations
 
 initial_condition_sine(x, t, equation) = sinpi(x[1])
 source_zero(u, x, t, equations) = SVector(0, 0, 0)
-params = DiscSimulationsjl.Parameters(512, 0, 2π, initial_condition_sine, source_zero)
+params = DiscSimulations.Parameters(512, 0.0, 2π, initial_condition_sine, source_zero)
 
-solution = DiscSimulationsjl.main(params, DiscSimulationsjl.BurgersSimulation())
-DiscSimulationsjl.plotgif(solution, 0.0, 2.0)
+solution = DiscSimulations.main(params, DiscSimulations.BurgersSimulation())
+s = DiscSimulations.DiscSolution(solution.sol, solution.semi, DiscSimulations.OneDimension())
+kwargs = ["ylims = (-1, 1.1)", "legend = :topright"]
+DiscSimulations.plotgif(s, 0.0, 2.0)#, kwargs)

test.jl

@@ -0,0 +1,9 @@
+include("src/DiscSimulations.jl")
+using .DiscSimulations
+
+initial_condition_sine(x, t, equation) = sinpi(x[1])
+source_zero(u, x, t, equations) = SVector(0, 0, 0)
+params = DiscSimulations.Parameters(512, 0, 2π, initial_condition_sine, source_zero)
+
+solution = DiscSimulations.main(params, DiscSimulations.BurgersSimulation())
+DiscSimulations.plotgif(solution, 0.0, 2.0)