From 0a5ce9f1533af6ec16ae3bdbf5fa91540c89012d Mon Sep 17 00:00:00 2001 From: "Documenter.jl" Date: Sat, 11 Jan 2025 09:48:09 +0000 Subject: [PATCH] build based on 1a464fd --- dev/.documenter-siteinfo.json | 2 +- dev/assets/documenter.js | 302 +++++++++++++++-------------- dev/background/index.html | 2 +- dev/diffusion_1D/index.html | 8 +- dev/diffusion_2D/index.html | 7 +- dev/diffusion_3D/index.html | 2 +- dev/diffusion_spherical/index.html | 8 +- dev/index.html | 2 +- dev/reference/index.html | 2 +- dev/saving_output/index.html | 2 +- dev/search_index.js | 2 +- dev/updating_PT/index.html | 2 +- 12 files changed, 182 insertions(+), 159 deletions(-) diff --git a/dev/.documenter-siteinfo.json b/dev/.documenter-siteinfo.json index e266061..4aebbcd 100644 --- a/dev/.documenter-siteinfo.json +++ b/dev/.documenter-siteinfo.json @@ -1 +1 @@ -{"documenter":{"julia_version":"1.11.1","generation_timestamp":"2024-10-23T16:21:09","documenter_version":"1.7.0"}} \ No newline at end of file +{"documenter":{"julia_version":"1.11.2","generation_timestamp":"2025-01-11T09:48:02","documenter_version":"1.8.0"}} \ No newline at end of file diff --git a/dev/assets/documenter.js b/dev/assets/documenter.js index d4d5683..a7b7cf7 100644 --- a/dev/assets/documenter.js +++ b/dev/assets/documenter.js @@ -612,176 +612,194 @@ function worker_function(documenterSearchIndex, documenterBaseURL, filters) { }; } -// `worker = Threads.@spawn worker_function(documenterSearchIndex)`, but in JavaScript! -const filters = [ - ...new Set(documenterSearchIndex["docs"].map((x) => x.category)), -]; -const worker_str = - "(" + - worker_function.toString() + - ")(" + - JSON.stringify(documenterSearchIndex["docs"]) + - "," + - JSON.stringify(documenterBaseURL) + - "," + - JSON.stringify(filters) + - ")"; -const worker_blob = new Blob([worker_str], { type: "text/javascript" }); -const worker = new Worker(URL.createObjectURL(worker_blob)); - /////// SEARCH MAIN /////// -// Whether the worker is currently handling a search. This is a boolean -// as the worker only ever handles 1 or 0 searches at a time. -var worker_is_running = false; - -// The last search text that was sent to the worker. This is used to determine -// if the worker should be launched again when it reports back results. -var last_search_text = ""; - -// The results of the last search. This, in combination with the state of the filters -// in the DOM, is used compute the results to display on calls to update_search. -var unfiltered_results = []; - -// Which filter is currently selected -var selected_filter = ""; - -$(document).on("input", ".documenter-search-input", function (event) { - if (!worker_is_running) { - launch_search(); - } -}); - -function launch_search() { - worker_is_running = true; - last_search_text = $(".documenter-search-input").val(); - worker.postMessage(last_search_text); -} - -worker.onmessage = function (e) { - if (last_search_text !== $(".documenter-search-input").val()) { - launch_search(); - } else { - worker_is_running = false; - } - - unfiltered_results = e.data; - update_search(); -}; +function runSearchMainCode() { + // `worker = Threads.@spawn worker_function(documenterSearchIndex)`, but in JavaScript! + const filters = [ + ...new Set(documenterSearchIndex["docs"].map((x) => x.category)), + ]; + const worker_str = + "(" + + worker_function.toString() + + ")(" + + JSON.stringify(documenterSearchIndex["docs"]) + + "," + + JSON.stringify(documenterBaseURL) + + "," + + JSON.stringify(filters) + + ")"; + const worker_blob = new Blob([worker_str], { type: "text/javascript" }); + const worker = new Worker(URL.createObjectURL(worker_blob)); + + // Whether the worker is currently handling a search. This is a boolean + // as the worker only ever handles 1 or 0 searches at a time. + var worker_is_running = false; + + // The last search text that was sent to the worker. This is used to determine + // if the worker should be launched again when it reports back results. + var last_search_text = ""; + + // The results of the last search. This, in combination with the state of the filters + // in the DOM, is used compute the results to display on calls to update_search. + var unfiltered_results = []; + + // Which filter is currently selected + var selected_filter = ""; + + $(document).on("input", ".documenter-search-input", function (event) { + if (!worker_is_running) { + launch_search(); + } + }); -$(document).on("click", ".search-filter", function () { - if ($(this).hasClass("search-filter-selected")) { - selected_filter = ""; - } else { - selected_filter = $(this).text().toLowerCase(); + function launch_search() { + worker_is_running = true; + last_search_text = $(".documenter-search-input").val(); + worker.postMessage(last_search_text); } - // This updates search results and toggles classes for UI: - update_search(); -}); + worker.onmessage = function (e) { + if (last_search_text !== $(".documenter-search-input").val()) { + launch_search(); + } else { + worker_is_running = false; + } -/** - * Make/Update the search component - */ -function update_search() { - let querystring = $(".documenter-search-input").val(); + unfiltered_results = e.data; + update_search(); + }; - if (querystring.trim()) { - if (selected_filter == "") { - results = unfiltered_results; + $(document).on("click", ".search-filter", function () { + if ($(this).hasClass("search-filter-selected")) { + selected_filter = ""; } else { - results = unfiltered_results.filter((result) => { - return selected_filter == result.category.toLowerCase(); - }); + selected_filter = $(this).text().toLowerCase(); } - let search_result_container = ``; - let modal_filters = make_modal_body_filters(); - let search_divider = `
`; + // This updates search results and toggles classes for UI: + update_search(); + }); - if (results.length) { - let links = []; - let count = 0; - let search_results = ""; - - for (var i = 0, n = results.length; i < n && count < 200; ++i) { - let result = results[i]; - if (result.location && !links.includes(result.location)) { - search_results += result.div; - count++; - links.push(result.location); - } - } + /** + * Make/Update the search component + */ + function update_search() { + let querystring = $(".documenter-search-input").val(); - if (count == 1) { - count_str = "1 result"; - } else if (count == 200) { - count_str = "200+ results"; + if (querystring.trim()) { + if (selected_filter == "") { + results = unfiltered_results; } else { - count_str = count + " results"; + results = unfiltered_results.filter((result) => { + return selected_filter == result.category.toLowerCase(); + }); } - let result_count = `
${count_str}
`; - search_result_container = ` + let search_result_container = ``; + let modal_filters = make_modal_body_filters(); + let search_divider = `
`; + + if (results.length) { + let links = []; + let count = 0; + let search_results = ""; + + for (var i = 0, n = results.length; i < n && count < 200; ++i) { + let result = results[i]; + if (result.location && !links.includes(result.location)) { + search_results += result.div; + count++; + links.push(result.location); + } + } + + if (count == 1) { + count_str = "1 result"; + } else if (count == 200) { + count_str = "200+ results"; + } else { + count_str = count + " results"; + } + let result_count = `
${count_str}
`; + + search_result_container = ` +
+ ${modal_filters} + ${search_divider} + ${result_count} +
+ ${search_results} +
+
+ `; + } else { + search_result_container = `
${modal_filters} ${search_divider} - ${result_count} -
- ${search_results} -
-
+
0 result(s)
+ +
No result found!
`; - } else { - search_result_container = ` -
- ${modal_filters} - ${search_divider} -
0 result(s)
-
-
No result found!
- `; - } + } - if ($(".search-modal-card-body").hasClass("is-justify-content-center")) { - $(".search-modal-card-body").removeClass("is-justify-content-center"); - } + if ($(".search-modal-card-body").hasClass("is-justify-content-center")) { + $(".search-modal-card-body").removeClass("is-justify-content-center"); + } - $(".search-modal-card-body").html(search_result_container); - } else { - if (!$(".search-modal-card-body").hasClass("is-justify-content-center")) { - $(".search-modal-card-body").addClass("is-justify-content-center"); + $(".search-modal-card-body").html(search_result_container); + } else { + if (!$(".search-modal-card-body").hasClass("is-justify-content-center")) { + $(".search-modal-card-body").addClass("is-justify-content-center"); + } + + $(".search-modal-card-body").html(` +
Type something to get started!
+ `); } + } - $(".search-modal-card-body").html(` -
Type something to get started!
- `); + /** + * Make the modal filter html + * + * @returns string + */ + function make_modal_body_filters() { + let str = filters + .map((val) => { + if (selected_filter == val.toLowerCase()) { + return `${val}`; + } else { + return `${val}`; + } + }) + .join(""); + + return ` +
+ Filters: + ${str} +
`; } } -/** - * Make the modal filter html - * - * @returns string - */ -function make_modal_body_filters() { - let str = filters - .map((val) => { - if (selected_filter == val.toLowerCase()) { - return `${val}`; - } else { - return `${val}`; - } - }) - .join(""); - - return ` -
- Filters: - ${str} -
`; +function waitUntilSearchIndexAvailable() { + // It is possible that the documenter.js script runs before the page + // has finished loading and documenterSearchIndex gets defined. + // So we need to wait until the search index actually loads before setting + // up all the search-related stuff. + if (typeof documenterSearchIndex !== "undefined") { + runSearchMainCode(); + } else { + console.warn("Search Index not available, waiting"); + setTimeout(waitUntilSearchIndexAvailable, 1000); + } } +// The actual entry point to the search code +waitUntilSearchIndexAvailable(); + }) //////////////////////////////////////////////////////////////////////////////// require(['jquery'], function($) { diff --git a/dev/background/index.html b/dev/background/index.html index cbad57d..1caff06 100644 --- a/dev/background/index.html +++ b/dev/background/index.html @@ -17,4 +17,4 @@ \end{equation}\]

with $D_i^*$ the tracer diffusion coefficient of the $i^{th}$ component, $\delta_{ij}$ the Kronecker delta, $z$ the charge of the species, and $D_n^*$ the tracer diffusion coefficient of the dependent molar fraction (here Ca).

The tracer diffusion coefficient $D_i^*$ can be calculated from an Arrhenius relationship:

\[\begin{equation} \label{arrhenius} D_i^* = D_{0,i} \exp \left( -\frac{E_{a,i} - (P-1)\Delta V^+_i}{RT} \right), -\end{equation}\]

with $D_{0,i}$ the pre-exponential constant, $E_{a,i}$ the activation energy of diffusion, $\Delta V^+_i$ the activation volume of diffusion at 1 bar, $R$ the universal gas constant, $T$ the temperature, and $P$ the pressure.

In DiffusionGarnet.jl, $D_{0,i}$, $E_{a,i}$, and $\Delta V^+_i$ are those of Chakraborty & Ganguly (1992) [2]. The tracer diffusion coefficient of Ca is defined as $0.5D_{Fe}$, following the approach of Loomis et al. (1985) [3].

Numerical approach

By defining the PT conditions of the metamorphic event of interest, (3) can be solved for each component, and the diffusion coefficient tensor can be calculated using (2) from the initial major composition data. In DiffusionGarnet.jl, (1) is then discretised in space using finite differences, and the resulting system of ordinary differential equations is solved with ROCK2, a stabilised explicit method (Abdulle & Medovikov, 2001 [4]) using the DifferentialEquations.jl ecosystem.

References

[1] Lasaga, A. C. (1979). Multicomponent exchange and diffusion in silicates. Geochimica et Cosmochimica Acta, 43(4), 455-469.

[2] Chakraborty, S., & Ganguly, J. (1992). Cation diffusion in aluminosilicate garnets: experimental determination in spessartine-almandine diffusion couples, evaluation of effective binary diffusion coefficients, and applications. Contributions to Mineralogy and petrology, 111(1), 74-86.

[3] Loomis, T. P., Ganguly, J., Elphick, S. C., 1985. Experimental determinations of cation diffusitivities in aluminosilicate garnets. II. Multicomponent simulation and tracer diffusion coefficients. Contributions to Mineralogy and Petrology 90, 45–51.

[4] Abdulle, A., & Medovikov, A. A. (2001). Second order Chebyshev methods based on orthogonal polynomials. Numerische Mathematik, 90, 1-18.

+\end{equation}\]

with $D_{0,i}$ the pre-exponential constant, $E_{a,i}$ the activation energy of diffusion, $\Delta V^+_i$ the activation volume of diffusion at 1 bar, $R$ the universal gas constant, $T$ the temperature, and $P$ the pressure.

In DiffusionGarnet.jl, $D_{0,i}$, $E_{a,i}$, and $\Delta V^+_i$ are those of Chakraborty & Ganguly (1992) [2]. The tracer diffusion coefficient of Ca is defined as $0.5D_{Fe}$, following the approach of Loomis et al. (1985) [3].

Numerical approach

By defining the PT conditions of the metamorphic event of interest, (3) can be solved for each component, and the diffusion coefficient tensor can be calculated using (2) from the initial major composition data. In DiffusionGarnet.jl, (1) is then discretised in space using finite differences, and the resulting system of ordinary differential equations is solved with ROCK2, a stabilised explicit method (Abdulle & Medovikov, 2001 [4]) using the DifferentialEquations.jl ecosystem.

References

[1] Lasaga, A. C. (1979). Multicomponent exchange and diffusion in silicates. Geochimica et Cosmochimica Acta, 43(4), 455-469.

[2] Chakraborty, S., & Ganguly, J. (1992). Cation diffusion in aluminosilicate garnets: experimental determination in spessartine-almandine diffusion couples, evaluation of effective binary diffusion coefficients, and applications. Contributions to Mineralogy and petrology, 111(1), 74-86.

[3] Loomis, T. P., Ganguly, J., Elphick, S. C., 1985. Experimental determinations of cation diffusitivities in aluminosilicate garnets. II. Multicomponent simulation and tracer diffusion coefficients. Contributions to Mineralogy and Petrology 90, 45–51.

[4] Abdulle, A., & Medovikov, A. A. (2001). Second order Chebyshev methods based on orthogonal polynomials. Numerische Mathematik, 90, 1-18.

diff --git a/dev/diffusion_1D/index.html b/dev/diffusion_1D/index.html index 8fc7266..50c8b56 100644 --- a/dev/diffusion_1D/index.html +++ b/dev/diffusion_1D/index.html @@ -30,11 +30,13 @@ # define a Domain struct containing the definition of our problem and nondimensionalised variables domain1D = Domain(IC1D, T, P)
Note

Lx, tfinal, T and P need to contain units, following the syntax of the package Unitful. This allows the user to specify the units that suit their problem.

Domain1D contains all the information that DiffusionGarnet needs to solve our coupled diffusion problem, at 900 °C and 0.6 GPa for a duration of 15 Myr.

This can be achieved with the function simulate():

# solve the problem using DifferentialEquations.jl
-sol = simulate(domain1D)

which outputs the time spent on the solver, for example, on the second run:

  0.399870 seconds (31.93 k allocations: 18.212 MiB)

simulate() uses the DifferentialEquations.jl 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. It basically holds all the information from our simulation.

We can now plot the solution to our problem. By default, DiffusionGarnet saves the time in Myr.


+sol = simulate(domain1D)

which outputs the time spent on the solver, for example, on the second run:

  0.399870 seconds (31.93 k allocations: 18.212 MiB)

simulate() uses the DifferentialEquations.jl 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. It basically holds all the information from our simulation.

We can now plot the solution to our problem.

# extract characteristic time to convert back to dimensional time
+@unpack t_charact = domain1D
+
 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;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*t_charact;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)")
@@ -51,4 +53,4 @@
 
 println("Now, generating the gif...")
 gif(anim, "Grt_1D.gif", fps = 7)
-println("...Done!")

Here is the resulting gif obtained:

1D diffusion profil of a garnet

It shows the compositional evolution of a 1D profile through a garnet grain with homogeneous Dirichlet boundaries on both sides.

+println("...Done!")

Here is the resulting gif obtained:

1D diffusion profil of a garnet

It shows the compositional evolution of a 1D profile through a garnet grain with homogeneous Dirichlet boundaries on both sides.

diff --git a/dev/diffusion_2D/index.html b/dev/diffusion_2D/index.html index 81320ca..a9b7ac3 100644 --- a/dev/diffusion_2D/index.html +++ b/dev/diffusion_2D/index.html @@ -35,11 +35,12 @@ plot(p1, p2, p3, p4, layout = l , plot_title="Initial conditions")

which outputs:

Initial conditions 2D.

Like with spherical or 1D Cartesian coordinates, we need to define our InitialConditions and Domain structures, that will hold all the information of our model.

IC2D = InitialConditions2D(Mg0, Fe0, Mn0, Lx, Ly, tfinal; grt_boundary = grt_boundary)
 domain2D = Domain(IC2D, T, P)
Note

The argument grt_boundary is a binary array specifiying the coordinates where the homogeneous Dirichlet boundaries are to be enforced, here on the contour of the grain. If not provided, the homogeneous Dirichlet boundaries are applied to the sides of the domain.

Using the function simulate() to solve our system:

# solve the problem using DifferentialEquations.jl
-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.

We can now plot the results:


+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.

We can now plot the results:

@unpack t_charact = domain2D
+
 println("Plotting...")
 anim = @animate for i = tqdm(LinRange(0, sol.t[end], 20))
 
-    time = round(i;digits=2)
+    time = round(i*t_charact;digits=2)
 
     l = @layout [a b ; c d ]
     Ca = 1 .- sol(i)[:,:,1] .- sol(i)[:,:,2] .- sol(i)[:,:,3]
@@ -58,4 +59,4 @@
 println("Now, generating the gif...")
 gif(anim, "Grt_2D_900.gif", fps = 3)
 println("...Done!")
-

Here is the resulting gif:

2D diffusion of a garnet

Note

2D modelling ignores the effect of the third dimension on diffusion and is equivalent to considering the garnet grain to be cylindrical rather than spherical.

+

Here is the resulting gif:

2D diffusion of a garnet

Note

2D modelling ignores the effect of the third dimension on diffusion and is equivalent to considering the garnet grain to be cylindrical rather than spherical.

diff --git a/dev/diffusion_3D/index.html b/dev/diffusion_3D/index.html index 883fcb4..be9aafd 100644 --- a/dev/diffusion_3D/index.html +++ b/dev/diffusion_3D/index.html @@ -26,4 +26,4 @@ save_data_callback = PresetTimeCallback(time_save_ad, save_data_paraview) path_save = "Grt_3D.h5" # choose the name and the path of the HDF5 output file (make sure to add .h5 or .hdf5 at the end)
Note

We used here save_data_paraview() instead of save_data() to save the data in a format that can be read by Paraview.

We can now use the function simulate() to solve our system:

# solve the problem using DifferentialEquations.jl
-sol = simulate(Domain3D; callback=save_data_callback, path_save=path_save, save_everystep=false, progress=true, progress_steps=1)
Note

The calculation can take some time, about 4 minutes on my machine with 16 threads.

Two files should have been created in the same folder as your current session: Grt_3D.h5 and Grt_3D.xdmf. The first one is the HDF5 file containing the results of the simulation, and the second one is an XDMF file describing the HDF5 file. This last file can be opened with visualisation software programs, such as Paraview.

Warning

For any visualisation software, make sure you open the XDMF file and not the HDF5 file. For Paraview, select the XDMF Reader as the reader when you open your data. Only Paraview has been tested with this package, but other software should work as well.

+sol = simulate(Domain3D; callback=save_data_callback, path_save=path_save, save_everystep=false, progress=true, progress_steps=1)
Note

The calculation can take some time, about 4 minutes on my machine with 16 threads.

Two files should have been created in the same folder as your current session: Grt_3D.h5 and Grt_3D.xdmf. The first one is the HDF5 file containing the results of the simulation, and the second one is an XDMF file describing the HDF5 file. This last file can be opened with visualisation software programs, such as Paraview.

Warning

For any visualisation software, make sure you open the XDMF file and not the HDF5 file. For Paraview, select the XDMF Reader as the reader when you open your data. Only Paraview has been tested with this package, but other software should work as well.

diff --git a/dev/diffusion_spherical/index.html b/dev/diffusion_spherical/index.html index 158607f..d86a6ad 100644 --- a/dev/diffusion_spherical/index.html +++ b/dev/diffusion_spherical/index.html @@ -30,10 +30,12 @@ DomainSph = Domain(ICSph, T, P) Domain1D = Domain(IC1D, T, P; bc_neumann = (true, false))
Tip

InitialConditionsSpherical always assumes that the core of the profile is on the left and the edge is on the right side.

Note

We use bc_neumann in Domain to indicate that a homogeneous Neumann boundary must be ensured on the left side of the 1D profile to simulate the core of the garnet grain.

We can then use the function simulate() to solve our system:

# solve the problem using DifferentialEquations.jl
 sol_sph = simulate(DomainSph)
-sol_1D = simulate(Domain1D)

and compare the 2 solutions:

anim = @animate for i = LinRange(0, sol_sph.t[end], 100)
+sol_1D = simulate(Domain1D)

and compare the 2 solutions:

@unpack t_charact = DomainSph
+
+anim = @animate for i = LinRange(0, sol_sph.t[end], 100)
     l = @layout [a ; b]
 
-    p1 = plot(distance, Fe0, label="Fe initial", linestyle = :dash, linewidth=1, dpi=200, title = "Total Time = $(round(i;digits=2)) Ma", legend=:outerbottomright, linecolor=1,xlabel = "Distance (µm)")
+    p1 = plot(distance, Fe0, label="Fe initial", linestyle = :dash, linewidth=1, dpi=200, title = "Total Time = $(round(i*t_charact;digits=2)) Ma", legend=:outerbottomright, linecolor=1,xlabel = "Distance (µm)")
     p1 = plot!(distance, sol_sph(i)[:,2], label="Fe Sph",linecolor=1, linewidth=1)
     p1 = plot!(distance, sol_1D(i)[:,2], label="Fe 1D",linecolor=5, linewidth=1)
 
@@ -55,4 +57,4 @@
 
 println("Now, generating the gif...")
 gif(anim, "Grt_Spherical+1D.gif", fps = 7)
-println("...Done!")

With the resulting gif:

Spherical diffusion profil of a garnet

It shows that using spherical coordinates makes the garnet diffuse faster. Using 1D Cartesian coordinates can overestimate the equilibration time.

+println("...Done!")

With the resulting gif:

Spherical diffusion profil of a garnet

It shows that using spherical coordinates makes the garnet diffuse faster. Using 1D Cartesian coordinates can overestimate the equilibration time.

diff --git a/dev/index.html b/dev/index.html index cd7478a..bd73854 100644 --- a/dev/index.html +++ b/dev/index.html @@ -1,4 +1,4 @@ Home · DiffusionGarnet.jl
GitHub

DiffusionGarnet.jl

This is the documentation page for DiffusionGarnet.jl, a high-level Julia package to do coupled diffusion modelling of major elements on real garnet data. It currently supports 1D, spherical, 2D and 3D coordinates for evenly spaced data.

It is built on top of the DifferentialEquations.jl package ecosystem and uses Unitful.jl to allow the user to define appropriate units for their problems. For 2D and 3D models, it uses ParallelStencil.jl to support multithreading on CPU and parallel computing on GPU.

This package can be used as a teaching tool or for research purposes.

Installation

DiffusionGarnet.jl is a registered package and may be installed directly from the REPL:

julia>]
   pkg> add DiffusionGarnet
-  pkg> test DiffusionGarnet
+ pkg> test DiffusionGarnet diff --git a/dev/reference/index.html b/dev/reference/index.html index 860d758..dd739a0 100644 --- a/dev/reference/index.html +++ b/dev/reference/index.html @@ -1,2 +1,2 @@ -List of functions · DiffusionGarnet.jl
GitHub
DiffusionGarnet.Domain — Type
Domain(IC::InitialConditions3D, T::Union{Unitful.Temperature,Array{<:Unitful.Temperature{<:Real}, 1}}, P::Union{Unitful.Pressure,Array{<:Unitful.Pressure{<:Real}, 1}}, time_update::Union{Unitful.Time,Array{<:Unitful.Time{<:Real}, 1}}=0u"Myr")

When applied to 3D initial conditions, define corresponding 3D domain.

source
DiffusionGarnet.Domain — Type
Domain(IC::InitialConditions2D, T::Union{Unitful.Temperature,Array{<:Unitful.Temperature{<:Real}, 1}}, P::Union{Unitful.Pressure,Array{<:Unitful.Pressure{<:Real}, 1}}, time_update::Union{Unitful.Time,Array{<:Unitful.Time{<:Real}, 1}}=0u"Myr")

When applied to 2D initial conditions, define corresponding 2D domain.

source
DiffusionGarnet.Domain — Type
Domain(IC::InitialConditionsSpherical, T::Union{Unitful.Temperature,Array{<:Unitful.Temperature{<:Real}, 1}}, P::Union{Unitful.Pressure,Array{<:Unitful.Pressure{<:Real}, 1}}, time_update::Union{Unitful.Time,Array{<:Unitful.Time{<:Real}, 1}}=0u"Myr")

When applied to spherical initial conditions, define corresponding spherical domain. Assume that the center of the grain is on the left side.

source
DiffusionGarnet.Domain — Type
Domain(IC::InitialConditions1D, T::Union{Unitful.Temperature,Array{<:Unitful.Temperature{<:Real}, 1}}, P::Union{Unitful.Pressure,Array{<:Unitful.Pressure{<:Real}, 1}}, time_update::Union{Unitful.Time,Array{<:Unitful.Time{<:Real}, 1}}=0u"Myr"; bc_neumann::Tuple=(false, false))

When applied to 1D initial conditions, define corresponding 1D domain. bc_neumann can be used to define Neumann boundary conditions on the left or right side of the domain if set to true.

source
DiffusionGarnet.Domain — Type
Domain(IC, T, P, time_update=0u"Myr")

Return a struct containing the struct IC, the PT conditions T and P and the nondimensionalised parameters of the problem. time_update can be used to update the PT conditions during the simulation. If so, T and P have to be of the same size as time_update and correspond to the values of PT to update to.

source
DiffusionGarnet.InitialConditions1D — Method
InitialConditions1D(CMg0::Array{<:Real, 1}, CFe0::Array{<:Real, 1}, CMn0::Array{<:Real, 1}, Lx::Unitful.Length, tfinal::Unitful.Time)

Return a structure containing the initial conditions for a 1D diffusion problem. CMg0, CFe0 and CMn0 need to be in molar fraction. Convert the Lxandtfinal`` to µm and Myr respectively.

source
DiffusionGarnet.InitialConditions2D — Method
InitialConditions2D(CMg0::AbstractArray{<:Real, 2}, CFe0::AbstractArray{<:Real, 2}, CMn0::AbstractArray{<:Real, 2}, Lx::Unitful.Length, Ly::Unitful.Length, tfinal::Unitful.Time; grt_boundary::AbstractArray{<:Real, 2}=zeros(eltype(CMg0), size(CMg0)...))

Return a structure containing the initial conditions for a 2D diffusion problem. CMg0, CFe0 and CMn0 need to be in molar fraction. Convert Lx, Ly and tfinal to µm, µm and Myr respectively. grt_boundary is a matrix of the same size as CMg0, CFe0 and CMn0. It is a binary matrix with 1 where the contour of the garnet is present and 0 otherwise. It is used to apply the Dirichlet boundary condition in the model. If not provided, it is equal to zero everywhere.

source
DiffusionGarnet.InitialConditions3D — Method
InitialConditions3D(CMg0::AbstractArray{<:Real, 3}, CFe0::AbstractArray{<:Real, 3}, CMn0::AbstractArray{<:Real, 3}, Lx::Unitful.Length, Ly::Unitful.Length, Lz::Unitful.Length, tfinal::Unitful.Time; grt_boundary::Union{AbstractArray{<:Real, 3}, AbstractArray{<:Bool, 3}}=zeros(Bool, size(CMg0)...))

Return a structure containing the initial conditions for a 3D diffusion problem. CMg0, CFe0 and CMn0 need to be in molar fraction. Convert Lx, Ly, Lz and tfinal to µm, µm, µm and Myr respectively. grt_boundary is a matrix of the same size as CMg0, CFe0 and CMn0. It is a binary matrix with 1 where the contour of the garnet is present and 0 otherwise. It is used to apply the Dirichlet boundary condition in the model. If not provided, it is equal to zero everywhere.

source
DiffusionGarnet.InitialConditionsSpherical — Method
InitialConditionsSpherical(CMg0::AbstractArray{<:Real, 1}, CFe0::AbstractArray{<:Real, 1}, CMn0::AbstractArray{<:Real, 1}, Lr::Unitful.Length, tfinal::Unitful.Time)

Return a structure containing the initial conditions for a spherical diffusion problem. CMg0, CFe0 and CMn0 need to be in molar fraction. Convert Lr and tfinal to µm and Myr respectively.

source
DiffusionGarnet.save_data — Method
save_data(integrator)

Callback function used to save major element compositions to an HDF5 file at a specific timestep.

source
DiffusionGarnet.save_data_paraview — Method
save_data_paraview(integrator)

Callback function used to save data to an HDF5 file and produce an XMDF file. XMDF files can be used to read data on software like Paraview.

Note that data are saved in row major format. For column major, use save_data() instead.

source
DiffusionGarnet.simulate — Function
simulate(domain; path_save=nothing, solver=ROCK2(), 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. The time of the solution is in Myr.

path_save is an optional argument, which can be used to define the path of the HDF5 output file. Default is to nothing.

solver is an optional argument, which can be used to define the solver to use for the time discretisation. Default is the ROCK2 method. All other ODE solvers accepted as the ones from DifferentialEquations (see https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/).

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/commonsolveropts/).

source
DiffusionGarnet.simulate — Method
simulate(domain::Domain1D; path_save=nothing, solver=ROCK2(), abstol=1e-8, reltol=1e-6, kwargs...)

Solve the coupled major element diffusion equations in 1D. Save all timesteps in the output solution type variable by default.

source
DiffusionGarnet.simulate — Method
simulate(domain::Domain2D; path_save=nothing, solver=ROCK2(), kwargs...)

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.

source
DiffusionGarnet.simulate — Method
simulate(domain::Domain3D; path_save=nothing, solver=ROCK2(), kwargs...)

Solve the coupled major element diffusion equations in 3D. Save only the first and last timestep in the output solution type variable by default.

source
DiffusionGarnet.simulate — Method
simulate(domain::DomainSpherical; path_save=nothing, solver=ROCK2(), abstol=1e-8, reltol=1e-6, kwargs...)

Solve the coupled major element diffusion equations in spherical coordinates. Save all timesteps in the output solution type variable by default.

source
DiffusionGarnet.update_diffusion_coef — Method
update_diffusion_coef(integrator)

Callback function to update the diffusion coefficients at a given time from a new pressure and temperature. To use with the callback PresetTimeCallback (https://docs.sciml.ai/stable/basics/callbacks/#PresetTimeCallback-1).

Follows the syntax of callback functions defined by DiffEqCallbacks.jl (https://docs.sciml.ai/DiffEqCallbacks/stable/).

source
+List of functions · DiffusionGarnet.jl
GitHub
DiffusionGarnet.Domain — Type
Domain(IC::InitialConditions2D, T::Union{Unitful.Temperature,Array{<:Unitful.Temperature{<:Real}, 1}}, P::Union{Unitful.Pressure,Array{<:Unitful.Pressure{<:Real}, 1}}, time_update::Union{Unitful.Time,Array{<:Unitful.Time{<:Real}, 1}}=0u"Myr")

When applied to 2D initial conditions, define corresponding 2D domain.

source
DiffusionGarnet.Domain — Type
Domain(IC::InitialConditions1D, T::Union{Unitful.Temperature,Array{<:Unitful.Temperature{<:Real}, 1}}, P::Union{Unitful.Pressure,Array{<:Unitful.Pressure{<:Real}, 1}}, time_update::Union{Unitful.Time,Array{<:Unitful.Time{<:Real}, 1}}=0u"Myr"; bc_neumann::Tuple=(false, false))

When applied to 1D initial conditions, define corresponding 1D domain. bc_neumann can be used to define Neumann boundary conditions on the left or right side of the domain if set to true.

source
DiffusionGarnet.Domain — Type
Domain(IC::InitialConditionsSpherical, T::Union{Unitful.Temperature,Array{<:Unitful.Temperature{<:Real}, 1}}, P::Union{Unitful.Pressure,Array{<:Unitful.Pressure{<:Real}, 1}}, time_update::Union{Unitful.Time,Array{<:Unitful.Time{<:Real}, 1}}=0u"Myr")

When applied to spherical initial conditions, define corresponding spherical domain. Assume that the center of the grain is on the left side.

source
DiffusionGarnet.Domain — Type
Domain(IC::InitialConditions3D, T::Union{Unitful.Temperature,Array{<:Unitful.Temperature{<:Real}, 1}}, P::Union{Unitful.Pressure,Array{<:Unitful.Pressure{<:Real}, 1}}, time_update::Union{Unitful.Time,Array{<:Unitful.Time{<:Real}, 1}}=0u"Myr")

When applied to 3D initial conditions, define corresponding 3D domain.

source
DiffusionGarnet.Domain — Type
Domain(IC, T, P, time_update=0u"Myr")

Return a struct containing the struct IC, the PT conditions T and P and the nondimensionalised parameters of the problem. time_update can be used to update the PT conditions during the simulation. If so, T and P have to be of the same size as time_update and correspond to the values of PT to update to.

source
DiffusionGarnet.InitialConditions1D — Method
InitialConditions1D(CMg0::Array{<:Real, 1}, CFe0::Array{<:Real, 1}, CMn0::Array{<:Real, 1}, Lx::Unitful.Length, tfinal::Unitful.Time)

Return a structure containing the initial conditions for a 1D diffusion problem. CMg0, CFe0 and CMn0 need to be in molar fraction. Convert the Lxandtfinal`` to µm and Myr respectively.

source
DiffusionGarnet.InitialConditions2D — Method
InitialConditions2D(CMg0::AbstractArray{<:Real, 2}, CFe0::AbstractArray{<:Real, 2}, CMn0::AbstractArray{<:Real, 2}, Lx::Unitful.Length, Ly::Unitful.Length, tfinal::Unitful.Time; grt_boundary::AbstractArray{<:Real, 2}=zeros(eltype(CMg0), size(CMg0)...))

Return a structure containing the initial conditions for a 2D diffusion problem. CMg0, CFe0 and CMn0 need to be in molar fraction. Convert Lx, Ly and tfinal to µm, µm and Myr respectively. grt_boundary is a matrix of the same size as CMg0, CFe0 and CMn0. It is a binary matrix with 1 where the contour of the garnet is present and 0 otherwise. It is used to apply the Dirichlet boundary condition in the model. If not provided, it is equal to zero everywhere.

source
DiffusionGarnet.InitialConditions3D — Method
InitialConditions3D(CMg0::AbstractArray{<:Real, 3}, CFe0::AbstractArray{<:Real, 3}, CMn0::AbstractArray{<:Real, 3}, Lx::Unitful.Length, Ly::Unitful.Length, Lz::Unitful.Length, tfinal::Unitful.Time; grt_boundary::Union{AbstractArray{<:Real, 3}, AbstractArray{<:Bool, 3}}=zeros(Bool, size(CMg0)...))

Return a structure containing the initial conditions for a 3D diffusion problem. CMg0, CFe0 and CMn0 need to be in molar fraction. Convert Lx, Ly, Lz and tfinal to µm, µm, µm and Myr respectively. grt_boundary is a matrix of the same size as CMg0, CFe0 and CMn0. It is a binary matrix with 1 where the contour of the garnet is present and 0 otherwise. It is used to apply the Dirichlet boundary condition in the model. If not provided, it is equal to zero everywhere.

source
DiffusionGarnet.InitialConditionsSpherical — Method
InitialConditionsSpherical(CMg0::AbstractArray{<:Real, 1}, CFe0::AbstractArray{<:Real, 1}, CMn0::AbstractArray{<:Real, 1}, Lr::Unitful.Length, tfinal::Unitful.Time)

Return a structure containing the initial conditions for a spherical diffusion problem. CMg0, CFe0 and CMn0 need to be in molar fraction. Convert Lr and tfinal to µm and Myr respectively.

source
DiffusionGarnet.save_data — Method
save_data(integrator)

Callback function used to save major element compositions to an HDF5 file at a specific timestep.

source
DiffusionGarnet.save_data_paraview — Method
save_data_paraview(integrator)

Callback function used to save data to an HDF5 file and produce an XMDF file. XMDF files can be used to read data on software like Paraview.

Note that data are saved in row major format. For column major, use save_data() instead.

source
DiffusionGarnet.simulate — Function
simulate(domain; path_save=nothing, solver=ROCK2(), 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. The time of the solution is non-dimensional but can be converted back using the characteristic time (t_charact contained in the Domain structure).

path_save is an optional argument, which can be used to define the path of the HDF5 output file. Default is to nothing.

solver is an optional argument, which can be used to define the solver to use for the time discretisation. Default is the ROCK2 method. All other ODE solvers accepted as the ones from DifferentialEquations (see https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/).

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/commonsolveropts/).

source
DiffusionGarnet.simulate — Method
simulate(domain::Domain1D; path_save=nothing, solver=ROCK2(), abstol=1e-8, reltol=1e-6, kwargs...)

Solve the coupled major element diffusion equations in 1D. Save all timesteps in the output solution type variable by default.

source
DiffusionGarnet.simulate — Method
simulate(domain::Domain2D; path_save=nothing, solver=ROCK2(), kwargs...)

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.

source
DiffusionGarnet.simulate — Method
simulate(domain::Domain3D; path_save=nothing, solver=ROCK2(), kwargs...)

Solve the coupled major element diffusion equations in 3D. Save only the first and last timestep in the output solution type variable by default.

source
DiffusionGarnet.simulate — Method
simulate(domain::DomainSpherical; path_save=nothing, solver=ROCK2(), abstol=1e-8, reltol=1e-6, kwargs...)

Solve the coupled major element diffusion equations in spherical coordinates. Save all timesteps in the output solution type variable by default.

source
DiffusionGarnet.update_diffusion_coef — Method
update_diffusion_coef(integrator)

Callback function to update the diffusion coefficients at a given time from a new pressure and temperature. To use with the callback PresetTimeCallback (https://docs.sciml.ai/stable/basics/callbacks/#PresetTimeCallback-1).

Follows the syntax of callback functions defined by DiffEqCallbacks.jl (https://docs.sciml.ai/DiffEqCallbacks/stable/).

source
diff --git a/dev/saving_output/index.html b/dev/saving_output/index.html index df4a0aa..8aff0f9 100644 --- a/dev/saving_output/index.html +++ b/dev/saving_output/index.html @@ -96,4 +96,4 @@ ├─ 🏷️ Center ├─ 🏷️ DataType └─ 🔢 Mn -julia> close(hf5) # close the HDF5 file

We can see that the HDF5 file contains some general information about the model as attributes (🏷️), such as characteristic values or the dimensions and the total time of the model. In addition, the three timesteps are stored as groups (📂) with the composition in Ca, Fe, Mg and Mn.

+julia> close(hf5) # close the HDF5 file

We can see that the HDF5 file contains some general information about the model as attributes (🏷️), such as characteristic values or the dimensions and the total time of the model. In addition, the three timesteps are stored as groups (📂) with the composition in Ca, Fe, Mg and Mn.

diff --git a/dev/search_index.js b/dev/search_index.js index 176c887..8492a8a 100644 --- a/dev/search_index.js +++ b/dev/search_index.js @@ -1,3 +1,3 @@ var documenterSearchIndex = {"docs": -[{"location":"diffusion_2D/#2D_diffusion_CPU","page":"Diffusion in 2D Cartesian coordinates on CPU","title":"Diffusion in 2D Cartesian coordinates on CPU","text":"","category":"section"},{"location":"diffusion_2D/","page":"Diffusion in 2D Cartesian coordinates on CPU","title":"Diffusion in 2D Cartesian coordinates on CPU","text":"From chemical maps (e.g. obtained by microprobe or SEM), we can acquire 2D composition of garnet. Assuming an initial composition, we can diffuse a garnet grain major element composition.","category":"page"},{"location":"diffusion_2D/","page":"Diffusion in 2D Cartesian coordinates on CPU","title":"Diffusion in 2D Cartesian coordinates on CPU","text":"When dealing with 2D or 3D data, the question of performance comes into play. It is important to write performant code that scales well. To help us with this, DiffusionGarnet internally uses the package ParallelStencil.jl, which allows to write high-level code for parallel high-performance stencil computation on both GPUs and CPUs. We will focus on CPUs in this tutorial, as they are easier to set up, but significant performance gains can be expected on GPUs.","category":"page"},{"location":"diffusion_2D/","page":"Diffusion in 2D Cartesian coordinates on CPU","title":"Diffusion in 2D Cartesian coordinates on CPU","text":"To get the most out of this approach, it is important to start Julia with multiple threads as multithreading is the approach used to parallelise the code. See the official Julia documentation for more information.","category":"page"},{"location":"diffusion_2D/","page":"Diffusion in 2D Cartesian coordinates on CPU","title":"Diffusion in 2D Cartesian coordinates on CPU","text":"To achieve this, the simpler is to start Julia with:","category":"page"},{"location":"diffusion_2D/","page":"Diffusion in 2D Cartesian coordinates on CPU","title":"Diffusion in 2D Cartesian coordinates on CPU","text":"$ julia --threads auto","category":"page"},{"location":"diffusion_2D/","page":"Diffusion in 2D Cartesian coordinates on CPU","title":"Diffusion in 2D Cartesian coordinates on CPU","text":"which should select the maximum number of threads available on your machine.","category":"page"},{"location":"diffusion_2D/","page":"Diffusion in 2D Cartesian coordinates on CPU","title":"Diffusion in 2D Cartesian coordinates on CPU","text":"The check that, we can use:","category":"page"},{"location":"diffusion_2D/","page":"Diffusion in 2D Cartesian coordinates on CPU","title":"Diffusion in 2D Cartesian coordinates on CPU","text":"julia> println(\"Number of threads: $(Base.Threads.nthreads())\")\n16","category":"page"},{"location":"diffusion_2D/","page":"Diffusion in 2D Cartesian coordinates on CPU","title":"Diffusion in 2D Cartesian coordinates on CPU","text":"which outputs the current number of threads used in the Julia session (here 16).","category":"page"},{"location":"diffusion_2D/","page":"Diffusion in 2D Cartesian coordinates on CPU","title":"Diffusion in 2D Cartesian coordinates on CPU","text":"tip: Tip\nThe degree of parallelisation that can be achieved increases with the number of threads available.","category":"page"},{"location":"diffusion_2D/","page":"Diffusion in 2D Cartesian coordinates on CPU","title":"Diffusion in 2D Cartesian coordinates on CPU","text":"For this tutorial, we will use example data from the 2D examples section, which contain four composition matrices, and the contour of the garnet grain.","category":"page"},{"location":"diffusion_2D/","page":"Diffusion in 2D Cartesian coordinates on CPU","title":"Diffusion in 2D Cartesian coordinates on CPU","text":"First, we load the data, which should be in the same folder as your current session:","category":"page"},{"location":"diffusion_2D/","page":"Diffusion in 2D Cartesian coordinates on CPU","title":"Diffusion in 2D Cartesian coordinates on CPU","text":"using DiffusionGarnet # this can take a while\nusing DelimitedFiles\nusing Plots\nusing ProgressBars\n\nprintln(\"Number of threads: $(Threads.nthreads())\")\n\n# load the data of your choice\n# (here from the text files located in https://github.com/Iddingsite/DiffusionGarnet.jl/tree/main/examples/2D,\n# place it in the same folder as where you are running the code)\n\nMg0 = DelimitedFiles.readdlm(\"Xprp.txt\", '\\t', '\\n', header=false)\nFe0 = DelimitedFiles.readdlm(\"Xalm.txt\", '\\t', '\\n', header=false)\nMn0 = DelimitedFiles.readdlm(\"Xsps.txt\", '\\t', '\\n', header=false)\nCa0 = DelimitedFiles.readdlm(\"Xgrs.txt\", '\\t', '\\n', header=false)\ngrt_boundary = DelimitedFiles.readdlm(\"contour_Grt.txt\", '\\t', '\\n', header=false)\n\n# define total length in x and y\nLx = 9000.0u\"µm\"\nLy = 9000.0u\"µm\"\n# define total time for the model\ntfinal = 1.0u\"Myr\"\n# define the pressure and temperature conditions\nT = 900u\"°C\"\nP = 0.6u\"GPa\"","category":"page"},{"location":"diffusion_2D/","page":"Diffusion in 2D Cartesian coordinates on CPU","title":"Diffusion in 2D Cartesian coordinates on CPU","text":"Plotting the initial data:","category":"page"},{"location":"diffusion_2D/","page":"Diffusion in 2D Cartesian coordinates on CPU","title":"Diffusion in 2D Cartesian coordinates on CPU","text":"distance = LinRange(0, ustrip(u\"µm\", Lx), size(Mg0,1))\n\nl = @layout [a b ; c d ]\n\np1 = heatmap(distance, distance, Mg0, label=\"Mg\", dpi=200, title=\"Mg\", clim=(0, maximum(Mg0)), ylabel= \"Distance (µm)\")\np2 = heatmap(distance, distance, Fe0, label=\"Fe\", dpi=200, title=\"Fe\", clim=(0.7, maximum(Fe0)))\np3 = heatmap(distance, distance, Mn0, label=\"Mn\", dpi=200, title=\"Mn\", clim=(0, maximum(Mn0)), xlabel= \"Distance (µm)\", ylabel= \"Distance (µm)\")\np4 = heatmap(distance, distance, Ca0, label=\"Ca\", dpi=200, title=\"Ca\", clim=(0, 0.1), xlabel= \"Distance (µm)\")\n\nplot(p1, p2, p3, p4, layout = l , plot_title=\"Initial conditions\")","category":"page"},{"location":"diffusion_2D/","page":"Diffusion in 2D Cartesian coordinates on CPU","title":"Diffusion in 2D Cartesian coordinates on CPU","text":"which outputs:","category":"page"},{"location":"diffusion_2D/","page":"Diffusion in 2D Cartesian coordinates on CPU","title":"Diffusion in 2D Cartesian coordinates on CPU","text":"(Image: Initial conditions 2D.)","category":"page"},{"location":"diffusion_2D/","page":"Diffusion in 2D Cartesian coordinates on CPU","title":"Diffusion in 2D Cartesian coordinates on CPU","text":"Like with spherical or 1D Cartesian coordinates, we need to define our InitialConditions and Domain structures, that will hold all the information of our model.","category":"page"},{"location":"diffusion_2D/","page":"Diffusion in 2D Cartesian coordinates on CPU","title":"Diffusion in 2D Cartesian coordinates on CPU","text":"IC2D = InitialConditions2D(Mg0, Fe0, Mn0, Lx, Ly, tfinal; grt_boundary = grt_boundary)\ndomain2D = Domain(IC2D, T, P)","category":"page"},{"location":"diffusion_2D/","page":"Diffusion in 2D Cartesian coordinates on CPU","title":"Diffusion in 2D Cartesian coordinates on CPU","text":"note: Note\nThe argument grt_boundary is a binary array specifiying the coordinates where the homogeneous Dirichlet boundaries are to be enforced, here on the contour of the grain. If not provided, the homogeneous Dirichlet boundaries are applied to the sides of the domain.","category":"page"},{"location":"diffusion_2D/","page":"Diffusion in 2D Cartesian coordinates on CPU","title":"Diffusion in 2D Cartesian coordinates on CPU","text":"Using the function simulate() to solve our system:","category":"page"},{"location":"diffusion_2D/","page":"Diffusion in 2D Cartesian coordinates on CPU","title":"Diffusion in 2D Cartesian coordinates on CPU","text":"# solve the problem using DifferentialEquations.jl\nsol = simulate(Domain2D, progress=true, progress_steps=1)","category":"page"},{"location":"diffusion_2D/","page":"Diffusion in 2D Cartesian coordinates on CPU","title":"Diffusion in 2D Cartesian coordinates on CPU","text":"We use here the arguments progress=true and progress_steps=1 to display a progress bar during the calculation.","category":"page"},{"location":"diffusion_2D/","page":"Diffusion in 2D Cartesian coordinates on CPU","title":"Diffusion in 2D Cartesian coordinates on CPU","text":"note: Note\nThe calculation can take some time, about 6 minutes on my machine with 16 threads.","category":"page"},{"location":"diffusion_2D/","page":"Diffusion in 2D Cartesian coordinates on CPU","title":"Diffusion in 2D Cartesian coordinates on CPU","text":"We can now plot the results:","category":"page"},{"location":"diffusion_2D/","page":"Diffusion in 2D Cartesian coordinates on CPU","title":"Diffusion in 2D Cartesian coordinates on CPU","text":"\nprintln(\"Plotting...\")\nanim = @animate for i = tqdm(LinRange(0, sol.t[end], 20))\n\n time = round(i;digits=2)\n\n l = @layout [a b ; c d ]\n Ca = 1 .- sol(i)[:,:,1] .- sol(i)[:,:,2] .- sol(i)[:,:,3]\n replace!(Ca, 1=>0)\n\n p1 = heatmap(distance, distance, sol(i)[:,:,1], label=\"Mg\", dpi=100, title=\"Mg\", clim=(0, maximum(sol(0)[:,:,1])), ylabel= \"Distance (µm)\")\n p2 = heatmap(distance, distance, sol(i)[:,:,2], label=\"Fe\", dpi=100, title=\"Fe\", clim=(0.7, maximum(sol(0)[:,:,2])))\n p3 = heatmap(distance, distance, sol(i)[:,:,3], label=\"Mn\", dpi=100, title=\"Mn\", clim=(0, maximum(sol(0)[:,:,3])), xlabel= \"Distance (µm)\", ylabel= \"Distance (µm)\")\n p4 = heatmap(distance, distance, Ca, label=\"Ca\", dpi=100, title=\"Ca\", clim=(0, 0.1), xlabel= \"Distance (µm)\")\n\n plot(p1, p2, p3, p4, layout = l , plot_title=\"Total Time = $(time) Ma, T=$(round(ustrip.(u\"°C\", T); digits=2)) °C\")\nend every 1\n\nprintln(\"...Done!\")\n\nprintln(\"Now, generating the gif...\")\ngif(anim, \"Grt_2D_900.gif\", fps = 3)\nprintln(\"...Done!\")\n","category":"page"},{"location":"diffusion_2D/","page":"Diffusion in 2D Cartesian coordinates on CPU","title":"Diffusion in 2D Cartesian coordinates on CPU","text":"Here is the resulting gif:","category":"page"},{"location":"diffusion_2D/","page":"Diffusion in 2D Cartesian coordinates on CPU","title":"Diffusion in 2D Cartesian coordinates on CPU","text":"(Image: 2D diffusion of a garnet)","category":"page"},{"location":"diffusion_2D/","page":"Diffusion in 2D Cartesian coordinates on CPU","title":"Diffusion in 2D Cartesian coordinates on CPU","text":"note: Note\n2D modelling ignores the effect of the third dimension on diffusion and is equivalent to considering the garnet grain to be cylindrical rather than spherical.","category":"page"},{"location":"diffusion_3D/#3D_diffusion_CPU","page":"Diffusion in 3D Cartesian coordinates on CPU","title":"Diffusion in 3D Cartesian coordinates on CPU","text":"","category":"section"},{"location":"diffusion_3D/","page":"Diffusion in 3D Cartesian coordinates on CPU","title":"Diffusion in 3D Cartesian coordinates on CPU","text":"note: Note\nThis tutorial is combining the knowledge acquired from the 2D modelling and the Saving output as HDF5 files tutorials. Make sure to understand them before starting this one.","category":"page"},{"location":"diffusion_3D/","page":"Diffusion in 3D Cartesian coordinates on CPU","title":"Diffusion in 3D Cartesian coordinates on CPU","text":"The 3D geometry of garnets can be obtained from computed tomography (CT) scans or other similar techniques. Combined with an initial composition, the full grain can be modelled with realistic geometry.","category":"page"},{"location":"diffusion_3D/","page":"Diffusion in 3D Cartesian coordinates on CPU","title":"Diffusion in 3D Cartesian coordinates on CPU","text":"In this tutorial, we will use fake geometry data and fake composition data with a low resolution to get a reasonable runtime. The acquisition and processing of original data is beyond the scope of this tutorial and this package, and is left to the user. The goal of this tutorial is only to illustrate the use of DiffusionGarnet.jl in 3D Cartesian coordinates and show how to visualise the results.","category":"page"},{"location":"diffusion_3D/","page":"Diffusion in 3D Cartesian coordinates on CPU","title":"Diffusion in 3D Cartesian coordinates on CPU","text":"To do so, we will use a callback function to save the results of the simulation to disk at regular intervals to be able to visualize the results using the software Paraview.","category":"page"},{"location":"diffusion_3D/","page":"Diffusion in 3D Cartesian coordinates on CPU","title":"Diffusion in 3D Cartesian coordinates on CPU","text":"As mentioned in the tutorial for 2D modelling DiffusionGarnet internally uses the package ParallelStencil.jl. Make sure to start with multiple threads to get the most out of this approach if you run the model on CPU.","category":"page"},{"location":"diffusion_3D/","page":"Diffusion in 3D Cartesian coordinates on CPU","title":"Diffusion in 3D Cartesian coordinates on CPU","text":"For this tutorial, we will use example data from the 3D examples section, which contains four composition matrices, and the contour of a fake spherical grain. The total resolution is of 128x128x128 voxels.","category":"page"},{"location":"diffusion_3D/","page":"Diffusion in 3D Cartesian coordinates on CPU","title":"Diffusion in 3D Cartesian coordinates on CPU","text":"First, we load the data, which should be in the same folder as your current session:","category":"page"},{"location":"diffusion_3D/","page":"Diffusion in 3D Cartesian coordinates on CPU","title":"Diffusion in 3D Cartesian coordinates on CPU","text":"using DiffusionGarnet # this can take a while\nusing JLD2\nusing Plots\n\nprintln(\"Number of threads: $(Threads.nthreads())\")\n# should ideally output more than 1 thread\n\n# use JLD2 to load data\nfile = jldopen(\"3D_data.jld2\", \"r\")\n@unpack Mg0, Fe0, Mn0, Ca0, grt_boundary = file\nclose(file)","category":"page"},{"location":"diffusion_3D/","page":"Diffusion in 3D Cartesian coordinates on CPU","title":"Diffusion in 3D Cartesian coordinates on CPU","text":"We will use the same conditions as in the 2D modelling tutorial:","category":"page"},{"location":"diffusion_3D/","page":"Diffusion in 3D Cartesian coordinates on CPU","title":"Diffusion in 3D Cartesian coordinates on CPU","text":"# define total length in x and y\nLx = 9000.0u\"µm\"\nLy = 9000.0u\"µm\"\nLz = 9000.0u\"µm\"\n# define total time for the model\ntfinal = 1.0u\"Myr\"\n# define the pressure and temperature conditions\nT = 900u\"°C\"\nP = 0.6u\"GPa\"\n\nIC3D = InitialConditions3D(Mg0, Fe0, Mn0, Lx, Ly, Lz, tfinal; grt_boundary = grt_boundary)\ndomain3D = Domain(IC3D, T, P)","category":"page"},{"location":"diffusion_3D/","page":"Diffusion in 3D Cartesian coordinates on CPU","title":"Diffusion in 3D Cartesian coordinates on CPU","text":"Now, let's define a callback function to save the results of the simulation in a HDF5 file, similar to the tutorial Saving output as HDF5 files:","category":"page"},{"location":"diffusion_3D/","page":"Diffusion in 3D Cartesian coordinates on CPU","title":"Diffusion in 3D Cartesian coordinates on CPU","text":"@unpack t_charact = domain3D # unpack characteristic time to nondimensionalise the time for the simulation\ntime_save_ad = ustrip.(u\"Myr\", time_save) ./ t_charact # convert to Myr, remove units, and convert to nondimensional time\n\nsave_data_callback = PresetTimeCallback(time_save_ad, save_data_paraview)\n\npath_save = \"Grt_3D.h5\" # choose the name and the path of the HDF5 output file (make sure to add .h5 or .hdf5 at the end)","category":"page"},{"location":"diffusion_3D/","page":"Diffusion in 3D Cartesian coordinates on CPU","title":"Diffusion in 3D Cartesian coordinates on CPU","text":"note: Note\nWe used here save_data_paraview() instead of save_data() to save the data in a format that can be read by Paraview.","category":"page"},{"location":"diffusion_3D/","page":"Diffusion in 3D Cartesian coordinates on CPU","title":"Diffusion in 3D Cartesian coordinates on CPU","text":"We can now use the function simulate() to solve our system:","category":"page"},{"location":"diffusion_3D/","page":"Diffusion in 3D Cartesian coordinates on CPU","title":"Diffusion in 3D Cartesian coordinates on CPU","text":"# solve the problem using DifferentialEquations.jl\nsol = simulate(Domain3D; callback=save_data_callback, path_save=path_save, save_everystep=false, progress=true, progress_steps=1)","category":"page"},{"location":"diffusion_3D/","page":"Diffusion in 3D Cartesian coordinates on CPU","title":"Diffusion in 3D Cartesian coordinates on CPU","text":"note: Note\nThe calculation can take some time, about 4 minutes on my machine with 16 threads.","category":"page"},{"location":"diffusion_3D/","page":"Diffusion in 3D Cartesian coordinates on CPU","title":"Diffusion in 3D Cartesian coordinates on CPU","text":"Two files should have been created in the same folder as your current session: Grt_3D.h5 and Grt_3D.xdmf. The first one is the HDF5 file containing the results of the simulation, and the second one is an XDMF file describing the HDF5 file. This last file can be opened with visualisation software programs, such as Paraview.","category":"page"},{"location":"diffusion_3D/","page":"Diffusion in 3D Cartesian coordinates on CPU","title":"Diffusion in 3D Cartesian coordinates on CPU","text":"warning: Warning\nFor any visualisation software, make sure you open the XDMF file and not the HDF5 file. For Paraview, select the XDMF Reader as the reader when you open your data. Only Paraview has been tested with this package, but other software should work as well.","category":"page"},{"location":"diffusion_spherical/#sph_diffusion","page":"Diffusion in spherical coordinates","title":"Diffusion in spherical coordinates","text":"","category":"section"},{"location":"diffusion_spherical/","page":"Diffusion in spherical coordinates","title":"Diffusion in spherical coordinates","text":"When working with 1D profiles of major elements in garnet, it is often better to consider spherical coordinates. This approximates a garnet grain as a sphere and allows the change in volume relative to the distance from the core to be taken into account.","category":"page"},{"location":"diffusion_spherical/","page":"Diffusion in spherical coordinates","title":"Diffusion in spherical coordinates","text":"In DiffusionGarnet, spherical coordinates are as easy to use as 1D Cartesian coordinates. To illustrate this, we will compare the diffusion of a core-rim profile using the 2 coordinate systems. We will follow the same procedure as in the 1D diffusion tutorial.","category":"page"},{"location":"diffusion_spherical/","page":"Diffusion in spherical coordinates","title":"Diffusion in spherical coordinates","text":"First we will load the data of the core-rim profile from the Spherical examples section, which should be in the same folder as your current session:","category":"page"},{"location":"diffusion_spherical/","page":"Diffusion in spherical coordinates","title":"Diffusion in spherical coordinates","text":"using DiffusionGarnet\nusing DelimitedFiles\n\ndata = DelimitedFiles.readdlm(\"Data_Grt_Sph.txt\", '\\t', '\\n', header=true)[1]\n\nMg0 = data[:, 4]\nFe0 = data[:, 2]\nMn0 = data[:, 3]\nCa0 = data[:, 5]\ndistance = data[:, 1]\nLx = Lr = (data[end,1] - data[1,1])u\"µm\"\ntfinal = 15u\"Myr\"","category":"page"},{"location":"diffusion_spherical/","page":"Diffusion in spherical coordinates","title":"Diffusion in spherical coordinates","text":"We can visualize our data:","category":"page"},{"location":"diffusion_spherical/","page":"Diffusion in spherical coordinates","title":"Diffusion in spherical coordinates","text":"using Plots\n\nl = @layout [a ; b]\n\np1 = plot(distance, Fe0, label=\"Fe initial\", linestyle = :dash, linewidth=1, dpi=200, title = \"Initial conditions\", legend=:outerbottomright, linecolor=1,xlabel = \"Distance (µm)\", ylabel=\"Molar fraction\")\n\np2 = plot(distance, Mg0, label=\"Mg initial\", linestyle = :dash, linewidth=1, dpi=200,legend=:outerbottomright,linecolor=2,xlabel = \"Distance (µm)\")\np2 = plot!(distance, Mn0, label=\"Mn initial\", linestyle = :dash, linewidth=1, linecolor=3)\np2 = plot!(distance, Ca0, label=\"Ca initial\", linestyle = :dash, linewidth=1, linecolor=4, ylabel=\"Molar fraction\")\n\nplot(p1, p2, layout = l)","category":"page"},{"location":"diffusion_spherical/","page":"Diffusion in spherical coordinates","title":"Diffusion in spherical coordinates","text":"which outputs:","category":"page"},{"location":"diffusion_spherical/","page":"Diffusion in spherical coordinates","title":"Diffusion in spherical coordinates","text":"(Image: Initial conditions.)","category":"page"},{"location":"diffusion_spherical/","page":"Diffusion in spherical coordinates","title":"Diffusion in spherical coordinates","text":"Then, we can define our initial conditions, both for the spherical and 1D Cartesian coordinates.","category":"page"},{"location":"diffusion_spherical/","page":"Diffusion in spherical coordinates","title":"Diffusion in spherical coordinates","text":"ICSph = InitialConditionsSpherical(Mg0, Fe0, Mn0, Lr, tfinal)\nIC1D = InitialConditions1D(Mg0, Fe0, Mn0, Lx, tfinal)\n\n# define the pressure and temperature conditions of diffusion\nT = 900u\"°C\"\nP = 0.6u\"GPa\"\n\nDomainSph = Domain(ICSph, T, P)\nDomain1D = Domain(IC1D, T, P; bc_neumann = (true, false))","category":"page"},{"location":"diffusion_spherical/","page":"Diffusion in spherical coordinates","title":"Diffusion in spherical coordinates","text":"tip: Tip\nInitialConditionsSpherical always assumes that the core of the profile is on the left and the edge is on the right side.","category":"page"},{"location":"diffusion_spherical/","page":"Diffusion in spherical coordinates","title":"Diffusion in spherical coordinates","text":"note: Note\nWe use bc_neumann in Domain to indicate that a homogeneous Neumann boundary must be ensured on the left side of the 1D profile to simulate the core of the garnet grain.","category":"page"},{"location":"diffusion_spherical/","page":"Diffusion in spherical coordinates","title":"Diffusion in spherical coordinates","text":"We can then use the function simulate() to solve our system:","category":"page"},{"location":"diffusion_spherical/","page":"Diffusion in spherical coordinates","title":"Diffusion in spherical coordinates","text":"# solve the problem using DifferentialEquations.jl\nsol_sph = simulate(DomainSph)\nsol_1D = simulate(Domain1D)","category":"page"},{"location":"diffusion_spherical/","page":"Diffusion in spherical coordinates","title":"Diffusion in spherical coordinates","text":"and compare the 2 solutions:","category":"page"},{"location":"diffusion_spherical/","page":"Diffusion in spherical coordinates","title":"Diffusion in spherical coordinates","text":"anim = @animate for i = LinRange(0, sol_sph.t[end], 100)\n l = @layout [a ; b]\n\n p1 = plot(distance, Fe0, label=\"Fe initial\", linestyle = :dash, linewidth=1, dpi=200, title = \"Total Time = $(round(i;digits=2)) Ma\", legend=:outerbottomright, linecolor=1,xlabel = \"Distance (µm)\")\n p1 = plot!(distance, sol_sph(i)[:,2], label=\"Fe Sph\",linecolor=1, linewidth=1)\n p1 = plot!(distance, sol_1D(i)[:,2], label=\"Fe 1D\",linecolor=5, linewidth=1)\n\n\n p2 = plot(distance, Mg0, label=\"Mg initial\", linestyle = :dash, linewidth=1, dpi=200,legend=:outerbottomright,linecolor=2,xlabel = \"Distance (µm)\")\n p2 = plot!(distance, Mn0, label=\"Mn initial\", linestyle = :dash, linewidth=1, linecolor=3)\n p2 = plot!(distance, Ca0, label=\"Ca initial\", linestyle = :dash, linewidth=1, linecolor=4)\n p2 = plot!(distance, sol_sph(i)[:,1], label=\"Mg Sph\",linecolor=2, linewidth=1)\n p2 = plot!(distance, sol_1D(i)[:,1], label=\"Mg 1D\",linecolor=6, linewidth=1)\n\n p2 = plot!(distance, sol_sph(i)[:,3], label=\"Mn Sph\", linecolor=3, linewidth=1)\n p2 = plot!(distance, sol_1D(i)[:,3], label=\"Mn 1D\", linecolor=7, linewidth=1)\n\n p2 = plot!(distance, 1 .- sol_sph(i)[:,1] .- sol_sph(i)[:,2] .- sol_sph(i)[:,3], label=\"Ca Sph\", linecolor=4, linewidth=1)\n p2 = plot!(distance, 1 .- sol_sph(i)[:,1] .- sol_sph(i)[:,2] .- sol_sph(i)[:,3], label=\"Ca 1D\", linecolor=8, linewidth=1)\n\n plot(p1, p2, layout = l)\nend every 1\n\nprintln(\"Now, generating the gif...\")\ngif(anim, \"Grt_Spherical+1D.gif\", fps = 7)\nprintln(\"...Done!\")","category":"page"},{"location":"diffusion_spherical/","page":"Diffusion in spherical coordinates","title":"Diffusion in spherical coordinates","text":"With the resulting gif:","category":"page"},{"location":"diffusion_spherical/","page":"Diffusion in spherical coordinates","title":"Diffusion in spherical coordinates","text":"(Image: Spherical diffusion profil of a garnet)","category":"page"},{"location":"diffusion_spherical/","page":"Diffusion in spherical coordinates","title":"Diffusion in spherical coordinates","text":"It shows that using spherical coordinates makes the garnet diffuse faster. Using 1D Cartesian coordinates can overestimate the equilibration time.","category":"page"},{"location":"saving_output/#saving_output","page":"Saving output as HDF5 files","title":"Saving output as HDF5 files","text":"","category":"section"},{"location":"saving_output/","page":"Saving output as HDF5 files","title":"Saving output as HDF5 files","text":"When dealing with 2D or especially 3D data, it can be impractical or even impossible to keep every timestep in memory, as the RAM on a single machine can quickly become saturated. It may then be relevant to save certain timesteps of interest to disk for later post-processing. DiffusionGarnet has a built-in function for this purpose, using a callback function based on the DiffEqCallbacks package to produce HDF5 files.","category":"page"},{"location":"saving_output/","page":"Saving output as HDF5 files","title":"Saving output as HDF5 files","text":"HDF5 is a data format designed to store and organise large amount of data and is the data format chosen for DiffusionGarnet.","category":"page"},{"location":"saving_output/","page":"Saving output as HDF5 files","title":"Saving output as HDF5 files","text":"As an example, we will use the data from the Diffusion in 2D Cartesian coordinates on CPU tutorial but this is also applicable to 1D Cartesian or spherical coordinates in the same manner.","category":"page"},{"location":"saving_output/","page":"Saving output as HDF5 files","title":"Saving output as HDF5 files","text":"tip: Tip\nMake sure to start Julia with multiple threads of execution to speed up the simulation.","category":"page"},{"location":"saving_output/","page":"Saving output as HDF5 files","title":"Saving output as HDF5 files","text":"using DiffusionGarnet\nusing DelimitedFiles\nusing Plots\n\nMg0 = DelimitedFiles.readdlm(\"Xprp.txt\", '\\t', '\\n', header=false)\nFe0 = DelimitedFiles.readdlm(\"Xalm.txt\", '\\t', '\\n', header=false)\nMn0 = DelimitedFiles.readdlm(\"Xsps.txt\", '\\t', '\\n', header=false)\nCa0 = DelimitedFiles.readdlm(\"Xgrs.txt\", '\\t', '\\n', header=false)\ngrt_boundary = DelimitedFiles.readdlm(\"contour_Grt.txt\", '\\t', '\\n', header=false)\n\nLx = 9000.0u\"µm\"\nLy = 9000.0u\"µm\"\ntfinal = 1.0u\"Myr\"\nT = 900u\"°C\"\nP = 0.6u\"GPa\"\n\nIC2D = InitialConditions2D(Mg0, Fe0, Mn0, Lx, Ly, tfinal; grt_boundary = grt_boundary)\ndomain2D = Domain(IC2D, T, P)","category":"page"},{"location":"saving_output/","page":"Saving output as HDF5 files","title":"Saving output as HDF5 files","text":"We then need to decide at which timesteps we want to save our data, here at 0, 0.5 and 1 Myr:","category":"page"},{"location":"saving_output/","page":"Saving output as HDF5 files","title":"Saving output as HDF5 files","text":"time_save = [0, 0.5, 1]u\"Myr\" # define times at which to save","category":"page"},{"location":"saving_output/","page":"Saving output as HDF5 files","title":"Saving output as HDF5 files","text":"note: Note\nMake sure to always include 0 in your timesteps to save the initial conditions.","category":"page"},{"location":"saving_output/","page":"Saving output as HDF5 files","title":"Saving output as HDF5 files","text":"Two additional steps are then required, we need to convert this time to dimensionless time so that the solver can use it, and we need to define our callback function:","category":"page"},{"location":"saving_output/","page":"Saving output as HDF5 files","title":"Saving output as HDF5 files","text":"@unpack t_charact = domain2D # unpack characteristic time to nondimensionalise the time for the simulation\ntime_save_ad = ustrip.(u\"Myr\", time_save) ./ t_charact # convert to Myr, remove units, and convert to nondimensional time","category":"page"},{"location":"saving_output/","page":"Saving output as HDF5 files","title":"Saving output as HDF5 files","text":"time_save_ad is the equivalent dimensionless time to time_save:","category":"page"},{"location":"saving_output/","page":"Saving output as HDF5 files","title":"Saving output as HDF5 files","text":"julia > time_save_ad\n3-element Vector{Float64}:\n 0.0\n 0.011476148087304199\n 0.022952296174608398","category":"page"},{"location":"saving_output/","page":"Saving output as HDF5 files","title":"Saving output as HDF5 files","text":"warning: Warning\nThe time provided to the callback function should always be dimensionless.","category":"page"},{"location":"saving_output/","page":"Saving output as HDF5 files","title":"Saving output as HDF5 files","text":"We can now define our callback function that will save the data at these timesteps, using PresetTimeCallback():","category":"page"},{"location":"saving_output/","page":"Saving output as HDF5 files","title":"Saving output as HDF5 files","text":"save_data_callback = PresetTimeCallback(time_save_ad, save_data)","category":"page"},{"location":"saving_output/","page":"Saving output as HDF5 files","title":"Saving output as HDF5 files","text":"save_data is the function that will be called at each timestep in our callback function. save_data_callback can then be passed to simulate after defining the path of the output file:","category":"page"},{"location":"saving_output/","page":"Saving output as HDF5 files","title":"Saving output as HDF5 files","text":"path_save = \"Grt_2D.h5\" # chose the name and the path of the HDF5 output file (make sure to add .h5 or .hdf5 at the end)\nsol = simulate(domain2D; callback=save_data_callback, path_save=path_save, save_everystep=false);","category":"page"},{"location":"saving_output/","page":"Saving output as HDF5 files","title":"Saving output as HDF5 files","text":"note: Note\nWe use save_everystep=false in simulation() to prevent saving every timestep in the model as we save it to disk. That should speed up the computation and reduce the burden on the RAM.","category":"page"},{"location":"saving_output/","page":"Saving output as HDF5 files","title":"Saving output as HDF5 files","text":"This output:","category":"page"},{"location":"saving_output/","page":"Saving output as HDF5 files","title":"Saving output as HDF5 files","text":"Data saved at 0.0 Myr.\nData saved at 0.5 Myr.\nData saved at 1.0 Myr.\n352.464431 seconds (17.47 M allocations: 1.789 GiB, 0.15% gc time, 4.01% compilation time)","category":"page"},{"location":"saving_output/","page":"Saving output as HDF5 files","title":"Saving output as HDF5 files","text":"which indicates the time at which we saved our data and the total run time of the simulation.","category":"page"},{"location":"saving_output/","page":"Saving output as HDF5 files","title":"Saving output as HDF5 files","text":"The output file produced can be read with the official viewer or using HDF5.jl in Julia or other programming languages and external programs.","category":"page"},{"location":"saving_output/","page":"Saving output as HDF5 files","title":"Saving output as HDF5 files","text":"Using HDF5.jl, we can look at the structure of the HDF5 file produced:","category":"page"},{"location":"saving_output/","page":"Saving output as HDF5 files","title":"Saving output as HDF5 files","text":"julia> using HDF5\njulia> hf5 = h5open(\"Grt_2D.h5\")\n🗂️ HDF5.File: (read-only) Grt_2D.h5\n└─ 📂 Diffusion_Grt\n ├─ 🏷️ CharacteristicDiffusionCoefficient\n ├─ 🏷️ CharacteristicLength\n ├─ 🏷️ CharacteristicTime\n ├─ 🏷️ Coordinates\n ├─ 🏷️ Dx(µm)\n ├─ 🏷️ Dy(µm)\n ├─ 🏷️ LengthX(µm)\n ├─ 🏷️ LengthY(µm)\n ├─ 🏷️ Nx\n ├─ 🏷️ Ny\n ├─ 🏷️ TotalTime(Myr)\n ├─ 📂 t0000\n │ ├─ 🏷️ Time(Myr)\n │ ├─ 📂 Ca\n │ │ ├─ 🏷️ Center\n │ │ ├─ 🏷️ DataType\n │ │ └─ 🔢 Ca\n │ ├─ 📂 Fe\n │ │ ├─ 🏷️ Center\n │ │ ├─ 🏷️ DataType\n │ │ └─ 🔢 Fe\n │ ├─ 📂 Mg\n │ │ ├─ 🏷️ Center\n │ │ ├─ 🏷️ DataType\n │ │ └─ 🔢 Mg\n │ └─ 📂 Mn\n │ ├─ 🏷️ Center\n │ ├─ 🏷️ DataType\n │ └─ 🔢 Mn\n ├─ 📂 t0001\n │ ├─ 🏷️ CurrentDt(Myr)\n │ ├─ 🏷️ Time(Myr)\n │ ├─ 📂 Ca\n │ │ ├─ 🏷️ Center\n │ │ ├─ 🏷️ DataType\n │ │ └─ 🔢 Ca\n │ ├─ 📂 Fe\n │ │ ├─ 🏷️ Center\n │ │ ├─ 🏷️ DataType\n │ │ └─ 🔢 Fe\n │ ├─ 📂 Mg\n │ │ ├─ 🏷️ Center\n │ │ ├─ 🏷️ DataType\n │ │ └─ 🔢 Mg\n │ └─ 📂 Mn\n │ ├─ 🏷️ Center\n │ ├─ 🏷️ DataType\n │ └─ 🔢 Mn\n └─ 📂 t0002\n ├─ 🏷️ CurrentDt(Myr)\n ├─ 🏷️ Time(Myr)\n ├─ 📂 Ca\n │ ├─ 🏷️ Center\n │ ├─ 🏷️ DataType\n │ └─ 🔢 Ca\n ├─ 📂 Fe\n │ ├─ 🏷️ Center\n │ ├─ 🏷️ DataType\n │ └─ 🔢 Fe\n ├─ 📂 Mg\n │ ├─ 🏷️ Center\n │ ├─ 🏷️ DataType\n │ └─ 🔢 Mg\n └─ 📂 Mn\n ├─ 🏷️ Center\n ├─ 🏷️ DataType\n └─ 🔢 Mn\njulia> close(hf5) # close the HDF5 file","category":"page"},{"location":"saving_output/","page":"Saving output as HDF5 files","title":"Saving output as HDF5 files","text":"We can see that the HDF5 file contains some general information about the model as attributes (🏷️), such as characteristic values or the dimensions and the total time of the model. In addition, the three timesteps are stored as groups (📂) with the composition in Ca, Fe, Mg and Mn.","category":"page"},{"location":"updating_PT/#PT_callback","page":"Updating pressure and temperature conditions","title":"Updating pressure and temperature conditions","text":"","category":"section"},{"location":"updating_PT/","page":"Updating pressure and temperature conditions","title":"Updating pressure and temperature conditions","text":"When modelling diffusion in garnet, it may be relevant to update the pressure and temperature (PT) conditions during the simulation, i.e. to model a decompression event in the retrograde path.","category":"page"},{"location":"updating_PT/","page":"Updating pressure and temperature conditions","title":"Updating pressure and temperature conditions","text":"To do this, we can use a callback function. A callback function is a function that can be called in our solver when a certain condition is met, i.e. when a certain time is reached. It is based on the DiffEqCallbacks package from the DifferentialEquations.jl ecosystem.","category":"page"},{"location":"updating_PT/","page":"Updating pressure and temperature conditions","title":"Updating pressure and temperature conditions","text":"For this tutorial, we will use the same data from the tutorial in 1D Cartesian coordinates:","category":"page"},{"location":"updating_PT/","page":"Updating pressure and temperature conditions","title":"Updating pressure and temperature conditions","text":"using DiffusionGarnet # this can take a while\nusing DelimitedFiles\n# load the data of your choice (here from the text file located in https://github.com/Iddingsite/DiffusionGarnet.jl/tree/main/examples/1D, place it in the same folder as where you are running the code)\ndata = DelimitedFiles.readdlm(\"Data_Grt_1D.txt\", '\\t', '\\n', header=true)[1]\n\nMg0 = data[:, 4] # load initial Mg mole fraction\nFe0 = data[:, 2] # load initial Fe mole fraction\nMn0 = data[:, 3] # load initial Mn mole fraction\nCa0 = data[:, 5] # load initial Ca mole fraction\ndistance = data[:, 1]\n\nLx = (data[end,1] - data[1,1])u\"µm\" # length in x of the model, here in µm\ntfinal = 15u\"Myr\" # total time of the model, here in Myr\n\n# define the initial conditions in 1D of your problem in that order.\nIC1D = InitialConditions1D(Mg0, Fe0, Mn0, Lx, tfinal)","category":"page"},{"location":"updating_PT/","page":"Updating pressure and temperature conditions","title":"Updating pressure and temperature conditions","text":"We will now define the PT evolution of our model. To do this, we define 3 arrays containing the pressure, temperature and the time at which the conditions are updated:","category":"page"},{"location":"updating_PT/","page":"Updating pressure and temperature conditions","title":"Updating pressure and temperature conditions","text":"T = [900, 850, 800, 750, 700, 650, 600, 550]u\"°C\"\nP = [0.6, 0.5, 0.4, 0.3, 0.3, 0.3, 0.3, 0.3]u\"GPa\"\ntime_update = [0, 2, 4, 6, 8, 10, 12, 14]u\"Myr\"","category":"page"},{"location":"updating_PT/","page":"Updating pressure and temperature conditions","title":"Updating pressure and temperature conditions","text":"note: Note\nAll three arrays must have the same size and with units. Make sure to also define the conditions at time 0.","category":"page"},{"location":"updating_PT/","page":"Updating pressure and temperature conditions","title":"Updating pressure and temperature conditions","text":"This can be supplied to Domain:","category":"page"},{"location":"updating_PT/","page":"Updating pressure and temperature conditions","title":"Updating pressure and temperature conditions","text":"domain1D = Domain(IC1D, T, P, time_update)","category":"page"},{"location":"updating_PT/","page":"Updating pressure and temperature conditions","title":"Updating pressure and temperature conditions","text":"Now that we have the time and PT conditions to update, we need to create our callback function. We will use the PresetTimeCallback() function which takes as input the times at which the callback is to be triggered and the function to be called at that time.","category":"page"},{"location":"updating_PT/","page":"Updating pressure and temperature conditions","title":"Updating pressure and temperature conditions","text":"# extract time_update_ad from domain1D\n@unpack time_update_ad = domain1D\n\n# define our callback function\nupdate_diffusion_coef_call = PresetTimeCallback(time_update_ad, update_diffusion_coef)","category":"page"},{"location":"updating_PT/","page":"Updating pressure and temperature conditions","title":"Updating pressure and temperature conditions","text":"note: Note\ntime_update_ad is the equivalent of time_update in nondimensional time. update_diffusion_coef is a function that updates the diffusion coefficients according to the new PT conditions.","category":"page"},{"location":"updating_PT/","page":"Updating pressure and temperature conditions","title":"Updating pressure and temperature conditions","text":"We can now provide our callback function to simulate and run our model:","category":"page"},{"location":"updating_PT/","page":"Updating pressure and temperature conditions","title":"Updating pressure and temperature conditions","text":"sol = simulate(domain1D; callback=update_diffusion_coef_call);","category":"page"},{"location":"updating_PT/","page":"Updating pressure and temperature conditions","title":"Updating pressure and temperature conditions","text":"This outputs the time at which the callback function was called and the new PT conditions, and the total solver runtime:","category":"page"},{"location":"updating_PT/","page":"Updating pressure and temperature conditions","title":"Updating pressure and temperature conditions","text":"New temperature and pressure: 850.0 °C and 5.0 kbar, updated at 2.0 Myr.\nNew temperature and pressure: 800.0 °C and 4.0 kbar, updated at 4.0 Myr.\nNew temperature and pressure: 750.0 °C and 3.0 kbar, updated at 6.0 Myr.\nNew temperature and pressure: 700.0 °C and 3.0 kbar, updated at 8.0 Myr.\nNew temperature and pressure: 650.0 °C and 3.0 kbar, updated at 10.0 Myr.\nNew temperature and pressure: 600.0 °C and 3.0 kbar, updated at 12.0 Myr.\nNew temperature and pressure: 550.0 °C and 3.0 kbar, updated at 14.0 Myr.\n 0.752661 seconds (44.65 k allocations: 25.676 MiB)","category":"page"},{"location":"updating_PT/","page":"Updating pressure and temperature conditions","title":"Updating pressure and temperature conditions","text":"We can now plot our results and see how the change in PT has affected our diffusion profiles:","category":"page"},{"location":"updating_PT/","page":"Updating pressure and temperature conditions","title":"Updating pressure and temperature conditions","text":"@unpack tfinal_ad, t_charact = domain1D\n\nanim = @animate for i = LinRange(0, tfinal_ad, 100)\n l = @layout [a ; b]\n\n p1 = plot(distance, Fe0, label=\"Fe initial\", linestyle = :dash, linewidth=1, dpi=200, title = \"Total Time = $(round(((i)* t_charact);digits=2)) Ma\", legend=:outerbottomright, linecolor=1,xlabel = \"Distance (µm)\")\n p1 = plot!(distance, sol(i)[:,2], label=\"Fe\",linecolor=1, linewidth=1)\n\n\n p2 = plot(distance, Mg0, label=\"Mg initial\", linestyle = :dash, linewidth=1, dpi=200,legend=:outerbottomright,linecolor=2,xlabel = \"Distance (µm)\")\n p2 = plot!(distance, Mn0, label=\"Mn initial\", linestyle = :dash, linewidth=1, linecolor=3)\n p2 = plot!(distance, Ca0, label=\"Ca initial\", linestyle = :dash, linewidth=1, linecolor=4)\n p2 = plot!(distance, sol(i)[:,1], label=\"Mg\",linecolor=2, linewidth=1)\n\n p2 = plot!(distance, sol(i)[:,3], label=\"Mn\", linecolor=3, linewidth=1)\n\n p2 = plot!(distance, 1 .- sol(i)[:,1] .- sol(i)[:,2] .- sol(i)[:,3], label=\"Ca\", linecolor=4, linewidth=1)\n\n plot(p1, p2, layout = l)\nend every 1\n\nprintln(\"Now, generating the gif...\")\ngif(anim, \"Grt_1D_var_TP.gif\", fps = 7)\nprintln(\"...Done!\")\n","category":"page"},{"location":"updating_PT/","page":"Updating pressure and temperature conditions","title":"Updating pressure and temperature conditions","text":"which outputs:","category":"page"},{"location":"updating_PT/","page":"Updating pressure and temperature conditions","title":"Updating pressure and temperature conditions","text":"(Image: 1D diffusion profil of a garnet with PT)","category":"page"},{"location":"updating_PT/","page":"Updating pressure and temperature conditions","title":"Updating pressure and temperature conditions","text":"We can clearly observe the slowing down of the diffusion processes as the PT conditions decrease with time.","category":"page"},{"location":"background/#background","page":"Background","title":"Background","text":"","category":"section"},{"location":"background/#Theory","page":"Background","title":"Theory","text":"","category":"section"},{"location":"background/","page":"Background","title":"Background","text":"Garnet is a mineral commonly used in metamorphic petrology to better understand geological processes, as it occurs in a variety of rock types. This mineral often exhibits a wide range of compositional zoning, which can be interpreted as recording ranges of pressure and temperature (PT) conditions. Modelling diffusion processes can help to better understand this zoning and better constrain the pressure-temperature-time (PTt) conditions of the metamorphic event of interest.","category":"page"},{"location":"background/","page":"Background","title":"Background","text":"Major element diffusion in garnet is described by a coupled multicomponent system between four poles: Mg, Fe, Mn and Ca. This can be expressed by a system of parabolic partial differential equations (PDEs) following Fick's first law of diffusion:","category":"page"},{"location":"background/","page":"Background","title":"Background","text":"beginequation\n labelficks_law\nbeginaligned\n fracpartial C_ipartial t = mathbfnabla cdot sum_j=1^n-1 textbfD_ij mathbfnabla C_j\nendaligned\nendequation","category":"page"},{"location":"background/","page":"Background","title":"Background","text":"with C_i the molar fraction of the i^th pole of garnet, n the total number of poles, and textbfD_ij the interdiffusion coefficient tensor.","category":"page"},{"location":"background/","page":"Background","title":"Background","text":"For a system of n components, only n-1 equations are needed, as the sum of the molar fractions of all components is equal to 1. Thus, Ca is made dependent on the other concentrations and (1) is made up of three coupled PDEs.","category":"page"},{"location":"background/","page":"Background","title":"Background","text":"According to Lasaga (1979) [1], textbfD_ij can be defined, assuming an ideal system for garnet, with:","category":"page"},{"location":"background/","page":"Background","title":"Background","text":"beginequation\n labeldiffusion_coefficient\nbeginaligned\n textbfD_ij = D_i^* delta_ij - fracC_i z_i z_j D_i^*sum_k=1^n C_k z_k^2 D_k^* (D_j^* - D_n^*) =\nbeginbmatrix\n D_MgMg D_MgFe D_MgMn \n D_FeMg D_FeFe D_FeMn \n D_MnMg D_MnFe D_MnMn \nendbmatrix\nendaligned\nendequation","category":"page"},{"location":"background/","page":"Background","title":"Background","text":"with D_i^* the tracer diffusion coefficient of the i^th component, delta_ij the Kronecker delta, z the charge of the species, and D_n^* the tracer diffusion coefficient of the dependent molar fraction (here Ca).","category":"page"},{"location":"background/","page":"Background","title":"Background","text":"The tracer diffusion coefficient D_i^* can be calculated from an Arrhenius relationship:","category":"page"},{"location":"background/","page":"Background","title":"Background","text":"beginequation\n labelarrhenius\nD_i^* = D_0i exp left( -fracE_ai - (P-1)Delta V^+_iRT right)\nendequation","category":"page"},{"location":"background/","page":"Background","title":"Background","text":"with D_0i the pre-exponential constant, E_ai the activation energy of diffusion, Delta V^+_i the activation volume of diffusion at 1 bar, R the universal gas constant, T the temperature, and P the pressure.","category":"page"},{"location":"background/","page":"Background","title":"Background","text":"In DiffusionGarnet.jl, D_0i, E_ai, and Delta V^+_i are those of Chakraborty & Ganguly (1992) [2]. The tracer diffusion coefficient of Ca is defined as 05D_Fe, following the approach of Loomis et al. (1985) [3].","category":"page"},{"location":"background/#Numerical-approach","page":"Background","title":"Numerical approach","text":"","category":"section"},{"location":"background/","page":"Background","title":"Background","text":"By defining the PT conditions of the metamorphic event of interest, (3) can be solved for each component, and the diffusion coefficient tensor can be calculated using (2) from the initial major composition data. In DiffusionGarnet.jl, (1) is then discretised in space using finite differences, and the resulting system of ordinary differential equations is solved with ROCK2, a stabilised explicit method (Abdulle & Medovikov, 2001 [4]) using the DifferentialEquations.jl ecosystem.","category":"page"},{"location":"background/#refs","page":"Background","title":"References","text":"","category":"section"},{"location":"background/","page":"Background","title":"Background","text":"[1] Lasaga, A. C. (1979). Multicomponent exchange and diffusion in silicates. Geochimica et Cosmochimica Acta, 43(4), 455-469.","category":"page"},{"location":"background/","page":"Background","title":"Background","text":"[2] Chakraborty, S., & Ganguly, J. (1992). Cation diffusion in aluminosilicate garnets: experimental determination in spessartine-almandine diffusion couples, evaluation of effective binary diffusion coefficients, and applications. Contributions to Mineralogy and petrology, 111(1), 74-86.","category":"page"},{"location":"background/","page":"Background","title":"Background","text":"[3] Loomis, T. P., Ganguly, J., Elphick, S. C., 1985. Experimental determinations of cation diffusitivities in aluminosilicate garnets. II. Multicomponent simulation and tracer diffusion coefficients. Contributions to Mineralogy and Petrology 90, 45–51.","category":"page"},{"location":"background/","page":"Background","title":"Background","text":"[4] Abdulle, A., & Medovikov, A. A. (2001). Second order Chebyshev methods based on orthogonal polynomials. Numerische Mathematik, 90, 1-18.","category":"page"},{"location":"reference/","page":"List of functions","title":"List of functions","text":"","category":"page"},{"location":"reference/","page":"List of functions","title":"List of functions","text":"Modules = [DiffusionGarnet]\n","category":"page"},{"location":"reference/#DiffusionGarnet.Domain","page":"List of functions","title":"DiffusionGarnet.Domain","text":"Domain(IC::InitialConditions3D, T::Union{Unitful.Temperature,Array{<:Unitful.Temperature{<:Real}, 1}}, P::Union{Unitful.Pressure,Array{<:Unitful.Pressure{<:Real}, 1}}, time_update::Union{Unitful.Time,Array{<:Unitful.Time{<:Real}, 1}}=0u\"Myr\")\n\nWhen applied to 3D initial conditions, define corresponding 3D domain.\n\n\n\n\n\n","category":"type"},{"location":"reference/#DiffusionGarnet.Domain-2","page":"List of functions","title":"DiffusionGarnet.Domain","text":"Domain(IC::InitialConditions2D, T::Union{Unitful.Temperature,Array{<:Unitful.Temperature{<:Real}, 1}}, P::Union{Unitful.Pressure,Array{<:Unitful.Pressure{<:Real}, 1}}, time_update::Union{Unitful.Time,Array{<:Unitful.Time{<:Real}, 1}}=0u\"Myr\")\n\nWhen applied to 2D initial conditions, define corresponding 2D domain.\n\n\n\n\n\n","category":"type"},{"location":"reference/#DiffusionGarnet.Domain-3","page":"List of functions","title":"DiffusionGarnet.Domain","text":"Domain(IC::InitialConditionsSpherical, T::Union{Unitful.Temperature,Array{<:Unitful.Temperature{<:Real}, 1}}, P::Union{Unitful.Pressure,Array{<:Unitful.Pressure{<:Real}, 1}}, time_update::Union{Unitful.Time,Array{<:Unitful.Time{<:Real}, 1}}=0u\"Myr\")\n\nWhen applied to spherical initial conditions, define corresponding spherical domain. Assume that the center of the grain is on the left side.\n\n\n\n\n\n","category":"type"},{"location":"reference/#DiffusionGarnet.Domain-4","page":"List of functions","title":"DiffusionGarnet.Domain","text":"Domain(IC::InitialConditions1D, T::Union{Unitful.Temperature,Array{<:Unitful.Temperature{<:Real}, 1}}, P::Union{Unitful.Pressure,Array{<:Unitful.Pressure{<:Real}, 1}}, time_update::Union{Unitful.Time,Array{<:Unitful.Time{<:Real}, 1}}=0u\"Myr\"; bc_neumann::Tuple=(false, false))\n\nWhen applied to 1D initial conditions, define corresponding 1D domain. bc_neumann can be used to define Neumann boundary conditions on the left or right side of the domain if set to true.\n\n\n\n\n\n","category":"type"},{"location":"reference/#DiffusionGarnet.Domain-5","page":"List of functions","title":"DiffusionGarnet.Domain","text":"Domain(IC, T, P, time_update=0u\"Myr\")\n\nReturn a struct containing the struct IC, the PT conditions T and P and the nondimensionalised parameters of the problem. time_update can be used to update the PT conditions during the simulation. If so, T and P have to be of the same size as time_update and correspond to the values of PT to update to.\n\n\n\n\n\n","category":"type"},{"location":"reference/#DiffusionGarnet.InitialConditions1D-Tuple{AbstractVector{<:Real}, AbstractVector{<:Real}, AbstractVector{<:Real}, Union{Quantity{T, 𝐋, U}, Level{L, S, Quantity{T, 𝐋, U}} where {L, S}} where {T, U}, Union{Quantity{T, 𝐓, U}, Level{L, S, Quantity{T, 𝐓, U}} where {L, S}} where {T, U}}","page":"List of functions","title":"DiffusionGarnet.InitialConditions1D","text":"InitialConditions1D(CMg0::Array{<:Real, 1}, CFe0::Array{<:Real, 1}, CMn0::Array{<:Real, 1}, Lx::Unitful.Length, tfinal::Unitful.Time)\n\nReturn a structure containing the initial conditions for a 1D diffusion problem. CMg0, CFe0 and CMn0 need to be in molar fraction. Convert the Lxandtfinal`` to µm and Myr respectively.\n\n\n\n\n\n","category":"method"},{"location":"reference/#DiffusionGarnet.InitialConditions2D-Tuple{AbstractMatrix{<:Real}, AbstractMatrix{<:Real}, AbstractMatrix{<:Real}, Union{Quantity{T, 𝐋, U}, Level{L, S, Quantity{T, 𝐋, U}} where {L, S}} where {T, U}, Union{Quantity{T, 𝐋, U}, Level{L, S, Quantity{T, 𝐋, U}} where {L, S}} where {T, U}, Union{Quantity{T, 𝐓, U}, Level{L, S, Quantity{T, 𝐓, U}} where {L, S}} where {T, U}}","page":"List of functions","title":"DiffusionGarnet.InitialConditions2D","text":"InitialConditions2D(CMg0::AbstractArray{<:Real, 2}, CFe0::AbstractArray{<:Real, 2}, CMn0::AbstractArray{<:Real, 2}, Lx::Unitful.Length, Ly::Unitful.Length, tfinal::Unitful.Time; grt_boundary::AbstractArray{<:Real, 2}=zeros(eltype(CMg0), size(CMg0)...))\n\nReturn a structure containing the initial conditions for a 2D diffusion problem. CMg0, CFe0 and CMn0 need to be in molar fraction. Convert Lx, Ly and tfinal to µm, µm and Myr respectively. grt_boundary is a matrix of the same size as CMg0, CFe0 and CMn0. It is a binary matrix with 1 where the contour of the garnet is present and 0 otherwise. It is used to apply the Dirichlet boundary condition in the model. If not provided, it is equal to zero everywhere.\n\n\n\n\n\n","category":"method"},{"location":"reference/#DiffusionGarnet.InitialConditions3D-Tuple{AbstractArray{<:Real, 3}, AbstractArray{<:Real, 3}, AbstractArray{<:Real, 3}, Union{Quantity{T, 𝐋, U}, Level{L, S, Quantity{T, 𝐋, U}} where {L, S}} where {T, U}, Union{Quantity{T, 𝐋, U}, Level{L, S, Quantity{T, 𝐋, U}} where {L, S}} where {T, U}, Union{Quantity{T, 𝐋, U}, Level{L, S, Quantity{T, 𝐋, U}} where {L, S}} where {T, U}, Union{Quantity{T, 𝐓, U}, Level{L, S, Quantity{T, 𝐓, U}} where {L, S}} where {T, U}}","page":"List of functions","title":"DiffusionGarnet.InitialConditions3D","text":"InitialConditions3D(CMg0::AbstractArray{<:Real, 3}, CFe0::AbstractArray{<:Real, 3}, CMn0::AbstractArray{<:Real, 3}, Lx::Unitful.Length, Ly::Unitful.Length, Lz::Unitful.Length, tfinal::Unitful.Time; grt_boundary::Union{AbstractArray{<:Real, 3}, AbstractArray{<:Bool, 3}}=zeros(Bool, size(CMg0)...))\n\nReturn a structure containing the initial conditions for a 3D diffusion problem. CMg0, CFe0 and CMn0 need to be in molar fraction. Convert Lx, Ly, Lz and tfinal to µm, µm, µm and Myr respectively. grt_boundary is a matrix of the same size as CMg0, CFe0 and CMn0. It is a binary matrix with 1 where the contour of the garnet is present and 0 otherwise. It is used to apply the Dirichlet boundary condition in the model. If not provided, it is equal to zero everywhere.\n\n\n\n\n\n","category":"method"},{"location":"reference/#DiffusionGarnet.InitialConditionsSpherical-Tuple{AbstractVector{<:Real}, AbstractVector{<:Real}, AbstractVector{<:Real}, Union{Quantity{T, 𝐋, U}, Level{L, S, Quantity{T, 𝐋, U}} where {L, S}} where {T, U}, Union{Quantity{T, 𝐓, U}, Level{L, S, Quantity{T, 𝐓, U}} where {L, S}} where {T, U}}","page":"List of functions","title":"DiffusionGarnet.InitialConditionsSpherical","text":"InitialConditionsSpherical(CMg0::AbstractArray{<:Real, 1}, CFe0::AbstractArray{<:Real, 1}, CMn0::AbstractArray{<:Real, 1}, Lr::Unitful.Length, tfinal::Unitful.Time)\n\nReturn a structure containing the initial conditions for a spherical diffusion problem. CMg0, CFe0 and CMn0 need to be in molar fraction. Convert Lr and tfinal to µm and Myr respectively.\n\n\n\n\n\n","category":"method"},{"location":"reference/#DiffusionGarnet.save_data-Tuple{Any}","page":"List of functions","title":"DiffusionGarnet.save_data","text":"save_data(integrator)\n\nCallback function used to save major element compositions to an HDF5 file at a specific timestep.\n\n\n\n\n\n","category":"method"},{"location":"reference/#DiffusionGarnet.save_data_paraview-Tuple{Any}","page":"List of functions","title":"DiffusionGarnet.save_data_paraview","text":"save_data_paraview(integrator)\n\nCallback function used to save data to an HDF5 file and produce an XMDF file. XMDF files can be used to read data on software like Paraview.\n\nNote that data are saved in row major format. For column major, use save_data() instead.\n\n\n\n\n\n","category":"method"},{"location":"reference/#DiffusionGarnet.simulate","page":"List of functions","title":"DiffusionGarnet.simulate","text":"simulate(domain; path_save=nothing, solver=ROCK2(), kwargs...)\n\nSolve the coupled major element diffusion equations for a given domain using finite differences for the discretisation in space and return a solution type variable.\n\nThe 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.\n\nThe 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. The time of the solution is in Myr.\n\npath_save is an optional argument, which can be used to define the path of the HDF5 output file. Default is to nothing.\n\nsolver is an optional argument, which can be used to define the solver to use for the time discretisation. Default is the ROCK2 method. All other ODE solvers accepted as the ones from DifferentialEquations (see https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/).\n\nAll 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/commonsolveropts/).\n\n\n\n\n\n","category":"function"},{"location":"reference/#DiffusionGarnet.simulate-Tuple{DiffusionGarnet.Domain1D}","page":"List of functions","title":"DiffusionGarnet.simulate","text":"simulate(domain::Domain1D; path_save=nothing, solver=ROCK2(), abstol=1e-8, reltol=1e-6, kwargs...)\n\nSolve the coupled major element diffusion equations in 1D. Save all timesteps in the output solution type variable by default.\n\n\n\n\n\n","category":"method"},{"location":"reference/#DiffusionGarnet.simulate-Tuple{DiffusionGarnet.Domain2D}","page":"List of functions","title":"DiffusionGarnet.simulate","text":"simulate(domain::Domain2D; path_save=nothing, solver=ROCK2(), kwargs...)\n\nSolve 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.\n\n\n\n\n\n","category":"method"},{"location":"reference/#DiffusionGarnet.simulate-Tuple{DiffusionGarnet.Domain3D}","page":"List of functions","title":"DiffusionGarnet.simulate","text":"simulate(domain::Domain3D; path_save=nothing, solver=ROCK2(), kwargs...)\n\nSolve the coupled major element diffusion equations in 3D. Save only the first and last timestep in the output solution type variable by default.\n\n\n\n\n\n","category":"method"},{"location":"reference/#DiffusionGarnet.simulate-Tuple{DiffusionGarnet.DomainSpherical}","page":"List of functions","title":"DiffusionGarnet.simulate","text":"simulate(domain::DomainSpherical; path_save=nothing, solver=ROCK2(), abstol=1e-8, reltol=1e-6, kwargs...)\n\nSolve the coupled major element diffusion equations in spherical coordinates. Save all timesteps in the output solution type variable by default.\n\n\n\n\n\n","category":"method"},{"location":"reference/#DiffusionGarnet.update_diffusion_coef-Tuple{Any}","page":"List of functions","title":"DiffusionGarnet.update_diffusion_coef","text":"update_diffusion_coef(integrator)\n\nCallback function to update the diffusion coefficients at a given time from a new pressure and temperature. To use with the callback PresetTimeCallback (https://docs.sciml.ai/stable/basics/callbacks/#PresetTimeCallback-1).\n\nFollows the syntax of callback functions defined by DiffEqCallbacks.jl (https://docs.sciml.ai/DiffEqCallbacks/stable/).\n\n\n\n\n\n","category":"method"},{"location":"#DiffusionGarnet.jl","page":"Home","title":"DiffusionGarnet.jl","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"This is the documentation page for DiffusionGarnet.jl, a high-level Julia package to do coupled diffusion modelling of major elements on real garnet data. It currently supports 1D, spherical, 2D and 3D coordinates for evenly spaced data.","category":"page"},{"location":"","page":"Home","title":"Home","text":"It is built on top of the DifferentialEquations.jl package ecosystem and uses Unitful.jl to allow the user to define appropriate units for their problems. For 2D and 3D models, it uses ParallelStencil.jl to support multithreading on CPU and parallel computing on GPU.","category":"page"},{"location":"","page":"Home","title":"Home","text":"This package can be used as a teaching tool or for research purposes.","category":"page"},{"location":"#Installation","page":"Home","title":"Installation","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"DiffusionGarnet.jl is a registered package and may be installed directly from the REPL:","category":"page"},{"location":"","page":"Home","title":"Home","text":"julia>]\n pkg> add DiffusionGarnet\n pkg> test DiffusionGarnet","category":"page"},{"location":"diffusion_1D/#1D_diffusion","page":"Diffusion in 1D Cartesian coordinates","title":"Diffusion in 1D Cartesian coordinates","text":"","category":"section"},{"location":"diffusion_1D/","page":"Diffusion in 1D Cartesian coordinates","title":"Diffusion in 1D Cartesian coordinates","text":"DiffusionGarnet expects the user to provide real natural data for modelling major element diffusion in garnet. Note that the profiles must be evenly spaced. A set of example data can be found in the repository of the package in the 1D examples section for 1D profile called Data_grt_1D.txt. This is what we will use for this tutorial.","category":"page"},{"location":"diffusion_1D/","page":"Diffusion in 1D Cartesian coordinates","title":"Diffusion in 1D Cartesian coordinates","text":"First, we will load the data, which should be in the same folder as your current session:","category":"page"},{"location":"diffusion_1D/","page":"Diffusion in 1D Cartesian coordinates","title":"Diffusion in 1D Cartesian coordinates","text":"using DiffusionGarnet # this can take a while\nusing DelimitedFiles\n# load the data of your choice (here from the text file located in https://github.com/Iddingsite/DiffusionGarnet.jl/tree/main/examples/1D, place it in the same folder as where you are running the code)\ndata = DelimitedFiles.readdlm(\"Data_Grt_1D.txt\", '\\t', '\\n', header=true)[1]\n\nMg0 = data[:, 4] # load initial Mg mole fraction\nFe0 = data[:, 2] # load initial Fe mole fraction\nMn0 = data[:, 3] # load initial Mn mole fraction\nCa0 = data[:, 5] # load initial Ca mole fraction\ndistance = data[:, 1]","category":"page"},{"location":"diffusion_1D/","page":"Diffusion in 1D Cartesian coordinates","title":"Diffusion in 1D Cartesian coordinates","text":"We can visualize our data:","category":"page"},{"location":"diffusion_1D/","page":"Diffusion in 1D Cartesian coordinates","title":"Diffusion in 1D Cartesian coordinates","text":"using Plots\n\nl = @layout [a ; b]\n\np1 = plot(distance, Fe0, label=\"Fe initial\", linestyle = :dash, linewidth=1, dpi=200, title = \"Initial conditions\", legend=:outerbottomright, linecolor=1,xlabel = \"Distance (µm)\", ylabel=\"Molar fraction\")\n\np2 = plot(distance, Mg0, label=\"Mg initial\", linestyle = :dash, linewidth=1, dpi=200,legend=:outerbottomright,linecolor=2,xlabel = \"Distance (µm)\")\np2 = plot!(distance, Mn0, label=\"Mn initial\", linestyle = :dash, linewidth=1, linecolor=3)\np2 = plot!(distance, Ca0, label=\"Ca initial\", linestyle = :dash, linewidth=1, linecolor=4, ylabel=\"Molar fraction\")\n\nplot(p1, p2, layout = l)","category":"page"},{"location":"diffusion_1D/","page":"Diffusion in 1D Cartesian coordinates","title":"Diffusion in 1D Cartesian coordinates","text":"which outputs:","category":"page"},{"location":"diffusion_1D/","page":"Diffusion in 1D Cartesian coordinates","title":"Diffusion in 1D Cartesian coordinates","text":"(Image: Initial conditions.)","category":"page"},{"location":"diffusion_1D/","page":"Diffusion in 1D Cartesian coordinates","title":"Diffusion in 1D Cartesian coordinates","text":"Then, we will define 2 structures from the constructors InitialConditions1D and Domain, which will contain all the information DiffusionGarnet needs to run a simulation.","category":"page"},{"location":"diffusion_1D/","page":"Diffusion in 1D Cartesian coordinates","title":"Diffusion in 1D Cartesian coordinates","text":"Lx = (data[end,1] - data[1,1])u\"µm\" # length in x of the model, here in µm\ntfinal = 15u\"Myr\" # total time of the model, here in Myr\n\n# define the initial conditions in 1D of your problem in that order.\nIC1D = InitialConditions1D(Mg0, Fe0, Mn0, Lx, tfinal)\n\n# define the pressure and temperature conditions of diffusion\nT = 900u\"°C\"\nP = 0.6u\"GPa\"\n\n# define a Domain struct containing the definition of our problem and nondimensionalised variables\ndomain1D = Domain(IC1D, T, P)","category":"page"},{"location":"diffusion_1D/","page":"Diffusion in 1D Cartesian coordinates","title":"Diffusion in 1D Cartesian coordinates","text":"note: Note\nLx, tfinal, T and P need to contain units, following the syntax of the package Unitful. This allows the user to specify the units that suit their problem.","category":"page"},{"location":"diffusion_1D/","page":"Diffusion in 1D Cartesian coordinates","title":"Diffusion in 1D Cartesian coordinates","text":"Domain1D contains all the information that DiffusionGarnet needs to solve our coupled diffusion problem, at 900 °C and 0.6 GPa for a duration of 15 Myr.","category":"page"},{"location":"diffusion_1D/","page":"Diffusion in 1D Cartesian coordinates","title":"Diffusion in 1D Cartesian coordinates","text":"This can be achieved with the function simulate():","category":"page"},{"location":"diffusion_1D/","page":"Diffusion in 1D Cartesian coordinates","title":"Diffusion in 1D Cartesian coordinates","text":"# solve the problem using DifferentialEquations.jl\nsol = simulate(domain1D)","category":"page"},{"location":"diffusion_1D/","page":"Diffusion in 1D Cartesian coordinates","title":"Diffusion in 1D Cartesian coordinates","text":"which outputs the time spent on the solver, for example, on the second run:","category":"page"},{"location":"diffusion_1D/","page":"Diffusion in 1D Cartesian coordinates","title":"Diffusion in 1D Cartesian coordinates","text":" 0.399870 seconds (31.93 k allocations: 18.212 MiB)","category":"page"},{"location":"diffusion_1D/","page":"Diffusion in 1D Cartesian coordinates","title":"Diffusion in 1D Cartesian coordinates","text":"simulate() uses the DifferentialEquations.jl 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. It basically holds all the information from our simulation.","category":"page"},{"location":"diffusion_1D/","page":"Diffusion in 1D Cartesian coordinates","title":"Diffusion in 1D Cartesian coordinates","text":"We can now plot the solution to our problem. By default, DiffusionGarnet saves the time in Myr.","category":"page"},{"location":"diffusion_1D/","page":"Diffusion in 1D Cartesian coordinates","title":"Diffusion in 1D Cartesian coordinates","text":"\nanim = @animate for i = LinRange(0, sol.t[end], 100)\n l = @layout [a ; b]\n\n 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)\")\n p1 = plot!(distance, sol(i)[:,2], label=\"Fe\",linecolor=1, linewidth=1)\n\n p2 = plot(distance, Mg0, label=\"Mg initial\", linestyle = :dash, linewidth=1, dpi=200,legend=:outerbottomright,linecolor=2,xlabel = \"Distance (µm)\")\n p2 = plot!(distance, Mn0, label=\"Mn initial\", linestyle = :dash, linewidth=1, linecolor=3)\n p2 = plot!(distance, Ca0, label=\"Ca initial\", linestyle = :dash, linewidth=1, linecolor=4)\n p2 = plot!(distance, sol(i)[:,1], label=\"Mg\",linecolor=2, linewidth=1)\n\n p2 = plot!(distance, sol(i)[:,3], label=\"Mn\", linecolor=3, linewidth=1)\n\n p2 = plot!(distance, 1 .- sol(i)[:,1] .- sol(i)[:,2] .- sol(i)[:,3], label=\"Ca\", linecolor=4, linewidth=1)\n\n plot(p1, p2, layout = l)\nend every 1\n\nprintln(\"Now, generating the gif...\")\ngif(anim, \"Grt_1D.gif\", fps = 7)\nprintln(\"...Done!\")","category":"page"},{"location":"diffusion_1D/","page":"Diffusion in 1D Cartesian coordinates","title":"Diffusion in 1D Cartesian coordinates","text":"Here is the resulting gif obtained:","category":"page"},{"location":"diffusion_1D/","page":"Diffusion in 1D Cartesian coordinates","title":"Diffusion in 1D Cartesian coordinates","text":"(Image: 1D diffusion profil of a garnet)","category":"page"},{"location":"diffusion_1D/","page":"Diffusion in 1D Cartesian coordinates","title":"Diffusion in 1D Cartesian coordinates","text":"It shows the compositional evolution of a 1D profile through a garnet grain with homogeneous Dirichlet boundaries on both sides.","category":"page"}] +[{"location":"diffusion_2D/#2D_diffusion_CPU","page":"Diffusion in 2D Cartesian coordinates on CPU","title":"Diffusion in 2D Cartesian coordinates on CPU","text":"","category":"section"},{"location":"diffusion_2D/","page":"Diffusion in 2D Cartesian coordinates on CPU","title":"Diffusion in 2D Cartesian coordinates on CPU","text":"From chemical maps (e.g. obtained by microprobe or SEM), we can acquire 2D composition of garnet. Assuming an initial composition, we can diffuse a garnet grain major element composition.","category":"page"},{"location":"diffusion_2D/","page":"Diffusion in 2D Cartesian coordinates on CPU","title":"Diffusion in 2D Cartesian coordinates on CPU","text":"When dealing with 2D or 3D data, the question of performance comes into play. It is important to write performant code that scales well. To help us with this, DiffusionGarnet internally uses the package ParallelStencil.jl, which allows to write high-level code for parallel high-performance stencil computation on both GPUs and CPUs. We will focus on CPUs in this tutorial, as they are easier to set up, but significant performance gains can be expected on GPUs.","category":"page"},{"location":"diffusion_2D/","page":"Diffusion in 2D Cartesian coordinates on CPU","title":"Diffusion in 2D Cartesian coordinates on CPU","text":"To get the most out of this approach, it is important to start Julia with multiple threads as multithreading is the approach used to parallelise the code. See the official Julia documentation for more information.","category":"page"},{"location":"diffusion_2D/","page":"Diffusion in 2D Cartesian coordinates on CPU","title":"Diffusion in 2D Cartesian coordinates on CPU","text":"To achieve this, the simpler is to start Julia with:","category":"page"},{"location":"diffusion_2D/","page":"Diffusion in 2D Cartesian coordinates on CPU","title":"Diffusion in 2D Cartesian coordinates on CPU","text":"$ julia --threads auto","category":"page"},{"location":"diffusion_2D/","page":"Diffusion in 2D Cartesian coordinates on CPU","title":"Diffusion in 2D Cartesian coordinates on CPU","text":"which should select the maximum number of threads available on your machine.","category":"page"},{"location":"diffusion_2D/","page":"Diffusion in 2D Cartesian coordinates on CPU","title":"Diffusion in 2D Cartesian coordinates on CPU","text":"The check that, we can use:","category":"page"},{"location":"diffusion_2D/","page":"Diffusion in 2D Cartesian coordinates on CPU","title":"Diffusion in 2D Cartesian coordinates on CPU","text":"julia> println(\"Number of threads: $(Base.Threads.nthreads())\")\n16","category":"page"},{"location":"diffusion_2D/","page":"Diffusion in 2D Cartesian coordinates on CPU","title":"Diffusion in 2D Cartesian coordinates on CPU","text":"which outputs the current number of threads used in the Julia session (here 16).","category":"page"},{"location":"diffusion_2D/","page":"Diffusion in 2D Cartesian coordinates on CPU","title":"Diffusion in 2D Cartesian coordinates on CPU","text":"tip: Tip\nThe degree of parallelisation that can be achieved increases with the number of threads available.","category":"page"},{"location":"diffusion_2D/","page":"Diffusion in 2D Cartesian coordinates on CPU","title":"Diffusion in 2D Cartesian coordinates on CPU","text":"For this tutorial, we will use example data from the 2D examples section, which contain four composition matrices, and the contour of the garnet grain.","category":"page"},{"location":"diffusion_2D/","page":"Diffusion in 2D Cartesian coordinates on CPU","title":"Diffusion in 2D Cartesian coordinates on CPU","text":"First, we load the data, which should be in the same folder as your current session:","category":"page"},{"location":"diffusion_2D/","page":"Diffusion in 2D Cartesian coordinates on CPU","title":"Diffusion in 2D Cartesian coordinates on CPU","text":"using DiffusionGarnet # this can take a while\nusing DelimitedFiles\nusing Plots\nusing ProgressBars\n\nprintln(\"Number of threads: $(Threads.nthreads())\")\n\n# load the data of your choice\n# (here from the text files located in https://github.com/Iddingsite/DiffusionGarnet.jl/tree/main/examples/2D,\n# place it in the same folder as where you are running the code)\n\nMg0 = DelimitedFiles.readdlm(\"Xprp.txt\", '\\t', '\\n', header=false)\nFe0 = DelimitedFiles.readdlm(\"Xalm.txt\", '\\t', '\\n', header=false)\nMn0 = DelimitedFiles.readdlm(\"Xsps.txt\", '\\t', '\\n', header=false)\nCa0 = DelimitedFiles.readdlm(\"Xgrs.txt\", '\\t', '\\n', header=false)\ngrt_boundary = DelimitedFiles.readdlm(\"contour_Grt.txt\", '\\t', '\\n', header=false)\n\n# define total length in x and y\nLx = 9000.0u\"µm\"\nLy = 9000.0u\"µm\"\n# define total time for the model\ntfinal = 1.0u\"Myr\"\n# define the pressure and temperature conditions\nT = 900u\"°C\"\nP = 0.6u\"GPa\"","category":"page"},{"location":"diffusion_2D/","page":"Diffusion in 2D Cartesian coordinates on CPU","title":"Diffusion in 2D Cartesian coordinates on CPU","text":"Plotting the initial data:","category":"page"},{"location":"diffusion_2D/","page":"Diffusion in 2D Cartesian coordinates on CPU","title":"Diffusion in 2D Cartesian coordinates on CPU","text":"distance = LinRange(0, ustrip(u\"µm\", Lx), size(Mg0,1))\n\nl = @layout [a b ; c d ]\n\np1 = heatmap(distance, distance, Mg0, label=\"Mg\", dpi=200, title=\"Mg\", clim=(0, maximum(Mg0)), ylabel= \"Distance (µm)\")\np2 = heatmap(distance, distance, Fe0, label=\"Fe\", dpi=200, title=\"Fe\", clim=(0.7, maximum(Fe0)))\np3 = heatmap(distance, distance, Mn0, label=\"Mn\", dpi=200, title=\"Mn\", clim=(0, maximum(Mn0)), xlabel= \"Distance (µm)\", ylabel= \"Distance (µm)\")\np4 = heatmap(distance, distance, Ca0, label=\"Ca\", dpi=200, title=\"Ca\", clim=(0, 0.1), xlabel= \"Distance (µm)\")\n\nplot(p1, p2, p3, p4, layout = l , plot_title=\"Initial conditions\")","category":"page"},{"location":"diffusion_2D/","page":"Diffusion in 2D Cartesian coordinates on CPU","title":"Diffusion in 2D Cartesian coordinates on CPU","text":"which outputs:","category":"page"},{"location":"diffusion_2D/","page":"Diffusion in 2D Cartesian coordinates on CPU","title":"Diffusion in 2D Cartesian coordinates on CPU","text":"(Image: Initial conditions 2D.)","category":"page"},{"location":"diffusion_2D/","page":"Diffusion in 2D Cartesian coordinates on CPU","title":"Diffusion in 2D Cartesian coordinates on CPU","text":"Like with spherical or 1D Cartesian coordinates, we need to define our InitialConditions and Domain structures, that will hold all the information of our model.","category":"page"},{"location":"diffusion_2D/","page":"Diffusion in 2D Cartesian coordinates on CPU","title":"Diffusion in 2D Cartesian coordinates on CPU","text":"IC2D = InitialConditions2D(Mg0, Fe0, Mn0, Lx, Ly, tfinal; grt_boundary = grt_boundary)\ndomain2D = Domain(IC2D, T, P)","category":"page"},{"location":"diffusion_2D/","page":"Diffusion in 2D Cartesian coordinates on CPU","title":"Diffusion in 2D Cartesian coordinates on CPU","text":"note: Note\nThe argument grt_boundary is a binary array specifiying the coordinates where the homogeneous Dirichlet boundaries are to be enforced, here on the contour of the grain. If not provided, the homogeneous Dirichlet boundaries are applied to the sides of the domain.","category":"page"},{"location":"diffusion_2D/","page":"Diffusion in 2D Cartesian coordinates on CPU","title":"Diffusion in 2D Cartesian coordinates on CPU","text":"Using the function simulate() to solve our system:","category":"page"},{"location":"diffusion_2D/","page":"Diffusion in 2D Cartesian coordinates on CPU","title":"Diffusion in 2D Cartesian coordinates on CPU","text":"# solve the problem using DifferentialEquations.jl\nsol = simulate(Domain2D, progress=true, progress_steps=1)","category":"page"},{"location":"diffusion_2D/","page":"Diffusion in 2D Cartesian coordinates on CPU","title":"Diffusion in 2D Cartesian coordinates on CPU","text":"We use here the arguments progress=true and progress_steps=1 to display a progress bar during the calculation.","category":"page"},{"location":"diffusion_2D/","page":"Diffusion in 2D Cartesian coordinates on CPU","title":"Diffusion in 2D Cartesian coordinates on CPU","text":"note: Note\nThe calculation can take some time, about 6 minutes on my machine with 16 threads.","category":"page"},{"location":"diffusion_2D/","page":"Diffusion in 2D Cartesian coordinates on CPU","title":"Diffusion in 2D Cartesian coordinates on CPU","text":"We can now plot the results:","category":"page"},{"location":"diffusion_2D/","page":"Diffusion in 2D Cartesian coordinates on CPU","title":"Diffusion in 2D Cartesian coordinates on CPU","text":"@unpack t_charact = domain2D\n\nprintln(\"Plotting...\")\nanim = @animate for i = tqdm(LinRange(0, sol.t[end], 20))\n\n time = round(i*t_charact;digits=2)\n\n l = @layout [a b ; c d ]\n Ca = 1 .- sol(i)[:,:,1] .- sol(i)[:,:,2] .- sol(i)[:,:,3]\n replace!(Ca, 1=>0)\n\n p1 = heatmap(distance, distance, sol(i)[:,:,1], label=\"Mg\", dpi=100, title=\"Mg\", clim=(0, maximum(sol(0)[:,:,1])), ylabel= \"Distance (µm)\")\n p2 = heatmap(distance, distance, sol(i)[:,:,2], label=\"Fe\", dpi=100, title=\"Fe\", clim=(0.7, maximum(sol(0)[:,:,2])))\n p3 = heatmap(distance, distance, sol(i)[:,:,3], label=\"Mn\", dpi=100, title=\"Mn\", clim=(0, maximum(sol(0)[:,:,3])), xlabel= \"Distance (µm)\", ylabel= \"Distance (µm)\")\n p4 = heatmap(distance, distance, Ca, label=\"Ca\", dpi=100, title=\"Ca\", clim=(0, 0.1), xlabel= \"Distance (µm)\")\n\n plot(p1, p2, p3, p4, layout = l , plot_title=\"Total Time = $(time) Ma, T=$(round(ustrip.(u\"°C\", T); digits=2)) °C\")\nend every 1\n\nprintln(\"...Done!\")\n\nprintln(\"Now, generating the gif...\")\ngif(anim, \"Grt_2D_900.gif\", fps = 3)\nprintln(\"...Done!\")\n","category":"page"},{"location":"diffusion_2D/","page":"Diffusion in 2D Cartesian coordinates on CPU","title":"Diffusion in 2D Cartesian coordinates on CPU","text":"Here is the resulting gif:","category":"page"},{"location":"diffusion_2D/","page":"Diffusion in 2D Cartesian coordinates on CPU","title":"Diffusion in 2D Cartesian coordinates on CPU","text":"(Image: 2D diffusion of a garnet)","category":"page"},{"location":"diffusion_2D/","page":"Diffusion in 2D Cartesian coordinates on CPU","title":"Diffusion in 2D Cartesian coordinates on CPU","text":"note: Note\n2D modelling ignores the effect of the third dimension on diffusion and is equivalent to considering the garnet grain to be cylindrical rather than spherical.","category":"page"},{"location":"diffusion_3D/#3D_diffusion_CPU","page":"Diffusion in 3D Cartesian coordinates on CPU","title":"Diffusion in 3D Cartesian coordinates on CPU","text":"","category":"section"},{"location":"diffusion_3D/","page":"Diffusion in 3D Cartesian coordinates on CPU","title":"Diffusion in 3D Cartesian coordinates on CPU","text":"note: Note\nThis tutorial is combining the knowledge acquired from the 2D modelling and the Saving output as HDF5 files tutorials. Make sure to understand them before starting this one.","category":"page"},{"location":"diffusion_3D/","page":"Diffusion in 3D Cartesian coordinates on CPU","title":"Diffusion in 3D Cartesian coordinates on CPU","text":"The 3D geometry of garnets can be obtained from computed tomography (CT) scans or other similar techniques. Combined with an initial composition, the full grain can be modelled with realistic geometry.","category":"page"},{"location":"diffusion_3D/","page":"Diffusion in 3D Cartesian coordinates on CPU","title":"Diffusion in 3D Cartesian coordinates on CPU","text":"In this tutorial, we will use fake geometry data and fake composition data with a low resolution to get a reasonable runtime. The acquisition and processing of original data is beyond the scope of this tutorial and this package, and is left to the user. The goal of this tutorial is only to illustrate the use of DiffusionGarnet.jl in 3D Cartesian coordinates and show how to visualise the results.","category":"page"},{"location":"diffusion_3D/","page":"Diffusion in 3D Cartesian coordinates on CPU","title":"Diffusion in 3D Cartesian coordinates on CPU","text":"To do so, we will use a callback function to save the results of the simulation to disk at regular intervals to be able to visualize the results using the software Paraview.","category":"page"},{"location":"diffusion_3D/","page":"Diffusion in 3D Cartesian coordinates on CPU","title":"Diffusion in 3D Cartesian coordinates on CPU","text":"As mentioned in the tutorial for 2D modelling DiffusionGarnet internally uses the package ParallelStencil.jl. Make sure to start with multiple threads to get the most out of this approach if you run the model on CPU.","category":"page"},{"location":"diffusion_3D/","page":"Diffusion in 3D Cartesian coordinates on CPU","title":"Diffusion in 3D Cartesian coordinates on CPU","text":"For this tutorial, we will use example data from the 3D examples section, which contains four composition matrices, and the contour of a fake spherical grain. The total resolution is of 128x128x128 voxels.","category":"page"},{"location":"diffusion_3D/","page":"Diffusion in 3D Cartesian coordinates on CPU","title":"Diffusion in 3D Cartesian coordinates on CPU","text":"First, we load the data, which should be in the same folder as your current session:","category":"page"},{"location":"diffusion_3D/","page":"Diffusion in 3D Cartesian coordinates on CPU","title":"Diffusion in 3D Cartesian coordinates on CPU","text":"using DiffusionGarnet # this can take a while\nusing JLD2\nusing Plots\n\nprintln(\"Number of threads: $(Threads.nthreads())\")\n# should ideally output more than 1 thread\n\n# use JLD2 to load data\nfile = jldopen(\"3D_data.jld2\", \"r\")\n@unpack Mg0, Fe0, Mn0, Ca0, grt_boundary = file\nclose(file)","category":"page"},{"location":"diffusion_3D/","page":"Diffusion in 3D Cartesian coordinates on CPU","title":"Diffusion in 3D Cartesian coordinates on CPU","text":"We will use the same conditions as in the 2D modelling tutorial:","category":"page"},{"location":"diffusion_3D/","page":"Diffusion in 3D Cartesian coordinates on CPU","title":"Diffusion in 3D Cartesian coordinates on CPU","text":"# define total length in x and y\nLx = 9000.0u\"µm\"\nLy = 9000.0u\"µm\"\nLz = 9000.0u\"µm\"\n# define total time for the model\ntfinal = 1.0u\"Myr\"\n# define the pressure and temperature conditions\nT = 900u\"°C\"\nP = 0.6u\"GPa\"\n\nIC3D = InitialConditions3D(Mg0, Fe0, Mn0, Lx, Ly, Lz, tfinal; grt_boundary = grt_boundary)\ndomain3D = Domain(IC3D, T, P)","category":"page"},{"location":"diffusion_3D/","page":"Diffusion in 3D Cartesian coordinates on CPU","title":"Diffusion in 3D Cartesian coordinates on CPU","text":"Now, let's define a callback function to save the results of the simulation in a HDF5 file, similar to the tutorial Saving output as HDF5 files:","category":"page"},{"location":"diffusion_3D/","page":"Diffusion in 3D Cartesian coordinates on CPU","title":"Diffusion in 3D Cartesian coordinates on CPU","text":"@unpack t_charact = domain3D # unpack characteristic time to nondimensionalise the time for the simulation\ntime_save_ad = ustrip.(u\"Myr\", time_save) ./ t_charact # convert to Myr, remove units, and convert to nondimensional time\n\nsave_data_callback = PresetTimeCallback(time_save_ad, save_data_paraview)\n\npath_save = \"Grt_3D.h5\" # choose the name and the path of the HDF5 output file (make sure to add .h5 or .hdf5 at the end)","category":"page"},{"location":"diffusion_3D/","page":"Diffusion in 3D Cartesian coordinates on CPU","title":"Diffusion in 3D Cartesian coordinates on CPU","text":"note: Note\nWe used here save_data_paraview() instead of save_data() to save the data in a format that can be read by Paraview.","category":"page"},{"location":"diffusion_3D/","page":"Diffusion in 3D Cartesian coordinates on CPU","title":"Diffusion in 3D Cartesian coordinates on CPU","text":"We can now use the function simulate() to solve our system:","category":"page"},{"location":"diffusion_3D/","page":"Diffusion in 3D Cartesian coordinates on CPU","title":"Diffusion in 3D Cartesian coordinates on CPU","text":"# solve the problem using DifferentialEquations.jl\nsol = simulate(Domain3D; callback=save_data_callback, path_save=path_save, save_everystep=false, progress=true, progress_steps=1)","category":"page"},{"location":"diffusion_3D/","page":"Diffusion in 3D Cartesian coordinates on CPU","title":"Diffusion in 3D Cartesian coordinates on CPU","text":"note: Note\nThe calculation can take some time, about 4 minutes on my machine with 16 threads.","category":"page"},{"location":"diffusion_3D/","page":"Diffusion in 3D Cartesian coordinates on CPU","title":"Diffusion in 3D Cartesian coordinates on CPU","text":"Two files should have been created in the same folder as your current session: Grt_3D.h5 and Grt_3D.xdmf. The first one is the HDF5 file containing the results of the simulation, and the second one is an XDMF file describing the HDF5 file. This last file can be opened with visualisation software programs, such as Paraview.","category":"page"},{"location":"diffusion_3D/","page":"Diffusion in 3D Cartesian coordinates on CPU","title":"Diffusion in 3D Cartesian coordinates on CPU","text":"warning: Warning\nFor any visualisation software, make sure you open the XDMF file and not the HDF5 file. For Paraview, select the XDMF Reader as the reader when you open your data. Only Paraview has been tested with this package, but other software should work as well.","category":"page"},{"location":"diffusion_spherical/#sph_diffusion","page":"Diffusion in spherical coordinates","title":"Diffusion in spherical coordinates","text":"","category":"section"},{"location":"diffusion_spherical/","page":"Diffusion in spherical coordinates","title":"Diffusion in spherical coordinates","text":"When working with 1D profiles of major elements in garnet, it is often better to consider spherical coordinates. This approximates a garnet grain as a sphere and allows the change in volume relative to the distance from the core to be taken into account.","category":"page"},{"location":"diffusion_spherical/","page":"Diffusion in spherical coordinates","title":"Diffusion in spherical coordinates","text":"In DiffusionGarnet, spherical coordinates are as easy to use as 1D Cartesian coordinates. To illustrate this, we will compare the diffusion of a core-rim profile using the 2 coordinate systems. We will follow the same procedure as in the 1D diffusion tutorial.","category":"page"},{"location":"diffusion_spherical/","page":"Diffusion in spherical coordinates","title":"Diffusion in spherical coordinates","text":"First we will load the data of the core-rim profile from the Spherical examples section, which should be in the same folder as your current session:","category":"page"},{"location":"diffusion_spherical/","page":"Diffusion in spherical coordinates","title":"Diffusion in spherical coordinates","text":"using DiffusionGarnet\nusing DelimitedFiles\n\ndata = DelimitedFiles.readdlm(\"Data_Grt_Sph.txt\", '\\t', '\\n', header=true)[1]\n\nMg0 = data[:, 4]\nFe0 = data[:, 2]\nMn0 = data[:, 3]\nCa0 = data[:, 5]\ndistance = data[:, 1]\nLx = Lr = (data[end,1] - data[1,1])u\"µm\"\ntfinal = 15u\"Myr\"","category":"page"},{"location":"diffusion_spherical/","page":"Diffusion in spherical coordinates","title":"Diffusion in spherical coordinates","text":"We can visualize our data:","category":"page"},{"location":"diffusion_spherical/","page":"Diffusion in spherical coordinates","title":"Diffusion in spherical coordinates","text":"using Plots\n\nl = @layout [a ; b]\n\np1 = plot(distance, Fe0, label=\"Fe initial\", linestyle = :dash, linewidth=1, dpi=200, title = \"Initial conditions\", legend=:outerbottomright, linecolor=1,xlabel = \"Distance (µm)\", ylabel=\"Molar fraction\")\n\np2 = plot(distance, Mg0, label=\"Mg initial\", linestyle = :dash, linewidth=1, dpi=200,legend=:outerbottomright,linecolor=2,xlabel = \"Distance (µm)\")\np2 = plot!(distance, Mn0, label=\"Mn initial\", linestyle = :dash, linewidth=1, linecolor=3)\np2 = plot!(distance, Ca0, label=\"Ca initial\", linestyle = :dash, linewidth=1, linecolor=4, ylabel=\"Molar fraction\")\n\nplot(p1, p2, layout = l)","category":"page"},{"location":"diffusion_spherical/","page":"Diffusion in spherical coordinates","title":"Diffusion in spherical coordinates","text":"which outputs:","category":"page"},{"location":"diffusion_spherical/","page":"Diffusion in spherical coordinates","title":"Diffusion in spherical coordinates","text":"(Image: Initial conditions.)","category":"page"},{"location":"diffusion_spherical/","page":"Diffusion in spherical coordinates","title":"Diffusion in spherical coordinates","text":"Then, we can define our initial conditions, both for the spherical and 1D Cartesian coordinates.","category":"page"},{"location":"diffusion_spherical/","page":"Diffusion in spherical coordinates","title":"Diffusion in spherical coordinates","text":"ICSph = InitialConditionsSpherical(Mg0, Fe0, Mn0, Lr, tfinal)\nIC1D = InitialConditions1D(Mg0, Fe0, Mn0, Lx, tfinal)\n\n# define the pressure and temperature conditions of diffusion\nT = 900u\"°C\"\nP = 0.6u\"GPa\"\n\nDomainSph = Domain(ICSph, T, P)\nDomain1D = Domain(IC1D, T, P; bc_neumann = (true, false))","category":"page"},{"location":"diffusion_spherical/","page":"Diffusion in spherical coordinates","title":"Diffusion in spherical coordinates","text":"tip: Tip\nInitialConditionsSpherical always assumes that the core of the profile is on the left and the edge is on the right side.","category":"page"},{"location":"diffusion_spherical/","page":"Diffusion in spherical coordinates","title":"Diffusion in spherical coordinates","text":"note: Note\nWe use bc_neumann in Domain to indicate that a homogeneous Neumann boundary must be ensured on the left side of the 1D profile to simulate the core of the garnet grain.","category":"page"},{"location":"diffusion_spherical/","page":"Diffusion in spherical coordinates","title":"Diffusion in spherical coordinates","text":"We can then use the function simulate() to solve our system:","category":"page"},{"location":"diffusion_spherical/","page":"Diffusion in spherical coordinates","title":"Diffusion in spherical coordinates","text":"# solve the problem using DifferentialEquations.jl\nsol_sph = simulate(DomainSph)\nsol_1D = simulate(Domain1D)","category":"page"},{"location":"diffusion_spherical/","page":"Diffusion in spherical coordinates","title":"Diffusion in spherical coordinates","text":"and compare the 2 solutions:","category":"page"},{"location":"diffusion_spherical/","page":"Diffusion in spherical coordinates","title":"Diffusion in spherical coordinates","text":"@unpack t_charact = DomainSph\n\nanim = @animate for i = LinRange(0, sol_sph.t[end], 100)\n l = @layout [a ; b]\n\n p1 = plot(distance, Fe0, label=\"Fe initial\", linestyle = :dash, linewidth=1, dpi=200, title = \"Total Time = $(round(i*t_charact;digits=2)) Ma\", legend=:outerbottomright, linecolor=1,xlabel = \"Distance (µm)\")\n p1 = plot!(distance, sol_sph(i)[:,2], label=\"Fe Sph\",linecolor=1, linewidth=1)\n p1 = plot!(distance, sol_1D(i)[:,2], label=\"Fe 1D\",linecolor=5, linewidth=1)\n\n\n p2 = plot(distance, Mg0, label=\"Mg initial\", linestyle = :dash, linewidth=1, dpi=200,legend=:outerbottomright,linecolor=2,xlabel = \"Distance (µm)\")\n p2 = plot!(distance, Mn0, label=\"Mn initial\", linestyle = :dash, linewidth=1, linecolor=3)\n p2 = plot!(distance, Ca0, label=\"Ca initial\", linestyle = :dash, linewidth=1, linecolor=4)\n p2 = plot!(distance, sol_sph(i)[:,1], label=\"Mg Sph\",linecolor=2, linewidth=1)\n p2 = plot!(distance, sol_1D(i)[:,1], label=\"Mg 1D\",linecolor=6, linewidth=1)\n\n p2 = plot!(distance, sol_sph(i)[:,3], label=\"Mn Sph\", linecolor=3, linewidth=1)\n p2 = plot!(distance, sol_1D(i)[:,3], label=\"Mn 1D\", linecolor=7, linewidth=1)\n\n p2 = plot!(distance, 1 .- sol_sph(i)[:,1] .- sol_sph(i)[:,2] .- sol_sph(i)[:,3], label=\"Ca Sph\", linecolor=4, linewidth=1)\n p2 = plot!(distance, 1 .- sol_sph(i)[:,1] .- sol_sph(i)[:,2] .- sol_sph(i)[:,3], label=\"Ca 1D\", linecolor=8, linewidth=1)\n\n plot(p1, p2, layout = l)\nend every 1\n\nprintln(\"Now, generating the gif...\")\ngif(anim, \"Grt_Spherical+1D.gif\", fps = 7)\nprintln(\"...Done!\")","category":"page"},{"location":"diffusion_spherical/","page":"Diffusion in spherical coordinates","title":"Diffusion in spherical coordinates","text":"With the resulting gif:","category":"page"},{"location":"diffusion_spherical/","page":"Diffusion in spherical coordinates","title":"Diffusion in spherical coordinates","text":"(Image: Spherical diffusion profil of a garnet)","category":"page"},{"location":"diffusion_spherical/","page":"Diffusion in spherical coordinates","title":"Diffusion in spherical coordinates","text":"It shows that using spherical coordinates makes the garnet diffuse faster. Using 1D Cartesian coordinates can overestimate the equilibration time.","category":"page"},{"location":"saving_output/#saving_output","page":"Saving output as HDF5 files","title":"Saving output as HDF5 files","text":"","category":"section"},{"location":"saving_output/","page":"Saving output as HDF5 files","title":"Saving output as HDF5 files","text":"When dealing with 2D or especially 3D data, it can be impractical or even impossible to keep every timestep in memory, as the RAM on a single machine can quickly become saturated. It may then be relevant to save certain timesteps of interest to disk for later post-processing. DiffusionGarnet has a built-in function for this purpose, using a callback function based on the DiffEqCallbacks package to produce HDF5 files.","category":"page"},{"location":"saving_output/","page":"Saving output as HDF5 files","title":"Saving output as HDF5 files","text":"HDF5 is a data format designed to store and organise large amount of data and is the data format chosen for DiffusionGarnet.","category":"page"},{"location":"saving_output/","page":"Saving output as HDF5 files","title":"Saving output as HDF5 files","text":"As an example, we will use the data from the Diffusion in 2D Cartesian coordinates on CPU tutorial but this is also applicable to 1D Cartesian or spherical coordinates in the same manner.","category":"page"},{"location":"saving_output/","page":"Saving output as HDF5 files","title":"Saving output as HDF5 files","text":"tip: Tip\nMake sure to start Julia with multiple threads of execution to speed up the simulation.","category":"page"},{"location":"saving_output/","page":"Saving output as HDF5 files","title":"Saving output as HDF5 files","text":"using DiffusionGarnet\nusing DelimitedFiles\nusing Plots\n\nMg0 = DelimitedFiles.readdlm(\"Xprp.txt\", '\\t', '\\n', header=false)\nFe0 = DelimitedFiles.readdlm(\"Xalm.txt\", '\\t', '\\n', header=false)\nMn0 = DelimitedFiles.readdlm(\"Xsps.txt\", '\\t', '\\n', header=false)\nCa0 = DelimitedFiles.readdlm(\"Xgrs.txt\", '\\t', '\\n', header=false)\ngrt_boundary = DelimitedFiles.readdlm(\"contour_Grt.txt\", '\\t', '\\n', header=false)\n\nLx = 9000.0u\"µm\"\nLy = 9000.0u\"µm\"\ntfinal = 1.0u\"Myr\"\nT = 900u\"°C\"\nP = 0.6u\"GPa\"\n\nIC2D = InitialConditions2D(Mg0, Fe0, Mn0, Lx, Ly, tfinal; grt_boundary = grt_boundary)\ndomain2D = Domain(IC2D, T, P)","category":"page"},{"location":"saving_output/","page":"Saving output as HDF5 files","title":"Saving output as HDF5 files","text":"We then need to decide at which timesteps we want to save our data, here at 0, 0.5 and 1 Myr:","category":"page"},{"location":"saving_output/","page":"Saving output as HDF5 files","title":"Saving output as HDF5 files","text":"time_save = [0, 0.5, 1]u\"Myr\" # define times at which to save","category":"page"},{"location":"saving_output/","page":"Saving output as HDF5 files","title":"Saving output as HDF5 files","text":"note: Note\nMake sure to always include 0 in your timesteps to save the initial conditions.","category":"page"},{"location":"saving_output/","page":"Saving output as HDF5 files","title":"Saving output as HDF5 files","text":"Two additional steps are then required, we need to convert this time to dimensionless time so that the solver can use it, and we need to define our callback function:","category":"page"},{"location":"saving_output/","page":"Saving output as HDF5 files","title":"Saving output as HDF5 files","text":"@unpack t_charact = domain2D # unpack characteristic time to nondimensionalise the time for the simulation\ntime_save_ad = ustrip.(u\"Myr\", time_save) ./ t_charact # convert to Myr, remove units, and convert to nondimensional time","category":"page"},{"location":"saving_output/","page":"Saving output as HDF5 files","title":"Saving output as HDF5 files","text":"time_save_ad is the equivalent dimensionless time to time_save:","category":"page"},{"location":"saving_output/","page":"Saving output as HDF5 files","title":"Saving output as HDF5 files","text":"julia > time_save_ad\n3-element Vector{Float64}:\n 0.0\n 0.011476148087304199\n 0.022952296174608398","category":"page"},{"location":"saving_output/","page":"Saving output as HDF5 files","title":"Saving output as HDF5 files","text":"warning: Warning\nThe time provided to the callback function should always be dimensionless.","category":"page"},{"location":"saving_output/","page":"Saving output as HDF5 files","title":"Saving output as HDF5 files","text":"We can now define our callback function that will save the data at these timesteps, using PresetTimeCallback():","category":"page"},{"location":"saving_output/","page":"Saving output as HDF5 files","title":"Saving output as HDF5 files","text":"save_data_callback = PresetTimeCallback(time_save_ad, save_data)","category":"page"},{"location":"saving_output/","page":"Saving output as HDF5 files","title":"Saving output as HDF5 files","text":"save_data is the function that will be called at each timestep in our callback function. save_data_callback can then be passed to simulate after defining the path of the output file:","category":"page"},{"location":"saving_output/","page":"Saving output as HDF5 files","title":"Saving output as HDF5 files","text":"path_save = \"Grt_2D.h5\" # chose the name and the path of the HDF5 output file (make sure to add .h5 or .hdf5 at the end)\nsol = simulate(domain2D; callback=save_data_callback, path_save=path_save, save_everystep=false);","category":"page"},{"location":"saving_output/","page":"Saving output as HDF5 files","title":"Saving output as HDF5 files","text":"note: Note\nWe use save_everystep=false in simulation() to prevent saving every timestep in the model as we save it to disk. That should speed up the computation and reduce the burden on the RAM.","category":"page"},{"location":"saving_output/","page":"Saving output as HDF5 files","title":"Saving output as HDF5 files","text":"This output:","category":"page"},{"location":"saving_output/","page":"Saving output as HDF5 files","title":"Saving output as HDF5 files","text":"Data saved at 0.0 Myr.\nData saved at 0.5 Myr.\nData saved at 1.0 Myr.\n352.464431 seconds (17.47 M allocations: 1.789 GiB, 0.15% gc time, 4.01% compilation time)","category":"page"},{"location":"saving_output/","page":"Saving output as HDF5 files","title":"Saving output as HDF5 files","text":"which indicates the time at which we saved our data and the total run time of the simulation.","category":"page"},{"location":"saving_output/","page":"Saving output as HDF5 files","title":"Saving output as HDF5 files","text":"The output file produced can be read with the official viewer or using HDF5.jl in Julia or other programming languages and external programs.","category":"page"},{"location":"saving_output/","page":"Saving output as HDF5 files","title":"Saving output as HDF5 files","text":"Using HDF5.jl, we can look at the structure of the HDF5 file produced:","category":"page"},{"location":"saving_output/","page":"Saving output as HDF5 files","title":"Saving output as HDF5 files","text":"julia> using HDF5\njulia> hf5 = h5open(\"Grt_2D.h5\")\n🗂️ HDF5.File: (read-only) Grt_2D.h5\n└─ 📂 Diffusion_Grt\n ├─ 🏷️ CharacteristicDiffusionCoefficient\n ├─ 🏷️ CharacteristicLength\n ├─ 🏷️ CharacteristicTime\n ├─ 🏷️ Coordinates\n ├─ 🏷️ Dx(µm)\n ├─ 🏷️ Dy(µm)\n ├─ 🏷️ LengthX(µm)\n ├─ 🏷️ LengthY(µm)\n ├─ 🏷️ Nx\n ├─ 🏷️ Ny\n ├─ 🏷️ TotalTime(Myr)\n ├─ 📂 t0000\n │ ├─ 🏷️ Time(Myr)\n │ ├─ 📂 Ca\n │ │ ├─ 🏷️ Center\n │ │ ├─ 🏷️ DataType\n │ │ └─ 🔢 Ca\n │ ├─ 📂 Fe\n │ │ ├─ 🏷️ Center\n │ │ ├─ 🏷️ DataType\n │ │ └─ 🔢 Fe\n │ ├─ 📂 Mg\n │ │ ├─ 🏷️ Center\n │ │ ├─ 🏷️ DataType\n │ │ └─ 🔢 Mg\n │ └─ 📂 Mn\n │ ├─ 🏷️ Center\n │ ├─ 🏷️ DataType\n │ └─ 🔢 Mn\n ├─ 📂 t0001\n │ ├─ 🏷️ CurrentDt(Myr)\n │ ├─ 🏷️ Time(Myr)\n │ ├─ 📂 Ca\n │ │ ├─ 🏷️ Center\n │ │ ├─ 🏷️ DataType\n │ │ └─ 🔢 Ca\n │ ├─ 📂 Fe\n │ │ ├─ 🏷️ Center\n │ │ ├─ 🏷️ DataType\n │ │ └─ 🔢 Fe\n │ ├─ 📂 Mg\n │ │ ├─ 🏷️ Center\n │ │ ├─ 🏷️ DataType\n │ │ └─ 🔢 Mg\n │ └─ 📂 Mn\n │ ├─ 🏷️ Center\n │ ├─ 🏷️ DataType\n │ └─ 🔢 Mn\n └─ 📂 t0002\n ├─ 🏷️ CurrentDt(Myr)\n ├─ 🏷️ Time(Myr)\n ├─ 📂 Ca\n │ ├─ 🏷️ Center\n │ ├─ 🏷️ DataType\n │ └─ 🔢 Ca\n ├─ 📂 Fe\n │ ├─ 🏷️ Center\n │ ├─ 🏷️ DataType\n │ └─ 🔢 Fe\n ├─ 📂 Mg\n │ ├─ 🏷️ Center\n │ ├─ 🏷️ DataType\n │ └─ 🔢 Mg\n └─ 📂 Mn\n ├─ 🏷️ Center\n ├─ 🏷️ DataType\n └─ 🔢 Mn\njulia> close(hf5) # close the HDF5 file","category":"page"},{"location":"saving_output/","page":"Saving output as HDF5 files","title":"Saving output as HDF5 files","text":"We can see that the HDF5 file contains some general information about the model as attributes (🏷️), such as characteristic values or the dimensions and the total time of the model. In addition, the three timesteps are stored as groups (📂) with the composition in Ca, Fe, Mg and Mn.","category":"page"},{"location":"updating_PT/#PT_callback","page":"Updating pressure and temperature conditions","title":"Updating pressure and temperature conditions","text":"","category":"section"},{"location":"updating_PT/","page":"Updating pressure and temperature conditions","title":"Updating pressure and temperature conditions","text":"When modelling diffusion in garnet, it may be relevant to update the pressure and temperature (PT) conditions during the simulation, i.e. to model a decompression event in the retrograde path.","category":"page"},{"location":"updating_PT/","page":"Updating pressure and temperature conditions","title":"Updating pressure and temperature conditions","text":"To do this, we can use a callback function. A callback function is a function that can be called in our solver when a certain condition is met, i.e. when a certain time is reached. It is based on the DiffEqCallbacks package from the DifferentialEquations.jl ecosystem.","category":"page"},{"location":"updating_PT/","page":"Updating pressure and temperature conditions","title":"Updating pressure and temperature conditions","text":"For this tutorial, we will use the same data from the tutorial in 1D Cartesian coordinates:","category":"page"},{"location":"updating_PT/","page":"Updating pressure and temperature conditions","title":"Updating pressure and temperature conditions","text":"using DiffusionGarnet # this can take a while\nusing DelimitedFiles\n# load the data of your choice (here from the text file located in https://github.com/Iddingsite/DiffusionGarnet.jl/tree/main/examples/1D, place it in the same folder as where you are running the code)\ndata = DelimitedFiles.readdlm(\"Data_Grt_1D.txt\", '\\t', '\\n', header=true)[1]\n\nMg0 = data[:, 4] # load initial Mg mole fraction\nFe0 = data[:, 2] # load initial Fe mole fraction\nMn0 = data[:, 3] # load initial Mn mole fraction\nCa0 = data[:, 5] # load initial Ca mole fraction\ndistance = data[:, 1]\n\nLx = (data[end,1] - data[1,1])u\"µm\" # length in x of the model, here in µm\ntfinal = 15u\"Myr\" # total time of the model, here in Myr\n\n# define the initial conditions in 1D of your problem in that order.\nIC1D = InitialConditions1D(Mg0, Fe0, Mn0, Lx, tfinal)","category":"page"},{"location":"updating_PT/","page":"Updating pressure and temperature conditions","title":"Updating pressure and temperature conditions","text":"We will now define the PT evolution of our model. To do this, we define 3 arrays containing the pressure, temperature and the time at which the conditions are updated:","category":"page"},{"location":"updating_PT/","page":"Updating pressure and temperature conditions","title":"Updating pressure and temperature conditions","text":"T = [900, 850, 800, 750, 700, 650, 600, 550]u\"°C\"\nP = [0.6, 0.5, 0.4, 0.3, 0.3, 0.3, 0.3, 0.3]u\"GPa\"\ntime_update = [0, 2, 4, 6, 8, 10, 12, 14]u\"Myr\"","category":"page"},{"location":"updating_PT/","page":"Updating pressure and temperature conditions","title":"Updating pressure and temperature conditions","text":"note: Note\nAll three arrays must have the same size and with units. Make sure to also define the conditions at time 0.","category":"page"},{"location":"updating_PT/","page":"Updating pressure and temperature conditions","title":"Updating pressure and temperature conditions","text":"This can be supplied to Domain:","category":"page"},{"location":"updating_PT/","page":"Updating pressure and temperature conditions","title":"Updating pressure and temperature conditions","text":"domain1D = Domain(IC1D, T, P, time_update)","category":"page"},{"location":"updating_PT/","page":"Updating pressure and temperature conditions","title":"Updating pressure and temperature conditions","text":"Now that we have the time and PT conditions to update, we need to create our callback function. We will use the PresetTimeCallback() function which takes as input the times at which the callback is to be triggered and the function to be called at that time.","category":"page"},{"location":"updating_PT/","page":"Updating pressure and temperature conditions","title":"Updating pressure and temperature conditions","text":"# extract time_update_ad from domain1D\n@unpack time_update_ad = domain1D\n\n# define our callback function\nupdate_diffusion_coef_call = PresetTimeCallback(time_update_ad, update_diffusion_coef)","category":"page"},{"location":"updating_PT/","page":"Updating pressure and temperature conditions","title":"Updating pressure and temperature conditions","text":"note: Note\ntime_update_ad is the equivalent of time_update in nondimensional time. update_diffusion_coef is a function that updates the diffusion coefficients according to the new PT conditions.","category":"page"},{"location":"updating_PT/","page":"Updating pressure and temperature conditions","title":"Updating pressure and temperature conditions","text":"We can now provide our callback function to simulate and run our model:","category":"page"},{"location":"updating_PT/","page":"Updating pressure and temperature conditions","title":"Updating pressure and temperature conditions","text":"sol = simulate(domain1D; callback=update_diffusion_coef_call);","category":"page"},{"location":"updating_PT/","page":"Updating pressure and temperature conditions","title":"Updating pressure and temperature conditions","text":"This outputs the time at which the callback function was called and the new PT conditions, and the total solver runtime:","category":"page"},{"location":"updating_PT/","page":"Updating pressure and temperature conditions","title":"Updating pressure and temperature conditions","text":"New temperature and pressure: 850.0 °C and 5.0 kbar, updated at 2.0 Myr.\nNew temperature and pressure: 800.0 °C and 4.0 kbar, updated at 4.0 Myr.\nNew temperature and pressure: 750.0 °C and 3.0 kbar, updated at 6.0 Myr.\nNew temperature and pressure: 700.0 °C and 3.0 kbar, updated at 8.0 Myr.\nNew temperature and pressure: 650.0 °C and 3.0 kbar, updated at 10.0 Myr.\nNew temperature and pressure: 600.0 °C and 3.0 kbar, updated at 12.0 Myr.\nNew temperature and pressure: 550.0 °C and 3.0 kbar, updated at 14.0 Myr.\n 0.752661 seconds (44.65 k allocations: 25.676 MiB)","category":"page"},{"location":"updating_PT/","page":"Updating pressure and temperature conditions","title":"Updating pressure and temperature conditions","text":"We can now plot our results and see how the change in PT has affected our diffusion profiles:","category":"page"},{"location":"updating_PT/","page":"Updating pressure and temperature conditions","title":"Updating pressure and temperature conditions","text":"@unpack tfinal_ad, t_charact = domain1D\n\nanim = @animate for i = LinRange(0, tfinal_ad, 100)\n l = @layout [a ; b]\n\n p1 = plot(distance, Fe0, label=\"Fe initial\", linestyle = :dash, linewidth=1, dpi=200, title = \"Total Time = $(round(((i)* t_charact);digits=2)) Ma\", legend=:outerbottomright, linecolor=1,xlabel = \"Distance (µm)\")\n p1 = plot!(distance, sol(i)[:,2], label=\"Fe\",linecolor=1, linewidth=1)\n\n\n p2 = plot(distance, Mg0, label=\"Mg initial\", linestyle = :dash, linewidth=1, dpi=200,legend=:outerbottomright,linecolor=2,xlabel = \"Distance (µm)\")\n p2 = plot!(distance, Mn0, label=\"Mn initial\", linestyle = :dash, linewidth=1, linecolor=3)\n p2 = plot!(distance, Ca0, label=\"Ca initial\", linestyle = :dash, linewidth=1, linecolor=4)\n p2 = plot!(distance, sol(i)[:,1], label=\"Mg\",linecolor=2, linewidth=1)\n\n p2 = plot!(distance, sol(i)[:,3], label=\"Mn\", linecolor=3, linewidth=1)\n\n p2 = plot!(distance, 1 .- sol(i)[:,1] .- sol(i)[:,2] .- sol(i)[:,3], label=\"Ca\", linecolor=4, linewidth=1)\n\n plot(p1, p2, layout = l)\nend every 1\n\nprintln(\"Now, generating the gif...\")\ngif(anim, \"Grt_1D_var_TP.gif\", fps = 7)\nprintln(\"...Done!\")\n","category":"page"},{"location":"updating_PT/","page":"Updating pressure and temperature conditions","title":"Updating pressure and temperature conditions","text":"which outputs:","category":"page"},{"location":"updating_PT/","page":"Updating pressure and temperature conditions","title":"Updating pressure and temperature conditions","text":"(Image: 1D diffusion profil of a garnet with PT)","category":"page"},{"location":"updating_PT/","page":"Updating pressure and temperature conditions","title":"Updating pressure and temperature conditions","text":"We can clearly observe the slowing down of the diffusion processes as the PT conditions decrease with time.","category":"page"},{"location":"background/#background","page":"Background","title":"Background","text":"","category":"section"},{"location":"background/#Theory","page":"Background","title":"Theory","text":"","category":"section"},{"location":"background/","page":"Background","title":"Background","text":"Garnet is a mineral commonly used in metamorphic petrology to better understand geological processes, as it occurs in a variety of rock types. This mineral often exhibits a wide range of compositional zoning, which can be interpreted as recording ranges of pressure and temperature (PT) conditions. Modelling diffusion processes can help to better understand this zoning and better constrain the pressure-temperature-time (PTt) conditions of the metamorphic event of interest.","category":"page"},{"location":"background/","page":"Background","title":"Background","text":"Major element diffusion in garnet is described by a coupled multicomponent system between four poles: Mg, Fe, Mn and Ca. This can be expressed by a system of parabolic partial differential equations (PDEs) following Fick's first law of diffusion:","category":"page"},{"location":"background/","page":"Background","title":"Background","text":"beginequation\n labelficks_law\nbeginaligned\n fracpartial C_ipartial t = mathbfnabla cdot sum_j=1^n-1 textbfD_ij mathbfnabla C_j\nendaligned\nendequation","category":"page"},{"location":"background/","page":"Background","title":"Background","text":"with C_i the molar fraction of the i^th pole of garnet, n the total number of poles, and textbfD_ij the interdiffusion coefficient tensor.","category":"page"},{"location":"background/","page":"Background","title":"Background","text":"For a system of n components, only n-1 equations are needed, as the sum of the molar fractions of all components is equal to 1. Thus, Ca is made dependent on the other concentrations and (1) is made up of three coupled PDEs.","category":"page"},{"location":"background/","page":"Background","title":"Background","text":"According to Lasaga (1979) [1], textbfD_ij can be defined, assuming an ideal system for garnet, with:","category":"page"},{"location":"background/","page":"Background","title":"Background","text":"beginequation\n labeldiffusion_coefficient\nbeginaligned\n textbfD_ij = D_i^* delta_ij - fracC_i z_i z_j D_i^*sum_k=1^n C_k z_k^2 D_k^* (D_j^* - D_n^*) =\nbeginbmatrix\n D_MgMg D_MgFe D_MgMn \n D_FeMg D_FeFe D_FeMn \n D_MnMg D_MnFe D_MnMn \nendbmatrix\nendaligned\nendequation","category":"page"},{"location":"background/","page":"Background","title":"Background","text":"with D_i^* the tracer diffusion coefficient of the i^th component, delta_ij the Kronecker delta, z the charge of the species, and D_n^* the tracer diffusion coefficient of the dependent molar fraction (here Ca).","category":"page"},{"location":"background/","page":"Background","title":"Background","text":"The tracer diffusion coefficient D_i^* can be calculated from an Arrhenius relationship:","category":"page"},{"location":"background/","page":"Background","title":"Background","text":"beginequation\n labelarrhenius\nD_i^* = D_0i exp left( -fracE_ai - (P-1)Delta V^+_iRT right)\nendequation","category":"page"},{"location":"background/","page":"Background","title":"Background","text":"with D_0i the pre-exponential constant, E_ai the activation energy of diffusion, Delta V^+_i the activation volume of diffusion at 1 bar, R the universal gas constant, T the temperature, and P the pressure.","category":"page"},{"location":"background/","page":"Background","title":"Background","text":"In DiffusionGarnet.jl, D_0i, E_ai, and Delta V^+_i are those of Chakraborty & Ganguly (1992) [2]. The tracer diffusion coefficient of Ca is defined as 05D_Fe, following the approach of Loomis et al. (1985) [3].","category":"page"},{"location":"background/#Numerical-approach","page":"Background","title":"Numerical approach","text":"","category":"section"},{"location":"background/","page":"Background","title":"Background","text":"By defining the PT conditions of the metamorphic event of interest, (3) can be solved for each component, and the diffusion coefficient tensor can be calculated using (2) from the initial major composition data. In DiffusionGarnet.jl, (1) is then discretised in space using finite differences, and the resulting system of ordinary differential equations is solved with ROCK2, a stabilised explicit method (Abdulle & Medovikov, 2001 [4]) using the DifferentialEquations.jl ecosystem.","category":"page"},{"location":"background/#refs","page":"Background","title":"References","text":"","category":"section"},{"location":"background/","page":"Background","title":"Background","text":"[1] Lasaga, A. C. (1979). Multicomponent exchange and diffusion in silicates. Geochimica et Cosmochimica Acta, 43(4), 455-469.","category":"page"},{"location":"background/","page":"Background","title":"Background","text":"[2] Chakraborty, S., & Ganguly, J. (1992). Cation diffusion in aluminosilicate garnets: experimental determination in spessartine-almandine diffusion couples, evaluation of effective binary diffusion coefficients, and applications. Contributions to Mineralogy and petrology, 111(1), 74-86.","category":"page"},{"location":"background/","page":"Background","title":"Background","text":"[3] Loomis, T. P., Ganguly, J., Elphick, S. C., 1985. Experimental determinations of cation diffusitivities in aluminosilicate garnets. II. Multicomponent simulation and tracer diffusion coefficients. Contributions to Mineralogy and Petrology 90, 45–51.","category":"page"},{"location":"background/","page":"Background","title":"Background","text":"[4] Abdulle, A., & Medovikov, A. A. (2001). Second order Chebyshev methods based on orthogonal polynomials. Numerische Mathematik, 90, 1-18.","category":"page"},{"location":"reference/","page":"List of functions","title":"List of functions","text":"","category":"page"},{"location":"reference/","page":"List of functions","title":"List of functions","text":"Modules = [DiffusionGarnet]\n","category":"page"},{"location":"reference/#DiffusionGarnet.Domain","page":"List of functions","title":"DiffusionGarnet.Domain","text":"Domain(IC::InitialConditions2D, T::Union{Unitful.Temperature,Array{<:Unitful.Temperature{<:Real}, 1}}, P::Union{Unitful.Pressure,Array{<:Unitful.Pressure{<:Real}, 1}}, time_update::Union{Unitful.Time,Array{<:Unitful.Time{<:Real}, 1}}=0u\"Myr\")\n\nWhen applied to 2D initial conditions, define corresponding 2D domain.\n\n\n\n\n\n","category":"type"},{"location":"reference/#DiffusionGarnet.Domain-2","page":"List of functions","title":"DiffusionGarnet.Domain","text":"Domain(IC::InitialConditions1D, T::Union{Unitful.Temperature,Array{<:Unitful.Temperature{<:Real}, 1}}, P::Union{Unitful.Pressure,Array{<:Unitful.Pressure{<:Real}, 1}}, time_update::Union{Unitful.Time,Array{<:Unitful.Time{<:Real}, 1}}=0u\"Myr\"; bc_neumann::Tuple=(false, false))\n\nWhen applied to 1D initial conditions, define corresponding 1D domain. bc_neumann can be used to define Neumann boundary conditions on the left or right side of the domain if set to true.\n\n\n\n\n\n","category":"type"},{"location":"reference/#DiffusionGarnet.Domain-3","page":"List of functions","title":"DiffusionGarnet.Domain","text":"Domain(IC::InitialConditionsSpherical, T::Union{Unitful.Temperature,Array{<:Unitful.Temperature{<:Real}, 1}}, P::Union{Unitful.Pressure,Array{<:Unitful.Pressure{<:Real}, 1}}, time_update::Union{Unitful.Time,Array{<:Unitful.Time{<:Real}, 1}}=0u\"Myr\")\n\nWhen applied to spherical initial conditions, define corresponding spherical domain. Assume that the center of the grain is on the left side.\n\n\n\n\n\n","category":"type"},{"location":"reference/#DiffusionGarnet.Domain-4","page":"List of functions","title":"DiffusionGarnet.Domain","text":"Domain(IC::InitialConditions3D, T::Union{Unitful.Temperature,Array{<:Unitful.Temperature{<:Real}, 1}}, P::Union{Unitful.Pressure,Array{<:Unitful.Pressure{<:Real}, 1}}, time_update::Union{Unitful.Time,Array{<:Unitful.Time{<:Real}, 1}}=0u\"Myr\")\n\nWhen applied to 3D initial conditions, define corresponding 3D domain.\n\n\n\n\n\n","category":"type"},{"location":"reference/#DiffusionGarnet.Domain-5","page":"List of functions","title":"DiffusionGarnet.Domain","text":"Domain(IC, T, P, time_update=0u\"Myr\")\n\nReturn a struct containing the struct IC, the PT conditions T and P and the nondimensionalised parameters of the problem. time_update can be used to update the PT conditions during the simulation. If so, T and P have to be of the same size as time_update and correspond to the values of PT to update to.\n\n\n\n\n\n","category":"type"},{"location":"reference/#DiffusionGarnet.InitialConditions1D-Tuple{AbstractVector{<:Real}, AbstractVector{<:Real}, AbstractVector{<:Real}, Union{Quantity{T, 𝐋, U}, Level{L, S, Quantity{T, 𝐋, U}} where {L, S}} where {T, U}, Union{Quantity{T, 𝐓, U}, Level{L, S, Quantity{T, 𝐓, U}} where {L, S}} where {T, U}}","page":"List of functions","title":"DiffusionGarnet.InitialConditions1D","text":"InitialConditions1D(CMg0::Array{<:Real, 1}, CFe0::Array{<:Real, 1}, CMn0::Array{<:Real, 1}, Lx::Unitful.Length, tfinal::Unitful.Time)\n\nReturn a structure containing the initial conditions for a 1D diffusion problem. CMg0, CFe0 and CMn0 need to be in molar fraction. Convert the Lxandtfinal`` to µm and Myr respectively.\n\n\n\n\n\n","category":"method"},{"location":"reference/#DiffusionGarnet.InitialConditions2D-Tuple{AbstractMatrix{<:Real}, AbstractMatrix{<:Real}, AbstractMatrix{<:Real}, Union{Quantity{T, 𝐋, U}, Level{L, S, Quantity{T, 𝐋, U}} where {L, S}} where {T, U}, Union{Quantity{T, 𝐋, U}, Level{L, S, Quantity{T, 𝐋, U}} where {L, S}} where {T, U}, Union{Quantity{T, 𝐓, U}, Level{L, S, Quantity{T, 𝐓, U}} where {L, S}} where {T, U}}","page":"List of functions","title":"DiffusionGarnet.InitialConditions2D","text":"InitialConditions2D(CMg0::AbstractArray{<:Real, 2}, CFe0::AbstractArray{<:Real, 2}, CMn0::AbstractArray{<:Real, 2}, Lx::Unitful.Length, Ly::Unitful.Length, tfinal::Unitful.Time; grt_boundary::AbstractArray{<:Real, 2}=zeros(eltype(CMg0), size(CMg0)...))\n\nReturn a structure containing the initial conditions for a 2D diffusion problem. CMg0, CFe0 and CMn0 need to be in molar fraction. Convert Lx, Ly and tfinal to µm, µm and Myr respectively. grt_boundary is a matrix of the same size as CMg0, CFe0 and CMn0. It is a binary matrix with 1 where the contour of the garnet is present and 0 otherwise. It is used to apply the Dirichlet boundary condition in the model. If not provided, it is equal to zero everywhere.\n\n\n\n\n\n","category":"method"},{"location":"reference/#DiffusionGarnet.InitialConditions3D-Tuple{AbstractArray{<:Real, 3}, AbstractArray{<:Real, 3}, AbstractArray{<:Real, 3}, Union{Quantity{T, 𝐋, U}, Level{L, S, Quantity{T, 𝐋, U}} where {L, S}} where {T, U}, Union{Quantity{T, 𝐋, U}, Level{L, S, Quantity{T, 𝐋, U}} where {L, S}} where {T, U}, Union{Quantity{T, 𝐋, U}, Level{L, S, Quantity{T, 𝐋, U}} where {L, S}} where {T, U}, Union{Quantity{T, 𝐓, U}, Level{L, S, Quantity{T, 𝐓, U}} where {L, S}} where {T, U}}","page":"List of functions","title":"DiffusionGarnet.InitialConditions3D","text":"InitialConditions3D(CMg0::AbstractArray{<:Real, 3}, CFe0::AbstractArray{<:Real, 3}, CMn0::AbstractArray{<:Real, 3}, Lx::Unitful.Length, Ly::Unitful.Length, Lz::Unitful.Length, tfinal::Unitful.Time; grt_boundary::Union{AbstractArray{<:Real, 3}, AbstractArray{<:Bool, 3}}=zeros(Bool, size(CMg0)...))\n\nReturn a structure containing the initial conditions for a 3D diffusion problem. CMg0, CFe0 and CMn0 need to be in molar fraction. Convert Lx, Ly, Lz and tfinal to µm, µm, µm and Myr respectively. grt_boundary is a matrix of the same size as CMg0, CFe0 and CMn0. It is a binary matrix with 1 where the contour of the garnet is present and 0 otherwise. It is used to apply the Dirichlet boundary condition in the model. If not provided, it is equal to zero everywhere.\n\n\n\n\n\n","category":"method"},{"location":"reference/#DiffusionGarnet.InitialConditionsSpherical-Tuple{AbstractVector{<:Real}, AbstractVector{<:Real}, AbstractVector{<:Real}, Union{Quantity{T, 𝐋, U}, Level{L, S, Quantity{T, 𝐋, U}} where {L, S}} where {T, U}, Union{Quantity{T, 𝐓, U}, Level{L, S, Quantity{T, 𝐓, U}} where {L, S}} where {T, U}}","page":"List of functions","title":"DiffusionGarnet.InitialConditionsSpherical","text":"InitialConditionsSpherical(CMg0::AbstractArray{<:Real, 1}, CFe0::AbstractArray{<:Real, 1}, CMn0::AbstractArray{<:Real, 1}, Lr::Unitful.Length, tfinal::Unitful.Time)\n\nReturn a structure containing the initial conditions for a spherical diffusion problem. CMg0, CFe0 and CMn0 need to be in molar fraction. Convert Lr and tfinal to µm and Myr respectively.\n\n\n\n\n\n","category":"method"},{"location":"reference/#DiffusionGarnet.save_data-Tuple{Any}","page":"List of functions","title":"DiffusionGarnet.save_data","text":"save_data(integrator)\n\nCallback function used to save major element compositions to an HDF5 file at a specific timestep.\n\n\n\n\n\n","category":"method"},{"location":"reference/#DiffusionGarnet.save_data_paraview-Tuple{Any}","page":"List of functions","title":"DiffusionGarnet.save_data_paraview","text":"save_data_paraview(integrator)\n\nCallback function used to save data to an HDF5 file and produce an XMDF file. XMDF files can be used to read data on software like Paraview.\n\nNote that data are saved in row major format. For column major, use save_data() instead.\n\n\n\n\n\n","category":"method"},{"location":"reference/#DiffusionGarnet.simulate","page":"List of functions","title":"DiffusionGarnet.simulate","text":"simulate(domain; path_save=nothing, solver=ROCK2(), kwargs...)\n\nSolve the coupled major element diffusion equations for a given domain using finite differences for the discretisation in space and return a solution type variable.\n\nThe 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.\n\nThe 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. The time of the solution is non-dimensional but can be converted back using the characteristic time (t_charact contained in the Domain structure).\n\npath_save is an optional argument, which can be used to define the path of the HDF5 output file. Default is to nothing.\n\nsolver is an optional argument, which can be used to define the solver to use for the time discretisation. Default is the ROCK2 method. All other ODE solvers accepted as the ones from DifferentialEquations (see https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/).\n\nAll 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/commonsolveropts/).\n\n\n\n\n\n","category":"function"},{"location":"reference/#DiffusionGarnet.simulate-Tuple{DiffusionGarnet.Domain1D}","page":"List of functions","title":"DiffusionGarnet.simulate","text":"simulate(domain::Domain1D; path_save=nothing, solver=ROCK2(), abstol=1e-8, reltol=1e-6, kwargs...)\n\nSolve the coupled major element diffusion equations in 1D. Save all timesteps in the output solution type variable by default.\n\n\n\n\n\n","category":"method"},{"location":"reference/#DiffusionGarnet.simulate-Tuple{DiffusionGarnet.Domain2D}","page":"List of functions","title":"DiffusionGarnet.simulate","text":"simulate(domain::Domain2D; path_save=nothing, solver=ROCK2(), kwargs...)\n\nSolve 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.\n\n\n\n\n\n","category":"method"},{"location":"reference/#DiffusionGarnet.simulate-Tuple{DiffusionGarnet.Domain3D}","page":"List of functions","title":"DiffusionGarnet.simulate","text":"simulate(domain::Domain3D; path_save=nothing, solver=ROCK2(), kwargs...)\n\nSolve the coupled major element diffusion equations in 3D. Save only the first and last timestep in the output solution type variable by default.\n\n\n\n\n\n","category":"method"},{"location":"reference/#DiffusionGarnet.simulate-Tuple{DiffusionGarnet.DomainSpherical}","page":"List of functions","title":"DiffusionGarnet.simulate","text":"simulate(domain::DomainSpherical; path_save=nothing, solver=ROCK2(), abstol=1e-8, reltol=1e-6, kwargs...)\n\nSolve the coupled major element diffusion equations in spherical coordinates. Save all timesteps in the output solution type variable by default.\n\n\n\n\n\n","category":"method"},{"location":"reference/#DiffusionGarnet.update_diffusion_coef-Tuple{Any}","page":"List of functions","title":"DiffusionGarnet.update_diffusion_coef","text":"update_diffusion_coef(integrator)\n\nCallback function to update the diffusion coefficients at a given time from a new pressure and temperature. To use with the callback PresetTimeCallback (https://docs.sciml.ai/stable/basics/callbacks/#PresetTimeCallback-1).\n\nFollows the syntax of callback functions defined by DiffEqCallbacks.jl (https://docs.sciml.ai/DiffEqCallbacks/stable/).\n\n\n\n\n\n","category":"method"},{"location":"#DiffusionGarnet.jl","page":"Home","title":"DiffusionGarnet.jl","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"This is the documentation page for DiffusionGarnet.jl, a high-level Julia package to do coupled diffusion modelling of major elements on real garnet data. It currently supports 1D, spherical, 2D and 3D coordinates for evenly spaced data.","category":"page"},{"location":"","page":"Home","title":"Home","text":"It is built on top of the DifferentialEquations.jl package ecosystem and uses Unitful.jl to allow the user to define appropriate units for their problems. For 2D and 3D models, it uses ParallelStencil.jl to support multithreading on CPU and parallel computing on GPU.","category":"page"},{"location":"","page":"Home","title":"Home","text":"This package can be used as a teaching tool or for research purposes.","category":"page"},{"location":"#Installation","page":"Home","title":"Installation","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"DiffusionGarnet.jl is a registered package and may be installed directly from the REPL:","category":"page"},{"location":"","page":"Home","title":"Home","text":"julia>]\n pkg> add DiffusionGarnet\n pkg> test DiffusionGarnet","category":"page"},{"location":"diffusion_1D/#1D_diffusion","page":"Diffusion in 1D Cartesian coordinates","title":"Diffusion in 1D Cartesian coordinates","text":"","category":"section"},{"location":"diffusion_1D/","page":"Diffusion in 1D Cartesian coordinates","title":"Diffusion in 1D Cartesian coordinates","text":"DiffusionGarnet expects the user to provide real natural data for modelling major element diffusion in garnet. Note that the profiles must be evenly spaced. A set of example data can be found in the repository of the package in the 1D examples section for 1D profile called Data_grt_1D.txt. This is what we will use for this tutorial.","category":"page"},{"location":"diffusion_1D/","page":"Diffusion in 1D Cartesian coordinates","title":"Diffusion in 1D Cartesian coordinates","text":"First, we will load the data, which should be in the same folder as your current session:","category":"page"},{"location":"diffusion_1D/","page":"Diffusion in 1D Cartesian coordinates","title":"Diffusion in 1D Cartesian coordinates","text":"using DiffusionGarnet # this can take a while\nusing DelimitedFiles\n# load the data of your choice (here from the text file located in https://github.com/Iddingsite/DiffusionGarnet.jl/tree/main/examples/1D, place it in the same folder as where you are running the code)\ndata = DelimitedFiles.readdlm(\"Data_Grt_1D.txt\", '\\t', '\\n', header=true)[1]\n\nMg0 = data[:, 4] # load initial Mg mole fraction\nFe0 = data[:, 2] # load initial Fe mole fraction\nMn0 = data[:, 3] # load initial Mn mole fraction\nCa0 = data[:, 5] # load initial Ca mole fraction\ndistance = data[:, 1]","category":"page"},{"location":"diffusion_1D/","page":"Diffusion in 1D Cartesian coordinates","title":"Diffusion in 1D Cartesian coordinates","text":"We can visualize our data:","category":"page"},{"location":"diffusion_1D/","page":"Diffusion in 1D Cartesian coordinates","title":"Diffusion in 1D Cartesian coordinates","text":"using Plots\n\nl = @layout [a ; b]\n\np1 = plot(distance, Fe0, label=\"Fe initial\", linestyle = :dash, linewidth=1, dpi=200, title = \"Initial conditions\", legend=:outerbottomright, linecolor=1,xlabel = \"Distance (µm)\", ylabel=\"Molar fraction\")\n\np2 = plot(distance, Mg0, label=\"Mg initial\", linestyle = :dash, linewidth=1, dpi=200,legend=:outerbottomright,linecolor=2,xlabel = \"Distance (µm)\")\np2 = plot!(distance, Mn0, label=\"Mn initial\", linestyle = :dash, linewidth=1, linecolor=3)\np2 = plot!(distance, Ca0, label=\"Ca initial\", linestyle = :dash, linewidth=1, linecolor=4, ylabel=\"Molar fraction\")\n\nplot(p1, p2, layout = l)","category":"page"},{"location":"diffusion_1D/","page":"Diffusion in 1D Cartesian coordinates","title":"Diffusion in 1D Cartesian coordinates","text":"which outputs:","category":"page"},{"location":"diffusion_1D/","page":"Diffusion in 1D Cartesian coordinates","title":"Diffusion in 1D Cartesian coordinates","text":"(Image: Initial conditions.)","category":"page"},{"location":"diffusion_1D/","page":"Diffusion in 1D Cartesian coordinates","title":"Diffusion in 1D Cartesian coordinates","text":"Then, we will define 2 structures from the constructors InitialConditions1D and Domain, which will contain all the information DiffusionGarnet needs to run a simulation.","category":"page"},{"location":"diffusion_1D/","page":"Diffusion in 1D Cartesian coordinates","title":"Diffusion in 1D Cartesian coordinates","text":"Lx = (data[end,1] - data[1,1])u\"µm\" # length in x of the model, here in µm\ntfinal = 15u\"Myr\" # total time of the model, here in Myr\n\n# define the initial conditions in 1D of your problem in that order.\nIC1D = InitialConditions1D(Mg0, Fe0, Mn0, Lx, tfinal)\n\n# define the pressure and temperature conditions of diffusion\nT = 900u\"°C\"\nP = 0.6u\"GPa\"\n\n# define a Domain struct containing the definition of our problem and nondimensionalised variables\ndomain1D = Domain(IC1D, T, P)","category":"page"},{"location":"diffusion_1D/","page":"Diffusion in 1D Cartesian coordinates","title":"Diffusion in 1D Cartesian coordinates","text":"note: Note\nLx, tfinal, T and P need to contain units, following the syntax of the package Unitful. This allows the user to specify the units that suit their problem.","category":"page"},{"location":"diffusion_1D/","page":"Diffusion in 1D Cartesian coordinates","title":"Diffusion in 1D Cartesian coordinates","text":"Domain1D contains all the information that DiffusionGarnet needs to solve our coupled diffusion problem, at 900 °C and 0.6 GPa for a duration of 15 Myr.","category":"page"},{"location":"diffusion_1D/","page":"Diffusion in 1D Cartesian coordinates","title":"Diffusion in 1D Cartesian coordinates","text":"This can be achieved with the function simulate():","category":"page"},{"location":"diffusion_1D/","page":"Diffusion in 1D Cartesian coordinates","title":"Diffusion in 1D Cartesian coordinates","text":"# solve the problem using DifferentialEquations.jl\nsol = simulate(domain1D)","category":"page"},{"location":"diffusion_1D/","page":"Diffusion in 1D Cartesian coordinates","title":"Diffusion in 1D Cartesian coordinates","text":"which outputs the time spent on the solver, for example, on the second run:","category":"page"},{"location":"diffusion_1D/","page":"Diffusion in 1D Cartesian coordinates","title":"Diffusion in 1D Cartesian coordinates","text":" 0.399870 seconds (31.93 k allocations: 18.212 MiB)","category":"page"},{"location":"diffusion_1D/","page":"Diffusion in 1D Cartesian coordinates","title":"Diffusion in 1D Cartesian coordinates","text":"simulate() uses the DifferentialEquations.jl 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. It basically holds all the information from our simulation.","category":"page"},{"location":"diffusion_1D/","page":"Diffusion in 1D Cartesian coordinates","title":"Diffusion in 1D Cartesian coordinates","text":"We can now plot the solution to our problem.","category":"page"},{"location":"diffusion_1D/","page":"Diffusion in 1D Cartesian coordinates","title":"Diffusion in 1D Cartesian coordinates","text":"# extract characteristic time to convert back to dimensional time\n@unpack t_charact = domain1D\n\nanim = @animate for i = LinRange(0, sol.t[end], 100)\n l = @layout [a ; b]\n\n 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)\")\n p1 = plot!(distance, sol(i)[:,2], label=\"Fe\",linecolor=1, linewidth=1)\n\n p2 = plot(distance, Mg0, label=\"Mg initial\", linestyle = :dash, linewidth=1, dpi=200,legend=:outerbottomright,linecolor=2,xlabel = \"Distance (µm)\")\n p2 = plot!(distance, Mn0, label=\"Mn initial\", linestyle = :dash, linewidth=1, linecolor=3)\n p2 = plot!(distance, Ca0, label=\"Ca initial\", linestyle = :dash, linewidth=1, linecolor=4)\n p2 = plot!(distance, sol(i)[:,1], label=\"Mg\",linecolor=2, linewidth=1)\n\n p2 = plot!(distance, sol(i)[:,3], label=\"Mn\", linecolor=3, linewidth=1)\n\n p2 = plot!(distance, 1 .- sol(i)[:,1] .- sol(i)[:,2] .- sol(i)[:,3], label=\"Ca\", linecolor=4, linewidth=1)\n\n plot(p1, p2, layout = l)\nend every 1\n\nprintln(\"Now, generating the gif...\")\ngif(anim, \"Grt_1D.gif\", fps = 7)\nprintln(\"...Done!\")","category":"page"},{"location":"diffusion_1D/","page":"Diffusion in 1D Cartesian coordinates","title":"Diffusion in 1D Cartesian coordinates","text":"Here is the resulting gif obtained:","category":"page"},{"location":"diffusion_1D/","page":"Diffusion in 1D Cartesian coordinates","title":"Diffusion in 1D Cartesian coordinates","text":"(Image: 1D diffusion profil of a garnet)","category":"page"},{"location":"diffusion_1D/","page":"Diffusion in 1D Cartesian coordinates","title":"Diffusion in 1D Cartesian coordinates","text":"It shows the compositional evolution of a 1D profile through a garnet grain with homogeneous Dirichlet boundaries on both sides.","category":"page"}] } diff --git a/dev/updating_PT/index.html b/dev/updating_PT/index.html index 1ce6b20..b81749b 100644 --- a/dev/updating_PT/index.html +++ b/dev/updating_PT/index.html @@ -51,4 +51,4 @@ println("Now, generating the gif...") gif(anim, "Grt_1D_var_TP.gif", fps = 7) println("...Done!") -

which outputs:

1D diffusion profil of a garnet with PT

We can clearly observe the slowing down of the diffusion processes as the PT conditions decrease with time.

+

which outputs:

1D diffusion profil of a garnet with PT

We can clearly observe the slowing down of the diffusion processes as the PT conditions decrease with time.