From 059156f0fe82ec2d6f3f9747997b6ab8831ff956 Mon Sep 17 00:00:00 2001 From: Pascal Aellig Date: Wed, 3 Jul 2024 16:52:55 +0200 Subject: [PATCH 01/20] adding Displacement formulation for BC's --- src/JustRelax.jl | 6 ++++- src/JustRelax_CPU.jl | 12 ++++++++-- src/Utils.jl | 15 ++++++++++++ src/boundaryconditions/BoundaryConditions.jl | 15 ++++++++++-- src/boundaryconditions/free_surface.jl | 2 +- src/boundaryconditions/types.jl | 23 ++++++++++++++++--- src/common.jl | 13 ++++++++--- src/ext/AMDGPU/2D.jl | 18 ++++++++++++--- src/ext/AMDGPU/3D.jl | 18 ++++++++++++--- src/ext/CUDA/2D.jl | 18 ++++++++++++--- src/ext/CUDA/3D.jl | 20 +++++++++++++--- src/stokes/Stokes2D.jl | 24 ++++++++++++++++---- src/stokes/Stokes3D.jl | 20 ++++++++++++---- 13 files changed, 171 insertions(+), 33 deletions(-) diff --git a/src/JustRelax.jl b/src/JustRelax.jl index 1d39792b..af675bc0 100644 --- a/src/JustRelax.jl +++ b/src/JustRelax.jl @@ -32,7 +32,11 @@ include("types/phases.jl") # export PhaseRatio include("boundaryconditions/types.jl") -export TemperatureBoundaryConditions, FlowBoundaryConditions +export AbstractBoundaryConditions, + TemperatureBoundaryConditions, + AbstractFlowBoundaryConditions, + DisplacementBoundaryConditions, + VelocityBoundaryConditions include("types/traits.jl") export BackendTrait, CPUBackendTrait, NonCPUBackendTrait diff --git a/src/JustRelax_CPU.jl b/src/JustRelax_CPU.jl index a75d0970..f2d39455 100644 --- a/src/JustRelax_CPU.jl +++ b/src/JustRelax_CPU.jl @@ -11,7 +11,11 @@ using MPI import JustRelax: IGG, BackendTrait, CPUBackendTrait, backend, CPUBackend import JustRelax: PTStokesCoeffs import JustRelax: - AbstractBoundaryConditions, TemperatureBoundaryConditions, FlowBoundaryConditions + AbstractBoundaryConditions, + TemperatureBoundaryConditions, + AbstractFlowBoundaryConditions, + DisplacementBoundaryConditions, + VelocityBoundaryConditions @init_parallel_stencil(Threads, Float64, 2) @@ -34,7 +38,11 @@ using MPI import JustRelax: IGG, BackendTrait, CPUBackendTrait, backend, CPUBackend import JustRelax: PTStokesCoeffs import JustRelax: - AbstractBoundaryConditions, TemperatureBoundaryConditions, FlowBoundaryConditions + AbstractBoundaryConditions, + TemperatureBoundaryConditions, + AbstractFlowBoundaryConditions, + DisplacementBoundaryConditions, + VelocityBoundaryConditions @init_parallel_stencil(Threads, Float64, 3) diff --git a/src/Utils.jl b/src/Utils.jl index 2be50507..25f657ac 100644 --- a/src/Utils.jl +++ b/src/Utils.jl @@ -104,6 +104,21 @@ end @inline unpack_velocity(V::JustRelax.Velocity{<:AbstractArray{T,3}}) where {T} = V.Vx, V.Vy, V.Vz +""" + @displacement(U) + +Unpacks the displacement arrays `U` from the StokesArrays `A`. +""" +macro displacement(A) + return quote + unpack_displacement(($(esc(A))).U) + end +end + +@inline unpack_displacement(U::JustRelax.Displacement{<:AbstractArray{T,2}}) where {T} = U.Ux, U.Uy +@inline unpack_displacement(U::JustRelax.Displacement{<:AbstractArray{T,3}}) where {T} = + U.Ux, U.Uy, U.Uz + """ @qT(V) diff --git a/src/boundaryconditions/BoundaryConditions.jl b/src/boundaryconditions/BoundaryConditions.jl index 81968a6a..88116e00 100644 --- a/src/boundaryconditions/BoundaryConditions.jl +++ b/src/boundaryconditions/BoundaryConditions.jl @@ -38,14 +38,25 @@ end # @inline thermal_bcs!(thermal::ThermalArrays, bcs::TemperatureBoundaryConditions) = thermal_bcs!(thermal.T, bcs) """ - flow_bcs!(stokes, bcs::FlowBoundaryConditions, di) + flow_bcs!(stokes, bcs::VelocityBoundaryConditions) Apply the prescribed flow boundary conditions `bc` on the `stokes` """ flow_bcs!(stokes, bcs) = flow_bcs!(backend(stokes), stokes, bcs) -function flow_bcs!(::CPUBackendTrait, stokes, bcs) + +function flow_bcs!(::CPUBackendTrait, stokes, bcs::VelocityBoundaryConditions) return _flow_bcs!(bcs, @velocity(stokes)) end + +""" + flow_bcs!(stokes, bcs::DisplacementBoundaryConditions) + +Apply the prescribed flow boundary conditions `bc` on the `stokes` +""" +function flow_bcs!(::CPUBackendTrait, stokes, bcs::DisplacementBoundaryConditions) + return _flow_bcs!(bcs, @displacement(stokes)) +end + function flow_bcs!(bcs, V::Vararg{T,N}) where {T,N} return _flow_bcs!(bcs, tuple(V...)) end diff --git a/src/boundaryconditions/free_surface.jl b/src/boundaryconditions/free_surface.jl index 1e4ff780..bee99f69 100644 --- a/src/boundaryconditions/free_surface.jl +++ b/src/boundaryconditions/free_surface.jl @@ -1,5 +1,5 @@ function free_surface_bcs!( - stokes, bcs::FlowBoundaryConditions, η, rheology, phase_ratios, dt, di + stokes, bcs::AbstractFlowBoundaryConditions, η, rheology, phase_ratios, dt, di ) indices_range(::Any, Vy) = @idx (size(Vy, 1) - 2) indices_range(::Any, ::Any, Vz) = @idx (size(Vz, 1) - 2, size(Vz, 2) - 2) diff --git a/src/boundaryconditions/types.jl b/src/boundaryconditions/types.jl index 2f8bdab6..fbb2c0a5 100644 --- a/src/boundaryconditions/types.jl +++ b/src/boundaryconditions/types.jl @@ -1,5 +1,5 @@ abstract type AbstractBoundaryConditions end - +abstract type AbstractFlowBoundaryConditions <: AbstractBoundaryConditions end struct TemperatureBoundaryConditions{T,nD} <: AbstractBoundaryConditions no_flux::T @@ -11,12 +11,29 @@ struct TemperatureBoundaryConditions{T,nD} <: AbstractBoundaryConditions end end -struct FlowBoundaryConditions{T,nD} <: AbstractBoundaryConditions +struct DisplacementBoundaryConditions{T,nD} <: AbstractFlowBoundaryConditions + no_slip::T + free_slip::T + free_surface::Bool + + function DisplacementBoundaryConditions(; + no_slip::T=(left=false, right=false, top=false, bot=false), + free_slip::T=(left=true, right=true, top=true, bot=true), + free_surface::Bool=false, + ) where {T} + @assert length(no_slip) === length(free_slip) + check_flow_bcs(no_slip, free_slip) + + nD = length(no_slip) == 4 ? 2 : 3 + return new{T,nD}(no_slip, free_slip, free_surface) + end +end +struct VelocityBoundaryConditions{T,nD} <: AbstractFlowBoundaryConditions no_slip::T free_slip::T free_surface::Bool - function FlowBoundaryConditions(; + function VelocityBoundaryConditions(; no_slip::T=(left=false, right=false, top=false, bot=false), free_slip::T=(left=true, right=true, top=true, bot=true), free_surface::Bool=false, diff --git a/src/common.jl b/src/common.jl index 45117ee1..e8b68e25 100644 --- a/src/common.jl +++ b/src/common.jl @@ -32,8 +32,15 @@ include("types/displacement.jl") export velocity2displacement!, displacement2velocity! include("boundaryconditions/BoundaryConditions.jl") -export FlowBoundaryConditions, - TemperatureBoundaryConditions, flow_bcs!, thermal_bcs!, pureshear_bc!, apply_free_slip! +export AbstractBoundaryConditions, + TemperatureBoundaryConditions, + AbstractFlowBoundaryConditions, + DisplacementBoundaryConditions, + VelocityBoundaryConditions, + flow_bcs!, + thermal_bcs!, + pureshear_bc!, + apply_free_slip! include("MiniKernels.jl") @@ -58,7 +65,7 @@ export vertex2center!, center2vertex!, temperature2center!, velocity2vertex! include("advection/weno5.jl") export WENO5, WENO_advection! -# Stokes +# Stokes include("rheology/GeoParams.jl") include("rheology/StressUpdate.jl") diff --git a/src/ext/AMDGPU/2D.jl b/src/ext/AMDGPU/2D.jl index 29159d42..f68e20ed 100644 --- a/src/ext/AMDGPU/2D.jl +++ b/src/ext/AMDGPU/2D.jl @@ -22,7 +22,11 @@ import JustRelax: Geometry, @cell import JustRelax: - AbstractBoundaryConditions, TemperatureBoundaryConditions, FlowBoundaryConditions + AbstractBoundaryConditions, + TemperatureBoundaryConditions, + AbstractFlowBoundaryConditions, + DisplacementBoundaryConditions, + VelocityBoundaryConditions @init_parallel_stencil(AMDGPU, Float64, 2) @@ -94,14 +98,22 @@ function JR2D.PTThermalCoeffs( end # Boundary conditions -function JR2D.flow_bcs!(::AMDGPUBackendTrait, stokes::JustRelax.StokesArrays, bcs) +function JR2D.flow_bcs!(::AMDGPUBackendTrait, stokes::JustRelax.StokesArrays, bcs::VelocityBoundaryConditions) return _flow_bcs!(bcs, @velocity(stokes)) end -function flow_bcs!(::AMDGPUBackendTrait, stokes::JustRelax.StokesArrays, bcs) +function flow_bcs!(::AMDGPUBackendTrait, stokes::JustRelax.StokesArrays, bcs::VelocityBoundaryConditions) return _flow_bcs!(bcs, @velocity(stokes)) end +function JR2D.flow_bcs!(::AMDGPUBackendTrait, stokes::JustRelax.StokesArrays, bcs::DisplacementBoundaryConditions) + return _flow_bcs!(bcs, @displacement(stokes)) +end + +function flow_bcs!(::AMDGPUBackendTrait, stokes::JustRelax.StokesArrays, bcs::DisplacementBoundaryConditions) + return _flow_bcs!(bcs, @displacement(stokes)) +end + function JR2D.thermal_bcs!(::AMDGPUBackendTrait, thermal::JustRelax.ThermalArrays, bcs) return thermal_bcs!(thermal.T, bcs) end diff --git a/src/ext/AMDGPU/3D.jl b/src/ext/AMDGPU/3D.jl index 4f2422f3..8464ee84 100644 --- a/src/ext/AMDGPU/3D.jl +++ b/src/ext/AMDGPU/3D.jl @@ -22,7 +22,11 @@ import JustRelax: Geometry, @cell import JustRelax: - AbstractBoundaryConditions, TemperatureBoundaryConditions, FlowBoundaryConditions +AbstractBoundaryConditions, +TemperatureBoundaryConditions, +AbstractFlowBoundaryConditions, +DisplacementBoundaryConditions, +VelocityBoundaryConditions @init_parallel_stencil(AMDGPU, Float64, 3) @@ -94,14 +98,22 @@ function JR3D.PTThermalCoeffs( end # Boundary conditions -function JR3D.flow_bcs!(::AMDGPUBackendTrait, stokes::JustRelax.StokesArrays, bcs) +function JR3D.flow_bcs!(::AMDGPUBackendTrait, stokes::JustRelax.StokesArrays, bcs::VelocityBoundaryConditions) return _flow_bcs!(bcs, @velocity(stokes)) end -function flow_bcs!(::AMDGPUBackendTrait, stokes::JustRelax.StokesArrays, bcs) +function flow_bcs!(::AMDGPUBackendTrait, stokes::JustRelax.StokesArrays, bcs::VelocityBoundaryConditions) return _flow_bcs!(bcs, @velocity(stokes)) end +function JR3D.flow_bcs!(::AMDGPUBackendTrait, stokes::JustRelax.StokesArrays, bcs::DisplacementBoundaryConditions) + return _flow_bcs!(bcs, @displacement(stokes)) +end + +function flow_bcs!(::AMDGPUBackendTrait, stokes::JustRelax.StokesArrays, bcs::DisplacementBoundaryConditions) + return _flow_bcs!(bcs, @displacement(stokes)) +end + function JR3D.thermal_bcs!(::AMDGPUBackendTrait, thermal::JustRelax.ThermalArrays, bcs) return thermal_bcs!(thermal.T, bcs) end diff --git a/src/ext/CUDA/2D.jl b/src/ext/CUDA/2D.jl index 526a0b58..5d968b65 100644 --- a/src/ext/CUDA/2D.jl +++ b/src/ext/CUDA/2D.jl @@ -21,7 +21,11 @@ import JustRelax: Geometry, @cell import JustRelax: - AbstractBoundaryConditions, TemperatureBoundaryConditions, FlowBoundaryConditions + AbstractBoundaryConditions, + TemperatureBoundaryConditions, + AbstractFlowBoundaryConditions, + DisplacementBoundaryConditions, + VelocityBoundaryConditions @init_parallel_stencil(CUDA, Float64, 2) @@ -89,11 +93,19 @@ function JR2D.PTThermalCoeffs( end # Boundary conditions -function JR2D.flow_bcs!(::CUDABackendTrait, stokes::JustRelax.StokesArrays, bcs) +function JR2D.flow_bcs!(::CUDABackendTrait, stokes::JustRelax.StokesArrays, bcs::VelocityBoundaryConditions) return _flow_bcs!(bcs, @velocity(stokes)) end -function flow_bcs!(::CUDABackendTrait, stokes::JustRelax.StokesArrays, bcs) +function flow_bcs!(::CUDABackendTrait, stokes::JustRelax.StokesArrays, bcs::VelocityBoundaryConditions) + return _flow_bcs!(bcs, @velocity(stokes)) +end + +function JR2D.flow_bcs!(::CUDABackendTrait, stokes::JustRelax.StokesArrays, bcs::DisplacementBoundaryConditions) + return _flow_bcs!(bcs, @displacement(stokes)) +end + +function flow_bcs!(::CUDABackendTrait, stokes::JustRelax.StokesArrays, bcs::DisplacementBoundaryConditions) return _flow_bcs!(bcs, @velocity(stokes)) end diff --git a/src/ext/CUDA/3D.jl b/src/ext/CUDA/3D.jl index bc9579b4..1694a53f 100644 --- a/src/ext/CUDA/3D.jl +++ b/src/ext/CUDA/3D.jl @@ -21,7 +21,12 @@ import JustRelax: Geometry, @cell import JustRelax: - AbstractBoundaryConditions, TemperatureBoundaryConditions, FlowBoundaryConditions +AbstractBoundaryConditions, +TemperatureBoundaryConditions, +AbstractFlowBoundaryConditions, +DisplacementBoundaryConditions, +VelocityBoundaryConditions + @init_parallel_stencil(CUDA, Float64, 3) @@ -100,14 +105,23 @@ function JR3D.update_pt_thermal_arrays!( end # Boundary conditions -function JR3D.flow_bcs!(::CUDABackendTrait, stokes::JustRelax.StokesArrays, bcs) +function JR3D.flow_bcs!(::CUDABackendTrait, stokes::JustRelax.StokesArrays, bcs::VelocityBoundaryConditions) return _flow_bcs!(bcs, @velocity(stokes)) end -function flow_bcs!(::CUDABackendTrait, stokes::JustRelax.StokesArrays, bcs) +function flow_bcs!(::CUDABackendTrait, stokes::JustRelax.StokesArrays, bcs::VelocityBoundaryConditions) return _flow_bcs!(bcs, @velocity(stokes)) end +# Boundary conditions +function JR3D.flow_bcs!(::CUDABackendTrait, stokes::JustRelax.StokesArrays, bcs::DisplacementBoundaryConditions) + return _flow_bcs!(bcs, @displacement(stokes)) +end + +function flow_bcs!(::CUDABackendTrait, stokes::JustRelax.StokesArrays, bcs::DisplacementBoundaryConditions) + return _flow_bcs!(bcs, @displacement(stokes)) +end + function JR3D.thermal_bcs!(::CUDABackendTrait, thermal::JustRelax.ThermalArrays, bcs) return thermal_bcs!(thermal.T, bcs) end diff --git a/src/stokes/Stokes2D.jl b/src/stokes/Stokes2D.jl index 4551fbbb..e4c3d356 100644 --- a/src/stokes/Stokes2D.jl +++ b/src/stokes/Stokes2D.jl @@ -20,7 +20,7 @@ function _solve!( stokes::JustRelax.StokesArrays, pt_stokes, di::NTuple{2,T}, - flow_bcs::FlowBoundaryConditions, + flow_bcs::AbstractFlowBoundaryConditions, ρg, K, dt, @@ -54,6 +54,9 @@ function _solve!( norm_Ry = Float64[] norm_∇V = Float64[] + # convert displacement to velocity + displacement2velocity!(stokes, dt) + # solver loop wtime0 = 0.0 while iter < 2 || (err > ϵ && iter ≤ iterMax) @@ -82,6 +85,7 @@ function _solve!( dt, ) # apply boundary conditions + velocity2displacement!(stokes, dt) flow_bcs!(stokes, flow_bcs) update_halo!(@velocity(stokes)...) end @@ -186,6 +190,9 @@ function _solve!( norm_Ry = Float64[] norm_∇V = Float64[] + # convert displacement to velocity + displacement2velocity!(stokes, dt) + # solver loop wtime0 = 0.0 while iter < 2 || (err > ϵ && iter ≤ iterMax) @@ -217,8 +224,9 @@ function _solve!( _di..., ) # free slip boundary conditions + velocity2displacement!(stokes, dt) flow_bcs!(stokes, flow_bcs) - update_halo!(stokes.V.Vx, stokes.V.Vy) + update_halo!(@velocity(stokes)...) end end @@ -324,6 +332,9 @@ function _solve!( compute_ρg!(ρg[end], phase_ratios, rheology, args) compute_viscosity!(stokes, phase_ratios, args, rheology, viscosity_cutoff) + # convert displacement to velocity + displacement2velocity!(stokes, dt) + while iter < 2 || (err > ϵ && iter ≤ iterMax) wtime0 += @elapsed begin @parallel (@idx ni) compute_∇V!(stokes.∇V, @velocity(stokes)..., _di...) @@ -380,8 +391,9 @@ function _solve!( dt * free_surface, ) # apply boundary conditions + velocity2displacement!(stokes, dt) flow_bcs!(stokes, flow_bcs) - update_halo!(stokes.V.Vx, stokes.V.Vy) + update_halo!(@velocity(stokes)...) end end @@ -424,7 +436,7 @@ function _solve!( end end - stokes.P .= θ + stokes.P .= θ # θ = P + plastic_overpressure @parallel (@idx ni .+ 1) multi_copy!(@tensor(stokes.τ_o), @tensor(stokes.τ)) @parallel (@idx ni) multi_copy!(@tensor_center(stokes.τ_o), @tensor_center(stokes.τ)) @@ -510,6 +522,7 @@ function _solve!( # compute buoyancy forces and viscosity compute_ρg!(ρg[end], phase_ratios, rheology, args) compute_viscosity!(stokes, phase_ratios, args, rheology, viscosity_cutoff) + displacement2velocity!(stokes, dt) while iter ≤ iterMax iterMin < iter && err < ϵ && break @@ -597,6 +610,7 @@ function _solve!( dt * free_surface, ) # apply boundary conditions + velocity2displacement!(stokes, dt) free_surface_bcs!(stokes, flow_bcs, η, rheology, phase_ratios, dt, di) flow_bcs!(stokes, flow_bcs) update_halo!(@velocity(stokes)...) @@ -650,7 +664,7 @@ function _solve!( end end - stokes.P .= θ + stokes.P .= θ # θ = P + plastic_overpressure @parallel (@idx ni .+ 1) multi_copy!(@tensor(stokes.τ_o), @tensor(stokes.τ)) @parallel (@idx ni) multi_copy!(@tensor_center(stokes.τ_o), @tensor_center(stokes.τ)) diff --git a/src/stokes/Stokes3D.jl b/src/stokes/Stokes3D.jl index b4ff2ee3..5abb5dbd 100644 --- a/src/stokes/Stokes3D.jl +++ b/src/stokes/Stokes3D.jl @@ -63,6 +63,9 @@ function _solve!( norm_Rz = Float64[] norm_∇V = Float64[] + # convert displacement to velocity + displacement2velocity!(stokes, dt) + # solver loop wtime0 = 0.0 while iter < 2 || (err > ϵ && iter ≤ iterMax) @@ -105,8 +108,9 @@ function _solve!( _di..., ) # apply boundary conditions + velocity2displacement!(stokes, dt) flow_bcs!(stokes, flow_bcs) - update_halo!(stokes.V.Vx, stokes.V.Vy, stokes.V.Vz) + update_halo!(@velocity(stokes)...) end end @@ -161,7 +165,7 @@ function _solve!( stokes::JustRelax.StokesArrays, pt_stokes, di::NTuple{3,T}, - flow_bcs::FlowBoundaryConditions, + flow_bcs::AbstractFlowBoundaryConditions, ρg, rheology::MaterialParams, args, @@ -207,6 +211,9 @@ function _solve!( compute_ρg!(ρg[end], phase_ratios, rheology, args) compute_viscosity!(stokes, phase_ratios, args, rheology, viscosity_cutoff) + # convert displacement to velocity + displacement2velocity!(stokes, dt) + # solver loop wtime0 = 0.0 while iter < 2 || (err > ϵ && iter ≤ iterMax) @@ -284,8 +291,9 @@ function _solve!( _di..., ) # apply boundary conditions + velocity2displacement!(stokes, dt) flow_bcs!(stokes, flow_bcs) - update_halo!(stokes.V.Vx, stokes.V.Vy, stokes.V.Vz) + update_halo!(@velocity(stokes)...) end end @@ -344,7 +352,7 @@ function _solve!( stokes::JustRelax.StokesArrays, pt_stokes, di::NTuple{3,T}, - flow_bc::FlowBoundaryConditions, + flow_bc::AbstractFlowBoundaryConditions, ρg, phase_ratios::JustRelax.PhaseRatio, rheology::NTuple{N,AbstractMaterialParamsStruct}, @@ -392,6 +400,9 @@ function _solve!( compute_ρg!(ρg[end], phase_ratios, rheology, args) compute_viscosity!(stokes, phase_ratios, args, rheology, viscosity_cutoff) + # convert displacement to velocity + displacement2velocity!(stokes, dt) + while iter < 2 || (err > ϵ && iter ≤ iterMax) wtime0 += @elapsed begin # ~preconditioner @@ -476,6 +487,7 @@ function _solve!( _di..., ) # apply boundary conditions + velocity2displacement!(stokes, dt) free_surface_bcs!(stokes, flow_bc, η, rheology, phase_ratios, dt, di) flow_bcs!(stokes, flow_bc) update_halo!(@velocity(stokes)...) From 555cae101b999e6fb4b5ebbe5a854ace6595126c Mon Sep 17 00:00:00 2001 From: Pascal Aellig Date: Wed, 3 Jul 2024 17:15:43 +0200 Subject: [PATCH 02/20] fix bug --- src/types/displacement.jl | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/src/types/displacement.jl b/src/types/displacement.jl index 1988c02f..f01a9fa7 100644 --- a/src/types/displacement.jl +++ b/src/types/displacement.jl @@ -11,20 +11,20 @@ function _velocity2displacement!(stokes::JustRelax.StokesArrays, dt) ni = size(stokes.P) (; V, U) = stokes @parallel (@idx ni .+ 2) _velocity2displacement!( - V.Vx, V.Vy, V.Vz, U.Ux, U.Uy, U.Uz, 1 / dt + V.Vx, V.Vy, V.Vz, U.Ux, U.Uy, U.Uz, dt ) return nothing end -@parallel_indices (I...) function _velocity2displacement!(Vx, Vy, Vz, Ux, Uy, Uz, _dt) +@parallel_indices (I...) function _velocity2displacement!(Vx, Vy, Vz, Ux, Uy, Uz, dt) if all(I .≤ size(Ux)) - Ux[I...] = Vx[I...] * _dt + Ux[I...] = Vx[I...] * dt end if all(I .≤ size(Uy)) - Uy[I...] = Vy[I...] * _dt + Uy[I...] = Vy[I...] * dt end if !isnothing(Vz) && all(I .≤ size(Uz)) - Uz[I...] = Vz[I...] * _dt + Uz[I...] = Vz[I...] * dt end return nothing end @@ -41,19 +41,19 @@ end function _displacement2velocity!(stokes::JustRelax.StokesArrays, dt) ni = size(stokes.P) (; V, U) = stokes - @parallel (@idx ni .+ 2) _displacement2velocity!(U.Ux, U.Uy, U.Uz, V.Vx, V.Vy, V.Vz, dt) + @parallel (@idx ni .+ 2) _displacement2velocity!(U.Ux, U.Uy, U.Uz, V.Vx, V.Vy, V.Vz, 1 / dt) return nothing end -@parallel_indices (I...) function _displacement2velocity!(Ux, Uy, Uz, Vx, Vy, Vz, dt) +@parallel_indices (I...) function _displacement2velocity!(Ux, Uy, Uz, Vx, Vy, Vz, _dt) if all(I .≤ size(Ux)) - Vx[I...] = Ux[I...] * dt + Vx[I...] = Ux[I...] * _dt end if all(I .≤ size(Uy)) - Vy[I...] = Uy[I...] * dt + Vy[I...] = Uy[I...] * _dt end if !isnothing(Vz) && all(I .≤ size(Uz)) - Vz[I...] = Uz[I...] * dt + Vz[I...] = Uz[I...] * _dt end return nothing end From 914e7d2b6dd9b905cc60e3badb221b4d398b0917 Mon Sep 17 00:00:00 2001 From: Pascal Aellig Date: Thu, 4 Jul 2024 10:19:02 +0200 Subject: [PATCH 03/20] update docs --- README.md | 18 +++++++++++++++--- docs/src/man/Blankenbach.md | 2 +- docs/src/man/ShearBands.md | 4 ++-- docs/src/man/subduction2D/subduction2D.md | 2 +- 4 files changed, 19 insertions(+), 7 deletions(-) diff --git a/README.md b/README.md index 9f83b30f..d8a76ef7 100644 --- a/README.md +++ b/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). diff --git a/docs/src/man/Blankenbach.md b/docs/src/man/Blankenbach.md index 96060e0d..10a7f763 100644 --- a/docs/src/man/Blankenbach.md +++ b/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(; diff --git a/docs/src/man/ShearBands.md b/docs/src/man/ShearBands.md index 1f9af1a8..c90591c6 100644 --- a/docs/src/man/ShearBands.md +++ b/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)...) ``` diff --git a/docs/src/man/subduction2D/subduction2D.md b/docs/src/man/subduction2D/subduction2D.md index 193a87e6..913d2cf9 100644 --- a/docs/src/man/subduction2D/subduction2D.md +++ b/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), ) ``` From d10aad27b632f635016d808f32fa52e29e726afa Mon Sep 17 00:00:00 2001 From: Pascal Aellig Date: Thu, 4 Jul 2024 10:19:46 +0200 Subject: [PATCH 04/20] update miniapps --- .../Blankenbach2D/Benchmark2D_WENO5.jl | 4 +- .../stokes2D/Blankenbach2D/Benchmark2D_sgd.jl | 4 +- .../Blankenbach2D/Benchmark2D_sgd_scaled.jl | 4 +- .../stokes2D/VanKeken.jl/VanKeken.jl | 4 +- .../elastic_buildup/Elastic_BuildUp.jl | 4 +- .../Elastic_BuildUp_phases_incompressible.jl | 4 +- .../PlumeFreeSurface_2D.jl | 2 +- .../RayleighTaylor2D.jl | 2 +- .../stokes2D/shear_band/ShearBand2D.jl | 8 +- .../stokes2D/shear_band/ShearBand2D_MPI.jl | 4 +- .../shear_band/ShearBand2D_displacement.jl | 196 ++++++++++++++++++ .../shear_band/ShearBand2D_softening.jl | 6 +- .../stokes2D/shear_heating/Shearheating2D.jl | 4 +- .../stokes2D/sinking_block/SinkingBlock2D.jl | 4 +- miniapps/benchmarks/stokes2D/solcx/SolCx.jl | 6 +- miniapps/benchmarks/stokes2D/solkz/SolKz.jl | 4 +- miniapps/benchmarks/stokes2D/solvi/SolVi.jl | 4 +- miniapps/benchmarks/stokes2D/solvi/SolViEl.jl | 4 +- .../stokes3D/burstedde/Burstedde.jl | 4 +- .../stokes3D/shear_band/ShearBand3D.jl | 4 +- .../stokes3D/shear_band/ShearBand3D_MPI.jl | 2 +- .../stokes3D/shear_heating/Shearheating3D.jl | 4 +- miniapps/benchmarks/stokes3D/solvi/SolVi3D.jl | 30 +-- .../stokes3D/taylor_green/TaylorGreen.jl | 4 +- .../Thermal_Stress_Magma_Chamber_nondim.jl | 4 +- .../Thermal_Stress_Magma_Chamber_nondim3D.jl | 2 +- .../convection/GlobalConvection2D_Upwind.jl | 4 +- .../convection/GlobalConvection2D_WENO5.jl | 2 +- .../convection/GlobalConvection3D_Upwind.jl | 2 +- .../Particles2D/Layered_convection2D.jl | 2 +- .../Layered_convection2D.jl | 4 +- .../Particles3D/Layered_convection3D.jl | 4 +- miniapps/convection/RisingBlob3D/Blob3D.jl | 10 +- .../convection/WENO5/WENO_convection2D.jl | 2 +- miniapps/subduction/2D/Subd2D.jl | 2 +- 35 files changed, 275 insertions(+), 79 deletions(-) create mode 100644 miniapps/benchmarks/stokes2D/shear_band/ShearBand2D_displacement.jl diff --git a/miniapps/benchmarks/stokes2D/Blankenbach2D/Benchmark2D_WENO5.jl b/miniapps/benchmarks/stokes2D/Blankenbach2D/Benchmark2D_WENO5.jl index 9a35b8a3..e563bb6a 100755 --- a/miniapps/benchmarks/stokes2D/Blankenbach2D/Benchmark2D_WENO5.jl +++ b/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 diff --git a/miniapps/benchmarks/stokes2D/Blankenbach2D/Benchmark2D_sgd.jl b/miniapps/benchmarks/stokes2D/Blankenbach2D/Benchmark2D_sgd.jl index 92930bf5..1e383b54 100755 --- a/miniapps/benchmarks/stokes2D/Blankenbach2D/Benchmark2D_sgd.jl +++ b/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 diff --git a/miniapps/benchmarks/stokes2D/Blankenbach2D/Benchmark2D_sgd_scaled.jl b/miniapps/benchmarks/stokes2D/Blankenbach2D/Benchmark2D_sgd_scaled.jl index 93cc85ab..00ff9245 100755 --- a/miniapps/benchmarks/stokes2D/Blankenbach2D/Benchmark2D_sgd_scaled.jl +++ b/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 diff --git a/miniapps/benchmarks/stokes2D/VanKeken.jl/VanKeken.jl b/miniapps/benchmarks/stokes2D/VanKeken.jl/VanKeken.jl index 8236e35d..da29592d 100644 --- a/miniapps/benchmarks/stokes2D/VanKeken.jl/VanKeken.jl +++ b/miniapps/benchmarks/stokes2D/VanKeken.jl/VanKeken.jl @@ -102,12 +102,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 diff --git a/miniapps/benchmarks/stokes2D/elastic_buildup/Elastic_BuildUp.jl b/miniapps/benchmarks/stokes2D/elastic_buildup/Elastic_BuildUp.jl index 9ea1014b..887b393a 100644 --- a/miniapps/benchmarks/stokes2D/elastic_buildup/Elastic_BuildUp.jl +++ b/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 diff --git a/miniapps/benchmarks/stokes2D/elastic_buildup/Elastic_BuildUp_phases_incompressible.jl b/miniapps/benchmarks/stokes2D/elastic_buildup/Elastic_BuildUp_phases_incompressible.jl index 78009d4a..f8c694eb 100644 --- a/miniapps/benchmarks/stokes2D/elastic_buildup/Elastic_BuildUp_phases_incompressible.jl +++ b/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 diff --git a/miniapps/benchmarks/stokes2D/free_surface_stabilization/PlumeFreeSurface_2D.jl b/miniapps/benchmarks/stokes2D/free_surface_stabilization/PlumeFreeSurface_2D.jl index 393b7f59..65fe5254 100644 --- a/miniapps/benchmarks/stokes2D/free_surface_stabilization/PlumeFreeSurface_2D.jl +++ b/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 ) diff --git a/miniapps/benchmarks/stokes2D/free_surface_stabilization/RayleighTaylor2D.jl b/miniapps/benchmarks/stokes2D/free_surface_stabilization/RayleighTaylor2D.jl index a96b43da..33aa5901 100644 --- a/miniapps/benchmarks/stokes2D/free_surface_stabilization/RayleighTaylor2D.jl +++ b/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, diff --git a/miniapps/benchmarks/stokes2D/shear_band/ShearBand2D.jl b/miniapps/benchmarks/stokes2D/shear_band/ShearBand2D.jl index 1c0be034..32affbd6 100644 --- a/miniapps/benchmarks/stokes2D/shear_band/ShearBand2D.jl +++ b/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,7 +82,7 @@ function main(igg; nx=64, ny=64, figdir="model_figs") Elasticity = el_inc, ), ) - + # perturbation array for the cohesion perturbation_C = @rand(ni...) @@ -90,7 +90,7 @@ function main(igg; nx=64, ny=64, figdir="model_figs") radius = 0.1 phase_ratios = PhaseRatio(ni, length(rheology)) init_phases!(phase_ratios, xci, radius) - + # STOKES --------------------------------------------- # Allocate arrays needed for every Stokes problem stokes = StokesArrays(backend, ni) @@ -105,7 +105,7 @@ 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), ) diff --git a/miniapps/benchmarks/stokes2D/shear_band/ShearBand2D_MPI.jl b/miniapps/benchmarks/stokes2D/shear_band/ShearBand2D_MPI.jl index 8fa018cf..0e15a353 100644 --- a/miniapps/benchmarks/stokes2D/shear_band/ShearBand2D_MPI.jl +++ b/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 diff --git a/miniapps/benchmarks/stokes2D/shear_band/ShearBand2D_displacement.jl b/miniapps/benchmarks/stokes2D/shear_band/ShearBand2D_displacement.jl new file mode 100644 index 00000000..101283d8 --- /dev/null +++ b/miniapps/benchmarks/stokes2D/shear_band/ShearBand2D_displacement.jl @@ -0,0 +1,196 @@ +using GeoParams, CairoMakie, CellArrays +using JustRelax, JustRelax.JustRelax2D +using ParallelStencil +@init_parallel_stencil(Threads, Float64, 2) + +const backend = CPUBackend + +# HELPER FUNCTIONS ----------------------------------- ---------------------------- +solution(ε, t, G, η) = 2 * ε * η * (1 - exp(-G * t / η)) + +# 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 (@idx ni) init_phases!(phase_ratios.center, xci..., origin..., radius) +end + +# MAIN SCRIPT -------------------------------------------------------------------- +function main(igg; nx=64, ny=64, figdir="model_figs") + + # Physical domain ------------------------------------ + ly = 1e0 # 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 + + # 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 + ) + + rheology = ( + # Low density phase + SetMaterialParams(; + Phase = 1, + Density = ConstantDensity(; ρ = 0.0), + Gravity = ConstantGravity(; g = 0.0), + CompositeRheology = CompositeRheology((visc, el_bg, pl)), + Elasticity = el_bg, + + ), + # High density phase + SetMaterialParams(; + Phase = 1, + Density = ConstantDensity(; ρ = 0.0), + Gravity = ConstantGravity(; g = 0.0), + CompositeRheology = CompositeRheology((visc, el_inc, pl)), + Elasticity = el_inc, + ), + ) + + # perturbation array for the cohesion + perturbation_C = @rand(ni...) + + # Initialize phase ratios ------------------------------- + radius = 0.1 + phase_ratios = PhaseRatio(ni, length(rheology)) + init_phases!(phase_ratios, xci, radius) + + # STOKES --------------------------------------------- + # Allocate arrays needed for every Stokes problem + stokes = StokesArrays(backend, ni) + pt_stokes = PTStokesCoeffs(li, di; ϵ=1e-6, CFL = 0.75 / √2.1) + + # Buoyancy forces + ρg = @zeros(ni...), @zeros(ni...) + args = (; T = @zeros(ni...), P = stokes.P, dt = dt, perturbation_C = perturbation_C) + + # Rheology + compute_viscosity!( + stokes, phase_ratios, args, rheology, (-Inf, Inf) + ) + # 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) + update_halo!(@velocity(stokes)...) + + # IO ------------------------------------------------- + take(figdir) + + # Time loop + t, it = 0.0, 0 + tmax = 3.5 + τII = Float64[] + sol = Float64[] + ttot = Float64[] + + while t < tmax + + # Stokes solver ---------------- + iters = solve!( + stokes, + pt_stokes, + di, + flow_bcs, + ρg, + phase_ratios, + rheology, + args, + dt, + igg; + kwargs = ( + verbose = false, + iterMax = 50e3, + nout = 1e2, + viscosity_cutoff = (-Inf, Inf) + ) + ) + tensor_invariant!(stokes.ε) + push!(τII, maximum(stokes.τ.xx)) + + it += 1 + t += dt + + push!(sol, solution(εbg, t, G0, η0)) + push!(ttot, t) + + println("it = $it; t = $t \n") + + # visualisation + 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 = "η_vep") + 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.viscosity.η_vep)) , 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) + + end + + return nothing +end + +n = 128 +nx = n +ny = n +figdir = "ShearBands2D_displacement" +igg = if !(JustRelax.MPI.Initialized()) + IGG(init_global_grid(nx, ny, 1; init_MPI = true)...) +else + igg +end +main(igg; figdir = figdir, nx = nx, ny = ny); diff --git a/miniapps/benchmarks/stokes2D/shear_band/ShearBand2D_softening.jl b/miniapps/benchmarks/stokes2D/shear_band/ShearBand2D_softening.jl index 786f87e6..94a8b264 100644 --- a/miniapps/benchmarks/stokes2D/shear_band/ShearBand2D_softening.jl +++ b/miniapps/benchmarks/stokes2D/shear_band/ShearBand2D_softening.jl @@ -25,7 +25,7 @@ function init_phases!(phase_ratios, xci, radius) return nothing end - @parallel (JustRelax.@idx ni) init_phases!(phase_ratios.center, xci..., origin..., radius) + @parallel (@idx ni) init_phases!(phase_ratios.center, xci..., origin..., radius) end # MAIN SCRIPT -------------------------------------------------------------------- @@ -102,7 +102,7 @@ function main(igg; nx=64, ny=64, figdir="model_figs") stokes, phase_ratios, args, rheology, (-Inf, Inf) ) # Boundary conditions - flow_bcs = FlowBoundaryConditions(; + flow_bcs = DisplacementBoundaryConditions(; free_slip = (left = true, right = true, top = true, bot = true), no_slip = (left = false, right = false, top = false, bot=false), ) @@ -183,7 +183,7 @@ end n = 128 nx = n ny = n -figdir = "ShearBands2D_lin_softening" +figdir = "Shearband_Displacement" igg = if !(JustRelax.MPI.Initialized()) IGG(init_global_grid(nx, ny, 1; init_MPI = true)...) else diff --git a/miniapps/benchmarks/stokes2D/shear_heating/Shearheating2D.jl b/miniapps/benchmarks/stokes2D/shear_heating/Shearheating2D.jl index 8c2ac6bd..544b2071 100644 --- a/miniapps/benchmarks/stokes2D/shear_heating/Shearheating2D.jl +++ b/miniapps/benchmarks/stokes2D/shear_heating/Shearheating2D.jl @@ -117,7 +117,7 @@ function main2D(igg; ar=8, ny=16, nx=ny*8, figdir="figs2D", do_vtk =false) ) # Boundary conditions - flow_bcs = FlowBoundaryConditions(; + flow_bcs = VelocityBoundaryConditions(; free_slip = (left = true, right=true, top=true, bot=true), ) ## Compression and not extension - fix this @@ -125,7 +125,7 @@ function main2D(igg; ar=8, ny=16, nx=ny*8, figdir="figs2D", do_vtk =false) stokes.V.Vx .= PTArray(backend_JR)([ -(x - lx/2) * εbg for x in xvi[1], _ in 1:ny+2]) stokes.V.Vy .= PTArray(backend_JR)([ (ly - abs(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)...) T_buffer = @zeros(ni.+1) Told_buffer = similar(T_buffer) diff --git a/miniapps/benchmarks/stokes2D/sinking_block/SinkingBlock2D.jl b/miniapps/benchmarks/stokes2D/sinking_block/SinkingBlock2D.jl index 416bf645..3b3d52a6 100644 --- a/miniapps/benchmarks/stokes2D/sinking_block/SinkingBlock2D.jl +++ b/miniapps/benchmarks/stokes2D/sinking_block/SinkingBlock2D.jl @@ -134,11 +134,11 @@ function sinking_block2D(igg; ar=8, ny=16, nx=ny*8, figdir="figs2D", thermal_per # ---------------------------------------------------- # 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)...) # Stokes solver ---------------- args = (; T = @ones(ni...), P = stokes.P, dt=dt, ΔTc = @zeros(ni...)) diff --git a/miniapps/benchmarks/stokes2D/solcx/SolCx.jl b/miniapps/benchmarks/stokes2D/solcx/SolCx.jl index 8953ab60..0085721a 100644 --- a/miniapps/benchmarks/stokes2D/solcx/SolCx.jl +++ b/miniapps/benchmarks/stokes2D/solcx/SolCx.jl @@ -105,15 +105,15 @@ function solCx( @views η2[end, :] .= η2[end-1, :] @views η2[:, 1] .= η2[:, 2] @views η2[:, end] .= η2[:, end-1] - η, η2 = η2, η # swap + η, η2 = η2, η # swap end ## 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)...) # Physical time loop t = 0.0 diff --git a/miniapps/benchmarks/stokes2D/solkz/SolKz.jl b/miniapps/benchmarks/stokes2D/solkz/SolKz.jl index acd44666..da66bf74 100644 --- a/miniapps/benchmarks/stokes2D/solkz/SolKz.jl +++ b/miniapps/benchmarks/stokes2D/solkz/SolKz.jl @@ -83,11 +83,11 @@ function solKz(; Kb = @fill(Inf, ni...) ## 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)...) # Physical time loop t = 0.0 diff --git a/miniapps/benchmarks/stokes2D/solvi/SolVi.jl b/miniapps/benchmarks/stokes2D/solvi/SolVi.jl index 7bfb0267..05d06a77 100644 --- a/miniapps/benchmarks/stokes2D/solvi/SolVi.jl +++ b/miniapps/benchmarks/stokes2D/solvi/SolVi.jl @@ -75,11 +75,11 @@ function solVi(; ## Boundary conditions pureshear_bc!(stokes, xci, xvi, εbg) - 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)...) # Physical time loop t = 0.0 diff --git a/miniapps/benchmarks/stokes2D/solvi/SolViEl.jl b/miniapps/benchmarks/stokes2D/solvi/SolViEl.jl index 66cae48f..acf2079a 100644 --- a/miniapps/benchmarks/stokes2D/solvi/SolViEl.jl +++ b/miniapps/benchmarks/stokes2D/solvi/SolViEl.jl @@ -72,11 +72,11 @@ function solViEl(; ## Boundary conditions pureshear_bc!(stokes, xci, xvi, εbg) - 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)...) # Physical time loop local iters diff --git a/miniapps/benchmarks/stokes3D/burstedde/Burstedde.jl b/miniapps/benchmarks/stokes3D/burstedde/Burstedde.jl index 6b654e3b..8f18c9e9 100644 --- a/miniapps/benchmarks/stokes3D/burstedde/Burstedde.jl +++ b/miniapps/benchmarks/stokes3D/burstedde/Burstedde.jl @@ -208,14 +208,14 @@ function burstedde(; nx=16, ny=16, nz=16, init_MPI=true, finalize_MPI=false) K = @fill(Inf, ni...) ## Boundary conditions - flow_bcs = FlowBoundaryConditions(; + flow_bcs = VelocityBoundaryConditionsionsions(; free_slip = (left=false, right=false, top=false, bot=false, back=false, front=false), no_slip = (left=false, right=false, top=false, bot=false, back=false, front=false), ) # impose analytical velociity at the boundaries of the domain velocity!(stokes, xci, xvi, di) flow_bcs!(stokes, flow_bcs) # apply boundary conditions - update_halo!(stokes.V.Vx, stokes.V.Vy, stokes.V.Vz) + update_halo!(@velocity(stokes)...) # Physical time loop t = 0.0 diff --git a/miniapps/benchmarks/stokes3D/shear_band/ShearBand3D.jl b/miniapps/benchmarks/stokes3D/shear_band/ShearBand3D.jl index 173cf517..48256186 100644 --- a/miniapps/benchmarks/stokes3D/shear_band/ShearBand3D.jl +++ b/miniapps/benchmarks/stokes3D/shear_band/ShearBand3D.jl @@ -98,7 +98,7 @@ function main(igg; nx=64, ny=64, nz=64, figdir="model_figs") compute_viscosity!(stokes, phase_ratios, args, rheology, cutoff_visc) # Boundary conditions - flow_bcs = FlowBoundaryConditions(; + flow_bcs = VelocityBoundaryConditions(; free_slip = (left = true , right = true , top = true , bot = true , back = true , front = true), no_slip = (left = false, right = false, top = false, bot = false, back = false, front = false), ) @@ -107,7 +107,7 @@ function main(igg; nx=64, ny=64, nz=64, figdir="model_figs") stokes.V.Vz .= PTArray([-z*εbg for _ in 1:nx+2, _ in 1:nx+2, z in xvi[3]]) flow_bcs!(stokes, flow_bcs) # apply boundary conditions update_halo!(@velocity(stokes)...) - + # IO ------------------------------------------------ # if it does not exist, make folder where figures are stored !isdir(figdir) && mkpath(figdir) diff --git a/miniapps/benchmarks/stokes3D/shear_band/ShearBand3D_MPI.jl b/miniapps/benchmarks/stokes3D/shear_band/ShearBand3D_MPI.jl index e61a9965..45bae0a2 100644 --- a/miniapps/benchmarks/stokes3D/shear_band/ShearBand3D_MPI.jl +++ b/miniapps/benchmarks/stokes3D/shear_band/ShearBand3D_MPI.jl @@ -99,7 +99,7 @@ function main(igg; nx=64, ny=64, nz=64, figdir="model_figs") compute_viscosity!(stokes, phase_ratios, args, rheology, cutoff_visc) # Boundary conditions - flow_bcs = FlowBoundaryConditions(; + flow_bcs = VelocityBoundaryConditions(; free_slip = (left = true , right = true , top = true , bot = true , back = true , front = true), no_slip = (left = false, right = false, top = false, bot = false, back = false, front = false), ) diff --git a/miniapps/benchmarks/stokes3D/shear_heating/Shearheating3D.jl b/miniapps/benchmarks/stokes3D/shear_heating/Shearheating3D.jl index 26de9fcc..b3c7ef3c 100644 --- a/miniapps/benchmarks/stokes3D/shear_heating/Shearheating3D.jl +++ b/miniapps/benchmarks/stokes3D/shear_heating/Shearheating3D.jl @@ -111,7 +111,7 @@ function main3D(igg; ar=8, ny=16, nx=ny*8, nz=ny*8, figdir="figs3D", do_vtk =fal ) # Boundary conditions - flow_bcs = FlowBoundaryConditions(; + flow_bcs = VelocityBoundaryConditions(; free_slip = (left = true , right = true , top = true , bot = true , front = true , back = true ), no_slip = (left = false, right = false, top = false, bot = false, front = false, back = false), ) @@ -121,7 +121,7 @@ function main3D(igg; ar=8, ny=16, nx=ny*8, nz=ny*8, figdir="figs3D", do_vtk =fal stokes.V.Vy .= PTArray([ -(y - ly/2) * εbg for _ in 1:nx+2, y in xvi[2], _ in 1:nz+2]) stokes.V.Vz .= PTArray([ (lz - abs(z)) * εbg for _ in 1:nx+2, _ in 1:ny+2, z in xvi[3]]) flow_bcs!(stokes, flow_bcs) # apply boundary conditions - update_halo!(stokes.V.Vx, stokes.V.Vy, stokes.V.Vz) + update_halo!(@velocity(stokes)...) # IO ----- ------------------------------------------- # if it does not exist, make folder where figures are stored diff --git a/miniapps/benchmarks/stokes3D/solvi/SolVi3D.jl b/miniapps/benchmarks/stokes3D/solvi/SolVi3D.jl index 9bdad110..c636af80 100644 --- a/miniapps/benchmarks/stokes3D/solvi/SolVi3D.jl +++ b/miniapps/benchmarks/stokes3D/solvi/SolVi3D.jl @@ -91,12 +91,12 @@ function solVi3D(; ## 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, back=true, front=true), no_slip = (left=false, right=false, top=false, bot=false, back=false, front=false), ) flow_bcs!(stokes, flow_bcs) # apply boundary conditions - update_halo!(stokes.V.Vx, stokes.V.Vy, stokes.V.Vz) + update_halo!(@velocity(stokes)...) ## Body forces ρg = ntuple(_ -> @zeros(ni...), Val(3)) @@ -106,18 +106,18 @@ function solVi3D(; local iters while t < ttot iters = solve!( - stokes, - pt_stokes, - di, - flow_bcs, - ρg, - Kb, - Gc, - dt, - igg; - kwargs = (; - iterMax=5000, - nout=100, + stokes, + pt_stokes, + di, + flow_bcs, + ρg, + Kb, + Gc, + dt, + igg; + kwargs = (; + iterMax=5000, + nout=100, verbose=false ), ) @@ -127,4 +127,4 @@ function solVi3D(; finalize_global_grid(; finalize_MPI=finalize_MPI) return (ni=ni, xci=xci, xvi=xvi, li=li, di=di), stokes, iters -end \ No newline at end of file +end diff --git a/miniapps/benchmarks/stokes3D/taylor_green/TaylorGreen.jl b/miniapps/benchmarks/stokes3D/taylor_green/TaylorGreen.jl index ea943937..a22dd577 100644 --- a/miniapps/benchmarks/stokes3D/taylor_green/TaylorGreen.jl +++ b/miniapps/benchmarks/stokes3D/taylor_green/TaylorGreen.jl @@ -113,14 +113,14 @@ function taylorGreen(; nx=16, ny=16, nz=16, init_MPI=true, finalize_MPI=false) K = @fill(Inf, ni...) ## Boundary conditions - flow_bcs = FlowBoundaryConditions(; + flow_bcs = VelocityBoundaryConditions(; free_slip=(left=false, right=false, top=false, bot=false, back=false, front=false), no_slip=(left=false, right=false, top=false, bot=false, back=false, front=false), ) # impose analytical velocity at the boundaries of the domain velocity!(stokes, xci, xvi) flow_bcs!(stokes, flow_bcs) # apply boundary conditions - update_halo!(stokes.V.Vx, stokes.V.Vy, stokes.V.Vz) + update_halo!(@velocity(stokes)...) # Physical time loop t = 0.0 diff --git a/miniapps/benchmarks/thermal_stress/Thermal_Stress_Magma_Chamber_nondim.jl b/miniapps/benchmarks/thermal_stress/Thermal_Stress_Magma_Chamber_nondim.jl index b2dccf60..a78a52db 100644 --- a/miniapps/benchmarks/thermal_stress/Thermal_Stress_Magma_Chamber_nondim.jl +++ b/miniapps/benchmarks/thermal_stress/Thermal_Stress_Magma_Chamber_nondim.jl @@ -281,12 +281,12 @@ function main2D(igg; figdir=figdir, nx=nx, ny=ny, do_vtk=false) (abs(y) - sticky_air) * εbg * (abs(y) > sticky_air) for _ in 1:(nx + 2), y in xvi[2] ]) - flow_bcs = FlowBoundaryConditions(; + flow_bcs = VelocityBoundaryConditions(; free_slip = (left=true, right=true, top=true, bot=true), free_surface = true, ) flow_bcs!(stokes, flow_bcs) - update_halo!(stokes.V.Vx, stokes.V.Vy) + update_halo!(@velocity(stokes)...) η = @ones(ni...) # initialise viscosity diff --git a/miniapps/benchmarks/thermal_stress/Thermal_Stress_Magma_Chamber_nondim3D.jl b/miniapps/benchmarks/thermal_stress/Thermal_Stress_Magma_Chamber_nondim3D.jl index 2df62091..dafd7705 100644 --- a/miniapps/benchmarks/thermal_stress/Thermal_Stress_Magma_Chamber_nondim3D.jl +++ b/miniapps/benchmarks/thermal_stress/Thermal_Stress_Magma_Chamber_nondim3D.jl @@ -275,7 +275,7 @@ function main3D(igg; figdir = "output", nx = 64, ny = 64, nz = 64, do_vtk = fals backend_JR, rheology, phase_ratios, args, dt, ni, di, li; ϵ=1e-5, CFL=0.8 / √3.1 ) - flow_bcs = FlowBoundaryConditions(; + flow_bcs = VelocityBoundaryConditions(; free_slip = (left=true, right=true, front=true, back=true, top=true, bot=true), no_slip = (left=false, right=false, front=false, back=false, top=false, bot=false), free_surface = true, diff --git a/miniapps/convection/GlobalConvection2D_Upwind.jl b/miniapps/convection/GlobalConvection2D_Upwind.jl index af97a365..851f4393 100644 --- a/miniapps/convection/GlobalConvection2D_Upwind.jl +++ b/miniapps/convection/GlobalConvection2D_Upwind.jl @@ -181,11 +181,11 @@ function thermal_convection2D(; ar=8, ny=16, nx=ny*8, figdir="figs2D", thermal_p compute_viscosity!(stokes, args, rheology, viscosity_cutoff) # 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 ------------------------------------------------- diff --git a/miniapps/convection/GlobalConvection2D_WENO5.jl b/miniapps/convection/GlobalConvection2D_WENO5.jl index d7d5d3e2..a7dbb3c6 100644 --- a/miniapps/convection/GlobalConvection2D_WENO5.jl +++ b/miniapps/convection/GlobalConvection2D_WENO5.jl @@ -187,7 +187,7 @@ function thermal_convection2D(igg; ar=8, ny=16, nx=ny*8, figdir="figs2D", therma ) # 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 diff --git a/miniapps/convection/GlobalConvection3D_Upwind.jl b/miniapps/convection/GlobalConvection3D_Upwind.jl index 37e83044..41dd7544 100644 --- a/miniapps/convection/GlobalConvection3D_Upwind.jl +++ b/miniapps/convection/GlobalConvection3D_Upwind.jl @@ -184,7 +184,7 @@ function thermal_convection3D(; ar=8, nz=16, nx=ny*8, ny=nx, figdir="figs3D", th compute_viscosity!(stokes, args, rheology, (1e18, 1e24)) # Boundary conditions - flow_bcs = FlowBoundaryConditions(; + flow_bcs = VelocityBoundaryConditions(; free_slip = (left=true , right=true , top=true , bot=true , front=true , back=true ), no_slip = (left=false, right=false, top=false, bot=false, front=false, back=false), ) diff --git a/miniapps/convection/Particles2D/Layered_convection2D.jl b/miniapps/convection/Particles2D/Layered_convection2D.jl index db7ae633..178ee5f8 100644 --- a/miniapps/convection/Particles2D/Layered_convection2D.jl +++ b/miniapps/convection/Particles2D/Layered_convection2D.jl @@ -169,7 +169,7 @@ function main2D(igg; ar=8, ny=16, nx=ny*8, figdir="figs2D", do_vtk =false) ) # 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 diff --git a/miniapps/convection/Particles2D_nonDim/Layered_convection2D.jl b/miniapps/convection/Particles2D_nonDim/Layered_convection2D.jl index 15222b3b..f79a4b4a 100644 --- a/miniapps/convection/Particles2D_nonDim/Layered_convection2D.jl +++ b/miniapps/convection/Particles2D_nonDim/Layered_convection2D.jl @@ -175,11 +175,11 @@ function main2D(igg; ar=8, ny=16, nx=ny*8, figdir="figs2D", do_vtk =false) ) # 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 diff --git a/miniapps/convection/Particles3D/Layered_convection3D.jl b/miniapps/convection/Particles3D/Layered_convection3D.jl index 3ffc9702..0ffd8ce1 100644 --- a/miniapps/convection/Particles3D/Layered_convection3D.jl +++ b/miniapps/convection/Particles3D/Layered_convection3D.jl @@ -151,12 +151,12 @@ function main3D(igg; ar=1, nx=16, ny=16, nz=16, figdir="figs3D", do_vtk =false) ) # Boundary conditions - flow_bcs = FlowBoundaryConditions(; + flow_bcs = VelocityBoundaryConditions(; free_slip = (left = true , right = true , top = true , bot = true , front = true , back = true ), no_slip = (left = false, right = false, top = false, bot = false, front = false, back = false), ) flow_bcs!(stokes, flow_bcs) # apply boundary conditions - update_halo!(stokes.V.Vx, stokes.V.Vy, stokes.V.Vz) + update_halo!(@velocity(stokes)...) # IO ------------------------------------------------- # if it does not exist, make folder where figures are stored diff --git a/miniapps/convection/RisingBlob3D/Blob3D.jl b/miniapps/convection/RisingBlob3D/Blob3D.jl index e2d6301b..5c186be6 100644 --- a/miniapps/convection/RisingBlob3D/Blob3D.jl +++ b/miniapps/convection/RisingBlob3D/Blob3D.jl @@ -1,14 +1,14 @@ # const isCUDA = false const isCUDA = true -@static if isCUDA +@static if isCUDA using CUDA end using JustRelax, JustRelax.JustRelax3D, JustRelax.DataIO import JustRelax.@cell -const backend_JR = @static if isCUDA +const backend_JR = @static if isCUDA CUDABackend # Options: CPUBackend, CUDABackend, AMDGPUBackend else JustRelax.CPUBackend # Options: CPUBackend, CUDABackend, AMDGPUBackend @@ -16,7 +16,7 @@ end using ParallelStencil, ParallelStencil.FiniteDifferences3D -@static if isCUDA +@static if isCUDA @init_parallel_stencil(CUDA, Float64, 3) else @init_parallel_stencil(Threads, Float64, 3) @@ -26,7 +26,7 @@ using JustPIC, JustPIC._3D # Threads is the default backend, # to run on a CUDA GPU load CUDA.jl (i.e. "using CUDA") at the beginning of the script, # and to run on an AMD GPU load AMDGPU.jl (i.e. "using AMDGPU") at the beginning of the script. -const backend = @static if isCUDA +const backend = @static if isCUDA CUDABackend # Options: CPUBackend, CUDABackend, AMDGPUBackend else JustPIC.CPUBackend # Options: CPUBackend, CUDABackend, AMDGPUBackend @@ -286,7 +286,7 @@ function main3D(igg; figdir = "output", nx = 64, ny = 64, nz = 64, do_vtk = fals backend_JR, rheology, phase_ratios, args, dt, ni, di, li; ϵ=1e-5, CFL=0.8 / √3.1 ) - flow_bcs = FlowBoundaryConditions(; + flow_bcs = VelocityBoundaryConditions(; free_slip = (left=true, right=true, front=true, back=true, top=true, bot=true), no_slip = (left=false, right=false, front=false, back=false, top=false, bot=false), free_surface = true, diff --git a/miniapps/convection/WENO5/WENO_convection2D.jl b/miniapps/convection/WENO5/WENO_convection2D.jl index 74d5dfe9..814812fc 100644 --- a/miniapps/convection/WENO5/WENO_convection2D.jl +++ b/miniapps/convection/WENO5/WENO_convection2D.jl @@ -161,7 +161,7 @@ function main2D(igg; ar=8, ny=16, nx=ny*8, figdir="figs2D", do_vtk =false) ) # 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 diff --git a/miniapps/subduction/2D/Subd2D.jl b/miniapps/subduction/2D/Subd2D.jl index ed86966f..0d924324 100644 --- a/miniapps/subduction/2D/Subd2D.jl +++ b/miniapps/subduction/2D/Subd2D.jl @@ -135,7 +135,7 @@ function main(li, origin, phases_GMG, igg; nx=16, ny=16, figdir="figs2D", do_vtk ) # Boundary conditions - flow_bcs = FlowBoundaryConditions(; + flow_bcs = VelocityBoundaryConditions(; free_slip = (left = true , right = true , top = true , bot = true), free_surface = false, ) From bc645a4179a1f282888c7b733855aab81e135de0 Mon Sep 17 00:00:00 2001 From: Pascal Aellig Date: Thu, 4 Jul 2024 10:19:57 +0200 Subject: [PATCH 05/20] update Tests --- test/test_VanKeken.jl | 4 +- test/test_boundary_conditions2D.jl | 95 +++++++++++++++++++++++++++--- test/test_shearband2D.jl | 4 +- test/test_shearband2D_softening.jl | 4 +- test/test_shearheating2D.jl | 4 +- test/test_shearheating3D.jl | 4 +- test/test_sinking_block.jl | 4 +- test/test_types.jl | 12 ++-- 8 files changed, 105 insertions(+), 26 deletions(-) diff --git a/test/test_VanKeken.jl b/test/test_VanKeken.jl index 09dea6b4..cc762699 100644 --- a/test/test_VanKeken.jl +++ b/test/test_VanKeken.jl @@ -130,12 +130,12 @@ function VanKeken2D(ny=32, nx=32) 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)...) # Buffer arrays to compute velocity rms Vx_v = @zeros(ni.+1...) diff --git a/test/test_boundary_conditions2D.jl b/test/test_boundary_conditions2D.jl index c2d5fb2c..6fc169e7 100644 --- a/test/test_boundary_conditions2D.jl +++ b/test/test_boundary_conditions2D.jl @@ -15,16 +15,17 @@ else CPUBackend end -@testset "Boundary Conditions" begin +@testset "Boundary Conditions 2D" begin if backend === CPUBackend @suppress begin + @testset "VelocityBoundaryConditions" begin # test incompatible boundary conditions - @test_throws ErrorException FlowBoundaryConditions(; + @test_throws ErrorException VelocityBoundaryConditions(; no_slip = (left=false, right=false, top=false, bot=false), free_slip = (left=false, right=true, top=true, bot=true), ) - @test_throws ErrorException FlowBoundaryConditions(; + @test_throws ErrorException VelocityBoundaryConditions(; no_slip = (left=false, right=false, top=false, bot=false), free_slip = (left=true , right=true , top=true , bot=false), ) @@ -32,7 +33,7 @@ end n = 5 # number of elements Vx, Vy = PTArray(backend)(rand(n + 1, n + 2)), PTArray(backend)(rand(n + 2, n + 1)) # free-slip - bcs = FlowBoundaryConditions(; + bcs = VelocityBoundaryConditions(; no_slip = (left=false, right=false, top=false, bot=false), free_slip = (left=true, right=true, top=true, bot=true), ) @@ -42,10 +43,11 @@ end @test @views Vx[: , end] == Vx[:, end - 1] @test @views Vy[1 , :] == Vy[2, :] @test @views Vy[end, :] == Vy[end - 1, :] - + @test typeof(bcs) <: AbstractFlowBoundaryConditions + @test typeof(bcs) <: VelocityBoundaryConditions # no-slip Vx, Vy = PTArray(backend)(rand(n + 1, n + 2)), PTArray(backend)(rand(n + 2, n + 1)) - bcs = FlowBoundaryConditions(; + bcs = VelocityBoundaryConditions(; no_slip = (left=true, right=true, top=true, bot=true), free_slip = (left=false, right=false, top=false, bot=false), ) @@ -65,7 +67,7 @@ end stokes.V.Vx .= PTArray(backend)(rand(n + 1, n + 2)) stokes.V.Vy .= PTArray(backend)(rand(n + 2, n + 1)) # free-slip - flow_bcs = FlowBoundaryConditions(; + flow_bcs = VelocityBoundaryConditions(; no_slip = (left=false, right=false, top=false, bot=false), free_slip = (left=true, right=true, top=true, bot=true), ) @@ -77,7 +79,7 @@ end @test @views stokes.V.Vy[end, :] == stokes.V.Vy[end - 1, :] # no-slip - flow_bcs = FlowBoundaryConditions(; + flow_bcs = VelocityBoundaryConditions(; no_slip = (left=true, right=true, top=true, bot=true), free_slip = (left=false, right=false, top=false, bot=false), ) @@ -91,6 +93,83 @@ end @test @views stokes.V.Vy[end, :] == -stokes.V.Vy[end - 1, :] @test @views stokes.V.Vx[: , 1] == -stokes.V.Vx[: , 2] @test @views stokes.V.Vx[: , end] == -stokes.V.Vx[: , end - 1] + end + + @testset "DisplacementBoundaryConditions" begin + + # test incompatible boundary conditions + @test_throws ErrorException DisplacementBoundaryConditions(; + no_slip = (left=false, right=false, top=false, bot=false), + free_slip = (left=false, right=true, top=true, bot=true), + ) + @test_throws ErrorException DisplacementBoundaryConditions(; + no_slip = (left=false, right=false, top=false, bot=false), + free_slip = (left=true , right=true , top=true , bot=false), + ) + n = 5 # number of elements + Ux, Uy = PTArray(backend)(rand(n + 1, n + 2)), PTArray(backend)(rand(n + 2, n + 1)) + # free-slip + bcs1 = DisplacementBoundaryConditions(; + no_slip = (left=false, right=false, top=false, bot=false), + free_slip = (left=true, right=true, top=true, bot=true), + ) + flow_bcs!(bcs1, Ux, Uy) + @test @views Ux[: , 1] == Ux[:, 2] + @test @views Ux[: , end] == Ux[:, end - 1] + @test @views Uy[1 , :] == Uy[2, :] + @test @views Uy[end, :] == Uy[end - 1, :] + @test typeof(bcs) <: AbstractFlowBoundaryConditions + @test typeof(bcs) <: DisplacementBoundaryConditions + + # no-slip + Ux, Uy = PTArray(backend)(rand(n + 1, n + 2)), PTArray(backend)(rand(n + 2, n + 1)) + bcs2 = DisplacementBoundaryConditions(; + no_slip = (left=true, right=true, top=true, bot=true), + free_slip = (left=false, right=false, top=false, bot=false), + ) + flow_bcs!(bcs2, Ux, Uy) + @test sum(!iszero(Ux[1 , i]) for i in axes(Ux,2)) == 0 + @test sum(!iszero(Ux[end, i]) for i in axes(Ux,2)) == 0 + @test sum(!iszero(Uy[i, 1]) for i in axes(Uy,1)) == 0 + @test sum(!iszero(Uy[i, 1]) for i in axes(Uy,1)) == 0 + @test @views Uy[1 , :] == -Uy[2 , :] + @test @views Uy[end, :] == -Uy[end - 1, :] + @test @views Ux[: , 1] == -Ux[: , 2] + @test @views Ux[: , end] == -Ux[: , end - 1] + + # test with StokesArrays + ni = 5, 5 + stokes = StokesArrays(backend, ni) + stokes.U.Ux .= PTArray(backend)(rand(n + 1, n + 2)) + stokes.U.Uy .= PTArray(backend)(rand(n + 2, n + 1)) + # free-slip + flow_bcs = DisplacementBoundaryConditions(; + no_slip = (left=false, right=false, top=false, bot=false), + free_slip = (left=true, right=true, top=true, bot=true), + ) + flow_bcs!(stokes, flow_bcs) + + @test @views stokes.U.Ux[:, 1] == stokes.U.Ux[:, 2] + @test @views stokes.U.Ux[:, end] == stokes.U.Ux[:, end - 1] + @test @views stokes.U.Uy[1, :] == stokes.U.Uy[2, :] + @test @views stokes.U.Uy[end, :] == stokes.U.Uy[end - 1, :] + # no-slip + flow_bcs = DisplacementBoundaryConditions(; + no_slip = (left=true, right=true, top=true, bot=true), + free_slip = (left=false, right=false, top=false, bot=false), + ) + flow_bcs!(stokes, flow_bcs) + + @test sum(!iszero(stokes.U.Ux[1 , i]) for i in axes(Ux,2)) == 0 + @test sum(!iszero(stokes.U.Ux[end, i]) for i in axes(Ux,2)) == 0 + @test sum(!iszero(stokes.U.Uy[i, 1]) for i in axes(Uy,1)) == 0 + @test sum(!iszero(stokes.U.Uy[i, 1]) for i in axes(Uy,1)) == 0 + @test @views stokes.U.Uy[1 , :] == -stokes.U.Uy[2 , :] + @test @views stokes.U.Uy[end, :] == -stokes.U.Uy[end - 1, :] + @test @views stokes.U.Ux[: , 1] == -stokes.U.Ux[: , 2] + @test @views stokes.U.Ux[: , end] == -stokes.U.Ux[: , end - 1] + + end end else @test true === true diff --git a/test/test_shearband2D.jl b/test/test_shearband2D.jl index 1cea992c..c23d512b 100644 --- a/test/test_shearband2D.jl +++ b/test/test_shearband2D.jl @@ -122,14 +122,14 @@ function ShearBand2D() 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)...) # Time loop t, it = 0.0, 0 diff --git a/test/test_shearband2D_softening.jl b/test/test_shearband2D_softening.jl index 1010a489..ddd28188 100644 --- a/test/test_shearband2D_softening.jl +++ b/test/test_shearband2D_softening.jl @@ -129,14 +129,14 @@ function ShearBand2D() ) # 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)...) # Time loop t, it = 0.0, 0 diff --git a/test/test_shearheating2D.jl b/test/test_shearheating2D.jl index db6e1188..ef405da8 100644 --- a/test/test_shearheating2D.jl +++ b/test/test_shearheating2D.jl @@ -147,7 +147,7 @@ function Shearheating2D() ) # Boundary conditions - flow_bcs = FlowBoundaryConditions(; + flow_bcs = VelocityBoundaryConditions(; free_slip = (left = true, right=true, top=true, bot=true), ) ## Compression and not extension - fix this @@ -155,7 +155,7 @@ function Shearheating2D() stokes.V.Vx .= PTArray(backend_JR)([ -(x - lx/2) * εbg for x in xvi[1], _ in 1:ny+2]) stokes.V.Vy .= PTArray(backend_JR)([ (ly - abs(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)...) T_buffer = @zeros(ni.+1) Told_buffer = similar(T_buffer) diff --git a/test/test_shearheating3D.jl b/test/test_shearheating3D.jl index df614386..d6cf4675 100644 --- a/test/test_shearheating3D.jl +++ b/test/test_shearheating3D.jl @@ -141,7 +141,7 @@ function Shearheating3D(nx=16, ny=16, nz=16) ) # Boundary conditions - flow_bcs = FlowBoundaryConditions(; + flow_bcs = VelocityBoundaryConditions(; free_slip = (left = true , right = true , top = true , bot = true , front = true , back = true ), no_slip = (left = false, right = false, top = false, bot = false, front = false, back = false), ) @@ -151,7 +151,7 @@ function Shearheating3D(nx=16, ny=16, nz=16) stokes.V.Vy .= PTArray(backend_JR)([ -(y - ly/2) * εbg for _ in 1:nx+2, y in xvi[2], _ in 1:nz+2]) stokes.V.Vz .= PTArray(backend_JR)([ (lz - abs(z)) * εbg for _ in 1:nx+2, _ in 1:ny+2, z in xvi[3]]) flow_bcs!(stokes, flow_bcs) # apply boundary conditions - update_halo!(stokes.V.Vx, stokes.V.Vy, stokes.V.Vz) + update_halo!(@velocity(stokes)...) grid2particle!(pT, xvi, thermal.T, particles) diff --git a/test/test_sinking_block.jl b/test/test_sinking_block.jl index 6e8f156b..e381db90 100644 --- a/test/test_sinking_block.jl +++ b/test/test_sinking_block.jl @@ -166,11 +166,11 @@ function Sinking_Block2D() # ---------------------------------------------------- # 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)...) # Stokes solver ---------------- iters = solve!( diff --git a/test/test_types.jl b/test/test_types.jl index e0305b1f..cb091490 100644 --- a/test/test_types.jl +++ b/test/test_types.jl @@ -89,17 +89,17 @@ end stokes.V.Vy .= 1.0 JR2.velocity2displacement!(stokes, 10) - @test all(stokes.U.Ux.==0.1) + @test all(stokes.U.Ux.== 10) JR2.displacement2velocity!(stokes, 5) - @test all(stokes.V.Vx.==0.5) + @test all(stokes.V.Vx.==2.0) end @testset "3D allocators" begin ni = nx, ny, nz = (2, 2, 2) stokes = JR3.StokesArrays(backend, ni) - + @test size(stokes.P) == ni @test size(stokes.P0) == ni @test size(stokes.∇V) == ni @@ -174,8 +174,8 @@ end stokes.V.Vz .= 1.0 JR3.velocity2displacement!(stokes, 10) - @test all(stokes.U.Ux.==0.1) + @test all(stokes.U.Ux.==10.0) JR3.displacement2velocity!(stokes, 5) - @test all(stokes.V.Vx.==0.5) -end \ No newline at end of file + @test all(stokes.V.Vx.==2.0) +end From 36184c65904a1ae885e8d5fb68405e5c62a9817d Mon Sep 17 00:00:00 2001 From: Pascal Aellig Date: Thu, 4 Jul 2024 11:55:51 +0200 Subject: [PATCH 06/20] fix bug ofr VelocityBCs --- .../benchmarks/stokes2D/shear_band/ShearBand2D.jl | 2 +- src/stokes/Stokes2D.jl | 14 +++++++------- src/stokes/Stokes3D.jl | 14 +++++++------- src/types/displacement.jl | 9 +++++++++ 4 files changed, 24 insertions(+), 15 deletions(-) diff --git a/miniapps/benchmarks/stokes2D/shear_band/ShearBand2D.jl b/miniapps/benchmarks/stokes2D/shear_band/ShearBand2D.jl index 32affbd6..2b2ca50e 100644 --- a/miniapps/benchmarks/stokes2D/shear_band/ShearBand2D.jl +++ b/miniapps/benchmarks/stokes2D/shear_band/ShearBand2D.jl @@ -112,7 +112,7 @@ function main(igg; nx=64, ny=64, figdir="model_figs") 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) diff --git a/src/stokes/Stokes2D.jl b/src/stokes/Stokes2D.jl index e4c3d356..3f0ed2c5 100644 --- a/src/stokes/Stokes2D.jl +++ b/src/stokes/Stokes2D.jl @@ -55,7 +55,7 @@ function _solve!( norm_∇V = Float64[] # convert displacement to velocity - displacement2velocity!(stokes, dt) + displacement2velocity!(stokes, dt, flow_bcs) # solver loop wtime0 = 0.0 @@ -155,7 +155,7 @@ function _solve!( stokes::JustRelax.StokesArrays, pt_stokes, di::NTuple{2,T}, - flow_bcs, + flow_bcs::AbstractFlowBoundaryConditions, ρg, G, K, @@ -191,7 +191,7 @@ function _solve!( norm_∇V = Float64[] # convert displacement to velocity - displacement2velocity!(stokes, dt) + displacement2velocity!(stokes, dt, flow_bcs) # solver loop wtime0 = 0.0 @@ -278,7 +278,7 @@ function _solve!( stokes::JustRelax.StokesArrays, pt_stokes, di::NTuple{2,T}, - flow_bcs, + flow_bcs::AbstractFlowBoundaryConditions, ρg, rheology::MaterialParams, args, @@ -333,7 +333,7 @@ function _solve!( compute_viscosity!(stokes, phase_ratios, args, rheology, viscosity_cutoff) # convert displacement to velocity - displacement2velocity!(stokes, dt) + displacement2velocity!(stokes, dt, flow_bcs) while iter < 2 || (err > ϵ && iter ≤ iterMax) wtime0 += @elapsed begin @@ -460,7 +460,7 @@ function _solve!( stokes::JustRelax.StokesArrays, pt_stokes, di::NTuple{2,T}, - flow_bcs, + flow_bcs::AbstractFlowBoundaryConditions, ρg, phase_ratios::JustRelax.PhaseRatio, rheology, @@ -522,7 +522,7 @@ function _solve!( # compute buoyancy forces and viscosity compute_ρg!(ρg[end], phase_ratios, rheology, args) compute_viscosity!(stokes, phase_ratios, args, rheology, viscosity_cutoff) - displacement2velocity!(stokes, dt) + displacement2velocity!(stokes, dt, flow_bcs) while iter ≤ iterMax iterMin < iter && err < ϵ && break diff --git a/src/stokes/Stokes3D.jl b/src/stokes/Stokes3D.jl index 5abb5dbd..6ff82faa 100644 --- a/src/stokes/Stokes3D.jl +++ b/src/stokes/Stokes3D.jl @@ -26,7 +26,7 @@ function _solve!( stokes::JustRelax.StokesArrays, pt_stokes, di::NTuple{3,T}, - flow_bcs, + flwo_bcs::AbstractFlowBoundaryConditions, ρg, K, G, @@ -64,7 +64,7 @@ function _solve!( norm_∇V = Float64[] # convert displacement to velocity - displacement2velocity!(stokes, dt) + displacement2velocity!(stokes, dt, flow_bcs) # solver loop wtime0 = 0.0 @@ -212,7 +212,7 @@ function _solve!( compute_viscosity!(stokes, phase_ratios, args, rheology, viscosity_cutoff) # convert displacement to velocity - displacement2velocity!(stokes, dt) + displacement2velocity!(stokes, dt, flow_bcs) # solver loop wtime0 = 0.0 @@ -352,7 +352,7 @@ function _solve!( stokes::JustRelax.StokesArrays, pt_stokes, di::NTuple{3,T}, - flow_bc::AbstractFlowBoundaryConditions, + flow_bcs::AbstractFlowBoundaryConditions, ρg, phase_ratios::JustRelax.PhaseRatio, rheology::NTuple{N,AbstractMaterialParamsStruct}, @@ -401,7 +401,7 @@ function _solve!( compute_viscosity!(stokes, phase_ratios, args, rheology, viscosity_cutoff) # convert displacement to velocity - displacement2velocity!(stokes, dt) + displacement2velocity!(stokes, dt, flow_bcs) while iter < 2 || (err > ϵ && iter ≤ iterMax) wtime0 += @elapsed begin @@ -488,8 +488,8 @@ function _solve!( ) # apply boundary conditions velocity2displacement!(stokes, dt) - free_surface_bcs!(stokes, flow_bc, η, rheology, phase_ratios, dt, di) - flow_bcs!(stokes, flow_bc) + free_surface_bcs!(stokes, flow_bcs, η, rheology, phase_ratios, dt, di) + flow_bcs!(stokes, flow_bcs) update_halo!(@velocity(stokes)...) end end diff --git a/src/types/displacement.jl b/src/types/displacement.jl index f01a9fa7..1aa92d6e 100644 --- a/src/types/displacement.jl +++ b/src/types/displacement.jl @@ -29,6 +29,15 @@ end return nothing end +function displacement2velocity!(stokes,dt, flow_bcs::AbstractFlowBoundaryConditions) + if typeof(flow_bcs) <: DisplacementBoundaryConditions + displacement2velocity!(stokes, backend(stokes), dt) + return nothing + elseif typeof(flow_bcs) <: VelocityBoundaryConditions + return nothing + end +end + function displacement2velocity!(stokes::JustRelax.StokesArrays, dt) displacement2velocity!(stokes, backend(stokes), dt) return nothing From e39ad94ab0aeeaaefdc52734f8251f7f304bcc1e Mon Sep 17 00:00:00 2001 From: aelligp Date: Thu, 4 Jul 2024 15:08:31 +0200 Subject: [PATCH 07/20] update --- .../stokes2D/shear_band/ShearBand2D.jl | 2 +- .../shear_band/ShearBand2D_displacement.jl | 2 +- src/ext/AMDGPU/2D.jl | 24 +- src/ext/AMDGPU/3D.jl | 23 +- src/ext/CUDA/2D.jl | 16 +- src/ext/CUDA/3D.jl | 26 +- src/types/displacement.jl | 42 ++- test/test_boundary_conditions2D.jl | 310 +++++++++--------- 8 files changed, 220 insertions(+), 225 deletions(-) diff --git a/miniapps/benchmarks/stokes2D/shear_band/ShearBand2D.jl b/miniapps/benchmarks/stokes2D/shear_band/ShearBand2D.jl index 2b2ca50e..61728b78 100644 --- a/miniapps/benchmarks/stokes2D/shear_band/ShearBand2D.jl +++ b/miniapps/benchmarks/stokes2D/shear_band/ShearBand2D.jl @@ -88,7 +88,7 @@ function main(igg; nx=64, ny=64, figdir="model_figs") # 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 --------------------------------------------- diff --git a/miniapps/benchmarks/stokes2D/shear_band/ShearBand2D_displacement.jl b/miniapps/benchmarks/stokes2D/shear_band/ShearBand2D_displacement.jl index 101283d8..486d8bc8 100644 --- a/miniapps/benchmarks/stokes2D/shear_band/ShearBand2D_displacement.jl +++ b/miniapps/benchmarks/stokes2D/shear_band/ShearBand2D_displacement.jl @@ -88,7 +88,7 @@ function main(igg; nx=64, ny=64, figdir="model_figs") # 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 --------------------------------------------- diff --git a/src/ext/AMDGPU/2D.jl b/src/ext/AMDGPU/2D.jl index f68e20ed..4eb2e970 100644 --- a/src/ext/AMDGPU/2D.jl +++ b/src/ext/AMDGPU/2D.jl @@ -38,18 +38,6 @@ function JR2D.StokesArrays(::Type{AMDGPUBackend}, ni::NTuple{N,Integer}) where { return StokesArrays(ni) end -function JR2D.velocity2displacement!( - stokes::JustRelax.StokesArrays, ::AMDGPUBackendTrait, dt -) - return _velocity2displacement!(stokes, dt) -end - -function JR2D.displacement2velocity!( - stokes::JustRelax.StokesArrays, ::AMDGPUBackendTrait, dt -) - return _displacement2velocity!(stokes, dt) -end - function JR2D.ThermalArrays(::Type{AMDGPUBackend}, ni::NTuple{N,Number}) where {N} return ThermalArrays(ni...) end @@ -207,6 +195,18 @@ function JR2D.velocity2vertex!( return nothing end +function JR2D.velocity2displacement!( + ::AMDGPUBackendTrait, stokes::JustRelax.StokesArrays, dt +) + return _velocity2displacement!(stokes, dt) +end + +function JR2D.displacement2velocity!( + ::AMDGPUBackendTrait, stokes::JustRelax.StokesArrays, dt +) + return _displacement2velocity!(stokes, dt) +end + # Solvers function JR2D.solve!(::AMDGPUBackendTrait, stokes, args...; kwargs) return _solve!(stokes, args...; kwargs...) diff --git a/src/ext/AMDGPU/3D.jl b/src/ext/AMDGPU/3D.jl index 8464ee84..e75f1695 100644 --- a/src/ext/AMDGPU/3D.jl +++ b/src/ext/AMDGPU/3D.jl @@ -38,18 +38,6 @@ function JR3D.StokesArrays(::Type{AMDGPUBackend}, ni::NTuple{N,Integer}) where { return StokesArrays(ni) end -function JR3D.velocity2displacement!( - stokes::JustRelax.StokesArrays, ::AMDGPUBackendTrait, dt -) - return _velocity2displacement!(stokes, dt) -end - -function JR3D.displacement2velocity!( - stokes::JustRelax.StokesArrays, ::AMDGPUBackendTrait, dt -) - return _displacement2velocity!(stokes, dt) -end - function JR3D.ThermalArrays(::Type{AMDGPUBackend}, ni::NTuple{N,Number}) where {N} return ThermalArrays(ni...) end @@ -207,6 +195,17 @@ function JR3D.velocity2vertex!( return nothing end +function JR3D.velocity2displacement!( + ::AMDGPUBackendTrait, stokes::JustRelax.StokesArrays, dt +) + return _velocity2displacement!(stokes, dt) +end + +function JR3D.displacement2velocity!( + ::AMDGPUBackendTrait, stokes::JustRelax.StokesArrays, dt +) + return _displacement2velocity!(stokes, dt) +end # Solvers function JR3D.solve!(::AMDGPUBackendTrait, stokes, args...; kwargs) return _solve!(stokes, args...; kwargs...) diff --git a/src/ext/CUDA/2D.jl b/src/ext/CUDA/2D.jl index 5d968b65..f674bf67 100644 --- a/src/ext/CUDA/2D.jl +++ b/src/ext/CUDA/2D.jl @@ -37,14 +37,6 @@ function JR2D.StokesArrays(::Type{CUDABackend}, ni::NTuple{N,Integer}) where {N} return StokesArrays(ni) end -function JR2D.velocity2displacement!(stokes::JustRelax.StokesArrays, ::CUDABackendTrait, dt) - return _velocity2displacement!(stokes, dt) -end - -function JR2D.displacement2velocity!(stokes::JustRelax.StokesArrays, ::CUDABackendTrait, dt) - return _displacement2velocity!(stokes, dt) -end - function JR2D.ThermalArrays(::Type{CUDABackend}, ni::NTuple{N,Number}) where {N} return ThermalArrays(ni...) end @@ -203,6 +195,14 @@ function JR2D.velocity2vertex!( return nothing end +function JR2D.velocity2displacement!(::CUDABackendTrait, stokes::JustRelax.StokesArrays, dt) + return _velocity2displacement!(stokes, dt) +end + +function JR2D.displacement2velocity!(::CUDABackendTrait, stokes::JustRelax.StokesArrays, dt) + return _displacement2velocity!(stokes, dt) +end + # Solvers function JR2D.solve!(::CUDABackendTrait, stokes, args...; kwargs) return _solve!(stokes, args...; kwargs...) diff --git a/src/ext/CUDA/3D.jl b/src/ext/CUDA/3D.jl index 1694a53f..3d93486d 100644 --- a/src/ext/CUDA/3D.jl +++ b/src/ext/CUDA/3D.jl @@ -21,11 +21,11 @@ import JustRelax: Geometry, @cell import JustRelax: -AbstractBoundaryConditions, -TemperatureBoundaryConditions, -AbstractFlowBoundaryConditions, -DisplacementBoundaryConditions, -VelocityBoundaryConditions + AbstractBoundaryConditions, + TemperatureBoundaryConditions, + AbstractFlowBoundaryConditions, + DisplacementBoundaryConditions, + VelocityBoundaryConditions @init_parallel_stencil(CUDA, Float64, 3) @@ -38,14 +38,6 @@ function JR3D.StokesArrays(::Type{CUDABackend}, ni::NTuple{N,Integer}) where {N} return StokesArrays(ni) end -function JR3D.velocity2displacement!(stokes::JustRelax.StokesArrays, ::CUDABackendTrait, dt) - return _velocity2displacement!(stokes, dt) -end - -function JR3D.displacement2velocity!(stokes::JustRelax.StokesArrays, ::CUDABackendTrait, dt) - return _displacement2velocity!(stokes, dt) -end - function JR3D.ThermalArrays(::Type{CUDABackend}, ni::NTuple{N,Number}) where {N} return ThermalArrays(ni...) end @@ -215,6 +207,14 @@ function JR3D.velocity2vertex!( return nothing end +function JR3D.velocity2displacement!(::CUDABackendTrait, stokes::JustRelax.StokesArrays, dt) + return _velocity2displacement!(stokes, dt) +end + +function JR3D.displacement2velocity!(::CUDABackendTrait, stokes::JustRelax.StokesArrays, dt) + return _displacement2velocity!(stokes, dt) +end + # Solvers function JR3D.solve!(::CUDABackendTrait, stokes, args...; kwargs) return _solve!(stokes, args...; kwargs...) diff --git a/src/types/displacement.jl b/src/types/displacement.jl index 1aa92d6e..332d2fe9 100644 --- a/src/types/displacement.jl +++ b/src/types/displacement.jl @@ -1,22 +1,21 @@ -function velocity2displacement!(stokes::JustRelax.StokesArrays, dt) - velocity2displacement!(stokes, backend(stokes), dt) - return nothing -end +# Velocity to displacement interpolation +velocity2displacement!(stokes, dt) = velocity2displacement!(backend(stokes), stokes, dt) + -function velocity2displacement!(stokes::JustRelax.StokesArrays, ::CPUBackendTrait, dt) +function velocity2displacement!(::CPUBackendTrait, stokes::JustRelax.StokesArrays, dt) return _velocity2displacement!(stokes, dt) end function _velocity2displacement!(stokes::JustRelax.StokesArrays, dt) ni = size(stokes.P) (; V, U) = stokes - @parallel (@idx ni .+ 2) _velocity2displacement!( + @parallel (@idx ni .+ 2) _velocity2displacement_kernel!( V.Vx, V.Vy, V.Vz, U.Ux, U.Uy, U.Uz, dt ) return nothing end -@parallel_indices (I...) function _velocity2displacement!(Vx, Vy, Vz, Ux, Uy, Uz, dt) +@parallel_indices (I...) function _velocity2displacement_kernel!(Vx, Vy, Vz, Ux, Uy, Uz, dt) if all(I .≤ size(Ux)) Ux[I...] = Vx[I...] * dt end @@ -29,32 +28,22 @@ end return nothing end -function displacement2velocity!(stokes,dt, flow_bcs::AbstractFlowBoundaryConditions) - if typeof(flow_bcs) <: DisplacementBoundaryConditions - displacement2velocity!(stokes, backend(stokes), dt) - return nothing - elseif typeof(flow_bcs) <: VelocityBoundaryConditions - return nothing - end -end +# Displacement to velocity interpolation -function displacement2velocity!(stokes::JustRelax.StokesArrays, dt) - displacement2velocity!(stokes, backend(stokes), dt) - return nothing -end +displacement2velocity!(stokes, dt) = displacement2velocity!(backend(stokes), stokes::JustRelax.StokesArrays, dt) -function displacement2velocity!(stokes::JustRelax.StokesArrays, ::CPUBackendTrait, dt) +function displacement2velocity!(::CPUBackendTrait, stokes::JustRelax.StokesArrays, dt) return _displacement2velocity!(stokes, dt) end function _displacement2velocity!(stokes::JustRelax.StokesArrays, dt) ni = size(stokes.P) (; V, U) = stokes - @parallel (@idx ni .+ 2) _displacement2velocity!(U.Ux, U.Uy, U.Uz, V.Vx, V.Vy, V.Vz, 1 / dt) + @parallel (@idx ni .+ 2) _displacement2velocity_kernel!(U.Ux, U.Uy, U.Uz, V.Vx, V.Vy, V.Vz, 1 / dt) return nothing end -@parallel_indices (I...) function _displacement2velocity!(Ux, Uy, Uz, Vx, Vy, Vz, _dt) +@parallel_indices (I...) function _displacement2velocity_kernel!(Ux, Uy, Uz, Vx, Vy, Vz, _dt) if all(I .≤ size(Ux)) Vx[I...] = Ux[I...] * _dt end @@ -66,3 +55,12 @@ end end return nothing end + +function displacement2velocity!(stokes,dt, flow_bcs::AbstractFlowBoundaryConditions) + if typeof(flow_bcs) <: DisplacementBoundaryConditions + displacement2velocity!(backend(stokes), stokes, dt) + return nothing + elseif typeof(flow_bcs) <: VelocityBoundaryConditions + return nothing + end +end diff --git a/test/test_boundary_conditions2D.jl b/test/test_boundary_conditions2D.jl index 6fc169e7..e56c454d 100644 --- a/test/test_boundary_conditions2D.jl +++ b/test/test_boundary_conditions2D.jl @@ -1,7 +1,9 @@ @static if ENV["JULIA_JUSTRELAX_BACKEND"] === "AMDGPU" using AMDGPU + AMDGPU.allowscalar(true) elseif ENV["JULIA_JUSTRELAX_BACKEND"] === "CUDA" using CUDA + CUDA.allowscalar(true) end using JustRelax, JustRelax.JustRelax2D @@ -16,162 +18,158 @@ else end @testset "Boundary Conditions 2D" begin - if backend === CPUBackend - @suppress begin - @testset "VelocityBoundaryConditions" begin - - # test incompatible boundary conditions - @test_throws ErrorException VelocityBoundaryConditions(; - no_slip = (left=false, right=false, top=false, bot=false), - free_slip = (left=false, right=true, top=true, bot=true), - ) - @test_throws ErrorException VelocityBoundaryConditions(; - no_slip = (left=false, right=false, top=false, bot=false), - free_slip = (left=true , right=true , top=true , bot=false), - ) - - n = 5 # number of elements - Vx, Vy = PTArray(backend)(rand(n + 1, n + 2)), PTArray(backend)(rand(n + 2, n + 1)) - # free-slip - bcs = VelocityBoundaryConditions(; - no_slip = (left=false, right=false, top=false, bot=false), - free_slip = (left=true, right=true, top=true, bot=true), - ) - flow_bcs!(bcs, Vx, Vy) - - @test @views Vx[: , 1] == Vx[:, 2] - @test @views Vx[: , end] == Vx[:, end - 1] - @test @views Vy[1 , :] == Vy[2, :] - @test @views Vy[end, :] == Vy[end - 1, :] - @test typeof(bcs) <: AbstractFlowBoundaryConditions - @test typeof(bcs) <: VelocityBoundaryConditions - # no-slip - Vx, Vy = PTArray(backend)(rand(n + 1, n + 2)), PTArray(backend)(rand(n + 2, n + 1)) - bcs = VelocityBoundaryConditions(; - no_slip = (left=true, right=true, top=true, bot=true), - free_slip = (left=false, right=false, top=false, bot=false), - ) - flow_bcs!(bcs, Vx, Vy) - @test sum(!iszero(Vx[1 , i]) for i in axes(Vx,2)) == 0 - @test sum(!iszero(Vx[end, i]) for i in axes(Vx,2)) == 0 - @test sum(!iszero(Vy[i, 1]) for i in axes(Vy,1)) == 0 - @test sum(!iszero(Vy[i, 1]) for i in axes(Vy,1)) == 0 - @test @views Vy[1 , :] == -Vy[2 , :] - @test @views Vy[end, :] == -Vy[end - 1, :] - @test @views Vx[: , 1] == -Vx[: , 2] - @test @views Vx[: , end] == -Vx[: , end - 1] - - # test with StokesArrays - ni = 5, 5 - stokes = StokesArrays(backend, ni) - stokes.V.Vx .= PTArray(backend)(rand(n + 1, n + 2)) - stokes.V.Vy .= PTArray(backend)(rand(n + 2, n + 1)) - # free-slip - flow_bcs = VelocityBoundaryConditions(; - no_slip = (left=false, right=false, top=false, bot=false), - free_slip = (left=true, right=true, top=true, bot=true), - ) - flow_bcs!(stokes, flow_bcs) - - @test @views stokes.V.Vx[:, 1] == stokes.V.Vx[:, 2] - @test @views stokes.V.Vx[:, end] == stokes.V.Vx[:, end - 1] - @test @views stokes.V.Vy[1, :] == stokes.V.Vy[2, :] - @test @views stokes.V.Vy[end, :] == stokes.V.Vy[end - 1, :] - - # no-slip - flow_bcs = VelocityBoundaryConditions(; - no_slip = (left=true, right=true, top=true, bot=true), - free_slip = (left=false, right=false, top=false, bot=false), - ) - flow_bcs!(stokes, flow_bcs) - - @test sum(!iszero(stokes.V.Vx[1 , i]) for i in axes(Vx,2)) == 0 - @test sum(!iszero(stokes.V.Vx[end, i]) for i in axes(Vx,2)) == 0 - @test sum(!iszero(stokes.V.Vy[i, 1]) for i in axes(Vy,1)) == 0 - @test sum(!iszero(stokes.V.Vy[i, 1]) for i in axes(Vy,1)) == 0 - @test @views stokes.V.Vy[1 , :] == -stokes.V.Vy[2 , :] - @test @views stokes.V.Vy[end, :] == -stokes.V.Vy[end - 1, :] - @test @views stokes.V.Vx[: , 1] == -stokes.V.Vx[: , 2] - @test @views stokes.V.Vx[: , end] == -stokes.V.Vx[: , end - 1] - end - - @testset "DisplacementBoundaryConditions" begin - - # test incompatible boundary conditions - @test_throws ErrorException DisplacementBoundaryConditions(; - no_slip = (left=false, right=false, top=false, bot=false), - free_slip = (left=false, right=true, top=true, bot=true), - ) - @test_throws ErrorException DisplacementBoundaryConditions(; - no_slip = (left=false, right=false, top=false, bot=false), - free_slip = (left=true , right=true , top=true , bot=false), - ) - n = 5 # number of elements - Ux, Uy = PTArray(backend)(rand(n + 1, n + 2)), PTArray(backend)(rand(n + 2, n + 1)) - # free-slip - bcs1 = DisplacementBoundaryConditions(; - no_slip = (left=false, right=false, top=false, bot=false), - free_slip = (left=true, right=true, top=true, bot=true), - ) - flow_bcs!(bcs1, Ux, Uy) - @test @views Ux[: , 1] == Ux[:, 2] - @test @views Ux[: , end] == Ux[:, end - 1] - @test @views Uy[1 , :] == Uy[2, :] - @test @views Uy[end, :] == Uy[end - 1, :] - @test typeof(bcs) <: AbstractFlowBoundaryConditions - @test typeof(bcs) <: DisplacementBoundaryConditions - - # no-slip - Ux, Uy = PTArray(backend)(rand(n + 1, n + 2)), PTArray(backend)(rand(n + 2, n + 1)) - bcs2 = DisplacementBoundaryConditions(; - no_slip = (left=true, right=true, top=true, bot=true), - free_slip = (left=false, right=false, top=false, bot=false), - ) - flow_bcs!(bcs2, Ux, Uy) - @test sum(!iszero(Ux[1 , i]) for i in axes(Ux,2)) == 0 - @test sum(!iszero(Ux[end, i]) for i in axes(Ux,2)) == 0 - @test sum(!iszero(Uy[i, 1]) for i in axes(Uy,1)) == 0 - @test sum(!iszero(Uy[i, 1]) for i in axes(Uy,1)) == 0 - @test @views Uy[1 , :] == -Uy[2 , :] - @test @views Uy[end, :] == -Uy[end - 1, :] - @test @views Ux[: , 1] == -Ux[: , 2] - @test @views Ux[: , end] == -Ux[: , end - 1] - - # test with StokesArrays - ni = 5, 5 - stokes = StokesArrays(backend, ni) - stokes.U.Ux .= PTArray(backend)(rand(n + 1, n + 2)) - stokes.U.Uy .= PTArray(backend)(rand(n + 2, n + 1)) - # free-slip - flow_bcs = DisplacementBoundaryConditions(; - no_slip = (left=false, right=false, top=false, bot=false), - free_slip = (left=true, right=true, top=true, bot=true), - ) - flow_bcs!(stokes, flow_bcs) - - @test @views stokes.U.Ux[:, 1] == stokes.U.Ux[:, 2] - @test @views stokes.U.Ux[:, end] == stokes.U.Ux[:, end - 1] - @test @views stokes.U.Uy[1, :] == stokes.U.Uy[2, :] - @test @views stokes.U.Uy[end, :] == stokes.U.Uy[end - 1, :] - # no-slip - flow_bcs = DisplacementBoundaryConditions(; - no_slip = (left=true, right=true, top=true, bot=true), - free_slip = (left=false, right=false, top=false, bot=false), - ) - flow_bcs!(stokes, flow_bcs) - - @test sum(!iszero(stokes.U.Ux[1 , i]) for i in axes(Ux,2)) == 0 - @test sum(!iszero(stokes.U.Ux[end, i]) for i in axes(Ux,2)) == 0 - @test sum(!iszero(stokes.U.Uy[i, 1]) for i in axes(Uy,1)) == 0 - @test sum(!iszero(stokes.U.Uy[i, 1]) for i in axes(Uy,1)) == 0 - @test @views stokes.U.Uy[1 , :] == -stokes.U.Uy[2 , :] - @test @views stokes.U.Uy[end, :] == -stokes.U.Uy[end - 1, :] - @test @views stokes.U.Ux[: , 1] == -stokes.U.Ux[: , 2] - @test @views stokes.U.Ux[: , end] == -stokes.U.Ux[: , end - 1] - - end + @suppress begin + @testset "VelocityBoundaryConditions" begin + + # test incompatible boundary conditions + @test_throws ErrorException VelocityBoundaryConditions(; + no_slip = (left=false, right=false, top=false, bot=false), + free_slip = (left=false, right=true, top=true, bot=true), + ) + @test_throws ErrorException VelocityBoundaryConditions(; + no_slip = (left=false, right=false, top=false, bot=false), + free_slip = (left=true , right=true , top=true , bot=false), + ) + + n = 5 # number of elements + Vx, Vy = PTArray(backend)(rand(n + 1, n + 2)), PTArray(backend)(rand(n + 2, n + 1)) + # free-slip + bcs = VelocityBoundaryConditions(; + no_slip = (left=false, right=false, top=false, bot=false), + free_slip = (left=true, right=true, top=true, bot=true), + ) + flow_bcs!(bcs, Vx, Vy) + + @test @views Vx[: , 1] == Vx[:, 2] + @test @views Vx[: , end] == Vx[:, end - 1] + @test @views Vy[1 , :] == Vy[2, :] + @test @views Vy[end, :] == Vy[end - 1, :] + @test typeof(bcs) <: AbstractFlowBoundaryConditions + @test typeof(bcs) <: VelocityBoundaryConditions + # no-slip + Vx, Vy = PTArray(backend)(rand(n + 1, n + 2)), PTArray(backend)(rand(n + 2, n + 1)) + bcs = VelocityBoundaryConditions(; + no_slip = (left=true, right=true, top=true, bot=true), + free_slip = (left=false, right=false, top=false, bot=false), + ) + flow_bcs!(bcs, Vx, Vy) + @test sum(!iszero(Vx[1 , i]) for i in axes(Vx,2)) == 0 + @test sum(!iszero(Vx[end, i]) for i in axes(Vx,2)) == 0 + @test sum(!iszero(Vy[i, 1]) for i in axes(Vy,1)) == 0 + @test sum(!iszero(Vy[i, 1]) for i in axes(Vy,1)) == 0 + @test @views Vy[1 , :] == -Vy[2 , :] + @test @views Vy[end, :] == -Vy[end - 1, :] + @test @views Vx[: , 1] == -Vx[: , 2] + @test @views Vx[: , end] == -Vx[: , end - 1] + + # test with StokesArrays + ni = 5, 5 + stokes = StokesArrays(backend, ni) + stokes.V.Vx .= PTArray(backend)(rand(n + 1, n + 2)) + stokes.V.Vy .= PTArray(backend)(rand(n + 2, n + 1)) + # free-slip + flow_bcs = VelocityBoundaryConditions(; + no_slip = (left=false, right=false, top=false, bot=false), + free_slip = (left=true, right=true, top=true, bot=true), + ) + flow_bcs!(stokes, flow_bcs) + + @test @views stokes.V.Vx[:, 1] == stokes.V.Vx[:, 2] + @test @views stokes.V.Vx[:, end] == stokes.V.Vx[:, end - 1] + @test @views stokes.V.Vy[1, :] == stokes.V.Vy[2, :] + @test @views stokes.V.Vy[end, :] == stokes.V.Vy[end - 1, :] + + # no-slip + flow_bcs = VelocityBoundaryConditions(; + no_slip = (left=true, right=true, top=true, bot=true), + free_slip = (left=false, right=false, top=false, bot=false), + ) + flow_bcs!(stokes, flow_bcs) + + @test sum(!iszero(stokes.V.Vx[1 , i]) for i in axes(Vx,2)) == 0 + @test sum(!iszero(stokes.V.Vx[end, i]) for i in axes(Vx,2)) == 0 + @test sum(!iszero(stokes.V.Vy[i, 1]) for i in axes(Vy,1)) == 0 + @test sum(!iszero(stokes.V.Vy[i, 1]) for i in axes(Vy,1)) == 0 + @test @views stokes.V.Vy[1 , :] == -stokes.V.Vy[2 , :] + @test @views stokes.V.Vy[end, :] == -stokes.V.Vy[end - 1, :] + @test @views stokes.V.Vx[: , 1] == -stokes.V.Vx[: , 2] + @test @views stokes.V.Vx[: , end] == -stokes.V.Vx[: , end - 1] + end + + @testset "DisplacementBoundaryConditions" begin + + # test incompatible boundary conditions + @test_throws ErrorException DisplacementBoundaryConditions(; + no_slip = (left=false, right=false, top=false, bot=false), + free_slip = (left=false, right=true, top=true, bot=true), + ) + @test_throws ErrorException DisplacementBoundaryConditions(; + no_slip = (left=false, right=false, top=false, bot=false), + free_slip = (left=true , right=true , top=true , bot=false), + ) + n = 5 # number of elements + Ux, Uy = PTArray(backend)(rand(n + 1, n + 2)), PTArray(backend)(rand(n + 2, n + 1)) + # free-slip + bcs1 = DisplacementBoundaryConditions(; + no_slip = (left=false, right=false, top=false, bot=false), + free_slip = (left=true, right=true, top=true, bot=true), + ) + flow_bcs!(bcs1, Ux, Uy) + @test @views Ux[: , 1] == Ux[:, 2] + @test @views Ux[: , end] == Ux[:, end - 1] + @test @views Uy[1 , :] == Uy[2, :] + @test @views Uy[end, :] == Uy[end - 1, :] + @test typeof(bcs1) <: AbstractFlowBoundaryConditions + @test typeof(bcs1) <: DisplacementBoundaryConditions + + # no-slip + Ux, Uy = PTArray(backend)(rand(n + 1, n + 2)), PTArray(backend)(rand(n + 2, n + 1)) + bcs2 = DisplacementBoundaryConditions(; + no_slip = (left=true, right=true, top=true, bot=true), + free_slip = (left=false, right=false, top=false, bot=false), + ) + flow_bcs!(bcs2, Ux, Uy) + @test sum(!iszero(Ux[1 , i]) for i in axes(Ux,2)) == 0 + @test sum(!iszero(Ux[end, i]) for i in axes(Ux,2)) == 0 + @test sum(!iszero(Uy[i, 1]) for i in axes(Uy,1)) == 0 + @test sum(!iszero(Uy[i, 1]) for i in axes(Uy,1)) == 0 + @test @views Uy[1 , :] == -Uy[2 , :] + @test @views Uy[end, :] == -Uy[end - 1, :] + @test @views Ux[: , 1] == -Ux[: , 2] + @test @views Ux[: , end] == -Ux[: , end - 1] + + # test with StokesArrays + ni = 5, 5 + stokes = StokesArrays(backend, ni) + stokes.U.Ux .= PTArray(backend)(rand(n + 1, n + 2)) + stokes.U.Uy .= PTArray(backend)(rand(n + 2, n + 1)) + # free-slip + flow_bcs = DisplacementBoundaryConditions(; + no_slip = (left=false, right=false, top=false, bot=false), + free_slip = (left=true, right=true, top=true, bot=true), + ) + flow_bcs!(stokes, flow_bcs) + + @test @views stokes.U.Ux[:, 1] == stokes.U.Ux[:, 2] + @test @views stokes.U.Ux[:, end] == stokes.U.Ux[:, end - 1] + @test @views stokes.U.Uy[1, :] == stokes.U.Uy[2, :] + @test @views stokes.U.Uy[end, :] == stokes.U.Uy[end - 1, :] + # no-slip + flow_bcs = DisplacementBoundaryConditions(; + no_slip = (left=true, right=true, top=true, bot=true), + free_slip = (left=false, right=false, top=false, bot=false), + ) + flow_bcs!(stokes, flow_bcs) + + @test sum(!iszero(stokes.U.Ux[1 , i]) for i in axes(Ux,2)) == 0 + @test sum(!iszero(stokes.U.Ux[end, i]) for i in axes(Ux,2)) == 0 + @test sum(!iszero(stokes.U.Uy[i, 1]) for i in axes(Uy,1)) == 0 + @test sum(!iszero(stokes.U.Uy[i, 1]) for i in axes(Uy,1)) == 0 + @test @views stokes.U.Uy[1 , :] == -stokes.U.Uy[2 , :] + @test @views stokes.U.Uy[end, :] == -stokes.U.Uy[end - 1, :] + @test @views stokes.U.Ux[: , 1] == -stokes.U.Ux[: , 2] + @test @views stokes.U.Ux[: , end] == -stokes.U.Ux[: , end - 1] + end - else - @test true === true end end From 2d9f951a58f70660d7f6e491c5ea17911ff8ed55 Mon Sep 17 00:00:00 2001 From: aelligp Date: Fri, 5 Jul 2024 10:23:35 +0200 Subject: [PATCH 08/20] rm type --- src/types/displacement.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/types/displacement.jl b/src/types/displacement.jl index 332d2fe9..b1a3e7da 100644 --- a/src/types/displacement.jl +++ b/src/types/displacement.jl @@ -30,7 +30,7 @@ end # Displacement to velocity interpolation -displacement2velocity!(stokes, dt) = displacement2velocity!(backend(stokes), stokes::JustRelax.StokesArrays, dt) +displacement2velocity!(stokes, dt) = displacement2velocity!(backend(stokes), stokes, dt) function displacement2velocity!(::CPUBackendTrait, stokes::JustRelax.StokesArrays, dt) return _displacement2velocity!(stokes, dt) From f0e1aeb7db344448c3a3aa24edb2f038f9b2ea61 Mon Sep 17 00:00:00 2001 From: albert-de-montserrat Date: Fri, 5 Jul 2024 14:57:49 +0200 Subject: [PATCH 09/20] fix displacement conversion on CUDA --- src/ext/CUDA/2D.jl | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/src/ext/CUDA/2D.jl b/src/ext/CUDA/2D.jl index f674bf67..e814ad7c 100644 --- a/src/ext/CUDA/2D.jl +++ b/src/ext/CUDA/2D.jl @@ -199,10 +199,18 @@ function JR2D.velocity2displacement!(::CUDABackendTrait, stokes::JustRelax.Stoke return _velocity2displacement!(stokes, dt) end +function velocity2displacement!(::CUDABackendTrait, stokes::JustRelax.StokesArrays, dt) + return _velocity2displacement!(stokes, dt) +end + function JR2D.displacement2velocity!(::CUDABackendTrait, stokes::JustRelax.StokesArrays, dt) return _displacement2velocity!(stokes, dt) end +function displacement2velocity!(::CUDABackendTrait, stokes::JustRelax.StokesArrays, dt) + return _displacement2velocity!(stokes, dt) +end + # Solvers function JR2D.solve!(::CUDABackendTrait, stokes, args...; kwargs) return _solve!(stokes, args...; kwargs...) From 9599ffa30dd50cd3514479dfd33a1918afd2f02f Mon Sep 17 00:00:00 2001 From: aelligp Date: Fri, 5 Jul 2024 15:34:28 +0200 Subject: [PATCH 10/20] update docs --- docs/src/man/boundary_conditions.md | 54 +++++++++++++++++++++++++---- 1 file changed, 48 insertions(+), 6 deletions(-) diff --git a/docs/src/man/boundary_conditions.md b/docs/src/man/boundary_conditions.md index 98665708..4d414aaa 100644 --- a/docs/src/man/boundary_conditions.md +++ b/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 ) -``` \ No newline at end of file +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)...) +``` From d8a30a6d563a0c2c16d9b22355182d4b3661dd63 Mon Sep 17 00:00:00 2001 From: Pascal Aellig Date: Fri, 5 Jul 2024 15:47:35 +0200 Subject: [PATCH 11/20] format --- src/JustRelax.jl | 10 +++++----- src/Utils.jl | 3 ++- src/common.jl | 18 +++++++++--------- src/ext/AMDGPU/2D.jl | 20 ++++++++++++++++---- src/ext/AMDGPU/3D.jl | 30 +++++++++++++++++++++--------- src/ext/CUDA/2D.jl | 16 ++++++++++++---- src/ext/CUDA/3D.jl | 17 ++++++++++++----- src/types/displacement.jl | 13 ++++++++----- 8 files changed, 85 insertions(+), 42 deletions(-) diff --git a/src/JustRelax.jl b/src/JustRelax.jl index af675bc0..3730ac8f 100644 --- a/src/JustRelax.jl +++ b/src/JustRelax.jl @@ -32,11 +32,11 @@ include("types/phases.jl") # export PhaseRatio include("boundaryconditions/types.jl") -export AbstractBoundaryConditions, - TemperatureBoundaryConditions, - AbstractFlowBoundaryConditions, - DisplacementBoundaryConditions, - VelocityBoundaryConditions +export AbstractBoundaryConditions, + TemperatureBoundaryConditions, + AbstractFlowBoundaryConditions, + DisplacementBoundaryConditions, + VelocityBoundaryConditions include("types/traits.jl") export BackendTrait, CPUBackendTrait, NonCPUBackendTrait diff --git a/src/Utils.jl b/src/Utils.jl index 25f657ac..1d779b9d 100644 --- a/src/Utils.jl +++ b/src/Utils.jl @@ -115,7 +115,8 @@ macro displacement(A) end end -@inline unpack_displacement(U::JustRelax.Displacement{<:AbstractArray{T,2}}) where {T} = U.Ux, U.Uy +@inline unpack_displacement(U::JustRelax.Displacement{<:AbstractArray{T,2}}) where {T} = + U.Ux, U.Uy @inline unpack_displacement(U::JustRelax.Displacement{<:AbstractArray{T,3}}) where {T} = U.Ux, U.Uy, U.Uz diff --git a/src/common.jl b/src/common.jl index e8b68e25..d32e5a9a 100644 --- a/src/common.jl +++ b/src/common.jl @@ -32,15 +32,15 @@ include("types/displacement.jl") export velocity2displacement!, displacement2velocity! include("boundaryconditions/BoundaryConditions.jl") -export AbstractBoundaryConditions, - TemperatureBoundaryConditions, - AbstractFlowBoundaryConditions, - DisplacementBoundaryConditions, - VelocityBoundaryConditions, - flow_bcs!, - thermal_bcs!, - pureshear_bc!, - apply_free_slip! +export AbstractBoundaryConditions, + TemperatureBoundaryConditions, + AbstractFlowBoundaryConditions, + DisplacementBoundaryConditions, + VelocityBoundaryConditions, + flow_bcs!, + thermal_bcs!, + pureshear_bc!, + apply_free_slip! include("MiniKernels.jl") diff --git a/src/ext/AMDGPU/2D.jl b/src/ext/AMDGPU/2D.jl index 4eb2e970..806ca20e 100644 --- a/src/ext/AMDGPU/2D.jl +++ b/src/ext/AMDGPU/2D.jl @@ -86,19 +86,31 @@ function JR2D.PTThermalCoeffs( end # Boundary conditions -function JR2D.flow_bcs!(::AMDGPUBackendTrait, stokes::JustRelax.StokesArrays, bcs::VelocityBoundaryConditions) +function JR2D.flow_bcs!( + ::AMDGPUBackendTrait, stokes::JustRelax.StokesArrays, bcs::VelocityBoundaryConditions +) return _flow_bcs!(bcs, @velocity(stokes)) end -function flow_bcs!(::AMDGPUBackendTrait, stokes::JustRelax.StokesArrays, bcs::VelocityBoundaryConditions) +function flow_bcs!( + ::AMDGPUBackendTrait, stokes::JustRelax.StokesArrays, bcs::VelocityBoundaryConditions +) return _flow_bcs!(bcs, @velocity(stokes)) end -function JR2D.flow_bcs!(::AMDGPUBackendTrait, stokes::JustRelax.StokesArrays, bcs::DisplacementBoundaryConditions) +function JR2D.flow_bcs!( + ::AMDGPUBackendTrait, + stokes::JustRelax.StokesArrays, + bcs::DisplacementBoundaryConditions, +) return _flow_bcs!(bcs, @displacement(stokes)) end -function flow_bcs!(::AMDGPUBackendTrait, stokes::JustRelax.StokesArrays, bcs::DisplacementBoundaryConditions) +function flow_bcs!( + ::AMDGPUBackendTrait, + stokes::JustRelax.StokesArrays, + bcs::DisplacementBoundaryConditions, +) return _flow_bcs!(bcs, @displacement(stokes)) end diff --git a/src/ext/AMDGPU/3D.jl b/src/ext/AMDGPU/3D.jl index e75f1695..e045279b 100644 --- a/src/ext/AMDGPU/3D.jl +++ b/src/ext/AMDGPU/3D.jl @@ -22,11 +22,11 @@ import JustRelax: Geometry, @cell import JustRelax: -AbstractBoundaryConditions, -TemperatureBoundaryConditions, -AbstractFlowBoundaryConditions, -DisplacementBoundaryConditions, -VelocityBoundaryConditions + AbstractBoundaryConditions, + TemperatureBoundaryConditions, + AbstractFlowBoundaryConditions, + DisplacementBoundaryConditions, + VelocityBoundaryConditions @init_parallel_stencil(AMDGPU, Float64, 3) @@ -86,19 +86,31 @@ function JR3D.PTThermalCoeffs( end # Boundary conditions -function JR3D.flow_bcs!(::AMDGPUBackendTrait, stokes::JustRelax.StokesArrays, bcs::VelocityBoundaryConditions) +function JR3D.flow_bcs!( + ::AMDGPUBackendTrait, stokes::JustRelax.StokesArrays, bcs::VelocityBoundaryConditions +) return _flow_bcs!(bcs, @velocity(stokes)) end -function flow_bcs!(::AMDGPUBackendTrait, stokes::JustRelax.StokesArrays, bcs::VelocityBoundaryConditions) +function flow_bcs!( + ::AMDGPUBackendTrait, stokes::JustRelax.StokesArrays, bcs::VelocityBoundaryConditions +) return _flow_bcs!(bcs, @velocity(stokes)) end -function JR3D.flow_bcs!(::AMDGPUBackendTrait, stokes::JustRelax.StokesArrays, bcs::DisplacementBoundaryConditions) +function JR3D.flow_bcs!( + ::AMDGPUBackendTrait, + stokes::JustRelax.StokesArrays, + bcs::DisplacementBoundaryConditions, +) return _flow_bcs!(bcs, @displacement(stokes)) end -function flow_bcs!(::AMDGPUBackendTrait, stokes::JustRelax.StokesArrays, bcs::DisplacementBoundaryConditions) +function flow_bcs!( + ::AMDGPUBackendTrait, + stokes::JustRelax.StokesArrays, + bcs::DisplacementBoundaryConditions, +) return _flow_bcs!(bcs, @displacement(stokes)) end diff --git a/src/ext/CUDA/2D.jl b/src/ext/CUDA/2D.jl index e814ad7c..4c41ed61 100644 --- a/src/ext/CUDA/2D.jl +++ b/src/ext/CUDA/2D.jl @@ -85,19 +85,27 @@ function JR2D.PTThermalCoeffs( end # Boundary conditions -function JR2D.flow_bcs!(::CUDABackendTrait, stokes::JustRelax.StokesArrays, bcs::VelocityBoundaryConditions) +function JR2D.flow_bcs!( + ::CUDABackendTrait, stokes::JustRelax.StokesArrays, bcs::VelocityBoundaryConditions +) return _flow_bcs!(bcs, @velocity(stokes)) end -function flow_bcs!(::CUDABackendTrait, stokes::JustRelax.StokesArrays, bcs::VelocityBoundaryConditions) +function flow_bcs!( + ::CUDABackendTrait, stokes::JustRelax.StokesArrays, bcs::VelocityBoundaryConditions +) return _flow_bcs!(bcs, @velocity(stokes)) end -function JR2D.flow_bcs!(::CUDABackendTrait, stokes::JustRelax.StokesArrays, bcs::DisplacementBoundaryConditions) +function JR2D.flow_bcs!( + ::CUDABackendTrait, stokes::JustRelax.StokesArrays, bcs::DisplacementBoundaryConditions +) return _flow_bcs!(bcs, @displacement(stokes)) end -function flow_bcs!(::CUDABackendTrait, stokes::JustRelax.StokesArrays, bcs::DisplacementBoundaryConditions) +function flow_bcs!( + ::CUDABackendTrait, stokes::JustRelax.StokesArrays, bcs::DisplacementBoundaryConditions +) return _flow_bcs!(bcs, @velocity(stokes)) end diff --git a/src/ext/CUDA/3D.jl b/src/ext/CUDA/3D.jl index 3d93486d..984c45a4 100644 --- a/src/ext/CUDA/3D.jl +++ b/src/ext/CUDA/3D.jl @@ -27,7 +27,6 @@ import JustRelax: DisplacementBoundaryConditions, VelocityBoundaryConditions - @init_parallel_stencil(CUDA, Float64, 3) include("../../common.jl") @@ -97,20 +96,28 @@ function JR3D.update_pt_thermal_arrays!( end # Boundary conditions -function JR3D.flow_bcs!(::CUDABackendTrait, stokes::JustRelax.StokesArrays, bcs::VelocityBoundaryConditions) +function JR3D.flow_bcs!( + ::CUDABackendTrait, stokes::JustRelax.StokesArrays, bcs::VelocityBoundaryConditions +) return _flow_bcs!(bcs, @velocity(stokes)) end -function flow_bcs!(::CUDABackendTrait, stokes::JustRelax.StokesArrays, bcs::VelocityBoundaryConditions) +function flow_bcs!( + ::CUDABackendTrait, stokes::JustRelax.StokesArrays, bcs::VelocityBoundaryConditions +) return _flow_bcs!(bcs, @velocity(stokes)) end # Boundary conditions -function JR3D.flow_bcs!(::CUDABackendTrait, stokes::JustRelax.StokesArrays, bcs::DisplacementBoundaryConditions) +function JR3D.flow_bcs!( + ::CUDABackendTrait, stokes::JustRelax.StokesArrays, bcs::DisplacementBoundaryConditions +) return _flow_bcs!(bcs, @displacement(stokes)) end -function flow_bcs!(::CUDABackendTrait, stokes::JustRelax.StokesArrays, bcs::DisplacementBoundaryConditions) +function flow_bcs!( + ::CUDABackendTrait, stokes::JustRelax.StokesArrays, bcs::DisplacementBoundaryConditions +) return _flow_bcs!(bcs, @displacement(stokes)) end diff --git a/src/types/displacement.jl b/src/types/displacement.jl index b1a3e7da..6328ceda 100644 --- a/src/types/displacement.jl +++ b/src/types/displacement.jl @@ -1,7 +1,6 @@ # Velocity to displacement interpolation velocity2displacement!(stokes, dt) = velocity2displacement!(backend(stokes), stokes, dt) - function velocity2displacement!(::CPUBackendTrait, stokes::JustRelax.StokesArrays, dt) return _velocity2displacement!(stokes, dt) end @@ -30,7 +29,7 @@ end # Displacement to velocity interpolation -displacement2velocity!(stokes, dt) = displacement2velocity!(backend(stokes), stokes, dt) +displacement2velocity!(stokes, dt) = displacement2velocity!(backend(stokes), stokes, dt) function displacement2velocity!(::CPUBackendTrait, stokes::JustRelax.StokesArrays, dt) return _displacement2velocity!(stokes, dt) @@ -39,11 +38,15 @@ end function _displacement2velocity!(stokes::JustRelax.StokesArrays, dt) ni = size(stokes.P) (; V, U) = stokes - @parallel (@idx ni .+ 2) _displacement2velocity_kernel!(U.Ux, U.Uy, U.Uz, V.Vx, V.Vy, V.Vz, 1 / dt) + @parallel (@idx ni .+ 2) _displacement2velocity_kernel!( + U.Ux, U.Uy, U.Uz, V.Vx, V.Vy, V.Vz, 1 / dt + ) return nothing end -@parallel_indices (I...) function _displacement2velocity_kernel!(Ux, Uy, Uz, Vx, Vy, Vz, _dt) +@parallel_indices (I...) function _displacement2velocity_kernel!( + Ux, Uy, Uz, Vx, Vy, Vz, _dt +) if all(I .≤ size(Ux)) Vx[I...] = Ux[I...] * _dt end @@ -56,7 +59,7 @@ end return nothing end -function displacement2velocity!(stokes,dt, flow_bcs::AbstractFlowBoundaryConditions) +function displacement2velocity!(stokes, dt, flow_bcs::AbstractFlowBoundaryConditions) if typeof(flow_bcs) <: DisplacementBoundaryConditions displacement2velocity!(backend(stokes), stokes, dt) return nothing From 0f2c8cc1c30c83843d61ebae29ec5155c3091358 Mon Sep 17 00:00:00 2001 From: Pascal Aellig Date: Fri, 5 Jul 2024 16:15:01 +0200 Subject: [PATCH 12/20] fix ext --- src/ext/AMDGPU/2D.jl | 12 ++++++++++++ src/ext/AMDGPU/3D.jl | 12 ++++++++++++ src/ext/CUDA/3D.jl | 8 ++++++++ 3 files changed, 32 insertions(+) diff --git a/src/ext/AMDGPU/2D.jl b/src/ext/AMDGPU/2D.jl index 806ca20e..13e02b72 100644 --- a/src/ext/AMDGPU/2D.jl +++ b/src/ext/AMDGPU/2D.jl @@ -213,12 +213,24 @@ function JR2D.velocity2displacement!( return _velocity2displacement!(stokes, dt) end +function velocity2displacement!( + ::AMDGPUBackendTrait, stokes::JustRelax.StokesArrays, dt +) + return _velocity2displacement!(stokes, dt) +end + function JR2D.displacement2velocity!( ::AMDGPUBackendTrait, stokes::JustRelax.StokesArrays, dt ) return _displacement2velocity!(stokes, dt) end +function displacement2velocity!( + ::AMDGPUBackendTrait, stokes::JustRelax.StokesArrays, dt +) + return _displacement2velocity!(stokes, dt) +end + # Solvers function JR2D.solve!(::AMDGPUBackendTrait, stokes, args...; kwargs) return _solve!(stokes, args...; kwargs...) diff --git a/src/ext/AMDGPU/3D.jl b/src/ext/AMDGPU/3D.jl index e045279b..e51adc99 100644 --- a/src/ext/AMDGPU/3D.jl +++ b/src/ext/AMDGPU/3D.jl @@ -213,11 +213,23 @@ function JR3D.velocity2displacement!( return _velocity2displacement!(stokes, dt) end +function velocity2displacement!( + ::AMDGPUBackendTrait, stokes::JustRelax.StokesArrays, dt +) + return _velocity2displacement!(stokes, dt) +end + function JR3D.displacement2velocity!( ::AMDGPUBackendTrait, stokes::JustRelax.StokesArrays, dt ) return _displacement2velocity!(stokes, dt) end + +function displacement2velocity!( + ::AMDGPUBackendTrait, stokes::JustRelax.StokesArrays, dt +) + return _displacement2velocity!(stokes, dt) +end # Solvers function JR3D.solve!(::AMDGPUBackendTrait, stokes, args...; kwargs) return _solve!(stokes, args...; kwargs...) diff --git a/src/ext/CUDA/3D.jl b/src/ext/CUDA/3D.jl index 984c45a4..d22552b0 100644 --- a/src/ext/CUDA/3D.jl +++ b/src/ext/CUDA/3D.jl @@ -218,10 +218,18 @@ function JR3D.velocity2displacement!(::CUDABackendTrait, stokes::JustRelax.Stoke return _velocity2displacement!(stokes, dt) end +function velocity2displacement!(::CUDABackendTrait, stokes::JustRelax.StokesArrays, dt) + return _velocity2displacement!(stokes, dt) +end + function JR3D.displacement2velocity!(::CUDABackendTrait, stokes::JustRelax.StokesArrays, dt) return _displacement2velocity!(stokes, dt) end +function displacement2velocity!(::CUDABackendTrait, stokes::JustRelax.StokesArrays, dt) + return _displacement2velocity!(stokes, dt) +end + # Solvers function JR3D.solve!(::CUDABackendTrait, stokes, args...; kwargs) return _solve!(stokes, args...; kwargs...) From cb8e2c98eee26f0162bca8665d2453c0171eaae3 Mon Sep 17 00:00:00 2001 From: Pascal Aellig Date: Fri, 5 Jul 2024 16:17:37 +0200 Subject: [PATCH 13/20] format --- src/ext/AMDGPU/2D.jl | 8 ++------ src/ext/AMDGPU/3D.jl | 8 ++------ 2 files changed, 4 insertions(+), 12 deletions(-) diff --git a/src/ext/AMDGPU/2D.jl b/src/ext/AMDGPU/2D.jl index 13e02b72..2f494174 100644 --- a/src/ext/AMDGPU/2D.jl +++ b/src/ext/AMDGPU/2D.jl @@ -213,9 +213,7 @@ function JR2D.velocity2displacement!( return _velocity2displacement!(stokes, dt) end -function velocity2displacement!( - ::AMDGPUBackendTrait, stokes::JustRelax.StokesArrays, dt -) +function velocity2displacement!(::AMDGPUBackendTrait, stokes::JustRelax.StokesArrays, dt) return _velocity2displacement!(stokes, dt) end @@ -225,9 +223,7 @@ function JR2D.displacement2velocity!( return _displacement2velocity!(stokes, dt) end -function displacement2velocity!( - ::AMDGPUBackendTrait, stokes::JustRelax.StokesArrays, dt -) +function displacement2velocity!(::AMDGPUBackendTrait, stokes::JustRelax.StokesArrays, dt) return _displacement2velocity!(stokes, dt) end diff --git a/src/ext/AMDGPU/3D.jl b/src/ext/AMDGPU/3D.jl index e51adc99..2a29f01e 100644 --- a/src/ext/AMDGPU/3D.jl +++ b/src/ext/AMDGPU/3D.jl @@ -213,9 +213,7 @@ function JR3D.velocity2displacement!( return _velocity2displacement!(stokes, dt) end -function velocity2displacement!( - ::AMDGPUBackendTrait, stokes::JustRelax.StokesArrays, dt -) +function velocity2displacement!(::AMDGPUBackendTrait, stokes::JustRelax.StokesArrays, dt) return _velocity2displacement!(stokes, dt) end @@ -225,9 +223,7 @@ function JR3D.displacement2velocity!( return _displacement2velocity!(stokes, dt) end -function displacement2velocity!( - ::AMDGPUBackendTrait, stokes::JustRelax.StokesArrays, dt -) +function displacement2velocity!(::AMDGPUBackendTrait, stokes::JustRelax.StokesArrays, dt) return _displacement2velocity!(stokes, dt) end # Solvers From eeba3bc91e1671c91d81c15254aff3cde9fbac69 Mon Sep 17 00:00:00 2001 From: Pascal Aellig Date: Fri, 5 Jul 2024 16:38:53 +0200 Subject: [PATCH 14/20] fix test_BC3D --- test/test_boundary_conditions3D.jl | 91 +++++++++++++++++++++++++----- 1 file changed, 77 insertions(+), 14 deletions(-) diff --git a/test/test_boundary_conditions3D.jl b/test/test_boundary_conditions3D.jl index 9f960392..a089909f 100644 --- a/test/test_boundary_conditions3D.jl +++ b/test/test_boundary_conditions3D.jl @@ -1,7 +1,9 @@ @static if ENV["JULIA_JUSTRELAX_BACKEND"] === "AMDGPU" using AMDGPU + AMDGPU.allowscalar(true) elseif ENV["JULIA_JUSTRELAX_BACKEND"] === "CUDA" using CUDA + CUDA.allowscalar(true) end using JustRelax, JustRelax.JustRelax3D @@ -15,24 +17,23 @@ else CPUBackend end -@testset "Boundary Conditions" begin - if backend === CPUBackend - +@testset "Boundary Conditions 3D" begin + @testset "VelocityBoundaryConditions" begin # test incompatible boundary conditions - @test_throws ErrorException FlowBoundaryConditions(; + @test_throws ErrorException VelocityBoundaryConditions(; no_slip = (left=true, right=true, front=true, back=true, top=true, bot=true), free_slip = (left=false, right=true, front=true, back=true, top=true, bot=true), ) - + # test with StokesArrays ni = 5, 5, 5 stokes = StokesArrays(backend, ni) stokes.V.Vx .= PTArray(backend)(rand(size(stokes.V.Vx)...)) stokes.V.Vy .= PTArray(backend)(rand(size(stokes.V.Vy)...)) stokes.V.Vz .= PTArray(backend)(rand(size(stokes.V.Vz)...)) - + # free-slip - flow_bcs = FlowBoundaryConditions(; + flow_bcs = VelocityBoundaryConditions(; no_slip = (left=false, right=false, front=false, back=false, top=false, bot=false), free_slip = (left=true, right=true, front=true, back=true, top=true, bot=true), ) @@ -51,14 +52,14 @@ end @test @views stokes.V.Vz[end, :, :] == stokes.V.Vz[end - 1, :, :] @test @views stokes.V.Vz[ :, 1, :] == stokes.V.Vz[:, 2, :] @test @views stokes.V.Vz[ :, end, :] == stokes.V.Vz[:, end - 1, :] - + # no-slip - flow_bcs = FlowBoundaryConditions(; + flow_bcs = VelocityBoundaryConditions(; no_slip = (left=true, right=true, front=true, back=true, top=true, bot=true), free_slip = (left=false, right=false, front=false, back=false, top=false, bot=false), ) flow_bcs!(stokes, flow_bcs) - + (; Vx, Vy, Vz) = stokes.V @test sum(!iszero(Vx[1 , i, j]) for i in axes(Vx,2), j in axes(Vx,3)) == 0 @test sum(!iszero(Vx[end, i, j]) for i in axes(Vx,2), j in axes(Vx,3)) == 0 @@ -78,7 +79,69 @@ end @test @views Vz[ :, end, :] == -Vz[ :, end - 1, :] @test @views Vz[ 1, :, :] == -Vz[ 2, :, :] @test @views Vz[end, :, :] == -Vz[end - 1, :, :] - else - @test true === true - end -end \ No newline at end of file + end + + @testset "DisplacementBoundaryConditions" begin + + # test incompatible boundary conditions + @test_throws ErrorException DisplacementBoundaryConditions(; + no_slip = (left=true, right=true, front=true, back=true, top=true, bot=true), + free_slip = (left=false, right=true, front=true, back=true, top=true, bot=true), + ) + + # test with StokesArrays + ni = 5, 5, 5 + stokes = StokesArrays(backend, ni) + stokes.U.Ux .= PTArray(backend)(rand(size(stokes.U.Ux)...)) + stokes.U.Uy .= PTArray(backend)(rand(size(stokes.U.Uy)...)) + stokes.U.Uz .= PTArray(backend)(rand(size(stokes.U.Uz)...)) + + # free-slip + flow_bcs = DisplacementBoundaryConditions(; + no_slip = (left=false, right=false, front=false, back=false, top=false, bot=false), + free_slip = (left=true, right=true, front=true, back=true, top=true, bot=true), + ) + flow_bcs!(stokes, flow_bcs) + flow_bcs!(stokes, flow_bcs) # just a trick to pass the CI + + @test @views stokes.U.Ux[ :, :, 1] == stokes.U.Ux[:, :, 2] + @test @views stokes.U.Ux[ :, :, end] == stokes.U.Ux[:, :, end - 1] + @test @views stokes.U.Ux[ :, 1, :] == stokes.U.Ux[:, 2, :] + @test @views stokes.U.Ux[ :, end, :] == stokes.U.Ux[:, end - 1, :] + @test @views stokes.U.Uy[ :, :, 1] == stokes.U.Uy[:, :, 2] + @test @views stokes.U.Uy[ :, :, end] == stokes.U.Uy[:, :, end - 1] + @test @views stokes.U.Uy[ 1, :, :] == stokes.U.Uy[2, :, :] + @test @views stokes.U.Uy[end, :, :] == stokes.U.Uy[end - 1, :, :] + @test @views stokes.U.Uz[ 1, :, :] == stokes.U.Uz[2, :, :] + @test @views stokes.U.Uz[end, :, :] == stokes.U.Uz[end - 1, :, :] + @test @views stokes.U.Uz[ :, 1, :] == stokes.U.Uz[:, 2, :] + @test @views stokes.U.Uz[ :, end, :] == stokes.U.Uz[:, end - 1, :] + + # no-slip + flow_bcs = DisplacementBoundaryConditions(; + no_slip = (left=true, right=true, front=true, back=true, top=true, bot=true), + free_slip = (left=false, right=false, front=false, back=false, top=false, bot=false), + ) + flow_bcs!(stokes, flow_bcs) + + (; Ux, Uy, Uz) = stokes.U + @test sum(!iszero(Ux[1 , i, j]) for i in axes(Ux,2), j in axes(Ux,3)) == 0 + @test sum(!iszero(Ux[end, i, j]) for i in axes(Ux,2), j in axes(Ux,3)) == 0 + @test sum(!iszero(Uy[i, 1, j]) for i in axes(Uy,1), j in axes(Uy,2)) == 0 + @test sum(!iszero(Uy[i, end, j]) for i in axes(Uy,1), j in axes(Uy,2)) == 0 + @test sum(!iszero(Uz[i, j, 1]) for i in axes(Uz,1), j in axes(Uz,3)) == 0 + @test sum(!iszero(Uz[i, j, end]) for i in axes(Uz,1), j in axes(Uz,3)) == 0 + @test @views Ux[ :, 1, :] == -Ux[ :, 2, :] + @test @views Ux[ :, end, :] == -Ux[ :, end - 1, :] + @test @views Ux[ :, :, 1] == -Ux[ :, :, 2] + @test @views Ux[ :, :, end] == -Ux[ :, :, end - 1] + @test @views Uy[ 1, :, :] == -Uy[ 2, :, :] + @test @views Uy[end, :, :] == -Uy[end - 1, :, :] + @test @views Uy[ :, :, 1] == -Uy[ :, :, 2] + @test @views Uy[ :, :, end] == -Uy[ :, :, end - 1] + @test @views Uz[ :, 1, :] == -Uz[ :, 2, :] + @test @views Uz[ :, end, :] == -Uz[ :, end - 1, :] + @test @views Uz[ 1, :, :] == -Uz[ 2, :, :] + @test @views Uz[end, :, :] == -Uz[end - 1, :, :] + end +end From d5bbfc7cc9eec67fafd09e385ec8b05e91c562fd Mon Sep 17 00:00:00 2001 From: Pascal Aellig Date: Mon, 8 Jul 2024 10:30:32 +0200 Subject: [PATCH 15/20] fix typo; add security to docs --- docs/make.jl | 4 ++++ src/stokes/Stokes3D.jl | 2 +- 2 files changed, 5 insertions(+), 1 deletion(-) diff --git a/docs/make.jl b/docs/make.jl index 1163d67c..430c4a4f 100644 --- a/docs/make.jl +++ b/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" ], ) diff --git a/src/stokes/Stokes3D.jl b/src/stokes/Stokes3D.jl index 383e4969..4df97ea7 100644 --- a/src/stokes/Stokes3D.jl +++ b/src/stokes/Stokes3D.jl @@ -26,7 +26,7 @@ function _solve!( stokes::JustRelax.StokesArrays, pt_stokes, di::NTuple{3,T}, - flwo_bcs::AbstractFlowBoundaryConditions, + flow_bcs::AbstractFlowBoundaryConditions, ρg, K, G, From abbe64dfa36d2675ea2c9070e55ec4e0d6980863 Mon Sep 17 00:00:00 2001 From: Pascal Aellig Date: Mon, 8 Jul 2024 11:33:28 +0200 Subject: [PATCH 16/20] address requested changes --- src/JustRelax.jl | 4 +--- src/types/displacement.jl | 11 +++-------- 2 files changed, 4 insertions(+), 11 deletions(-) diff --git a/src/JustRelax.jl b/src/JustRelax.jl index 3730ac8f..71d40d03 100644 --- a/src/JustRelax.jl +++ b/src/JustRelax.jl @@ -32,9 +32,7 @@ include("types/phases.jl") # export PhaseRatio include("boundaryconditions/types.jl") -export AbstractBoundaryConditions, - TemperatureBoundaryConditions, - AbstractFlowBoundaryConditions, +export TemperatureBoundaryConditions, DisplacementBoundaryConditions, VelocityBoundaryConditions diff --git a/src/types/displacement.jl b/src/types/displacement.jl index 6328ceda..85decb6b 100644 --- a/src/types/displacement.jl +++ b/src/types/displacement.jl @@ -59,11 +59,6 @@ end return nothing end -function displacement2velocity!(stokes, dt, flow_bcs::AbstractFlowBoundaryConditions) - if typeof(flow_bcs) <: DisplacementBoundaryConditions - displacement2velocity!(backend(stokes), stokes, dt) - return nothing - elseif typeof(flow_bcs) <: VelocityBoundaryConditions - return nothing - end -end +displacement2velocity!(stokes, dt, <:DisplacementBoundaryConditions) = displacement2velocity!(backend(stokes), stokes, dt) +displacement2velocity!(::Any, ::Any, <:VelocityBoundaryConditions) = nothing +displacement2velocity!(::Any, ::Any, ::T) where T = throw(ArgumentError("Unknown boundary conditions type: $T")) From e339ab20c31504d4a8967c3d160e57e47c63d3a2 Mon Sep 17 00:00:00 2001 From: Pascal Aellig Date: Mon, 8 Jul 2024 11:33:47 +0200 Subject: [PATCH 17/20] revert backto CPU backend --- test/test_boundary_conditions2D.jl | 299 +++++++++++++++-------------- test/test_boundary_conditions3D.jl | 229 +++++++++++----------- 2 files changed, 270 insertions(+), 258 deletions(-) diff --git a/test/test_boundary_conditions2D.jl b/test/test_boundary_conditions2D.jl index e56c454d..dbb0d6b6 100644 --- a/test/test_boundary_conditions2D.jl +++ b/test/test_boundary_conditions2D.jl @@ -20,156 +20,161 @@ end @testset "Boundary Conditions 2D" begin @suppress begin @testset "VelocityBoundaryConditions" begin - - # test incompatible boundary conditions - @test_throws ErrorException VelocityBoundaryConditions(; - no_slip = (left=false, right=false, top=false, bot=false), - free_slip = (left=false, right=true, top=true, bot=true), - ) - @test_throws ErrorException VelocityBoundaryConditions(; - no_slip = (left=false, right=false, top=false, bot=false), - free_slip = (left=true , right=true , top=true , bot=false), - ) - - n = 5 # number of elements - Vx, Vy = PTArray(backend)(rand(n + 1, n + 2)), PTArray(backend)(rand(n + 2, n + 1)) - # free-slip - bcs = VelocityBoundaryConditions(; - no_slip = (left=false, right=false, top=false, bot=false), - free_slip = (left=true, right=true, top=true, bot=true), - ) - flow_bcs!(bcs, Vx, Vy) - - @test @views Vx[: , 1] == Vx[:, 2] - @test @views Vx[: , end] == Vx[:, end - 1] - @test @views Vy[1 , :] == Vy[2, :] - @test @views Vy[end, :] == Vy[end - 1, :] - @test typeof(bcs) <: AbstractFlowBoundaryConditions - @test typeof(bcs) <: VelocityBoundaryConditions - # no-slip - Vx, Vy = PTArray(backend)(rand(n + 1, n + 2)), PTArray(backend)(rand(n + 2, n + 1)) - bcs = VelocityBoundaryConditions(; - no_slip = (left=true, right=true, top=true, bot=true), - free_slip = (left=false, right=false, top=false, bot=false), - ) - flow_bcs!(bcs, Vx, Vy) - @test sum(!iszero(Vx[1 , i]) for i in axes(Vx,2)) == 0 - @test sum(!iszero(Vx[end, i]) for i in axes(Vx,2)) == 0 - @test sum(!iszero(Vy[i, 1]) for i in axes(Vy,1)) == 0 - @test sum(!iszero(Vy[i, 1]) for i in axes(Vy,1)) == 0 - @test @views Vy[1 , :] == -Vy[2 , :] - @test @views Vy[end, :] == -Vy[end - 1, :] - @test @views Vx[: , 1] == -Vx[: , 2] - @test @views Vx[: , end] == -Vx[: , end - 1] - - # test with StokesArrays - ni = 5, 5 - stokes = StokesArrays(backend, ni) - stokes.V.Vx .= PTArray(backend)(rand(n + 1, n + 2)) - stokes.V.Vy .= PTArray(backend)(rand(n + 2, n + 1)) - # free-slip - flow_bcs = VelocityBoundaryConditions(; - no_slip = (left=false, right=false, top=false, bot=false), - free_slip = (left=true, right=true, top=true, bot=true), - ) - flow_bcs!(stokes, flow_bcs) - - @test @views stokes.V.Vx[:, 1] == stokes.V.Vx[:, 2] - @test @views stokes.V.Vx[:, end] == stokes.V.Vx[:, end - 1] - @test @views stokes.V.Vy[1, :] == stokes.V.Vy[2, :] - @test @views stokes.V.Vy[end, :] == stokes.V.Vy[end - 1, :] - - # no-slip - flow_bcs = VelocityBoundaryConditions(; - no_slip = (left=true, right=true, top=true, bot=true), - free_slip = (left=false, right=false, top=false, bot=false), - ) - flow_bcs!(stokes, flow_bcs) - - @test sum(!iszero(stokes.V.Vx[1 , i]) for i in axes(Vx,2)) == 0 - @test sum(!iszero(stokes.V.Vx[end, i]) for i in axes(Vx,2)) == 0 - @test sum(!iszero(stokes.V.Vy[i, 1]) for i in axes(Vy,1)) == 0 - @test sum(!iszero(stokes.V.Vy[i, 1]) for i in axes(Vy,1)) == 0 - @test @views stokes.V.Vy[1 , :] == -stokes.V.Vy[2 , :] - @test @views stokes.V.Vy[end, :] == -stokes.V.Vy[end - 1, :] - @test @views stokes.V.Vx[: , 1] == -stokes.V.Vx[: , 2] - @test @views stokes.V.Vx[: , end] == -stokes.V.Vx[: , end - 1] + if backend === CPUBackend + # test incompatible boundary conditions + @test_throws ErrorException VelocityBoundaryConditions(; + no_slip = (left=false, right=false, top=false, bot=false), + free_slip = (left=false, right=true, top=true, bot=true), + ) + @test_throws ErrorException VelocityBoundaryConditions(; + no_slip = (left=false, right=false, top=false, bot=false), + free_slip = (left=true , right=true , top=true , bot=false), + ) + + n = 5 # number of elements + Vx, Vy = PTArray(backend)(rand(n + 1, n + 2)), PTArray(backend)(rand(n + 2, n + 1)) + # free-slip + bcs = VelocityBoundaryConditions(; + no_slip = (left=false, right=false, top=false, bot=false), + free_slip = (left=true, right=true, top=true, bot=true), + ) + flow_bcs!(bcs, Vx, Vy) + + @test @views Vx[: , 1] == Vx[:, 2] + @test @views Vx[: , end] == Vx[:, end - 1] + @test @views Vy[1 , :] == Vy[2, :] + @test @views Vy[end, :] == Vy[end - 1, :] + @test typeof(bcs) <: AbstractFlowBoundaryConditions + @test typeof(bcs) <: VelocityBoundaryConditions + # no-slip + Vx, Vy = PTArray(backend)(rand(n + 1, n + 2)), PTArray(backend)(rand(n + 2, n + 1)) + bcs = VelocityBoundaryConditions(; + no_slip = (left=true, right=true, top=true, bot=true), + free_slip = (left=false, right=false, top=false, bot=false), + ) + flow_bcs!(bcs, Vx, Vy) + @test sum(!iszero(Vx[1 , i]) for i in axes(Vx,2)) == 0 + @test sum(!iszero(Vx[end, i]) for i in axes(Vx,2)) == 0 + @test sum(!iszero(Vy[i, 1]) for i in axes(Vy,1)) == 0 + @test sum(!iszero(Vy[i, 1]) for i in axes(Vy,1)) == 0 + @test @views Vy[1 , :] == -Vy[2 , :] + @test @views Vy[end, :] == -Vy[end - 1, :] + @test @views Vx[: , 1] == -Vx[: , 2] + @test @views Vx[: , end] == -Vx[: , end - 1] + + # test with StokesArrays + ni = 5, 5 + stokes = StokesArrays(backend, ni) + stokes.V.Vx .= PTArray(backend)(rand(n + 1, n + 2)) + stokes.V.Vy .= PTArray(backend)(rand(n + 2, n + 1)) + # free-slip + flow_bcs = VelocityBoundaryConditions(; + no_slip = (left=false, right=false, top=false, bot=false), + free_slip = (left=true, right=true, top=true, bot=true), + ) + flow_bcs!(stokes, flow_bcs) + + @test @views stokes.V.Vx[:, 1] == stokes.V.Vx[:, 2] + @test @views stokes.V.Vx[:, end] == stokes.V.Vx[:, end - 1] + @test @views stokes.V.Vy[1, :] == stokes.V.Vy[2, :] + @test @views stokes.V.Vy[end, :] == stokes.V.Vy[end - 1, :] + + # no-slip + flow_bcs = VelocityBoundaryConditions(; + no_slip = (left=true, right=true, top=true, bot=true), + free_slip = (left=false, right=false, top=false, bot=false), + ) + flow_bcs!(stokes, flow_bcs) + + @test sum(!iszero(stokes.V.Vx[1 , i]) for i in axes(Vx,2)) == 0 + @test sum(!iszero(stokes.V.Vx[end, i]) for i in axes(Vx,2)) == 0 + @test sum(!iszero(stokes.V.Vy[i, 1]) for i in axes(Vy,1)) == 0 + @test sum(!iszero(stokes.V.Vy[i, 1]) for i in axes(Vy,1)) == 0 + @test @views stokes.V.Vy[1 , :] == -stokes.V.Vy[2 , :] + @test @views stokes.V.Vy[end, :] == -stokes.V.Vy[end - 1, :] + @test @views stokes.V.Vx[: , 1] == -stokes.V.Vx[: , 2] + @test @views stokes.V.Vx[: , end] == -stokes.V.Vx[: , end - 1] + else + @test true === true + end end @testset "DisplacementBoundaryConditions" begin - - # test incompatible boundary conditions - @test_throws ErrorException DisplacementBoundaryConditions(; - no_slip = (left=false, right=false, top=false, bot=false), - free_slip = (left=false, right=true, top=true, bot=true), - ) - @test_throws ErrorException DisplacementBoundaryConditions(; - no_slip = (left=false, right=false, top=false, bot=false), - free_slip = (left=true , right=true , top=true , bot=false), - ) - n = 5 # number of elements - Ux, Uy = PTArray(backend)(rand(n + 1, n + 2)), PTArray(backend)(rand(n + 2, n + 1)) - # free-slip - bcs1 = DisplacementBoundaryConditions(; - no_slip = (left=false, right=false, top=false, bot=false), - free_slip = (left=true, right=true, top=true, bot=true), - ) - flow_bcs!(bcs1, Ux, Uy) - @test @views Ux[: , 1] == Ux[:, 2] - @test @views Ux[: , end] == Ux[:, end - 1] - @test @views Uy[1 , :] == Uy[2, :] - @test @views Uy[end, :] == Uy[end - 1, :] - @test typeof(bcs1) <: AbstractFlowBoundaryConditions - @test typeof(bcs1) <: DisplacementBoundaryConditions - - # no-slip - Ux, Uy = PTArray(backend)(rand(n + 1, n + 2)), PTArray(backend)(rand(n + 2, n + 1)) - bcs2 = DisplacementBoundaryConditions(; - no_slip = (left=true, right=true, top=true, bot=true), - free_slip = (left=false, right=false, top=false, bot=false), - ) - flow_bcs!(bcs2, Ux, Uy) - @test sum(!iszero(Ux[1 , i]) for i in axes(Ux,2)) == 0 - @test sum(!iszero(Ux[end, i]) for i in axes(Ux,2)) == 0 - @test sum(!iszero(Uy[i, 1]) for i in axes(Uy,1)) == 0 - @test sum(!iszero(Uy[i, 1]) for i in axes(Uy,1)) == 0 - @test @views Uy[1 , :] == -Uy[2 , :] - @test @views Uy[end, :] == -Uy[end - 1, :] - @test @views Ux[: , 1] == -Ux[: , 2] - @test @views Ux[: , end] == -Ux[: , end - 1] - - # test with StokesArrays - ni = 5, 5 - stokes = StokesArrays(backend, ni) - stokes.U.Ux .= PTArray(backend)(rand(n + 1, n + 2)) - stokes.U.Uy .= PTArray(backend)(rand(n + 2, n + 1)) - # free-slip - flow_bcs = DisplacementBoundaryConditions(; - no_slip = (left=false, right=false, top=false, bot=false), - free_slip = (left=true, right=true, top=true, bot=true), - ) - flow_bcs!(stokes, flow_bcs) - - @test @views stokes.U.Ux[:, 1] == stokes.U.Ux[:, 2] - @test @views stokes.U.Ux[:, end] == stokes.U.Ux[:, end - 1] - @test @views stokes.U.Uy[1, :] == stokes.U.Uy[2, :] - @test @views stokes.U.Uy[end, :] == stokes.U.Uy[end - 1, :] - # no-slip - flow_bcs = DisplacementBoundaryConditions(; - no_slip = (left=true, right=true, top=true, bot=true), - free_slip = (left=false, right=false, top=false, bot=false), - ) - flow_bcs!(stokes, flow_bcs) - - @test sum(!iszero(stokes.U.Ux[1 , i]) for i in axes(Ux,2)) == 0 - @test sum(!iszero(stokes.U.Ux[end, i]) for i in axes(Ux,2)) == 0 - @test sum(!iszero(stokes.U.Uy[i, 1]) for i in axes(Uy,1)) == 0 - @test sum(!iszero(stokes.U.Uy[i, 1]) for i in axes(Uy,1)) == 0 - @test @views stokes.U.Uy[1 , :] == -stokes.U.Uy[2 , :] - @test @views stokes.U.Uy[end, :] == -stokes.U.Uy[end - 1, :] - @test @views stokes.U.Ux[: , 1] == -stokes.U.Ux[: , 2] - @test @views stokes.U.Ux[: , end] == -stokes.U.Ux[: , end - 1] - + if backend === CPUBackend + # test incompatible boundary conditions + @test_throws ErrorException DisplacementBoundaryConditions(; + no_slip = (left=false, right=false, top=false, bot=false), + free_slip = (left=false, right=true, top=true, bot=true), + ) + @test_throws ErrorException DisplacementBoundaryConditions(; + no_slip = (left=false, right=false, top=false, bot=false), + free_slip = (left=true , right=true , top=true , bot=false), + ) + n = 5 # number of elements + Ux, Uy = PTArray(backend)(rand(n + 1, n + 2)), PTArray(backend)(rand(n + 2, n + 1)) + # free-slip + bcs1 = DisplacementBoundaryConditions(; + no_slip = (left=false, right=false, top=false, bot=false), + free_slip = (left=true, right=true, top=true, bot=true), + ) + flow_bcs!(bcs1, Ux, Uy) + @test @views Ux[: , 1] == Ux[:, 2] + @test @views Ux[: , end] == Ux[:, end - 1] + @test @views Uy[1 , :] == Uy[2, :] + @test @views Uy[end, :] == Uy[end - 1, :] + @test typeof(bcs1) <: AbstractFlowBoundaryConditions + @test typeof(bcs1) <: DisplacementBoundaryConditions + + # no-slip + Ux, Uy = PTArray(backend)(rand(n + 1, n + 2)), PTArray(backend)(rand(n + 2, n + 1)) + bcs2 = DisplacementBoundaryConditions(; + no_slip = (left=true, right=true, top=true, bot=true), + free_slip = (left=false, right=false, top=false, bot=false), + ) + flow_bcs!(bcs2, Ux, Uy) + @test sum(!iszero(Ux[1 , i]) for i in axes(Ux,2)) == 0 + @test sum(!iszero(Ux[end, i]) for i in axes(Ux,2)) == 0 + @test sum(!iszero(Uy[i, 1]) for i in axes(Uy,1)) == 0 + @test sum(!iszero(Uy[i, 1]) for i in axes(Uy,1)) == 0 + @test @views Uy[1 , :] == -Uy[2 , :] + @test @views Uy[end, :] == -Uy[end - 1, :] + @test @views Ux[: , 1] == -Ux[: , 2] + @test @views Ux[: , end] == -Ux[: , end - 1] + + # test with StokesArrays + ni = 5, 5 + stokes = StokesArrays(backend, ni) + stokes.U.Ux .= PTArray(backend)(rand(n + 1, n + 2)) + stokes.U.Uy .= PTArray(backend)(rand(n + 2, n + 1)) + # free-slip + flow_bcs = DisplacementBoundaryConditions(; + no_slip = (left=false, right=false, top=false, bot=false), + free_slip = (left=true, right=true, top=true, bot=true), + ) + flow_bcs!(stokes, flow_bcs) + + @test @views stokes.U.Ux[:, 1] == stokes.U.Ux[:, 2] + @test @views stokes.U.Ux[:, end] == stokes.U.Ux[:, end - 1] + @test @views stokes.U.Uy[1, :] == stokes.U.Uy[2, :] + @test @views stokes.U.Uy[end, :] == stokes.U.Uy[end - 1, :] + # no-slip + flow_bcs = DisplacementBoundaryConditions(; + no_slip = (left=true, right=true, top=true, bot=true), + free_slip = (left=false, right=false, top=false, bot=false), + ) + flow_bcs!(stokes, flow_bcs) + + @test sum(!iszero(stokes.U.Ux[1 , i]) for i in axes(Ux,2)) == 0 + @test sum(!iszero(stokes.U.Ux[end, i]) for i in axes(Ux,2)) == 0 + @test sum(!iszero(stokes.U.Uy[i, 1]) for i in axes(Uy,1)) == 0 + @test sum(!iszero(stokes.U.Uy[i, 1]) for i in axes(Uy,1)) == 0 + @test @views stokes.U.Uy[1 , :] == -stokes.U.Uy[2 , :] + @test @views stokes.U.Uy[end, :] == -stokes.U.Uy[end - 1, :] + @test @views stokes.U.Ux[: , 1] == -stokes.U.Ux[: , 2] + @test @views stokes.U.Ux[: , end] == -stokes.U.Ux[: , end - 1] + else + @test true === true + end end end end diff --git a/test/test_boundary_conditions3D.jl b/test/test_boundary_conditions3D.jl index a089909f..88aa96d6 100644 --- a/test/test_boundary_conditions3D.jl +++ b/test/test_boundary_conditions3D.jl @@ -19,129 +19,136 @@ end @testset "Boundary Conditions 3D" begin @testset "VelocityBoundaryConditions" begin - # test incompatible boundary conditions - @test_throws ErrorException VelocityBoundaryConditions(; - no_slip = (left=true, right=true, front=true, back=true, top=true, bot=true), - free_slip = (left=false, right=true, front=true, back=true, top=true, bot=true), - ) + if backend === CPUBackend + # test incompatible boundary conditions + @test_throws ErrorException VelocityBoundaryConditions(; + no_slip = (left=true, right=true, front=true, back=true, top=true, bot=true), + free_slip = (left=false, right=true, front=true, back=true, top=true, bot=true), + ) - # test with StokesArrays - ni = 5, 5, 5 - stokes = StokesArrays(backend, ni) - stokes.V.Vx .= PTArray(backend)(rand(size(stokes.V.Vx)...)) - stokes.V.Vy .= PTArray(backend)(rand(size(stokes.V.Vy)...)) - stokes.V.Vz .= PTArray(backend)(rand(size(stokes.V.Vz)...)) + # test with StokesArrays + ni = 5, 5, 5 + stokes = StokesArrays(backend, ni) + stokes.V.Vx .= PTArray(backend)(rand(size(stokes.V.Vx)...)) + stokes.V.Vy .= PTArray(backend)(rand(size(stokes.V.Vy)...)) + stokes.V.Vz .= PTArray(backend)(rand(size(stokes.V.Vz)...)) - # free-slip - flow_bcs = VelocityBoundaryConditions(; - no_slip = (left=false, right=false, front=false, back=false, top=false, bot=false), - free_slip = (left=true, right=true, front=true, back=true, top=true, bot=true), - ) - flow_bcs!(stokes, flow_bcs) - flow_bcs!(stokes, flow_bcs) # just a trick to pass the CI + # free-slip + flow_bcs = VelocityBoundaryConditions(; + no_slip = (left=false, right=false, front=false, back=false, top=false, bot=false), + free_slip = (left=true, right=true, front=true, back=true, top=true, bot=true), + ) + flow_bcs!(stokes, flow_bcs) + flow_bcs!(stokes, flow_bcs) # just a trick to pass the CI - @test @views stokes.V.Vx[ :, :, 1] == stokes.V.Vx[:, :, 2] - @test @views stokes.V.Vx[ :, :, end] == stokes.V.Vx[:, :, end - 1] - @test @views stokes.V.Vx[ :, 1, :] == stokes.V.Vx[:, 2, :] - @test @views stokes.V.Vx[ :, end, :] == stokes.V.Vx[:, end - 1, :] - @test @views stokes.V.Vy[ :, :, 1] == stokes.V.Vy[:, :, 2] - @test @views stokes.V.Vy[ :, :, end] == stokes.V.Vy[:, :, end - 1] - @test @views stokes.V.Vy[ 1, :, :] == stokes.V.Vy[2, :, :] - @test @views stokes.V.Vy[end, :, :] == stokes.V.Vy[end - 1, :, :] - @test @views stokes.V.Vz[ 1, :, :] == stokes.V.Vz[2, :, :] - @test @views stokes.V.Vz[end, :, :] == stokes.V.Vz[end - 1, :, :] - @test @views stokes.V.Vz[ :, 1, :] == stokes.V.Vz[:, 2, :] - @test @views stokes.V.Vz[ :, end, :] == stokes.V.Vz[:, end - 1, :] + @test @views stokes.V.Vx[ :, :, 1] == stokes.V.Vx[:, :, 2] + @test @views stokes.V.Vx[ :, :, end] == stokes.V.Vx[:, :, end - 1] + @test @views stokes.V.Vx[ :, 1, :] == stokes.V.Vx[:, 2, :] + @test @views stokes.V.Vx[ :, end, :] == stokes.V.Vx[:, end - 1, :] + @test @views stokes.V.Vy[ :, :, 1] == stokes.V.Vy[:, :, 2] + @test @views stokes.V.Vy[ :, :, end] == stokes.V.Vy[:, :, end - 1] + @test @views stokes.V.Vy[ 1, :, :] == stokes.V.Vy[2, :, :] + @test @views stokes.V.Vy[end, :, :] == stokes.V.Vy[end - 1, :, :] + @test @views stokes.V.Vz[ 1, :, :] == stokes.V.Vz[2, :, :] + @test @views stokes.V.Vz[end, :, :] == stokes.V.Vz[end - 1, :, :] + @test @views stokes.V.Vz[ :, 1, :] == stokes.V.Vz[:, 2, :] + @test @views stokes.V.Vz[ :, end, :] == stokes.V.Vz[:, end - 1, :] - # no-slip - flow_bcs = VelocityBoundaryConditions(; - no_slip = (left=true, right=true, front=true, back=true, top=true, bot=true), - free_slip = (left=false, right=false, front=false, back=false, top=false, bot=false), - ) - flow_bcs!(stokes, flow_bcs) + # no-slip + flow_bcs = VelocityBoundaryConditions(; + no_slip = (left=true, right=true, front=true, back=true, top=true, bot=true), + free_slip = (left=false, right=false, front=false, back=false, top=false, bot=false), + ) + flow_bcs!(stokes, flow_bcs) - (; Vx, Vy, Vz) = stokes.V - @test sum(!iszero(Vx[1 , i, j]) for i in axes(Vx,2), j in axes(Vx,3)) == 0 - @test sum(!iszero(Vx[end, i, j]) for i in axes(Vx,2), j in axes(Vx,3)) == 0 - @test sum(!iszero(Vy[i, 1, j]) for i in axes(Vy,1), j in axes(Vy,2)) == 0 - @test sum(!iszero(Vy[i, end, j]) for i in axes(Vy,1), j in axes(Vy,2)) == 0 - @test sum(!iszero(Vz[i, j, 1]) for i in axes(Vz,1), j in axes(Vz,3)) == 0 - @test sum(!iszero(Vz[i, j, end]) for i in axes(Vz,1), j in axes(Vz,3)) == 0 - @test @views Vx[ :, 1, :] == -Vx[ :, 2, :] - @test @views Vx[ :, end, :] == -Vx[ :, end - 1, :] - @test @views Vx[ :, :, 1] == -Vx[ :, :, 2] - @test @views Vx[ :, :, end] == -Vx[ :, :, end - 1] - @test @views Vy[ 1, :, :] == -Vy[ 2, :, :] - @test @views Vy[end, :, :] == -Vy[end - 1, :, :] - @test @views Vy[ :, :, 1] == -Vy[ :, :, 2] - @test @views Vy[ :, :, end] == -Vy[ :, :, end - 1] - @test @views Vz[ :, 1, :] == -Vz[ :, 2, :] - @test @views Vz[ :, end, :] == -Vz[ :, end - 1, :] - @test @views Vz[ 1, :, :] == -Vz[ 2, :, :] - @test @views Vz[end, :, :] == -Vz[end - 1, :, :] + (; Vx, Vy, Vz) = stokes.V + @test sum(!iszero(Vx[1 , i, j]) for i in axes(Vx,2), j in axes(Vx,3)) == 0 + @test sum(!iszero(Vx[end, i, j]) for i in axes(Vx,2), j in axes(Vx,3)) == 0 + @test sum(!iszero(Vy[i, 1, j]) for i in axes(Vy,1), j in axes(Vy,2)) == 0 + @test sum(!iszero(Vy[i, end, j]) for i in axes(Vy,1), j in axes(Vy,2)) == 0 + @test sum(!iszero(Vz[i, j, 1]) for i in axes(Vz,1), j in axes(Vz,3)) == 0 + @test sum(!iszero(Vz[i, j, end]) for i in axes(Vz,1), j in axes(Vz,3)) == 0 + @test @views Vx[ :, 1, :] == -Vx[ :, 2, :] + @test @views Vx[ :, end, :] == -Vx[ :, end - 1, :] + @test @views Vx[ :, :, 1] == -Vx[ :, :, 2] + @test @views Vx[ :, :, end] == -Vx[ :, :, end - 1] + @test @views Vy[ 1, :, :] == -Vy[ 2, :, :] + @test @views Vy[end, :, :] == -Vy[end - 1, :, :] + @test @views Vy[ :, :, 1] == -Vy[ :, :, 2] + @test @views Vy[ :, :, end] == -Vy[ :, :, end - 1] + @test @views Vz[ :, 1, :] == -Vz[ :, 2, :] + @test @views Vz[ :, end, :] == -Vz[ :, end - 1, :] + @test @views Vz[ 1, :, :] == -Vz[ 2, :, :] + @test @views Vz[end, :, :] == -Vz[end - 1, :, :] + else + @test true === true + end end @testset "DisplacementBoundaryConditions" begin + if backend === CPUBackend + # test incompatible boundary conditions + @test_throws ErrorException DisplacementBoundaryConditions(; + no_slip = (left=true, right=true, front=true, back=true, top=true, bot=true), + free_slip = (left=false, right=true, front=true, back=true, top=true, bot=true), + ) - # test incompatible boundary conditions - @test_throws ErrorException DisplacementBoundaryConditions(; - no_slip = (left=true, right=true, front=true, back=true, top=true, bot=true), - free_slip = (left=false, right=true, front=true, back=true, top=true, bot=true), - ) + # test with StokesArrays + ni = 5, 5, 5 + stokes = StokesArrays(backend, ni) + stokes.U.Ux .= PTArray(backend)(rand(size(stokes.U.Ux)...)) + stokes.U.Uy .= PTArray(backend)(rand(size(stokes.U.Uy)...)) + stokes.U.Uz .= PTArray(backend)(rand(size(stokes.U.Uz)...)) - # test with StokesArrays - ni = 5, 5, 5 - stokes = StokesArrays(backend, ni) - stokes.U.Ux .= PTArray(backend)(rand(size(stokes.U.Ux)...)) - stokes.U.Uy .= PTArray(backend)(rand(size(stokes.U.Uy)...)) - stokes.U.Uz .= PTArray(backend)(rand(size(stokes.U.Uz)...)) + # free-slip + flow_bcs = DisplacementBoundaryConditions(; + no_slip = (left=false, right=false, front=false, back=false, top=false, bot=false), + free_slip = (left=true, right=true, front=true, back=true, top=true, bot=true), + ) + flow_bcs!(stokes, flow_bcs) + flow_bcs!(stokes, flow_bcs) # just a trick to pass the CI - # free-slip - flow_bcs = DisplacementBoundaryConditions(; - no_slip = (left=false, right=false, front=false, back=false, top=false, bot=false), - free_slip = (left=true, right=true, front=true, back=true, top=true, bot=true), - ) - flow_bcs!(stokes, flow_bcs) - flow_bcs!(stokes, flow_bcs) # just a trick to pass the CI + @test @views stokes.U.Ux[ :, :, 1] == stokes.U.Ux[:, :, 2] + @test @views stokes.U.Ux[ :, :, end] == stokes.U.Ux[:, :, end - 1] + @test @views stokes.U.Ux[ :, 1, :] == stokes.U.Ux[:, 2, :] + @test @views stokes.U.Ux[ :, end, :] == stokes.U.Ux[:, end - 1, :] + @test @views stokes.U.Uy[ :, :, 1] == stokes.U.Uy[:, :, 2] + @test @views stokes.U.Uy[ :, :, end] == stokes.U.Uy[:, :, end - 1] + @test @views stokes.U.Uy[ 1, :, :] == stokes.U.Uy[2, :, :] + @test @views stokes.U.Uy[end, :, :] == stokes.U.Uy[end - 1, :, :] + @test @views stokes.U.Uz[ 1, :, :] == stokes.U.Uz[2, :, :] + @test @views stokes.U.Uz[end, :, :] == stokes.U.Uz[end - 1, :, :] + @test @views stokes.U.Uz[ :, 1, :] == stokes.U.Uz[:, 2, :] + @test @views stokes.U.Uz[ :, end, :] == stokes.U.Uz[:, end - 1, :] - @test @views stokes.U.Ux[ :, :, 1] == stokes.U.Ux[:, :, 2] - @test @views stokes.U.Ux[ :, :, end] == stokes.U.Ux[:, :, end - 1] - @test @views stokes.U.Ux[ :, 1, :] == stokes.U.Ux[:, 2, :] - @test @views stokes.U.Ux[ :, end, :] == stokes.U.Ux[:, end - 1, :] - @test @views stokes.U.Uy[ :, :, 1] == stokes.U.Uy[:, :, 2] - @test @views stokes.U.Uy[ :, :, end] == stokes.U.Uy[:, :, end - 1] - @test @views stokes.U.Uy[ 1, :, :] == stokes.U.Uy[2, :, :] - @test @views stokes.U.Uy[end, :, :] == stokes.U.Uy[end - 1, :, :] - @test @views stokes.U.Uz[ 1, :, :] == stokes.U.Uz[2, :, :] - @test @views stokes.U.Uz[end, :, :] == stokes.U.Uz[end - 1, :, :] - @test @views stokes.U.Uz[ :, 1, :] == stokes.U.Uz[:, 2, :] - @test @views stokes.U.Uz[ :, end, :] == stokes.U.Uz[:, end - 1, :] + # no-slip + flow_bcs = DisplacementBoundaryConditions(; + no_slip = (left=true, right=true, front=true, back=true, top=true, bot=true), + free_slip = (left=false, right=false, front=false, back=false, top=false, bot=false), + ) + flow_bcs!(stokes, flow_bcs) - # no-slip - flow_bcs = DisplacementBoundaryConditions(; - no_slip = (left=true, right=true, front=true, back=true, top=true, bot=true), - free_slip = (left=false, right=false, front=false, back=false, top=false, bot=false), - ) - flow_bcs!(stokes, flow_bcs) - - (; Ux, Uy, Uz) = stokes.U - @test sum(!iszero(Ux[1 , i, j]) for i in axes(Ux,2), j in axes(Ux,3)) == 0 - @test sum(!iszero(Ux[end, i, j]) for i in axes(Ux,2), j in axes(Ux,3)) == 0 - @test sum(!iszero(Uy[i, 1, j]) for i in axes(Uy,1), j in axes(Uy,2)) == 0 - @test sum(!iszero(Uy[i, end, j]) for i in axes(Uy,1), j in axes(Uy,2)) == 0 - @test sum(!iszero(Uz[i, j, 1]) for i in axes(Uz,1), j in axes(Uz,3)) == 0 - @test sum(!iszero(Uz[i, j, end]) for i in axes(Uz,1), j in axes(Uz,3)) == 0 - @test @views Ux[ :, 1, :] == -Ux[ :, 2, :] - @test @views Ux[ :, end, :] == -Ux[ :, end - 1, :] - @test @views Ux[ :, :, 1] == -Ux[ :, :, 2] - @test @views Ux[ :, :, end] == -Ux[ :, :, end - 1] - @test @views Uy[ 1, :, :] == -Uy[ 2, :, :] - @test @views Uy[end, :, :] == -Uy[end - 1, :, :] - @test @views Uy[ :, :, 1] == -Uy[ :, :, 2] - @test @views Uy[ :, :, end] == -Uy[ :, :, end - 1] - @test @views Uz[ :, 1, :] == -Uz[ :, 2, :] - @test @views Uz[ :, end, :] == -Uz[ :, end - 1, :] - @test @views Uz[ 1, :, :] == -Uz[ 2, :, :] - @test @views Uz[end, :, :] == -Uz[end - 1, :, :] + (; Ux, Uy, Uz) = stokes.U + @test sum(!iszero(Ux[1 , i, j]) for i in axes(Ux,2), j in axes(Ux,3)) == 0 + @test sum(!iszero(Ux[end, i, j]) for i in axes(Ux,2), j in axes(Ux,3)) == 0 + @test sum(!iszero(Uy[i, 1, j]) for i in axes(Uy,1), j in axes(Uy,2)) == 0 + @test sum(!iszero(Uy[i, end, j]) for i in axes(Uy,1), j in axes(Uy,2)) == 0 + @test sum(!iszero(Uz[i, j, 1]) for i in axes(Uz,1), j in axes(Uz,3)) == 0 + @test sum(!iszero(Uz[i, j, end]) for i in axes(Uz,1), j in axes(Uz,3)) == 0 + @test @views Ux[ :, 1, :] == -Ux[ :, 2, :] + @test @views Ux[ :, end, :] == -Ux[ :, end - 1, :] + @test @views Ux[ :, :, 1] == -Ux[ :, :, 2] + @test @views Ux[ :, :, end] == -Ux[ :, :, end - 1] + @test @views Uy[ 1, :, :] == -Uy[ 2, :, :] + @test @views Uy[end, :, :] == -Uy[end - 1, :, :] + @test @views Uy[ :, :, 1] == -Uy[ :, :, 2] + @test @views Uy[ :, :, end] == -Uy[ :, :, end - 1] + @test @views Uz[ :, 1, :] == -Uz[ :, 2, :] + @test @views Uz[ :, end, :] == -Uz[ :, end - 1, :] + @test @views Uz[ 1, :, :] == -Uz[ 2, :, :] + @test @views Uz[end, :, :] == -Uz[end - 1, :, :] + else + @test true === true + end end end From 2114c2a741ad8fad641afae57dcbe09c826745c6 Mon Sep 17 00:00:00 2001 From: Pascal Aellig Date: Mon, 8 Jul 2024 11:52:38 +0200 Subject: [PATCH 18/20] typo --- src/types/displacement.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/types/displacement.jl b/src/types/displacement.jl index 85decb6b..03ce86e4 100644 --- a/src/types/displacement.jl +++ b/src/types/displacement.jl @@ -59,6 +59,6 @@ end return nothing end -displacement2velocity!(stokes, dt, <:DisplacementBoundaryConditions) = displacement2velocity!(backend(stokes), stokes, dt) -displacement2velocity!(::Any, ::Any, <:VelocityBoundaryConditions) = nothing +displacement2velocity!(stokes, dt, ::DisplacementBoundaryConditions) = displacement2velocity!(backend(stokes), stokes, dt) +displacement2velocity!(::Any, ::Any, ::VelocityBoundaryConditions) = nothing displacement2velocity!(::Any, ::Any, ::T) where T = throw(ArgumentError("Unknown boundary conditions type: $T")) From 62c3e3bf81353af09479c2719dcfa46039ad22e3 Mon Sep 17 00:00:00 2001 From: Pascal Aellig Date: Mon, 8 Jul 2024 12:00:57 +0200 Subject: [PATCH 19/20] format --- src/types/displacement.jl | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/src/types/displacement.jl b/src/types/displacement.jl index 03ce86e4..d50c9008 100644 --- a/src/types/displacement.jl +++ b/src/types/displacement.jl @@ -59,6 +59,10 @@ end return nothing end -displacement2velocity!(stokes, dt, ::DisplacementBoundaryConditions) = displacement2velocity!(backend(stokes), stokes, dt) +function displacement2velocity!(stokes, dt, ::DisplacementBoundaryConditions) + return displacement2velocity!(backend(stokes), stokes, dt) +end displacement2velocity!(::Any, ::Any, ::VelocityBoundaryConditions) = nothing -displacement2velocity!(::Any, ::Any, ::T) where T = throw(ArgumentError("Unknown boundary conditions type: $T")) +function displacement2velocity!(::Any, ::Any, ::T) where {T} + throw(ArgumentError("Unknown boundary conditions type: $T")) +end From 4786b6a0ff989b85c54a8668c5c63aafdc3344ba Mon Sep 17 00:00:00 2001 From: Pascal Aellig Date: Mon, 8 Jul 2024 12:03:09 +0200 Subject: [PATCH 20/20] format again --- src/JustRelax.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/JustRelax.jl b/src/JustRelax.jl index 71d40d03..7d2ae38d 100644 --- a/src/JustRelax.jl +++ b/src/JustRelax.jl @@ -33,8 +33,7 @@ include("types/phases.jl") include("boundaryconditions/types.jl") export TemperatureBoundaryConditions, - DisplacementBoundaryConditions, - VelocityBoundaryConditions + DisplacementBoundaryConditions, VelocityBoundaryConditions include("types/traits.jl") export BackendTrait, CPUBackendTrait, NonCPUBackendTrait