diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index bb45dff3..9a1a0ab9 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -6,6 +6,12 @@ on: pull_request: branches: - main + +# Cancel redundant CI tests automatically +concurrency: + group: ${{ github.workflow }}-${{ github.ref }} + cancel-in-progress: true + jobs: test: name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }} diff --git a/miniapps/benchmarks/stokes2D/Blankenbach2D/Benchmark2D_WENO5.jl b/miniapps/benchmarks/stokes2D/Blankenbach2D/Benchmark2D_WENO5.jl index ee7dade3..ae2a9ff0 100755 --- a/miniapps/benchmarks/stokes2D/Blankenbach2D/Benchmark2D_WENO5.jl +++ b/miniapps/benchmarks/stokes2D/Blankenbach2D/Benchmark2D_WENO5.jl @@ -73,7 +73,7 @@ function main2D(igg; ar=1, nx=32, ny=32, nit = 1e1, figdir="figs2D", do_vtk =fal # ---------------------------------------------------- # Weno model ----------------------------------------- - weno = WENO5(Val(2), (nx,ny).+1) # ni.+1 for Temp + weno = WENO5(backend, Val(2), (nx,ny).+1) # ni.+1 for Temp # ---------------------------------------------------- # Initialize particles ------------------------------- diff --git a/miniapps/convection/GlobalConvection2D_WENO5.jl b/miniapps/convection/GlobalConvection2D_WENO5.jl index 76c136c5..d8a60d65 100644 --- a/miniapps/convection/GlobalConvection2D_WENO5.jl +++ b/miniapps/convection/GlobalConvection2D_WENO5.jl @@ -94,7 +94,7 @@ function thermal_convection2D(igg; ar=8, ny=16, nx=ny*8, figdir="figs2D", therma # ---------------------------------------------------- # Weno model ----------------------------------------- - weno = WENO5(Val(2), ni.+1) # ni.+1 for Temp + weno = WENO5(backend, Val(2), ni.+1) # ni.+1 for Temp # ---------------------------------------------------- # create rheology struct diff --git a/miniapps/convection/WENO5/WENO_convection2D.jl b/miniapps/convection/WENO5/WENO_convection2D.jl index 283b4eb5..6e57a614 100644 --- a/miniapps/convection/WENO5/WENO_convection2D.jl +++ b/miniapps/convection/WENO5/WENO_convection2D.jl @@ -101,7 +101,7 @@ function main2D(igg; ar=8, ny=16, nx=ny*8, figdir="figs2D", do_vtk =false) # Weno model ----------------------------------------- weno = WENO5(backend_JR, Val(2), ni.+1) # ---------------------------------------------------- - + # Initialize particles ------------------------------- nxcell, max_xcell, min_xcell = 20, 40, 1 particles = init_particles(backend, nxcell, max_xcell, min_xcell, xvi...); @@ -285,8 +285,8 @@ end figdir = "Weno2D" ar = 1 # aspect ratio n = 64 -nx = n*ar -ny = n +nx = n*ar +ny = n igg = if !(JustRelax.MPI.Initialized()) # initialize (or not) MPI grid IGG(init_global_grid(nx, ny, 1; init_MPI= true)...) else diff --git a/src/ext/AMDGPU/2D.jl b/src/ext/AMDGPU/2D.jl index 8a99e040..0b5b9804 100644 --- a/src/ext/AMDGPU/2D.jl +++ b/src/ext/AMDGPU/2D.jl @@ -364,4 +364,8 @@ function JR2D.compute_shear_heating!( return nothing end +function JR2D.WENO_advection!(u::ROCArrayArray, Vxi::NTuple, weno, di, dt) + return WENO_advection!(u, Vxi, weno, di, dt) +end + end diff --git a/src/ext/AMDGPU/3D.jl b/src/ext/AMDGPU/3D.jl index 00b53a2a..a40c754c 100644 --- a/src/ext/AMDGPU/3D.jl +++ b/src/ext/AMDGPU/3D.jl @@ -236,10 +236,6 @@ function JR3D.compute_melt_fraction!( return compute_melt_fraction!(ϕ, phase_ratios, rheology, args) end -function JR3D.compute_melt_fraction!(ϕ::ROCArray, phase_ratios, rheology, args) - return compute_melt_fraction!(ϕ, phase_ratios, rheology, args) -end - # Interpolations function JR3D.temperature2center!(::AMDGPUBackendTrait, thermal::JustRelax.ThermalArrays) return _temperature2center!(thermal) @@ -385,4 +381,9 @@ function JR3D.compute_shear_heating!( ) return nothing end + +function JR3D.WENO_advection!(u::ROCArray, Vxi::NTuple, weno, di, dt) + return WENO_advection!(u, Vxi, weno, di, dt) +end + end diff --git a/src/ext/CUDA/2D.jl b/src/ext/CUDA/2D.jl index cd6a1602..5f448607 100644 --- a/src/ext/CUDA/2D.jl +++ b/src/ext/CUDA/2D.jl @@ -348,4 +348,9 @@ function JR2D.compute_shear_heating!( return nothing end +# WENO-5 advection +function JR2D.WENO_advection!(u::CuArray, Vxi::NTuple, weno, di, dt) + return WENO_advection!(u, Vxi, weno, di, dt) +end + end diff --git a/src/ext/CUDA/3D.jl b/src/ext/CUDA/3D.jl index aa5737db..da702989 100644 --- a/src/ext/CUDA/3D.jl +++ b/src/ext/CUDA/3D.jl @@ -180,20 +180,6 @@ end function thermal_bcs!(::CUDABackendTrait, thermal::JustRelax.ThermalArrays, bcs) return thermal_bcs!(thermal.T, bcs) end - -# # Phases -# function JR3D.phase_ratios_center!( -# ::CUDABackendTrait, phase_ratios::JustPIC.PhaseRatios, particles, grid::Geometry, phases -# ) -# return _phase_ratios_center!(phase_ratios, particles, grid, phases) -# end - -# function JR3D.phase_ratios_vertex!( -# ::CUDABackendTrait, phase_ratios::JustPIC.PhaseRatios, particles, grid::Geometry, phases -# ) -# return _phase_ratios_vertex!(phase_ratios, particles, grid, phases) -# end - # Rheology ## viscosity @@ -252,10 +238,6 @@ function JR3D.compute_melt_fraction!( return compute_melt_fraction!(ϕ, phase_ratios, rheology, args) end -function JR3D.compute_melt_fraction!(ϕ::CuArray, phase_ratios, rheology, args) - return compute_melt_fraction!(ϕ, phase_ratios, rheology, args) -end - # Interpolations function JR3D.temperature2center!(::CUDABackendTrait, thermal::JustRelax.ThermalArrays) return _temperature2center!(thermal) @@ -399,4 +381,8 @@ function JR3D.compute_shear_heating!( return nothing end +function JR3D.WENO_advection!(u::CuArray, Vxi::NTuple, weno, di, dt) + return WENO_advection!(u, Vxi, weno, di, dt) +end + end diff --git a/src/rheology/Melting.jl b/src/rheology/Melting.jl index 7a340641..04ed5c82 100644 --- a/src/rheology/Melting.jl +++ b/src/rheology/Melting.jl @@ -1,31 +1,26 @@ function compute_melt_fraction!(ϕ, rheology, args) ni = size(ϕ) @parallel (@idx ni) compute_melt_fraction_kernel!(ϕ, rheology, args) + return nothing end @parallel_indices (I...) function compute_melt_fraction_kernel!(ϕ, rheology, args) args_ijk = ntuple_idx(args, I...) - ϕ[I...] = compute_melt_frac(rheology, args_ijk) + ϕ[I...] = compute_meltfraction(rheology, args_ijk) return nothing end -@inline function compute_melt_frac(rheology, args) - return compute_meltfraction(rheology, args) -end - -function compute_melt_fraction!(ϕ, phase_ratios, rheology, args) +function compute_melt_fraction!(ϕ, phase_ratios::JustPIC.PhaseRatios, rheology, args) ni = size(ϕ) - @parallel (@idx ni) compute_melt_fraction_kernel!(ϕ, phase_ratios, rheology, args) + @parallel (@idx ni) compute_melt_fraction_kernel!( + ϕ, phase_ratios.center, rheology, args + ) end @parallel_indices (I...) function compute_melt_fraction_kernel!( ϕ, phase_ratios, rheology, args ) args_ijk = ntuple_idx(args, I...) - ϕ[I...] = compute_melt_frac(rheology, args_ijk, phase_ratios[I...]) + ϕ[I...] = fn_ratio(compute_meltfraction, rheology, @cell(phase_ratios[I...]), args_ijk) return nothing end - -@inline function compute_melt_frac(rheology, args, phase_ratios) - return compute_meltfraction_ratio(phase_ratios, rheology, args) -end diff --git a/src/types/weno.jl b/src/types/weno.jl index f980d47d..1c31bd89 100644 --- a/src/types/weno.jl +++ b/src/types/weno.jl @@ -1,3 +1,4 @@ +using Adapt ## Weno5 advection scheme. Implementation based on the repository from # https://gmd.copernicus.org/preprints/gmd-2023-189/ @@ -59,3 +60,5 @@ struct WENO5{T,N,A,M} <: AbstractWENO # method method::M # 1:JS, 2:Z end + +Adapt.@adapt_structure WENO5 diff --git a/test/test_Blankenbach.jl b/test/test_Blankenbach.jl index fd5ded93..bcea11d4 100644 --- a/test/test_Blankenbach.jl +++ b/test/test_Blankenbach.jl @@ -2,10 +2,8 @@ push!(LOAD_PATH, "..") @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 Test, Suppressor @@ -137,12 +135,6 @@ function main2D(igg; ar=1, nx=32, ny=32, nit = 10) temperature2center!(thermal) # ---------------------------------------------------- - # Rayleigh number - ΔT = thermal.T[1,1] - thermal.T[1,end] - Ra = (rheology[1].Density[1].ρ0 * rheology[1].Gravity[1].g * rheology[1].Density[1].α * ΔT * ly^3.0 ) / - (κ * rheology[1].CompositeRheology[1].elements[1].η ) - @show Ra - args = (; T = thermal.Tc, P = stokes.P, dt = Inf) # Buoyancy forces & viscosity ---------------------- @@ -264,7 +256,7 @@ function main2D(igg; ar=1, nx=32, ny=32, nit = 10) # Compute U rms ----------------------------- # U₍ᵣₘₛ₎ = H*ρ₀*c₍ₚ₎/k * √ 1/H/L * ∫∫ (vx²+vz²) dx dz Urms_it = let - JustRelax.JustRelax2D.velocity2vertex!(Vx_v, Vy_v, stokes.V.Vx, stokes.V.Vy; ghost_nodes=true) + velocity2vertex!(Vx_v, Vy_v, stokes.V.Vx, stokes.V.Vy; ghost_nodes=true) @. Vx_v .= hypot.(Vx_v, Vy_v) # we reuse Vx_v to store the velocity magnitude sqrt( sum( Vx_v.^2 .* prod(di)) / lx /ly ) * ((ly * rheology[1].Density[1].ρ0 * rheology[1].HeatCapacity[1].Cp) / rheology[1].Conductivity[1].k ) diff --git a/test/test_Interpolations.jl b/test/test_Interpolations.jl index 08a2fca8..3bc0b317 100644 --- a/test/test_Interpolations.jl +++ b/test/test_Interpolations.jl @@ -2,10 +2,8 @@ push!(LOAD_PATH, "..") @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 Test @@ -27,7 +25,7 @@ end @testset "Interpolations" begin - + if backend_JR == CPUBackend # Set up mock data # Physical domain ------------------------------------ ly = 1.0 # domain length in y @@ -75,4 +73,7 @@ end velocity2vertex!(Vx_v, Vy_v, stokes.V.Vx, stokes.V.Vy; ghost_nodes=false) @test iszero(Vx_v[1,1]) @test Vy_v[1,1] == 10 + else + @test true == true + end end diff --git a/test/test_WENO5.jl b/test/test_WENO5.jl index ddd2b008..ccd7c2d4 100644 --- a/test/test_WENO5.jl +++ b/test/test_WENO5.jl @@ -2,10 +2,8 @@ push!(LOAD_PATH, "..") @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 Test, Suppressor @@ -112,7 +110,7 @@ function thermal_convection2D(igg; ar=8, ny=16, nx=ny*8, thermal_perturbation = # ---------------------------------------------------- # Weno model ----------------------------------------- - weno = WENO5(Val(2), ni.+1) # ni.+1 for Temp + weno = WENO5(backend_JR, Val(2), ni.+1) # ni.+1 for Temp # ---------------------------------------------------- # create rheology struct diff --git a/test/test_boundary_conditions2D.jl b/test/test_boundary_conditions2D.jl index dbb0d6b6..285c8e60 100644 --- a/test/test_boundary_conditions2D.jl +++ b/test/test_boundary_conditions2D.jl @@ -1,9 +1,7 @@ @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 diff --git a/test/test_boundary_conditions3D.jl b/test/test_boundary_conditions3D.jl index 88aa96d6..2630dfc9 100644 --- a/test/test_boundary_conditions3D.jl +++ b/test/test_boundary_conditions3D.jl @@ -1,9 +1,7 @@ @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 diff --git a/test/test_shearheating2D.jl b/test/test_shearheating2D.jl index b3335052..e71f5578 100644 --- a/test/test_shearheating2D.jl +++ b/test/test_shearheating2D.jl @@ -3,10 +3,8 @@ using Test, Suppressor @static if ENV["JULIA_JUSTRELAX_BACKEND"] === "AMDGPU" using AMDGPU - AMDGPU.allowscalar(true) elseif ENV["JULIA_JUSTRELAX_BACKEND"] === "CUDA" using CUDA - CUDA.allowscalar(true) end # Benchmark of Duretz et al. 2014 diff --git a/test/test_shearheating3D.jl b/test/test_shearheating3D.jl index e79cf381..68b29d58 100644 --- a/test/test_shearheating3D.jl +++ b/test/test_shearheating3D.jl @@ -3,10 +3,8 @@ using Test, Suppressor @static if ENV["JULIA_JUSTRELAX_BACKEND"] === "AMDGPU" using AMDGPU - AMDGPU.allowscalar(true) elseif ENV["JULIA_JUSTRELAX_BACKEND"] === "CUDA" using CUDA - CUDA.allowscalar(true) end # Benchmark of Duretz et al. 2014 diff --git a/test/test_sinking_block.jl b/test/test_sinking_block.jl index 3f9f9be6..50a823fa 100644 --- a/test/test_sinking_block.jl +++ b/test/test_sinking_block.jl @@ -2,10 +2,8 @@ push!(LOAD_PATH, "..") @static if ENV["JULIA_JUSTRELAX_BACKEND"] === "AMDGPU" using AMDGPU - AMDGPU.allowscalar(true) #used because of velocity2vertex! elseif ENV["JULIA_JUSTRELAX_BACKEND"] === "CUDA" using CUDA - CUDA.allowscalar(true) #used because of velocity2vertex! end diff --git a/test/test_stokes_solvi3D.jl b/test/test_stokes_solvi3D.jl index 5e24415b..ea227fa8 100644 --- a/test/test_stokes_solvi3D.jl +++ b/test/test_stokes_solvi3D.jl @@ -2,10 +2,8 @@ push!(LOAD_PATH, "..") @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 Test, Suppressor diff --git a/test/test_thermalstresses.jl b/test/test_thermalstresses.jl index ca538c7e..9dc806e4 100644 --- a/test/test_thermalstresses.jl +++ b/test/test_thermalstresses.jl @@ -2,10 +2,8 @@ push!(LOAD_PATH, "..") @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 Test, Suppressor @@ -222,10 +220,7 @@ function init_rheology(CharDim; is_compressible = false, steady_state=true) end -function main2D(; nx=32, ny=32) - - init_mpi = JustRelax.MPI.Initialized() ? false : true - igg = IGG(init_global_grid(nx, ny, 1; init_MPI = init_mpi)...) +function main2D(igg; nx=32, ny=32) # Characteristic lengths CharDim = GEO_units(;length=12.5km, viscosity=1e21, temperature = 1e3C) @@ -317,7 +312,7 @@ function main2D(; nx=32, ny=32) ϕ = @zeros(ni...) compute_melt_fraction!( - ϕ, phase_ratios.center, rheology, (T=thermal.Tc, P=stokes.P) + ϕ, phase_ratios, rheology, (T=thermal.Tc, P=stokes.P) ) # Buoyancy force ρg = @zeros(ni...), @zeros(ni...) # ρg[1] is the buoyancy force in the x direction, ρg[2] is the buoyancy force in the y direction @@ -343,8 +338,8 @@ function main2D(; nx=32, ny=32) @copy stokes.P0 stokes.P thermal.Told .= thermal.T P_init = deepcopy(stokes.P) - Tsurf = thermal.T[1, end] - Tbot = thermal.T[1, 1] + Tsurf = nondimensionalize(273K,CharDim) + Tbot = nondimensionalize(648K,CharDim) local ϕ, stokes, thermal while it < 1 @@ -416,7 +411,7 @@ function main2D(; nx=32, ny=32) ) # ------------------------------ compute_melt_fraction!( - ϕ, phase_ratios.center, rheology, (T=thermal.Tc, P=stokes.P) + ϕ, phase_ratios, rheology, (T=thermal.Tc, P=stokes.P) ) # Advection -------------------- # advect particles in space @@ -448,9 +443,17 @@ end @testset "thermal stresses" begin @suppress begin - ϕ, stokes, thermal = main2D(; nx=32, ny=32) - nx_T, ny_T = size(thermal.T) .>>> 1 - @test Array(thermal.T)[nx_T + 1, ny_T + 1] ≈ 0.5369 rtol = 1e-2 - @test Array(ϕ)[nx_T + 1, ny_T + 1] ≤ 1e-8 + nx, ny = 32, 32 # number of cells + igg = if !(JustRelax.MPI.Initialized()) # initialize (or not) MPI grid + IGG(init_global_grid(nx, ny, 1; init_MPI= true)...) + else + igg + end + + ϕ, stokes, thermal = main2D(igg; nx=32, ny=32) + + nx_T, ny_T = size(thermal.T) + @test Array(thermal.T)[nx_T >>> 1 + 1, ny_T >>> 1 + 1] ≈ 0.5369 rtol = 1e-2 + @test Array(ϕ)[nx_T >>> 1 + 1, ny_T >>> 1 + 1] ≈ 9.351e-9 rtol = 1e-1 end end \ No newline at end of file