@testset "DynamicSimulation.jl" begin
    @testset "Constructor" begin
        grid = TUG.Grid{Float64}(5, ones(1, 5))
        boundary = TUG.Boundary{Float64}(grid)
        simulation = TUG.DynamicSimulation{Float64}(grid, boundary; timestep = 0.1)
        @test simulation.grid == grid
        @test simulation.bc == boundary
        @test simulation.approach == BTCS
        @test simulation.iterations == 1
        @test simulation.timestep == 0.1

        grid = TUG.Grid{Float64}(5, ones(1, 5))
        boundary = TUG.Boundary{Float64}(grid)
        simulation =
            TUG.DynamicSimulation{Float64}(grid, boundary; approach = FTCS, timestep = 0.01)
        @test simulation.grid == grid
        @test simulation.bc == boundary
        @test simulation.approach == FTCS
        @test simulation.iterations == 1
        @test simulation.timestep == 0.01

        grid = TUG.Grid{Float64}(5, ones(1, 5))
        boundary = TUG.Boundary{Float64}(grid)
        simulation =
            TUG.DynamicSimulation{Float64}(grid, boundary; approach = FTCS, timestep = 2.33)
        @test simulation.grid == grid
        @test simulation.bc == boundary
        @test simulation.approach == FTCS
        @test simulation.iterations == 10
        @test simulation.timestep == 0.233
    end
    @testset "1D-Run" begin
        grid = TUG.Grid{Float64}(5, ones(1, 5))
        TUG.setConcentrations!(grid, [1.0 1.0 20.0 1.0 1.0])
        boundary = TUG.Boundary{Float64}(grid)
        simulation = TUG.DynamicSimulation{Float64}(grid, boundary; timestep = 0.01)
        TUG.createGrid(simulation)
        for _ = 1:20
            TUG.next(simulation)
        end
        expected_concentrations =
            [1.281106278320615 3.5643693033301567 14.309048836698485 3.5643693033301598 1.281106278320616]
        @test isapprox(
            TUG.getConcentrations(simulation, 1),
            expected_concentrations,
            atol = 1e-6,
        )

        grid = TUG.Grid{Float64}(5, ones(1, 5))
        TUG.setConcentrations!(grid, [1.0 1.0 20.0 1.0 1.0])
        boundary = TUG.Boundary{Float64}(grid)
        TUG.setBoundarySideConstant!(boundary, LEFT, 5.0)
        simulation = TUG.DynamicSimulation{Float64}(grid, boundary; timestep = 0.01)
        TUG.createGrid(simulation)
        for _ = 1:20
            TUG.next(simulation)
        end
        expected_concentrations =
            [2.4416160635284823 3.6810808789967466 14.317333805802393 3.5648326408458035 1.2811288426376255]
        @test isapprox(
            TUG.getConcentrations(simulation, 1),
            expected_concentrations,
            atol = 1e-6,
        )
    end
    @testset "2D-Run" begin
        grid = TUG.Grid{Float64}(5, 5, ones(5, 5), ones(5, 5))
        TUG.setConcentrations!(
            grid,
            [
                1.0 1.0 20.0 1.0 1.0
                1.0 1.0 20.0 1.0 1.0
                1.0 1.0 20.0 1.0 1.0
                1.0 1.0 20.0 1.0 1.0
                1.0 1.0 20.0 1.0 1.0
            ],
        )
        boundary = TUG.Boundary{Float64}(grid)
        simulation = TUG.DynamicSimulation{Float64}(grid, boundary; timestep = 0.01)
        TUG.createGrid(simulation)
        for _ = 1:20
            TUG.next(simulation)
        end
        expected_concentrations = [
            1.141904802011076 3.591390417498421 14.249599956958917 3.5913904174984217 1.1419048020110782
            1.1419048020110738 3.5913904174984173 14.2495999569589 3.5913904174984177 1.1419048020110767
            1.1419048020110725 3.591390417498413 14.249599956958875 3.5913904174984137 1.1419048020110751
            1.1419048020110738 3.5913904174984164 14.249599956958901 3.5913904174984173 1.141904802011077
            1.1419048020110774 3.5913904174984297 14.24959995695894 3.5913904174984297 1.1419048020110796
        ]
        @test isapprox(
            TUG.getConcentrations(simulation, 1),
            expected_concentrations,
            atol = 1e-6,
        )

        grid = TUG.Grid{Float64}(5, 5, ones(5, 5), ones(5, 5))
        TUG.setConcentrations!(
            grid,
            [
                1.0 1.0 20.0 1.0 1.0
                1.0 1.0 20.0 1.0 1.0
                1.0 1.0 20.0 1.0 1.0
                1.0 1.0 20.0 1.0 1.0
                1.0 1.0 20.0 1.0 1.0
            ],
        )
        boundary = TUG.Boundary{Float64}(grid)
        TUG.setBoundarySideConstant!(boundary, LEFT, 5.0)
        simulation = TUG.DynamicSimulation{Float64}(grid, boundary; timestep = 0.01)
        TUG.createGrid(simulation)
        for _ = 1:20
            TUG.next(simulation)
        end
        expected_concentrations = [
            2.116922560959072 3.6843335107065727 14.255652775906785 3.5916901126508205 1.141911092414447
            2.1169225609590687 3.684333510706568 14.255652775906759 3.5916901126508134 1.1419110924144453
            2.1169225609590674 3.6843335107065607 14.255652775906748 3.591690112650811 1.1419110924144442
            2.1169225609590687 3.684333510706565 14.255652775906766 3.5916901126508147 1.1419110924144449
            2.116922560959073 3.684333510706576 14.25565277590681 3.5916901126508267 1.141911092414448
        ]
        @test isapprox(
            TUG.getConcentrations(simulation, 1),
            expected_concentrations,
            atol = 1e-6,
        )
    end
end