add kwargs argument to simulate to fully support solve function

.github/dependabot.yml

@@ -0,0 +1,11 @@
+
+# https://docs.github.com/github/administering-a-repository/configuration-options-for-dependency-updates
+version: 2
+updates:
+  - package-ecosystem: "github-actions"
+    directory: "/" # Location of package manifests
+    schedule:
+      interval: "weekly"
+    ignore:
+      - dependency-name: "crate-ci/typos"
+        update-types: ["version-update:semver-patch"]

.github/workflows/CI.yml

@@ -30,7 +30,7 @@ jobs:
         with:
           version: ${{ matrix.version }}
           arch: ${{ matrix.arch }}
-      - uses: julia-actions/cache@v1
+      - uses: julia-actions/cache
       - uses: julia-actions/julia-buildpkg@v1
       - uses: julia-actions/julia-runtest@v1
       - uses: julia-actions/julia-processcoverage@v1

docs/src/diffusion_2D.md

@@ -91,8 +91,9 @@ Using the function `simulate()` to solve our system:
 
 ```julia
 # solve the problem using DifferentialEquations.jl
-sol = simulate(Domain2D)
+sol = simulate(Domain2D, progress=true, progress_steps=1)
 ```
+We use here the arguments `progress=true` and `progress_steps=1` to display a progress bar during the calculation.
 
 !!! note
     The calculation can take some time, about 6 minutes on my machine with 16 threads.

docs/src/diffusion_3D.md

@@ -63,7 +63,7 @@ We can now use the function `simulate()` to solve our system:
 
 ```julia
 # solve the problem using DifferentialEquations.jl
-sol = simulate(Domain3D; callback=save_data_callback, path_save=path_save, save_everystep=false)
+sol = simulate(Domain3D; callback=save_data_callback, path_save=path_save, save_everystep=false, progress=true, progress_steps=1)
 ```
 
 !!! note

examples/2D/2D_diffusion.jl

@@ -21,8 +21,7 @@ P      = 0.6u"GPa"
 IC2D     = InitialConditions2D(Mg0, Fe0, Mn0, Lx, Ly, tfinal; grt_boundary = grt_boundary)
 domain2D = Domain(IC2D, T, P)
 
-sol = simulate(domain2D; save_everystep=true)
-# 377.768768 seconds (16.12 M allocations: 15.863 GiB, 10.96% gc time, 5.41% compilation time)
+sol = simulate(domain2D; save_everystep=true, progress=true, progress_steps=1)
 
 @unpack tfinal_ad, t_charact = domain2D
 

examples/3D/3D_diffusion.jl

@@ -28,4 +28,4 @@ save_data_callback = PresetTimeCallback(time_save_ad, save_data_paraview)
 
 path_save = "Grt_3D.h5"  # chose the name and the path of the HDF5 output file (make sure to add .h5 or .hdf5 at the end)
 
-sol = simulate(domain3D; callback=save_data_callback, path_save=path_save, save_everystep=false);
+sol = simulate(domain3D; callback=save_data_callback, path_save=path_save, save_everystep=false, progress=true, progress_steps=1);

src/simulate/simulate.jl

@@ -1,18 +1,16 @@
 
 """
-    simulate(domain; callback=nothing, path_save=nothing, progress=true)
+    simulate(domain; path_save=nothing, kwargs...)
 
 Solve the coupled major element diffusion equations for a given domain using finite differences for the discretisation in space and return a solution type variable.
 
 The time discretisation is based on the ROCK2 method, a stabilized explicit method (Adbdulle and Medovikov, 2001 ; https://doi.org/10.1007/s002110100292) using OrdinaryDiffEq.jl.
 
 The solution type variable is following the format of OrdinaryDiffEq.jl (see https://docs.sciml.ai/DiffEqDocs/stable/basics/solution/), and can be used to plot the solution, and to extract the solution at a given time. As the system is nondimensionalised, the time of the solution is in nondimensional time.
 
-`callback` is an optional argument, which can be used to pass a callback function to the solver. It follows the format of DiffEqCallbacks.jl (see https://docs.sciml.ai/DiffEqCallbacks/stable/).
-
 `path_save` is an optional argument, which can be used to define the path of the HDF5 output file. Default is to nothing.
 
-`progress` is an optional argument, which can be used to display a progressbar during the simulation. Default is to true.
+All other accepted arguments such as `callback` or `progress` are the same as those of the `solve` function from DifferentialEquations (see https://docs.sciml.ai/DiffEqDocs/stable/basics/common_solver_opts/).
 """
 function simulate end
 
@@ -22,7 +20,7 @@ function simulate end
 Solve the coupled major element diffusion equations in 1D. Save all timesteps in the output solution type variable by default.
 
 """
-function simulate(domain::Domain1D; callback=nothing, path_save=nothing, progress=true, save_everystep=true)
+function simulate(domain::Domain1D; path_save=nothing, kwargs...)
 
     p = (domain = domain, path_save = path_save)
 
@@ -32,7 +30,7 @@ function simulate(domain::Domain1D; callback=nothing, path_save=nothing, progres
 
     prob = ODEProblem(semi_discretisation_diffusion_1D, u0, t, p)
 
-    sol = @time solve(prob, ROCK2(), progress=progress, progress_steps=1, save_start=true, abstol=1e-6,reltol=1e-6, callback=callback, save_everystep=save_everystep)
+    sol = @time solve(prob, ROCK2(); kwargs...)
 
     return sol
 end
@@ -42,7 +40,7 @@ end
 
 Solve the coupled major element diffusion equations in spherical coordinates. Save all timesteps in the output solution type variable by default.
 """
-function simulate(domain::DomainSpherical; callback=nothing, path_save=nothing, progress=true, save_everystep=true)
+function simulate(domain::DomainSpherical; path_save=nothing, kwargs...)
 
     p = (domain = domain, path_save = path_save)
 
@@ -52,7 +50,7 @@ function simulate(domain::DomainSpherical; callback=nothing, path_save=nothing,
 
     prob = ODEProblem(semi_discretisation_diffusion_spherical, u0, t, p)
 
-    sol = @time solve(prob, ROCK2(), progress=progress, progress_steps=1, save_start=true, abstol=1e-6,reltol=1e-6, callback=callback, save_everystep=save_everystep)
+    sol = @time solve(prob, ROCK2(); kwargs...)
 
     return sol
 end
@@ -62,7 +60,7 @@ end
 
 Solve the coupled major element diffusion equations in 2D. By default, save only the first and last timestep in the output solution type variable by default.
 """
-function simulate(domain::Domain2D; callback=nothing, path_save=nothing, progress=true, save_everystep=false, save_start=true)
+function simulate(domain::Domain2D; path_save=nothing, kwargs...)
 
     p = (domain = domain, path_save = path_save)
 
@@ -72,7 +70,7 @@ function simulate(domain::Domain2D; callback=nothing, path_save=nothing, progres
 
     prob = ODEProblem(semi_discretisation_diffusion_2D, u0, t, p)
 
-    sol = @time solve(prob, ROCK2(), progress=progress, progress_steps=1, save_start=save_start, abstol=1e-6,reltol=1e-6, save_everystep=save_everystep, callback=callback)
+    sol = @time solve(prob, ROCK2(); kwargs...)
 
     return sol
 end
@@ -83,7 +81,7 @@ end
 Solve the coupled major element diffusion equations in 3D. Save only the first and last timestep in the output solution type variable by default.
 
 """
-function simulate(domain::Domain3D; callback=nothing, path_save=nothing, progress=true, save_everystep=false, save_start=true)
+function simulate(domain::Domain3D; path_save=nothing, kwargs...)
 
     p = (domain = domain, path_save = path_save)
 
@@ -93,7 +91,7 @@ function simulate(domain::Domain3D; callback=nothing, path_save=nothing, progres
 
     prob = ODEProblem(semi_discretisation_diffusion_3D, u0, t, p)
 
-    sol = @time solve(prob, ROCK2(), progress=progress, progress_steps=1, save_start=save_start, abstol=1e-6,reltol=1e-6, save_everystep=save_everystep, callback=callback)
+    sol = @time solve(prob, ROCK2(); kwargs...)
 
     return sol
 end

test/runtests.jl

@@ -102,7 +102,7 @@ end
 
     domain1D = Domain(IC1D, T, P)
 
-    sol = simulate(domain1D; progress=false)
+    sol = simulate(domain1D; progress=false, abstol=1e-6,reltol=1e-6)
 
     @test norm(sum.(sol.u[end][:,1] .+ sol.u[end][:,2] .+ sol.u[end][:,3])) ≈ 28.64886878627501
 end
@@ -152,7 +152,7 @@ end
     IC2D = InitialConditions2D(CMg, CFe, CMn, Lx, Ly, tfinal; grt_boundary = grt_boundary)
     domain2D = Domain(IC2D, T, P)
 
-    sol = simulate(domain2D; save_everystep=false, progress=false)
+    sol = simulate(domain2D; save_everystep=false, progress=false, abstol=1e-6,reltol=1e-6)
 
     @test norm(sol.u[end][:,:,1]) ≈ 12.783357041653609
 end
@@ -174,7 +174,7 @@ end
     IC3D = InitialConditions3D(Mg0, Fe0, Mn0, Lx, Ly, Lz, tfinal; grt_boundary = grt_boundary)
     domain3D = Domain(IC3D, T, P)
 
-    sol = simulate(domain3D; save_everystep=false, progress=false);
+    sol = simulate(domain3D; save_everystep=false, progress=false, abstol=1e-6,reltol=1e-6);
 
     @test norm(sol.u[end][:,:,:,1]) ≈ 371.1477084396848
 end
@@ -222,13 +222,13 @@ end
 
     update_diffusion_coef_call = PresetTimeCallback(time_update_ad, update_diffusion_coef)
 
-    sol_sph = simulate(domainSph; callback=update_diffusion_coef_call, progress=false)
-    sol_1D = simulate(domain1D; callback=update_diffusion_coef_call, progress=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)
 
     @unpack time_update_ad = domain2D
     update_diffusion_coef_call = PresetTimeCallback(time_update_ad, update_diffusion_coef)
 
-    sol_2D = simulate(domain2D; callback=update_diffusion_coef_call, progress=false, save_everystep=false)
+    sol_2D = simulate(domain2D; callback=update_diffusion_coef_call, save_everystep=false, abstol=1e-6,reltol=1e-6)
 
     T=600  # in °C
     P=3  # in kbar
@@ -276,9 +276,9 @@ end
 
     save_data_callback = PresetTimeCallback(ustrip.(time_save) ./ domain1D.t_charact, save_data)
 
-    sol_1D = simulate(domain1D; callback=save_data_callback, path_save=(@__DIR__) * "/Grt_1D.h5", progress=false)
+    sol_1D = simulate(domain1D; callback=save_data_callback, path_save=(@__DIR__) * "/Grt_1D.h5", abstol=1e-6,reltol=1e-6)
 
-    sol_sph = simulate(domainSph; callback=save_data_callback, path_save=(@__DIR__) * "/Grt_Sph.h5", progress=false)
+    sol_sph = simulate(domainSph; callback=save_data_callback, path_save=(@__DIR__) * "/Grt_Sph.h5", abstol=1e-6,reltol=1e-6)
 
     save_data_callback = PresetTimeCallback(ustrip.(time_save) ./ domain2D.t_charact, save_data)