Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Clean up docs and shorten README #206

Merged
merged 13 commits into from
Jul 31, 2024
4 changes: 2 additions & 2 deletions CONTRIBUTING.md
Original file line number Diff line number Diff line change
@@ -1,12 +1,12 @@
# Contributing

JustRelax.jl is an open-source project and we are very happy to accept contributions
`JustRelax.jl` is an open-source project and we are very happy to accept contributions
from the community. Please feel free to open issues or submit patches (preferably
as pull requests) any time. For planned larger contributions, it is often
beneficial to get in contact with one of the principal developers first (see
[AUTHORS.md](AUTHORS.md)).

JustRelax.jl and its contributions are licensed under the MIT license. As a contributor, you certify that all your
`JustRelax.jl` and its contributions are licensed under the MIT license. As a contributor, you certify that all your
contributions are in conformance with the *Developer Certificate of Origin
(Version 1.1)*, which is reproduced below.

Expand Down
263 changes: 3 additions & 260 deletions README.md
Original file line number Diff line number Diff line change
Expand Up @@ -3,8 +3,9 @@
[![Dev](https://img.shields.io/badge/docs-dev-blue.svg)](https://ptsolvers.github.io/JustRelax.jl/dev/)
[![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.10212422.svg)](https://doi.org/10.5281/zenodo.10212422)
![CI](https://github.com/PTSolvers/JustRelax.jl/actions/workflows/ci.yml/badge.svg)
[![Build status](https://badge.buildkite.com/6b970b1066dc828a56a75bccc65a8bc896a8bb76012a61fe96.svg)](https://buildkite.com/julialang/justrelax-dot-jl)
[![Build status](https://badge.buildkite.com/6b970b1066dc828a56a75bccc65a8bc896a8bb76012a61fe96.svg?branch=main)](https://buildkite.com/julialang/justrelax-dot-jl)
[![codecov](https://codecov.io/gh/PTsolvers/JustRelax.jl/graph/badge.svg?token=4ZJO7ZGT8H)](https://codecov.io/gh/PTsolvers/JustRelax.jl)
[![License: MIT](https://img.shields.io/badge/License-MIT-yellow.svg)](https://opensource.org/licenses/MIT)

<p align="center"><img src="./docs/src/assets/logo.png" alt="JustRelax.jl" width="200"></p>

Expand Down Expand Up @@ -59,268 +60,10 @@ julia> ]
```
The test will take a while, so grab a :coffee: or :tea:

## Example: shear band localisation (2D)

![ShearBand2D](docs/src/assets/movies/DP_nx2058_2D.gif)

This example displays how the package can be used to simulate shear band localisation. The example is based on the [ShearBands2D.jl](miniapps/benchmarks/stokes2D/shear_band/ShearBand2D.jl).


### Load packages

Load JustRelax necessary modules and define backend.
```julia
using JustRelax, JustRelax.JustRelax2D, JustRelax.DataIO
const backend_JR = CPUBackend
```

For this specific example we use particles to define the material phases, for which we rely on [JustPIC.jl](https://github.com/JuliaGeodynamics/JustPIC.jl). As in `JustRelax.jl`, we need to set up the environment of `JustPIC.jl`. This is done by running/including the following commands:

```julia
using JustPIC
using JustPIC._2D

const backend = CPUBackend # Threads is the default backend
const backend = CUDABackend # to run on a CUDA GPU load CUDA.jl (i.e. "using CUDA") at the beginning of the script
const backend = AMDGPUBackend # and to run on an AMD GPU load AMDGPU.jl (i.e. "using AMDGPU") at the beginning of the script.
```

We will also use `ParallelStencil.jl` to write some device-agnostic helper functions:
```julia
using ParallelStencil
@init_parallel_stencil(Threads, Float64, 2) #or (CUDA, Float64, 2) or (AMDGPU, Float64, 2)
```
and will use [GeoParams.jl](https://github.com/JuliaGeodynamics/GeoParams.jl/tree/main) to define and compute physical properties of the materials:
```julia
using GeoParams
```

## Model setup


For the initial setup, you will need to specify the number of nodes in x- and y- direction `nx` and `ny` as well as the directory where the figures are stored (`figdir`). The initialisation of the global grid and MPI environment is done with `igg = IGG(init_global_grid(nx, ny, 0; init_MPI = true)...)`:

```julia
N = 128
n = N + 2
nx = n - 2
ny = n - 2
figdir = "ShearBands2D"
igg = IGG(init_global_grid(nx, ny, 0; init_MPI = true)...)
```

Initialisation of the physical domain and the grid. As `JustRelax.jl` relies on [ImplicitGlobalGrid.jl](https://github.com/omlins/ImplicitGlobalGrid.jl), the grid can be `MPIAWARE`. This makes it a global grid and the grid steps are automatically distributed over the MPI processes.

```julia
# Physical domain ------------------------------------
ly = 1.0 # domain length in y
lx = ly # domain length in x
ni = nx, ny # number of cells
li = lx, ly # domain length in x- and y-
di = @. li / ni # grid step in x- and -y
origin = 0.0, 0.0 # origin coordinates
grid = Geometry(ni, li; origin = origin)
(; xci, xvi) = grid # nodes at the center and vertices of the cells
dt = Inf
```

### Initialisation of the rheology with [GeoParams.jl](https://github.com/JuliaGeodynamics/GeoParams.jl).

The rheology can be tailored to the specific problem with different creep laws and material parameters or the miniapps in the [convection folder](miniapps/convection).

```julia
# Physical properties using GeoParams ----------------
τ_y = 1.6 # yield stress. If do_DP=true, τ_y stand for the cohesion: c*cos(ϕ)
ϕ = 30 # friction angle
C = τ_y # Cohesion
η0 = 1.0 # viscosity
G0 = 1.0 # elastic shear modulus
Gi = G0/(6.0-4.0) # elastic shear modulus perturbation
εbg = 1.0 # background strain-rate
η_reg = 8e-3 # regularisation "viscosity"
dt = η0/G0/4.0 # assumes Maxwell time of 4
el_bg = ConstantElasticity(; G=G0, Kb=4)
el_inc = ConstantElasticity(; G=Gi, Kb=4)
visc = LinearViscous(; η=η0)
pl = DruckerPrager_regularised(; # non-regularized plasticity
C = C,
ϕ = ϕ,
η_vp = η_reg,
Ψ = 0,
) # viscoplasticity model from e.g. Duretz, T., Räss, L., de Borst, R., & Hageman, T. (2023). A comparison of plasticity regularization approaches for geodynamic modeling. Geochemistry, Geophysics, Geosystems, 24, e2022GC010675. https://doi.org/101029/2022GC010675
rheology = (
# Matrix phase
SetMaterialParams(;
Phase = 1,
Density = ConstantDensity(; ρ = 0.0),
Gravity = ConstantGravity(; g = 0.0),
CompositeRheology = CompositeRheology((visc, el_bg, pl)),
Elasticity = el_bg,
),
# Inclusion phase
SetMaterialParams(;
Density = ConstantDensity(; ρ = 0.0),
Gravity = ConstantGravity(; g = 0.0),
CompositeRheology = CompositeRheology((visc, el_inc, pl)),
Elasticity = el_inc,
),
)
```

Since we have two material phases, we need their phase ratio at the cell centers. For that purpose, we define a helper `init_phases!`. This function is parallelized with the `@parallel_indices` macro from [ParallelStencil.jl](https://github.com/omlins/ParallelStencil.jl). In this case, this is the only function that needs to be tailored to the specific problem, everything else is handled by `JustRelax.jl` itself.

```julia
# Initialize phases on the particles
function init_phases!(phase_ratios, xci, radius)
ni = size(phase_ratios.center)
origin = 0.5, 0.5

@parallel_indices (i, j) function init_phases!(phases, xc, yc, o_x, o_y, radius)
x, y = xc[i], yc[j]
if ((x-o_x)^2 + (y-o_y)^2) > radius^2
JustRelax.@cell phases[1, i, j] = 1.0
JustRelax.@cell phases[2, i, j] = 0.0

else
JustRelax.@cell phases[1, i, j] = 0.0
JustRelax.@cell phases[2, i, j] = 1.0
end
return nothing
end

@parallel (JustRelax.@idx ni) init_phases!(phase_ratios.center, xci..., origin..., radius)
end
```

### Initialisation of the Stokes arrays and buoyancy forces

The rheology is computed with `compute_viscosity!` which is a function from `JustRelax.jl` and computes the viscosity according to the strain rate and the phase ratios.

```julia
# Initialize phase ratios -------------------------------
radius = 0.1
phase_ratios = PhaseRatio(backend_JR, ni, length(rheology))
init_phases!(phase_ratios, xci, radius)
# STOKES ---------------------------------------------
# Allocate arrays needed for every Stokes problem
stokes = StokesArrays(backend_JR, ni)
pt_stokes = PTStokesCoeffs(li, di; ϵ=1e-6, CFL = 0.75 / √2.1)
# PT coefficients after Räss, L., Utkin, I., Duretz, T., Omlin, S., and Podladchikov, Y. Y.: Assessing the robustness and scalability of the accelerated pseudo-transient method, Geosci. Model Dev., 15, 5757–5786, https://doi.org/10.5194/gmd-15-5757-2022, 2022.
# Buoyancy forces
ρg = @zeros(ni...), @zeros(ni...)
η = @ones(ni...)
args = (; T = thermal.Tc, P = stokes.P, dt = Inf)
compute_ρg!(ρg[2], phase_ratios, rheology, args)
compute_viscosity!(stokes, 1.0, phase_ratios, args, rheology, (-Inf, Inf))
```

### Boundary conditions

The boundary conditions for this example are pure sher velocity conditions with free slip in al the boundaries. conditions. The boundary conditions are applied with `flow_bcs!` and the halo is updated with `update_halo!`.

```julia
flow_bcs = VelocityBoundaryConditions(;
free_slip = (left = true, right = true, top = true, bot = true),
no_slip = (left = false, right = false, top = false, bot=false),
)
stokes.V.Vx .= PTArray([ x*εbg for x in xvi[1], _ in 1:ny+2])
stokes.V.Vy .= PTArray([-y*εbg for _ in 1:nx+2, y in xvi[2]])
flow_bcs!(stokes, flow_bcs) # apply boundary conditions
```

Alternatively, we can also use boundary conditions in terms of displacements:
```julia
# Boundary conditions
flow_bcs = DisplacementBoundaryConditions(;
free_slip = (left = true, right = true, top = true, bot = true),
no_slip = (left = false, right = false, top = false, bot=false),
)
stokes.U.Ux .= PTArray(backend)([ x*εbg*lx*dt for x in xvi[1], _ in 1:ny+2])
stokes.U.Uy .= PTArray(backend)([-y*εbg*ly*dt for _ in 1:nx+2, y in xvi[2]])
flow_bcs!(stokes, flow_bcs) # apply boundary conditions
displacement2velocity!(stokes, dt) # convert displacement to velocity
update_halo!(@velocity(stokes)...)
```

Finally we can build out time-steping iterations where we solve the pseudo-transient stokes equations with `solve!`. The visualisation is done with [GLMakie.jl](https://github.com/MakieOrg/Makie.jl).

```julia
take(figdir) # if it does not exist, make a folder where figures are stored

# Time loop
t, it, tmax = 0.0, 0, 3.5
τII, sol, ttot = Float64[], Float64[], Float64[]

while t < tmax
# Stokes solver ----------------
solve!(
stokes,
pt_stokes,
di,
flow_bcs,
ρg,
phase_ratios,
rheology,
args,
dt,
igg;
kwargs = (;
iterMax = 150e3,
nout = 200,
viscosity_cutoff = (-Inf, Inf),
verbose = true
)
)
# Compute second invariant of the strain rate tensor
tensor_invariant!(stokes.ε)
push!(τII, maximum(stokes.τ.xx))

it += 1
t += dt

analytical_solution = εbg * η * (1 - exp(-G0 * t / η0))
push!(sol, analytical_solution)
push!(ttot, t)

println("it = $it; t = $t \n")

# visualisation of high density inclusion
th = 0:pi/50:3*pi;
xunit = @. radius * cos(th) + 0.5;
yunit = @. radius * sin(th) + 0.5;

fig = Figure(size = (1600, 1600), title = "t = $t")
ax1 = Axis(fig[1,1], aspect = 1, title = L"\tau_{II}", titlesize=35)
ax2 = Axis(fig[2,1], aspect = 1, title = L"E_{II}", titlesize=35)
ax3 = Axis(fig[1,2], aspect = 1, title = L"\log_{10}(\varepsilon_{II})", titlesize=35)
ax4 = Axis(fig[2,2], aspect = 1)
heatmap!(ax1, xci..., Array(stokes.τ.II) , colormap=:batlow)
heatmap!(ax2, xci..., Array(log10.(stokes.EII_pl)) , colormap=:batlow)
heatmap!(ax3, xci..., Array(log10.(stokes.ε.II)) , colormap=:batlow)
lines!(ax2, xunit, yunit, color = :black, linewidth = 5)
lines!(ax4, ttot, τII, color = :black)
lines!(ax4, ttot, sol, color = :red)
hidexdecorations!(ax1)
hidexdecorations!(ax3)
save(joinpath(figdir, "$(it).png"), fig)
fig
end
```

## Miniapps

Currently there are 4 convection miniapps with particles and 4 corresponding miniapps without. The miniapps with particles are:

* [Layered_convection2D.jl](miniapps/convection/Particles2D/Layered_convection2D.jl)
* [Layered_convection2D_nonDim.jl](miniapps/convection/Particles2D_nonDim/Layered_convection2D.jl)
* [Layered_convection3D.jl](miniapps/convection/Particles3D/Layered_convection3D.jl)
* [WENO_convection2D.jl](miniapps/convection/WENO5/WENO_convection2D.jl)
Available miniapps can be found [here](miniapps) and will be updated regularly. The miniapps are designed to be simple and easy to understand, while still providing a good basis for more complex applications. The miniapps are designed to be run on a single node, but can be easily extended to run on multiple nodes using [ImplicitGlobalGrid.jl](https://github.com/omlins/ImplicitGlobalGrid.jl) and [MPI.jl](https://github.com/JuliaParallel/MPI.jl).
aelligp marked this conversation as resolved.
Show resolved Hide resolved

The miniapps without particles are:
* [GlobalConvection2D_Upwind.jl](miniapps/convection/GlobalConvection2D_Upwind.jl)
* [GlobalConvection2D_WENO5.jl](miniapps/convection/GlobalConvection2D_WENO5.jl)
* [GlobalConvection2D_WENO5_MPI.jl](miniapps/convection/GlobalConvection2D_WENO5_MPI.jl)
* [GlobalConvection3D_Upwind.jl](miniapps/convection/GlobalConvection3D_Upwind.jl)

## Benchmarks

Expand Down
2 changes: 1 addition & 1 deletion docs/src/index.md
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,7 @@ CurrentModule = JustRelax

Need to solve a very large multi-physics problem on many GPUs in parallel? Just Relax!

`JustRelax` is a collection of accelerated iterative pseudo-transient solvers using MPI and multiple CPU or GPU backends. It's part of the [PTSolvers organisation](https://ptsolvers.github.io) and
`JustRelax.jl` is a collection of accelerated iterative pseudo-transient solvers using MPI and multiple CPU or GPU backends. It's part of the [PTSolvers organisation](https://ptsolvers.github.io) and
developed within the [GPU4GEO project](https://www.pasc-ch.org/projects/2021-2024/gpu4geo/). Current publications, outreach and news can be found on the [GPU4GEO website](https://ptsolvers.github.io/GPU4GEO/).

The package relies on other packages as building blocks and parallelisation tools:
Expand Down
8 changes: 4 additions & 4 deletions docs/src/man/Blankenbach.md
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,7 @@ Thermal convection benchmark from [Blankenbach et al., 1989](https://academic.o

# Initialize packages

Load JustRelax necessary modules and define backend.
Load `JustRelax.jl` necessary modules and define backend.
```julia
using JustRelax, JustRelax.JustRelax2D, JustRelax.DataIO
const backend_JR = CPUBackend
Expand All @@ -16,7 +16,7 @@ using JustPIC, JustPIC._2D
const backend = CPUBackend # Options: CPUBackend, CUDABackend, AMDGPUBackend
```

We will also use `ParallelStencil.jl` to write some device-agnostic helper functions:
We will also use [ParallelStencil.jl](https://github.com/omlins/ParallelStencil.jl) to write some device-agnostic helper functions:
```julia
using ParallelStencil
@init_parallel_stencil(Threads, Float64, 2) #or (CUDA, Float64, 2) or (AMDGPU, Float64, 2)
Expand Down Expand Up @@ -88,7 +88,7 @@ yc_anomaly = 1 / 3 # origin of thermal anomaly
r_anomaly = 0.1 / 2 # radius of perturbation
```

Helper function to initialize material phases with `ParallelStencil.jl`
Helper function to initialize material phases with [ParallelStencil.jl](https://github.com/omlins/ParallelStencil.jl)
```julia
function init_phases!(phases, particles)
ni = size(phases)
Expand Down Expand Up @@ -138,7 +138,7 @@ thermal = ThermalArrays(backend_JR, ni)

## Initialize thermal profile and viscosity fields

To initialize the thermal profile we use `ParallelStencil.jl` again
To initialize the thermal profile we use [ParallelStencil.jl](https://github.com/omlins/ParallelStencil.jl) again
```julia
@parallel_indices (i, j) function init_T!(T, y)
T[i, j] = 1 - y[j]
Expand Down
10 changes: 5 additions & 5 deletions docs/src/man/ShearBands.md
Original file line number Diff line number Diff line change
@@ -1,17 +1,17 @@
# ShearBand benchmark

Shear Band benchmark to test the visco-elasto-plastic rheology implementation in JustRelax.jl
Shear Band benchmark to test the visco-elasto-plastic rheology implementation in `JustRelax.jl`

# Initialize packages

Load JustRelax necessary modules and define backend.
Load `JustRelax.jl` necessary modules and define backend.
```julia
using JustRelax, JustRelax.JustRelax2D, JustRelax.DataIO
const backend_JR = CPUBackend
```


We will also use `ParallelStencil.jl` to write some device-agnostic helper functions:
We will also use [ParallelStencil.jl](https://github.com/omlins/ParallelStencil.jl) to write some device-agnostic helper functions:
```julia
using ParallelStencil
@init_parallel_stencil(Threads, Float64, 2) #or (CUDA, Float64, 2) or (AMDGPU, Float64, 2)
Expand Down Expand Up @@ -85,7 +85,7 @@ pl = DruckerPrager_regularised(; # non-regularized plasticity

# Phase anomaly

Helper function to initialize material phases with `ParallelStencil.jl`
Helper function to initialize material phases with [ParallelStencil.jl](https://github.com/omlins/ParallelStencil.jl)
```julia
function init_phases!(phase_ratios, xci, radius)
ni = size(phase_ratios.center)
Expand Down Expand Up @@ -200,7 +200,7 @@ push!(sol, solution(εbg, t, G0, η0))
push!(ttot, t)
```
# Visualization
We will use `Makie.jl` to visualize the results
We will use [Makie.jl](https://github.com/MakieOrg/Makie.jl) to visualize the results
```julia
using GLMakie
```
Expand Down
4 changes: 2 additions & 2 deletions docs/src/man/advection.md
Original file line number Diff line number Diff line change
@@ -1,8 +1,8 @@
# Field advection

## Particles-in-Cell
`JustRelax` relies on [JustPIC.jl](https://github.com/JuliaGeodynamics/JustPIC.jl) for advections of particles containing material information.
`JustRelax.jl` relies on [JustPIC.jl](https://github.com/JuliaGeodynamics/JustPIC.jl) for advections of particles containing material information.

## Upwind

## WENO5
## WENO5
2 changes: 1 addition & 1 deletion docs/src/man/backend.md
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
# Selecting the backend

JustRelax supports three backends: the default CPU backend, and two GPU backends for Nvidia and AMD GPUs. The default CPU backend is selected upon loading JustRelax:
`JustRelax.jl` supports three backends: the default CPU backend, and two GPU backends for Nvidia and AMD GPUs. The default CPU backend is selected upon loading JustRelax:

```julia
using JustRelax
Expand Down
Loading