Merge pull request #188 from PTsolvers/pa-BC

Add DisplacementBoundaryConditions

README.md

@@ -214,17 +214,29 @@ compute_ρg!(ρg[2], phase_ratios, rheology, args)
 compute_viscosity!(stokes, 1.0, phase_ratios, args, rheology, (-Inf, Inf))
 ```
 
-Define pure shear velocity boundary conditions
+Define pure shear velocity boundary conditions or displacement boundary conditions. The boundary conditions are applied with `flow_bcs!` and the halo is updated with `update_halo!`.
 ```julia
 # Boundary conditions
-flow_bcs     = FlowBoundaryConditions(;
+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
-update_halo!(stokes.V.Vx, stokes.V.Vy)
+update_halo!(@velocity(stokes)...)
+```
+```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)...)
 ```
 Pseudo-transient Stokes solver and visualisation of the results. The visualisation is done with [GLMakie.jl](https://github.com/MakieOrg/Makie.jl).
 

docs/make.jl

@@ -7,6 +7,9 @@ JR_root_dir = dirname(@__DIR__)
 license = read(joinpath(JR_root_dir, "LICENSE.md"), String)
 write(joinpath(@__DIR__, "src", "man", "license.md"), license)
 
+security = read(joinpath(JR_root_dir, "SECURITY.md"), String)
+write(joinpath(@__DIR__, "src", "man", "security.md"), security)
+
 # Copy list of authors to not need to synchronize it manually
 authors_text = read(joinpath(JR_root_dir, "AUTHORS.md"), String)
 # authors_text = replace(authors_text, "in the [LICENSE.md](LICENSE.md) file" => "under [License](@ref)")
@@ -94,6 +97,7 @@ makedocs(;
         "Authors" => "man/authors.md",
         "Contributing" => "man/contributing.md",
         "Code of Conduct" => "man/code_of_conduct.md",
+        "Security" => "man/security.md",
         "License" => "man/license.md"
     ],
 )

docs/src/man/Blankenbach.md

@@ -181,7 +181,7 @@ where `(-Inf, Inf)` is the viscosity cutoff.
 
 ## Boundary conditions
 ```julia
-flow_bcs      = FlowBoundaryConditions(;
+flow_bcs      = VelocityBoundaryConditions(;
     free_slip = (left = true, right=true, top=true, bot=true),
 )
 thermal_bc    = TemperatureBoundaryConditions(;

docs/src/man/ShearBands.md

@@ -136,14 +136,14 @@ where `(-Inf, Inf)` is the viscosity cutoff.
 
 ## Boundary conditions
 ```julia
-    flow_bcs     = FlowBoundaryConditions(;
+    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
-    update_halo!(stokes.V.Vx, stokes.V.Vy)
+    update_halo!(@velocity(stokes)...)
 
 ```
 

docs/src/man/boundary_conditions.md

@@ -14,12 +14,21 @@ Supported boundary conditions:
 
     $\sigma_z = 0 \rightarrow \tau_z = P$ at the top boundary
 
-## Defining the boundary contions
-Information regarding flow boundary conditions is defined in the `FlowBoundaryConditions` object. They can be switched on and off by setting them as `true` or `false` at the appropriate boundaries. Valid boundary names are `left` and `right`, `top` and `bot`, and for the 3D case, `front` and `back`. 
+## Defining the boundary conditions
+We have two ways of defining the boundary condition formulations:
+    - `VelocityBoundaryConditions`, and
+    - `DisplacementBoundaryConditions`.
+The first one is used for the velocity-pressure formulation, and the second one is used for the displacement-pressure formulation. The flow boundary conditions can be switched on and off by setting them as `true` or `false` at the appropriate boundaries. Valid boundary names are `left` and `right`, `top` and `bot`, and for the 3D case, `front` and `back`.
 
-For example, if we want to have free free-slip in every single boundary in a 2D simulation, we need to instantiate `FlowBoundaryConditions` as:
+
+For example, if we want to have free free-slip in every single boundary in a 2D simulation, we need to instantiate `VelocityBoundaryConditions` or `DisplacementBoundaryConditions` as:
 ```julia
-bcs = FlowBoundaryConditions(;
+bcs = VelocityBoundaryConditions(;
+    no_slip      = (left=false, right=false, top=false, bot=false),
+    free_slip    = (left=true, right=true, top=true, bot=true),
+    free_surface = false
+)
+bcs = DisplacementBoundaryConditions(;
     no_slip      = (left=false, right=false, top=false, bot=false),
     free_slip    = (left=true, right=true, top=true, bot=true),
     free_surface = false
@@ -28,9 +37,42 @@ bcs = FlowBoundaryConditions(;
 
 The equivalent for the 3D case would be:
 ```julia
-bcs = FlowBoundaryConditions(;
+bcs = VelocityBoundaryConditions(;
     no_slip      = (left=false, right=false, top=false, bot=false, front=false, back=false),
     free_slip    = (left=true, right=true, top=true, bot=true, front=true, back=true),
     free_surface = false
 )
-```
+bcs = DisplacementBoundaryConditions(;
+    no_slip      = (left=false, right=false, top=false, bot=false, front=false, back=false),
+    free_slip    = (left=true, right=true, top=true, bot=true, front=true, back=true),
+    free_surface = false
+)
+```
+## Prescribing the velocity/displacement boundary conditions
+Normally, one would prescribe the velocity/displacement boundary conditions by setting the velocity/displacement field at the boundary through the application of a background strain rate `εbg`.
+Depending on the formulation, the velocity/displacement field is set as follows for the 2D case:
+### Velocity formulation
+```julia
+stokes.V.Vx .= PTArray(backend)([ x*εbg for x in xvi[1], _ in 1:ny+2]) # Velocity in x direction
+stokes.V.Vy .= PTArray(backend)([-y*εbg for _ in 1:nx+2, y in xvi[2]]) # Velocity in y direction
+```
+Make sure to apply the set velocity to the boundary conditions. You do this by calling the `flow_bcs!` function,
+```julia
+flow_bcs!(stokes, flow_bcs)
+```
+and then applying the velocities to the halo
+```julia
+update_halo!(@velocity(stokes)...)
+```
+### Displacement formulation
+```julia
+stokes.U.Ux .= PTArray(backend)([ x*εbg*lx*dt for x in xvi[1], _ in 1:ny+2]) # Displacement in x direction
+stokes.U.Uy .= PTArray(backend)([-y*εbg*ly*dt for _ in 1:nx+2, y in xvi[2]]) # Displacement in y direction
+flow_bcs!(stokes, flow_bcs)
+```
+Make sure to initialize the displacement according to the extent of your domain. Here, lx and ly are the domain lengths in the x and y directions, respectively.
+Also for the displacement formulation it is important that the displacement is converted to velocity before updating the halo. This can be done by calling the `displacement2velocity!` function.
+```julia
+displacement2velocity!(stokes, dt) # convert displacement to velocity
+update_halo!(@velocity(stokes)...)
+```

docs/src/man/subduction2D/subduction2D.md

@@ -124,7 +124,7 @@ compute_viscosity!(stokes, phase_ratios, args0, rheology, viscosity_cutoff)
 We we will use free slip boundary conditions on all sides
 ```julia
 # Boundary conditions
-flow_bcs         = FlowBoundaryConditions(;
+flow_bcs         = VelocityBoundaryConditions(;
     free_slip    = (left = true , right = true , top = true , bot = true),
 )
 ```

miniapps/benchmarks/stokes2D/Blankenbach2D/Benchmark2D_WENO5.jl

@@ -133,11 +133,11 @@ function main2D(igg; ar=1, nx=32, ny=32, nit = 1e1, figdir="figs2D", do_vtk =fal
     )
 
     # Boundary conditions -------------------------------
-    flow_bcs         = FlowBoundaryConditions(;
+    flow_bcs         = VelocityBoundaryConditions(;
         free_slip    = (left = true, right=true, top=true, bot=true),
     )
     flow_bcs!(stokes, flow_bcs) # apply boundary conditions
-    update_halo!(stokes.V.Vx, stokes.V.Vy)
+    update_halo!(@velocity(stokes)...)
 
     # IO ------------------------------------------------
     # if it does not exist, make folder where figures are stored

miniapps/benchmarks/stokes2D/Blankenbach2D/Benchmark2D_sgd.jl

@@ -134,11 +134,11 @@ function main2D(igg; ar=1, nx=32, ny=32, nit = 1e1, figdir="figs2D", do_vtk =fal
     )
 
     # Boundary conditions -------------------------------
-    flow_bcs         = FlowBoundaryConditions(;
+    flow_bcs         = VelocityBoundaryConditions(;
         free_slip    = (left = true, right=true, top=true, bot=true),
     )
     flow_bcs!(stokes, flow_bcs) # apply boundary conditions
-    update_halo!(stokes.V.Vx, stokes.V.Vy)
+    update_halo!(@velocity(stokes)...)
 
     # IO ------------------------------------------------
     # if it does not exist, make folder where figures are stored

miniapps/benchmarks/stokes2D/Blankenbach2D/Benchmark2D_sgd_scaled.jl

@@ -124,11 +124,11 @@ function main2D(igg; ar=1, nx=32, ny=32, nit = 1e1, figdir="figs2D", do_vtk =fal
     )
 
     # Boundary conditions -------------------------------
-    flow_bcs         = FlowBoundaryConditions(;
+    flow_bcs         = VelocityBoundaryConditions(;
         free_slip    = (left = true, right=true, top=true, bot=true),
     )
     flow_bcs!(stokes, flow_bcs) # apply boundary conditions
-    update_halo!(stokes.V.Vx, stokes.V.Vy)
+    update_halo!(@velocity(stokes)...)
 
     # IO ------------------------------------------------
     # if it does not exist, make folder where figures are stored

miniapps/benchmarks/stokes2D/VanKeken.jl/VanKeken.jl

@@ -104,12 +104,12 @@ function main2D(igg; ny=64, nx=64, figdir="model_figs")
     compute_viscosity!(stokes, phase_ratios, args, rheology, (-Inf, Inf))
 
     # Boundary conditions
-    flow_bcs             = FlowBoundaryConditions(;
+    flow_bcs             = VelocityBoundaryConditions(;
         free_slip = (left =  true, right =  true, top = false, bot = false),
         no_slip   = (left = false, right = false, top =  true, bot =  true),
     )
     flow_bcs!(stokes, flow_bcs)
-    update_halo!(stokes.V.Vx, stokes.V.Vy)
+    update_halo!(@velocity(stokes)...)
 
     # IO ----- -------------------------------------------
     # if it does not exist, make folder where figures are stored

miniapps/benchmarks/stokes2D/elastic_buildup/Elastic_BuildUp.jl

@@ -57,11 +57,11 @@ function elastic_buildup(;
 
     ## Boundary conditions
     pureshear_bc!(stokes, xci, xvi, εbg, backend)
-    flow_bcs  = FlowBoundaryConditions(;
+    flow_bcs  = VelocityBoundaryConditions(;
         free_slip = (left = true, right = true, top = true, bot = true)
     )
     flow_bcs!(stokes, flow_bcs)
-    update_halo!(stokes.V.Vx, stokes.V.Vy)
+    update_halo!(@velocity(stokes)...)
 
     # Physical time loop
     t         = 0.0

miniapps/benchmarks/stokes2D/elastic_buildup/Elastic_BuildUp_phases_incompressible.jl

@@ -73,14 +73,14 @@ function main(igg; nx=64, ny=64, figdir="model_figs")
     compute_viscosity!(stokes, phase_ratios, args, rheology, (-Inf, Inf))
 
     # Boundary conditions
-    flow_bcs = FlowBoundaryConditions(;
+    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(backend_JR)([x * εbg for x in xvi[1], _ in 1:(ny + 2)])
     stokes.V.Vy .= PTArray(backend_JR)([-y * εbg for _ in 1:(nx + 2), y in xvi[2]])
     flow_bcs!(stokes, flow_bcs) # apply boundary conditions
-    update_halo!(stokes.V.Vx, stokes.V.Vy)
+    update_halo!(@velocity(stokes)...)
 
     # IO ------------------------------------------------
     # if it does not exist, make folder where figures are stored

miniapps/benchmarks/stokes2D/free_surface_stabilization/PlumeFreeSurface_2D.jl

@@ -141,7 +141,7 @@ function main(igg, nx, ny)
     compute_viscosity!(stokes, phase_ratios, args, rheology, (-Inf, Inf))
 
     # Boundary conditions
-    flow_bcs         = FlowBoundaryConditions(;
+    flow_bcs         = VelocityBoundaryConditions(;
         free_slip    = (left = true, right = true, top = true, bot = true),
         free_surface = true
     )

miniapps/benchmarks/stokes2D/free_surface_stabilization/RayleighTaylor2D.jl

@@ -139,7 +139,7 @@ function RT_2D(igg, nx, ny)
     compute_viscosity!(stokes, phase_ratios, args, rheology, (-Inf, Inf))
 
     # Boundary conditions
-    flow_bcs         = FlowBoundaryConditions(;
+    flow_bcs         = VelocityBoundaryConditions(;
         free_slip    = (left =  true, right =  true, top =  true, bot = false),
         no_slip      = (left = false, right = false, top = false, bot =  true),
         free_surface = true,

miniapps/benchmarks/stokes2D/shear_band/ShearBand2D.jl

@@ -1,4 +1,4 @@
-using GeoParams, GLMakie, CellArrays
+using GeoParams, CairoMakie, CellArrays
 using JustRelax, JustRelax.JustRelax2D
 using ParallelStencil
 @init_parallel_stencil(Threads, Float64, 2)
@@ -82,15 +82,15 @@ function main(igg; nx=64, ny=64, figdir="model_figs")
             Elasticity        = el_inc,
         ),
     )
-    
+
     # perturbation array for the cohesion
     perturbation_C = @rand(ni...)
 
     # Initialize phase ratios -------------------------------
     radius       = 0.1
-    phase_ratios = PhaseRatio(ni, length(rheology))
+    phase_ratios = PhaseRatio(backend, ni, length(rheology))
     init_phases!(phase_ratios, xci, radius)
-    
+
     # STOKES ---------------------------------------------
     # Allocate arrays needed for every Stokes problem
     stokes    = StokesArrays(backend, ni)
@@ -105,14 +105,14 @@ function main(igg; nx=64, ny=64, figdir="model_figs")
         stokes, phase_ratios, args, rheology, (-Inf, Inf)
     )
     # Boundary conditions
-    flow_bcs     = FlowBoundaryConditions(;
+    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(backend)([ x*εbg for x in xvi[1], _ in 1:ny+2])
     stokes.V.Vy .= PTArray(backend)([-y*εbg for _ in 1:nx+2, y in xvi[2]])
     flow_bcs!(stokes, flow_bcs) # apply boundary conditions
-    update_halo!(stokes.V.Vx, stokes.V.Vy)
+    update_halo!(@velocity(stokes)...)
 
     # IO -------------------------------------------------
     take(figdir)

miniapps/benchmarks/stokes2D/shear_band/ShearBand2D_MPI.jl

@@ -101,14 +101,14 @@ function main(igg; nx=64, ny=64, figdir="model_figs")
     )
 
     # Boundary conditions
-    flow_bcs      = FlowBoundaryConditions(;
+    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
-    update_halo!(stokes.V.Vx, stokes.V.Vy)
+    update_halo!(@velocity(stokes)...)
 
     # IO ------------------------------------------------
     # if it does not exist, make folder where figures are stored