diff --git a/README.md b/README.md
index 4e921e0..21b431c 100644
--- a/README.md
+++ b/README.md
@@ -46,16 +46,12 @@ sol = simulate(domain1D)
 
 # you can now plot the solutions from the sol variable
 
-# extract characteristic time to convert back to dimensional time
-@unpack tfinal_ad, t_charact = domain1D
-
-anim = @animate for i = LinRange(0, tfinal_ad, 100)
+anim = @animate for i = LinRange(0, sol.t[end], 100)
     l = @layout [a ; b]
 
-    p1 = plot(distance, Fe0, label="Fe initial", linestyle = :dash, linewidth=1, dpi=200, title = "Timestep = $(round(((i)* t_charact);digits=2)) Ma", legend=:outerbottomright, linecolor=1,xlabel = "Distance (µm)")
+    p1 = plot(distance, Fe0, label="Fe initial", linestyle = :dash, linewidth=1, dpi=200, title = "Timestep = $(round(i;digits=2)) Ma", legend=:outerbottomright, linecolor=1,xlabel = "Distance (µm)")
     p1 = plot!(distance, sol(i)[:,2], label="Fe",linecolor=1, linewidth=1)
 
-
     p2 = plot(distance, Mg0, label="Mg initial", linestyle = :dash, linewidth=1, dpi=200,legend=:outerbottomright,linecolor=2,xlabel = "Distance (µm)")
     p2 = plot!(distance, Mn0, label="Mn initial", linestyle = :dash, linewidth=1, linecolor=3)
     p2 = plot!(distance, Ca0, label="Ca initial", linestyle = :dash, linewidth=1, linecolor=4)
@@ -84,6 +80,12 @@ DiffusionGarnet may be installed directly from the REPL:
 ```julia-repl
 julia>]
   pkg> add DiffusionGarnet
+```
+
+And you can test the package with:
+
+```julia-repl
+julia>]
   pkg> test DiffusionGarnet
 ```
 
diff --git a/docs/src/diffusion_1D.md b/docs/src/diffusion_1D.md
index 02efc85..9cb7f00 100644
--- a/docs/src/diffusion_1D.md
+++ b/docs/src/diffusion_1D.md
@@ -74,7 +74,7 @@ which outputs the time spent on the solver, for example, on the second run:
 
 `simulate()` uses the [DifferentialEquations.jl](https://docs.sciml.ai/DiffEqDocs/stable/) package behind the hood to solve our problem efficiently. The returned variable `sol` is the common solution type from this package and more information can be found [here](https://docs.sciml.ai/DiffEqDocs/stable/basics/solution/). It basically holds all the information from our simulation.
 
-We can now plot the solution to our problem.
+We can now plot the solution to our problem. By default, DiffusionGarnet saves the time in Myr.
 
 ```julia
 
diff --git a/test/runtests.jl b/test/runtests.jl
index 5b1dcf4..cff1f56 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -4,351 +4,382 @@ using DelimitedFiles
 using JLD2
 using LinearAlgebra: norm
 
-@testset "initial conditions" begin
-    # 1D, test geometry.jl
-    CMg = ones(5) .* 0.1
-    CFe = ones(5) .* 0.1
-    CMn = ones(5) .* 0.1
-    Lx = 10.0u"µm"
-    tfinal = 1.0u"Myr"
-
-    IC1D = InitialConditions1D(CMg, CFe, CMn, Lx, tfinal)
-
-    @test IC1D.CMg0 == CMg
-    @test IC1D.Lx == ustrip(u"µm",Lx)
-    @test IC1D.tfinal == ustrip(u"Myr",tfinal)
-
-    ICSph = InitialConditionsSpherical(CMg, CFe, CMn, Lx, tfinal)
-
-    @test ICSph.CMg0 == CMg
-    @test ICSph.Lr == ustrip(u"µm",Lx)
-    @test ICSph.tfinal == ustrip(u"Myr",tfinal)
-
-    CMg = ones(5, 5) .* 0.1
-    CFe = ones(5, 5) .* 0.1
-    CMn = ones(5, 5) .* 0.1
-    Ly = 10.0u"µm"
-
-    IC2D = InitialConditions2D(CMg, CFe, CMn, Lx, Ly, tfinal)
-
-    @test IC2D.CMg0 == CMg
-    @test IC2D.Ly == ustrip(u"µm",Ly)
-    @test IC2D.tfinal == ustrip(u"Myr",tfinal)
-
-    CMg = ones(5, 5, 5) .* 0.1
-    CFe = ones(5, 5, 5) .* 0.1
-    CMn = ones(5, 5, 5) .* 0.1
-    Lz = 10.0u"µm"
-
-    IC3D = InitialConditions3D(CMg, CFe, CMn, Lx, Ly, Lz, tfinal)
-
-    @test IC3D.CMg0 == CMg
-    @test IC3D.Lz == ustrip(u"µm",Lz)
-    @test IC3D.tfinal == ustrip(u"Myr",tfinal)
-
-    T = 650u"°C"
-    P = 2u"GPa"
-    domain1D = Domain(IC1D, T, P)
-    @test domain1D.L_charact == ustrip(u"µm",Lx)
-    @test domain1D.t_charact ≈ 0.22518558662307234
-    @test domain1D.tfinal_ad ≈ 4.44078155709785
-    @test domain1D.Δxad_ == 4.0
-
-    domainSph = Domain(ICSph, T, P)
-    @test domainSph.L_charact == ustrip(u"µm",Lx)
-    @test domainSph.t_charact ≈ 0.22518558662307234
-    @test domainSph.tfinal_ad ≈ 4.44078155709785
-    @test domainSph.Δrad_ == 4.0
-    @test domainSph.r_ad[end] == 1.0
-
-    domain2D = Domain(IC2D, T, P)
-    @test domain2D.L_charact == ustrip(u"µm",Lx)
-    @test domain2D.Δyad_ == 4.0
-
-    domain3D = Domain(IC3D, T, P)
-    @test domain3D.L_charact == ustrip(u"µm",Lx)
-    @test domain3D.Δzad_ == 4.0
-end
+@testset "All tests" begin
+
+    @testset "initial conditions" begin
+        # 1D, test geometry.jl
+        CMg = ones(5) .* 0.1
+        CFe = ones(5) .* 0.1
+        CMn = ones(5) .* 0.1
+        Lx = 10.0u"µm"
+        tfinal = 1.0u"Myr"
+
+        IC1D = InitialConditions1D(CMg, CFe, CMn, Lx, tfinal)
+
+        @test IC1D.CMg0 == CMg
+        @test IC1D.Lx == ustrip(u"µm",Lx)
+        @test IC1D.tfinal == ustrip(u"Myr",tfinal)
+
+        ICSph = InitialConditionsSpherical(CMg, CFe, CMn, Lx, tfinal)
+
+        @test ICSph.CMg0 == CMg
+        @test ICSph.Lr == ustrip(u"µm",Lx)
+        @test ICSph.tfinal == ustrip(u"Myr",tfinal)
+
+        CMg = ones(5, 5) .* 0.1
+        CFe = ones(5, 5) .* 0.1
+        CMn = ones(5, 5) .* 0.1
+        Ly = 10.0u"µm"
+
+        IC2D = InitialConditions2D(CMg, CFe, CMn, Lx, Ly, tfinal)
+
+        @test IC2D.CMg0 == CMg
+        @test IC2D.Ly == ustrip(u"µm",Ly)
+        @test IC2D.tfinal == ustrip(u"Myr",tfinal)
+
+        CMg = ones(5, 5, 5) .* 0.1
+        CFe = ones(5, 5, 5) .* 0.1
+        CMn = ones(5, 5, 5) .* 0.1
+        Lz = 10.0u"µm"
+
+        IC3D = InitialConditions3D(CMg, CFe, CMn, Lx, Ly, Lz, tfinal)
+
+        @test IC3D.CMg0 == CMg
+        @test IC3D.Lz == ustrip(u"µm",Lz)
+        @test IC3D.tfinal == ustrip(u"Myr",tfinal)
+
+        T = 650u"°C"
+        P = 2u"GPa"
+        domain1D = Domain(IC1D, T, P)
+        @test domain1D.L_charact == ustrip(u"µm",Lx)
+        @test domain1D.t_charact ≈ 0.22518558662307234
+        @test domain1D.tfinal_ad ≈ 4.44078155709785
+        @test domain1D.Δxad_ == 4.0
+
+        domainSph = Domain(ICSph, T, P)
+        @test domainSph.L_charact == ustrip(u"µm",Lx)
+        @test domainSph.t_charact ≈ 0.22518558662307234
+        @test domainSph.tfinal_ad ≈ 4.44078155709785
+        @test domainSph.Δrad_ == 4.0
+        @test domainSph.r_ad[end] == 1.0
+
+        domain2D = Domain(IC2D, T, P)
+        @test domain2D.L_charact == ustrip(u"µm",Lx)
+        @test domain2D.Δyad_ == 4.0
+
+        domain3D = Domain(IC3D, T, P)
+        @test domain3D.L_charact == ustrip(u"µm",Lx)
+        @test domain3D.Δzad_ == 4.0
+    end
 
-@testset "diffusion coefficients" begin
-    # test Domain
-    T=650  # in °C
-    P=2  # in kbar
-    D0 = zeros(Float64, 4)
-    D_ini!(D0,T,P)
-
-    # benchmarked with TeriaG and 12.783356187311105
-    @test D0[1] ≈ 2.383676419323230e2
-    @test D0[2] ≈ 4.500484945403809e+02
-    @test D0[3] ≈ 6.232092221232668e+03
-    @test D0[4] ≈ 2.250242472701905e+02
-end
+    @testset "diffusion coefficients" begin
+        # test Domain
+        T=650  # in °C
+        P=2  # in kbar
+        D0 = zeros(Float64, 4)
+        D_ini!(D0,T,P)
+
+        # benchmarked with TeriaG and 12.783356187311105
+        @test D0[1] ≈ 2.383676419323230e2
+        @test D0[2] ≈ 4.500484945403809e+02
+        @test D0[3] ≈ 6.232092221232668e+03
+        @test D0[4] ≈ 2.250242472701905e+02
+    end
 
-@testset "1D diffusion" begin
+    @testset "1D diffusion" begin
 
-    data = DelimitedFiles.readdlm("./Data/1D/Data_Grt_1D.txt", '\t', '\n', header=true)[1]
+        path_1D = joinpath("Data", "1D", "Data_Grt_1D.txt")
 
-    Mg0 = data[:, 4]
-    Fe0 = data[:, 2]
-    Mn0 = data[:, 3]
-    Ca0 = data[:, 5]
-    distance = data[:, 1]
-    Lx = (data[end,1] - data[1,1])u"µm"
-    tfinal = 15u"Myr"
+        data = DelimitedFiles.readdlm(path_1D, '\t', '\n', header=true)[1]
 
-    IC1D = InitialConditions1D(Mg0, Fe0, Mn0, Lx, tfinal)
+        Mg0 = data[:, 4]
+        Fe0 = data[:, 2]
+        Mn0 = data[:, 3]
+        Ca0 = data[:, 5]
+        distance = data[:, 1]
+        Lx = (data[end,1] - data[1,1])u"µm"
+        tfinal = 15u"Myr"
 
-    T = 900u"°C"
-    P = 0.6u"GPa"
+        IC1D = InitialConditions1D(Mg0, Fe0, Mn0, Lx, tfinal)
 
-    domain1D = Domain(IC1D, T, P)
+        T = 900u"°C"
+        P = 0.6u"GPa"
 
-    sol = simulate(domain1D; progress=false, abstol=1e-6,reltol=1e-6)
+        domain1D = Domain(IC1D, T, P)
 
-    @test norm(sum.(sol.u[end][:,1] .+ sol.u[end][:,2] .+ sol.u[end][:,3])) ≈ 28.64886878627501
-end
+        sol = simulate(domain1D; progress=false, abstol=1e-6,reltol=1e-6)
 
-@testset "Spherical diffusion" begin
+        @test norm(sum.(sol.u[end][:,1] .+ sol.u[end][:,2] .+ sol.u[end][:,3])) ≈ 28.64886878627501
+    end
 
-    using LinearAlgebra: norm
+    @testset "Spherical diffusion" begin
 
+        using LinearAlgebra: norm
 
-    data = DelimitedFiles.readdlm("./Data/1D/Data_Grt_1D.txt", '\t', '\n', header=true)[1]
+        path_1D = joinpath("Data", "1D", "Data_Grt_1D.txt")
 
-    Mg0 = reverse(data[1:size(data,1)÷2, 4])
-    Fe0 = reverse(data[1:size(data,1)÷2, 2])
-    Mn0 = reverse(data[1:size(data,1)÷2, 3])
-    Ca0 = reverse(data[1:size(data,1)÷2, 5])
-    distance = data[1:size(data,1)÷2, 1]
-    Lr = (data[end,1] - data[1,1])u"µm"
-    tfinal = 15u"Myr"
+        data = DelimitedFiles.readdlm(path_1D, '\t', '\n', header=true)[1]
 
-    ICSph = InitialConditionsSpherical(Mg0, Fe0, Mn0, Lr, tfinal)
+        Mg0 = reverse(data[1:size(data,1)÷2, 4])
+        Fe0 = reverse(data[1:size(data,1)÷2, 2])
+        Mn0 = reverse(data[1:size(data,1)÷2, 3])
+        Ca0 = reverse(data[1:size(data,1)÷2, 5])
+        distance = data[1:size(data,1)÷2, 1]
+        Lr = (data[end,1] - data[1,1])u"µm"
+        tfinal = 15u"Myr"
 
-    T = 900u"°C"
-    P = 0.6u"GPa"
+        ICSph = InitialConditionsSpherical(Mg0, Fe0, Mn0, Lr, tfinal)
 
-    domainSph = Domain(ICSph, T, P)
+        T = 900u"°C"
+        P = 0.6u"GPa"
 
-    sol = simulate(domainSph; progress=false, abstol=1e-6,reltol=1e-6)
+        domainSph = Domain(ICSph, T, P)
 
-    @test norm(sum.(sol.u[end][:,1] .+ sol.u[end][:,2] .+ sol.u[end][:,3])) ≈ 20.268803083443927
-end
+        sol = simulate(domainSph; progress=false, abstol=1e-6,reltol=1e-6)
 
-@testset "2D Diffusion" begin
+        @test norm(sum.(sol.u[end][:,1] .+ sol.u[end][:,2] .+ sol.u[end][:,3])) ≈ 20.268803083443927
+    end
 
-    using LinearAlgebra: norm
+    @testset "2D Diffusion" begin
 
-    CMg = DelimitedFiles.readdlm("./Data/2D/Xprp_LR.txt", '\t', '\n', header=false)
-    CFe = DelimitedFiles.readdlm("./Data/2D/Xalm_LR.txt", '\t', '\n', header=false)
-    CMn = DelimitedFiles.readdlm("./Data/2D/Xsps_LR.txt", '\t', '\n', header=false)
-    grt_boundary = DelimitedFiles.readdlm("./Data/2D/Contour_LR.txt", '\t', '\n', header=false)
+        using LinearAlgebra: norm
 
-    Lx = 900.0u"µm"
-    Ly = 900.0u"µm"
-    tfinal = 1.0u"Myr"
-    T = 900u"°C"
-    P = 0.6u"GPa"
+        path_2D_Mg = joinpath("Data", "2D", "Xprp_LR.txt")
+        path_2D_Fe = joinpath("Data", "2D", "Xalm_LR.txt")
+        path_2D_Mn = joinpath("Data", "2D", "Xsps_LR.txt")
+        path_2D_grt = joinpath("Data", "2D", "Contour_LR.txt")
 
-    IC2D = InitialConditions2D(CMg, CFe, CMn, Lx, Ly, tfinal; grt_boundary = grt_boundary)
-    domain2D = Domain(IC2D, T, P)
+        CMg = DelimitedFiles.readdlm(path_2D_Mg, '\t', '\n', header=false)
+        CFe = DelimitedFiles.readdlm(path_2D_Fe, '\t', '\n', header=false)
+        CMn = DelimitedFiles.readdlm(path_2D_Mn, '\t', '\n', header=false)
+        grt_boundary = DelimitedFiles.readdlm(path_2D_grt, '\t', '\n', header=false)
 
-    sol = simulate(domain2D; save_everystep=false, save_start=false)
+        Lx = 900.0u"µm"
+        Ly = 900.0u"µm"
+        tfinal = 1.0u"Myr"
+        T = 900u"°C"
+        P = 0.6u"GPa"
 
-    @test norm(sol.u[end][:,:,1]) ≈ 12.783356187311105
-end
+        IC2D = InitialConditions2D(CMg, CFe, CMn, Lx, Ly, tfinal; grt_boundary = grt_boundary)
+        domain2D = Domain(IC2D, T, P)
 
-@testset "3D Diffusion" begin
+        sol = simulate(domain2D; save_everystep=false, save_start=false)
 
-    # use JLD2 to load data
-    file = jldopen("./Data/3D/3D_data.jld2", "r")
-    @unpack Mg0, Fe0, Mn0, Ca0, grt_boundary = file
-    close(file)
+        @test norm(sol.u[end][:,:,1]) ≈ 12.783356187311105
+    end
 
-    Lx = 9000.0u"µm"
-    Ly = 9000.0u"µm"
-    Lz = 9000.0u"µm"
-    tfinal = 0.2u"Myr"
-    T = 900u"°C"
-    P = 0.6u"GPa"
+    @testset "3D Diffusion" begin
 
-    IC3D = InitialConditions3D(Mg0, Fe0, Mn0, Lx, Ly, Lz, tfinal; grt_boundary = grt_boundary)
-    domain3D = Domain(IC3D, T, P)
+        # use JLD2 to load data
+        path_3D = joinpath("Data", "3D", "3D_data.jld2")
 
-    sol = simulate(domain3D; save_everystep=false, save_start=false);
+        file = jldopen(path_3D, "r")
+        @unpack Mg0, Fe0, Mn0, Ca0, grt_boundary = file
+        close(file)
 
-    @test norm(sol.u[end][:,:,:,1]) ≈ 371.1466260290486
-end
+        Lx = 9000.0u"µm"
+        Ly = 9000.0u"µm"
+        Lz = 9000.0u"µm"
+        tfinal = 0.2u"Myr"
+        T = 900u"°C"
+        P = 0.6u"GPa"
 
+        IC3D = InitialConditions3D(Mg0, Fe0, Mn0, Lx, Ly, Lz, tfinal; grt_boundary = grt_boundary)
+        domain3D = Domain(IC3D, T, P)
 
-@testset "Callback update D0" begin
+        sol = simulate(domain3D; save_everystep=false, save_start=false);
 
-    data = DelimitedFiles.readdlm("./Data/1D/Data_Grt_1D.txt", '\t', '\n', header=true)[1]
+        @test norm(sol.u[end][:,:,:,1]) ≈ 371.1466260290486
+    end
 
-    Mg0 = reverse(data[1:size(data,1)÷2, 4])
-    Fe0 = reverse(data[1:size(data,1)÷2, 2])
-    Mn0 = reverse(data[1:size(data,1)÷2, 3])
-    Ca0 = reverse(data[1:size(data,1)÷2, 5])
-    distance = data[1:size(data,1)÷2, 1]
-    Lx = Lr = (data[end,1] - data[1,1])u"µm"
-    tfinal = 1u"Myr"
-    T = [850u"°C", 600u"°C"]
-    P = [0.5u"GPa", 0.3u"GPa"]
 
-    ICSph = InitialConditionsSpherical(Mg0, Fe0, Mn0, Lr, tfinal)
-    IC1D = InitialConditions1D(Mg0, Fe0, Mn0, Lx, tfinal)
+    @testset "Callback update D0" begin
 
-    time_update = [0u"Myr", 1u"Myr"]
+        path_1D = joinpath("Data", "1D", "Data_Grt_1D.txt")
 
-    domainSph = Domain(ICSph, T, P, time_update)
-    domain1D = Domain(IC1D, T, P, time_update)
+        data = DelimitedFiles.readdlm(path_1D, '\t', '\n', header=true)[1]
 
-    @test domainSph.D0[1] ≈ 151880.41527919917
-    @test domain1D.D0[1] ≈ 151880.41527919917
+        Mg0 = reverse(data[1:size(data,1)÷2, 4])
+        Fe0 = reverse(data[1:size(data,1)÷2, 2])
+        Mn0 = reverse(data[1:size(data,1)÷2, 3])
+        Ca0 = reverse(data[1:size(data,1)÷2, 5])
+        distance = data[1:size(data,1)÷2, 1]
+        Lx = Lr = (data[end,1] - data[1,1])u"µm"
+        tfinal = 1u"Myr"
+        T = [850u"°C", 600u"°C"]
+        P = [0.5u"GPa", 0.3u"GPa"]
 
-    CMg = DelimitedFiles.readdlm("./Data/2D/Xprp_LR.txt", '\t', '\n', header=false)
-    CFe = DelimitedFiles.readdlm("./Data/2D/Xalm_LR.txt", '\t', '\n', header=false)
-    CMn = DelimitedFiles.readdlm("./Data/2D/Xsps_LR.txt", '\t', '\n', header=false)
-    grt_boundary = DelimitedFiles.readdlm("./Data/2D/Contour_LR.txt", '\t', '\n', header=false)
+        ICSph = InitialConditionsSpherical(Mg0, Fe0, Mn0, Lr, tfinal)
+        IC1D = InitialConditions1D(Mg0, Fe0, Mn0, Lx, tfinal)
 
-    Lx = 900.0u"µm"
-    Ly = 900.0u"µm"
+        time_update = [0u"Myr", 1u"Myr"]
 
-    IC2D = InitialConditions2D(CMg, CFe, CMn, Lx, Ly, tfinal; grt_boundary = grt_boundary)
-    domain2D = Domain(IC2D, T, P, time_update)
+        domainSph = Domain(ICSph, T, P, time_update)
+        domain1D = Domain(IC1D, T, P, time_update)
 
-    @test domain2D.D0[1] ≈ 151880.41527919917
+        @test domainSph.D0[1] ≈ 151880.41527919917
+        @test domain1D.D0[1] ≈ 151880.41527919917
 
-    @unpack time_update_ad = domainSph
+        path_2D_Mg = joinpath("Data", "2D", "Xprp_LR.txt")
+        path_2D_Fe = joinpath("Data", "2D", "Xalm_LR.txt")
+        path_2D_Mn = joinpath("Data", "2D", "Xsps_LR.txt")
+        path_2D_grt = joinpath("Data", "2D", "Contour_LR.txt")
 
-    update_diffusion_coef_call = PresetTimeCallback(time_update_ad, update_diffusion_coef)
+        CMg = DelimitedFiles.readdlm(path_2D_Mg, '\t', '\n', header=false)
+        CFe = DelimitedFiles.readdlm(path_2D_Fe, '\t', '\n', header=false)
+        CMn = DelimitedFiles.readdlm(path_2D_Mn, '\t', '\n', header=false)
+        grt_boundary = DelimitedFiles.readdlm(path_2D_grt, '\t', '\n', header=false)
 
-    sol_sph = simulate(domainSph; callback=update_diffusion_coef_call, progress=false, abstol=1e-6,reltol=1e-6)
-    sol_1D = simulate(domain1D; callback=update_diffusion_coef_call, progress=false, abstol=1e-6,reltol=1e-6)
+        Lx = 900.0u"µm"
+        Ly = 900.0u"µm"
 
-    @unpack time_update_ad = domain2D
-    update_diffusion_coef_call = PresetTimeCallback(time_update_ad, update_diffusion_coef)
+        IC2D = InitialConditions2D(CMg, CFe, CMn, Lx, Ly, tfinal; grt_boundary = grt_boundary)
+        domain2D = Domain(IC2D, T, P, time_update)
 
-    sol_2D = simulate(domain2D; callback=update_diffusion_coef_call,save_everystep=false, save_start=false)
+        @test domain2D.D0[1] ≈ 151880.41527919917
 
-    T=600  # in °C
-    P=3  # in kbar
-    D0 = zeros(Float64, 4)
-    D_ini!(D0,T,P)
+        @unpack time_update_ad = domainSph
 
-    @test sol_1D.prob.p.domain.D0[1] ≈ D0[1]
-    @test sol_sph.prob.p.domain.D0[1] ≈ D0[1]
-    @test sol_2D.prob.p.domain.D0[1] ≈ D0[1]
-end
+        update_diffusion_coef_call = PresetTimeCallback(time_update_ad, update_diffusion_coef)
 
-@testset "Callback output" begin
+        sol_sph = simulate(domainSph; callback=update_diffusion_coef_call, progress=false, abstol=1e-6,reltol=1e-6)
 
-    using HDF5
+        update_diffusion_coef_call = PresetTimeCallback(time_update_ad, update_diffusion_coef)
 
-    data_1D = DelimitedFiles.readdlm("./Data/1D/Data_Grt_1D.txt", '\t', '\n', header=true)[1]
+        sol_1D = simulate(domain1D; callback=update_diffusion_coef_call, progress=false, abstol=1e-6,reltol=1e-6)
 
-    Mg0 = data_1D[:, 4]
-    Fe0 = data_1D[:, 2]
-    Mn0 = data_1D[:, 3]
-    Ca0 = data_1D[:, 5]
-    distance = data_1D[:, 1]
-    Lx = (data_1D[end,1] - data_1D[1,1])u"µm"
-    tfinal = 15u"Myr"
+        @unpack time_update_ad = domain2D
+        update_diffusion_coef_call = PresetTimeCallback(time_update_ad, update_diffusion_coef)
 
-    IC1D = InitialConditions1D(Mg0, Fe0, Mn0, Lx, tfinal)
-    ICSph = InitialConditionsSpherical(Mg0, Fe0, Mn0, Lx, tfinal)
+        sol_2D = simulate(domain2D; callback=update_diffusion_coef_call,save_everystep=false, save_start=false)
 
-    T = 700u"°C"
-    P = 0.6u"GPa"
+        T=600  # in °C
+        P=3  # in kbar
+        D0 = zeros(Float64, 4)
+        D_ini!(D0,T,P)
 
-    domainSph = Domain(ICSph, T, P)
-    domain1D = Domain(IC1D, T, P)
+        @test sol_1D.prob.p.domain.D0[1] ≈ D0[1]
+        @test sol_sph.prob.p.domain.D0[1] ≈ D0[1]
+        @test sol_2D.prob.p.domain.D0[1] ≈ D0[1]
+    end
 
-    Mg_2D = DelimitedFiles.readdlm("./Data/2D/Xprp_LR.txt", '\t', '\n', header=false)
-    Fe_2D = DelimitedFiles.readdlm("./Data/2D/Xalm_LR.txt", '\t', '\n', header=false)
-    Mn_2D = DelimitedFiles.readdlm("./Data/2D/Xsps_LR.txt", '\t', '\n', header=false)
-    Lx = 900.0u"µm"
-    Ly = 900.0u"µm"
+    @testset "Callback output" begin
 
-    IC2D = InitialConditions2D(Mg_2D, Fe_2D, Mn_2D, Lx, Ly, tfinal)
-    domain2D = Domain(IC2D, T, P)
+        using HDF5
 
-    time_save = [0u"Myr", 2u"Myr", 5u"Myr", 15u"Myr"]
+        path_1D = joinpath("Data", "1D", "Data_Grt_1D.txt")
 
-    save_data_callback = PresetTimeCallback(ustrip.(time_save) ./ domain1D.t_charact, save_data)
+        data_1D = DelimitedFiles.readdlm(path_1D, '\t', '\n', header=true)[1]
 
-    sol_1D = simulate(domain1D; callback=save_data_callback, path_save=(@__DIR__) * "/Grt_1D.h5", abstol=1e-6,reltol=1e-6)
+        Mg0 = data_1D[:, 4]
+        Fe0 = data_1D[:, 2]
+        Mn0 = data_1D[:, 3]
+        Ca0 = data_1D[:, 5]
+        distance = data_1D[:, 1]
+        Lx = (data_1D[end,1] - data_1D[1,1])u"µm"
+        tfinal = 15u"Myr"
 
-    sol_sph = simulate(domainSph; callback=save_data_callback, path_save=(@__DIR__) * "/Grt_Sph.h5", abstol=1e-6,reltol=1e-6)
+        IC1D = InitialConditions1D(Mg0, Fe0, Mn0, Lx, tfinal)
+        ICSph = InitialConditionsSpherical(Mg0, Fe0, Mn0, Lx, tfinal)
 
-    save_data_callback = PresetTimeCallback(ustrip.(time_save) ./ domain2D.t_charact, save_data)
+        T = 700u"°C"
+        P = 0.6u"GPa"
 
-    sol_2D = simulate(domain2D; callback=save_data_callback, path_save=(@__DIR__) * "/Grt_2D.h5", progress=false, save_everystep=false)
+        domainSph = Domain(ICSph, T, P)
+        domain1D = Domain(IC1D, T, P)
 
-    h5open("./Grt_1D.h5", "r") do file
-        @test read(file["Diffusion_Grt"]["t0000"]["Mg"]["Mg"]) == convert(Array{Float32}, IC1D.CMg0)
-        @test read(file["Diffusion_Grt"]["t0000"]["Fe"]["Fe"]) == convert(Array{Float32}, IC1D.CFe0)
-        @test read(file["Diffusion_Grt"]["t0000"]["Mn"]["Mn"]) == convert(Array{Float32}, IC1D.CMn0)
-        @test read(file["Diffusion_Grt"]["t0000"]["Ca"]["Ca"]) == convert(Array{Float32}, replace!((1 .- IC1D.CMg0 .- IC1D.CFe0 .- IC1D.CMn0), 1=>0))
-        @test read(file["Diffusion_Grt"]["t0003"]["Mg"]["Mg"]) == convert(Array{Float32}, sol_1D.u[end][:,1])
-        @test read(file["Diffusion_Grt"]["t0003"]["Fe"]["Fe"]) == convert(Array{Float32}, sol_1D.u[end][:,2])
-        @test read(file["Diffusion_Grt"]["t0003"]["Mn"]["Mn"]) == convert(Array{Float32}, sol_1D.u[end][:,3])
-    end
+        path_2D_Mg = joinpath("Data", "2D", "Xprp_LR.txt")
+        path_2D_Fe = joinpath("Data", "2D", "Xalm_LR.txt")
+        path_2D_Mn = joinpath("Data", "2D", "Xsps_LR.txt")
 
-    h5open("./Grt_Sph.h5", "r") do file
-        @test read(file["Diffusion_Grt"]["t0000"]["Mg"]["Mg"]) == convert(Array{Float32}, ICSph.CMg0)
-        @test read(file["Diffusion_Grt"]["t0000"]["Fe"]["Fe"]) == convert(Array{Float32}, ICSph.CFe0)
-        @test read(file["Diffusion_Grt"]["t0000"]["Mn"]["Mn"]) == convert(Array{Float32}, ICSph.CMn0)
-        @test read(file["Diffusion_Grt"]["t0000"]["Ca"]["Ca"]) == convert(Array{Float32}, replace!((1 .- ICSph.CMg0 .- ICSph.CFe0 .- ICSph.CMn0), 1=>0))
-        @test read(file["Diffusion_Grt"]["t0003"]["Mg"]["Mg"]) == convert(Array{Float32}, sol_sph.u[end][:,1])
-        @test read(file["Diffusion_Grt"]["t0003"]["Fe"]["Fe"]) == convert(Array{Float32}, sol_sph.u[end][:,2])
-        @test read(file["Diffusion_Grt"]["t0003"]["Mn"]["Mn"]) == convert(Array{Float32}, sol_sph.u[end][:,3])
-    end
+        Mg_2D = DelimitedFiles.readdlm(path_2D_Mg, '\t', '\n', header=false)
+        Fe_2D = DelimitedFiles.readdlm(path_2D_Fe, '\t', '\n', header=false)
+        Mn_2D = DelimitedFiles.readdlm(path_2D_Mn, '\t', '\n', header=false)
+        Lx = 900.0u"µm"
+        Ly = 900.0u"µm"
 
-    h5open("./Grt_2D.h5", "r") do file
-        @test read(file["Diffusion_Grt"]["t0000"]["Mg"]["Mg"]) == convert(Array{Float32}, IC2D.CMg0)
-        @test read(file["Diffusion_Grt"]["t0000"]["Fe"]["Fe"]) == convert(Array{Float32}, IC2D.CFe0)
-        @test read(file["Diffusion_Grt"]["t0000"]["Mn"]["Mn"]) == convert(Array{Float32}, IC2D.CMn0)
-        @test read(file["Diffusion_Grt"]["t0000"]["Ca"]["Ca"]) == convert(Array{Float32}, replace!((1 .- IC2D.CMg0 .- IC2D.CFe0 .- IC2D.CMn0), 1=>0))
-        @test read(file["Diffusion_Grt"]["t0003"]["Mg"]["Mg"]) == convert(Array{Float32}, sol_2D.u[end][:,:,1])
-        @test read(file["Diffusion_Grt"]["t0003"]["Fe"]["Fe"]) == convert(Array{Float32}, sol_2D.u[end][:,:,2])
-        @test read(file["Diffusion_Grt"]["t0003"]["Mn"]["Mn"]) == convert(Array{Float32}, sol_2D.u[end][:,:,3])
-    end
+        IC2D = InitialConditions2D(Mg_2D, Fe_2D, Mn_2D, Lx, Ly, tfinal)
+        domain2D = Domain(IC2D, T, P)
 
-    # delete files "Grt_1D.h5" and "Grt_Sph.h5" if it exists
-    if isfile("Grt_1D.h5")
-        rm("Grt_1D.h5")
-    end
+        time_save = [0u"Myr", 2u"Myr", 5u"Myr", 15u"Myr"]
 
-    if isfile("Grt_Sph.h5")
-        rm("Grt_Sph.h5")
-    end
+        save_data_callback = PresetTimeCallback(ustrip.(time_save) ./ domain1D.t_charact, save_data)
 
-    if isfile("Grt_2D.h5")
-        rm("Grt_2D.h5")
-    end
+        sol_1D = simulate(domain1D; callback=save_data_callback, path_save="Grt_1D.h5", abstol=1e-6,reltol=1e-6)
 
-    save_data_callback = PresetTimeCallback(ustrip.(time_save) ./ domain2D.t_charact, save_data_paraview)
+        save_data_callback = PresetTimeCallback(ustrip.(time_save) ./ domain1D.t_charact, save_data)
 
-    sol_2D = simulate(domain2D; callback=save_data_callback, path_save=(@__DIR__) * "/Grt_2D.h5", progress=false, save_everystep=false, save_start=false)
+        sol_sph = simulate(domainSph; callback=save_data_callback, path_save="Grt_Sph.h5", abstol=1e-6,reltol=1e-6)
 
-    h5open("./Grt_2D.h5", "r") do file
-        @test read(file["Diffusion_Grt"]["t0000"]["Mg"]["Mg"])' == convert(Array{Float32}, IC2D.CMg0)
-        @test read(file["Diffusion_Grt"]["t0000"]["Fe"]["Fe"])' == convert(Array{Float32}, IC2D.CFe0)
-        @test read(file["Diffusion_Grt"]["t0000"]["Mn"]["Mn"])' == convert(Array{Float32}, IC2D.CMn0)
-        @test read(file["Diffusion_Grt"]["t0000"]["Ca"]["Ca"])' == convert(Array{Float32},replace!((1 .- IC2D.CMg0 .- IC2D.CFe0 .- IC2D.CMn0), 1=>0))
-        @test read(file["Diffusion_Grt"]["t0003"]["Mg"]["Mg"])' == convert(Array{Float32}, sol_2D.u[end][:,:,1])
-        @test read(file["Diffusion_Grt"]["t0003"]["Fe"]["Fe"])' == convert(Array{Float32}, sol_2D.u[end][:,:,2])
-        @test read(file["Diffusion_Grt"]["t0003"]["Mn"]["Mn"])' == convert(Array{Float32}, sol_2D.u[end][:,:,3])
-    end
+        save_data_callback = PresetTimeCallback(ustrip.(time_save) ./ domain2D.t_charact, save_data)
 
-    # delete files if it exists
-    if isfile("Grt_2D.h5")
-        rm("Grt_2D.h5")
-    end
+        sol_2D = simulate(domain2D; callback=save_data_callback, path_save="Grt_2D.h5", progress=false, save_everystep=false)
 
-    if isfile("Grt_2D.xdmf")
-        rm("Grt_2D.xdmf")
-    end
+        h5open("Grt_1D.h5", "r") do file
+            @test read(file["Diffusion_Grt"]["t0000"]["Mg"]["Mg"]) == convert(Array{Float32}, IC1D.CMg0)
+            @test read(file["Diffusion_Grt"]["t0000"]["Fe"]["Fe"]) == convert(Array{Float32}, IC1D.CFe0)
+            @test read(file["Diffusion_Grt"]["t0000"]["Mn"]["Mn"]) == convert(Array{Float32}, IC1D.CMn0)
+            @test read(file["Diffusion_Grt"]["t0000"]["Ca"]["Ca"]) == convert(Array{Float32}, replace!((1 .- IC1D.CMg0 .- IC1D.CFe0 .- IC1D.CMn0), 1=>0))
+            @test read(file["Diffusion_Grt"]["t0003"]["Mg"]["Mg"]) == convert(Array{Float32}, sol_1D.u[end][:,1])
+            @test read(file["Diffusion_Grt"]["t0003"]["Fe"]["Fe"]) == convert(Array{Float32}, sol_1D.u[end][:,2])
+            @test read(file["Diffusion_Grt"]["t0003"]["Mn"]["Mn"]) == convert(Array{Float32}, sol_1D.u[end][:,3])
+        end
+
+        h5open("Grt_Sph.h5", "r") do file
+            @test read(file["Diffusion_Grt"]["t0000"]["Mg"]["Mg"]) == convert(Array{Float32}, ICSph.CMg0)
+            @test read(file["Diffusion_Grt"]["t0000"]["Fe"]["Fe"]) == convert(Array{Float32}, ICSph.CFe0)
+            @test read(file["Diffusion_Grt"]["t0000"]["Mn"]["Mn"]) == convert(Array{Float32}, ICSph.CMn0)
+            @test read(file["Diffusion_Grt"]["t0000"]["Ca"]["Ca"]) == convert(Array{Float32}, replace!((1 .- ICSph.CMg0 .- ICSph.CFe0 .- ICSph.CMn0), 1=>0))
+            @test read(file["Diffusion_Grt"]["t0003"]["Mg"]["Mg"]) == convert(Array{Float32}, sol_sph.u[end][:,1])
+            @test read(file["Diffusion_Grt"]["t0003"]["Fe"]["Fe"]) == convert(Array{Float32}, sol_sph.u[end][:,2])
+            @test read(file["Diffusion_Grt"]["t0003"]["Mn"]["Mn"]) == convert(Array{Float32}, sol_sph.u[end][:,3])
+        end
 
+        h5open("Grt_2D.h5", "r") do file
+            @test read(file["Diffusion_Grt"]["t0000"]["Mg"]["Mg"]) == convert(Array{Float32}, IC2D.CMg0)
+            @test read(file["Diffusion_Grt"]["t0000"]["Fe"]["Fe"]) == convert(Array{Float32}, IC2D.CFe0)
+            @test read(file["Diffusion_Grt"]["t0000"]["Mn"]["Mn"]) == convert(Array{Float32}, IC2D.CMn0)
+            @test read(file["Diffusion_Grt"]["t0000"]["Ca"]["Ca"]) == convert(Array{Float32}, replace!((1 .- IC2D.CMg0 .- IC2D.CFe0 .- IC2D.CMn0), 1=>0))
+            @test read(file["Diffusion_Grt"]["t0003"]["Mg"]["Mg"]) == convert(Array{Float32}, sol_2D.u[end][:,:,1])
+            @test read(file["Diffusion_Grt"]["t0003"]["Fe"]["Fe"]) == convert(Array{Float32}, sol_2D.u[end][:,:,2])
+            @test read(file["Diffusion_Grt"]["t0003"]["Mn"]["Mn"]) == convert(Array{Float32}, sol_2D.u[end][:,:,3])
+        end
+
+        # delete files "Grt_1D.h5" and "Grt_Sph.h5" if it exists
+        if isfile("Grt_1D.h5")
+            rm("Grt_1D.h5")
+        end
+
+        if isfile("Grt_Sph.h5")
+            rm("Grt_Sph.h5")
+        end
+
+        if isfile("Grt_2D.h5")
+            rm("Grt_2D.h5")
+        end
+
+        save_data_callback = PresetTimeCallback(ustrip.(time_save) ./ domain2D.t_charact, save_data_paraview)
+
+        sol_2D = simulate(domain2D; callback=save_data_callback, path_save="Grt_2D.h5", progress=false, save_everystep=false, save_start=false)
+
+        h5open("Grt_2D.h5", "r") do file
+            @test read(file["Diffusion_Grt"]["t0000"]["Mg"]["Mg"])' == convert(Array{Float32}, IC2D.CMg0)
+            @test read(file["Diffusion_Grt"]["t0000"]["Fe"]["Fe"])' == convert(Array{Float32}, IC2D.CFe0)
+            @test read(file["Diffusion_Grt"]["t0000"]["Mn"]["Mn"])' == convert(Array{Float32}, IC2D.CMn0)
+            @test read(file["Diffusion_Grt"]["t0000"]["Ca"]["Ca"])' == convert(Array{Float32},replace!((1 .- IC2D.CMg0 .- IC2D.CFe0 .- IC2D.CMn0), 1=>0))
+            @test read(file["Diffusion_Grt"]["t0003"]["Mg"]["Mg"])' == convert(Array{Float32}, sol_2D.u[end][:,:,1])
+            @test read(file["Diffusion_Grt"]["t0003"]["Fe"]["Fe"])' == convert(Array{Float32}, sol_2D.u[end][:,:,2])
+            @test read(file["Diffusion_Grt"]["t0003"]["Mn"]["Mn"])' == convert(Array{Float32}, sol_2D.u[end][:,:,3])
+        end
+
+        # delete files if it exists
+        if isfile("Grt_2D.h5")
+            rm("Grt_2D.h5")
+        end
+
+        if isfile("Grt_2D.xdmf")
+            rm("Grt_2D.xdmf")
+        end
+
+    end
 end