From dffbaa58c24158b497899cd7294718ebd3b51e89 Mon Sep 17 00:00:00 2001 From: Pascal Aellig Date: Fri, 4 Oct 2024 11:29:44 +0200 Subject: [PATCH 01/16] test ext melt fraction --- src/ext/AMDGPU/2D.jl | 10 ++++++++++ src/ext/AMDGPU/3D.jl | 8 +++++++- src/ext/CUDA/2D.jl | 10 ++++++++++ src/ext/CUDA/3D.jl | 22 +++++++--------------- 4 files changed, 34 insertions(+), 16 deletions(-) diff --git a/src/ext/AMDGPU/2D.jl b/src/ext/AMDGPU/2D.jl index 8a99e040..82c62248 100644 --- a/src/ext/AMDGPU/2D.jl +++ b/src/ext/AMDGPU/2D.jl @@ -231,6 +231,16 @@ function JR2D.compute_melt_fraction!( return compute_melt_fraction!(ϕ, phase_ratios, rheology, args) end +function compute_melt_fraction!(ϕ::ROCArray, rheology, args) + return compute_melt_fraction!(ϕ, rheology, args) +end + +function compute_melt_fraction!( + ϕ::ROCArray, phase_ratios::JustPIC.PhaseRatios, rheology, args +) + return compute_melt_fraction!(ϕ, phase_ratios, rheology, args) +end + # Interpolations function JR2D.temperature2center!(::AMDGPUBackendTrait, thermal::JustRelax.ThermalArrays) return _temperature2center!(thermal) diff --git a/src/ext/AMDGPU/3D.jl b/src/ext/AMDGPU/3D.jl index 00b53a2a..651c6f6e 100644 --- a/src/ext/AMDGPU/3D.jl +++ b/src/ext/AMDGPU/3D.jl @@ -236,7 +236,13 @@ function JR3D.compute_melt_fraction!( return compute_melt_fraction!(ϕ, phase_ratios, rheology, args) end -function JR3D.compute_melt_fraction!(ϕ::ROCArray, phase_ratios, rheology, args) +function compute_melt_fraction!(ϕ::ROCArray, rheology, args) + return compute_melt_fraction!(ϕ, rheology, args) +end + +function compute_melt_fraction!( + ϕ::ROCArray, phase_ratios::JustPIC.PhaseRatios, rheology, args +) return compute_melt_fraction!(ϕ, phase_ratios, rheology, args) end diff --git a/src/ext/CUDA/2D.jl b/src/ext/CUDA/2D.jl index cd6a1602..aae5cf04 100644 --- a/src/ext/CUDA/2D.jl +++ b/src/ext/CUDA/2D.jl @@ -219,6 +219,16 @@ function JR2D.compute_melt_fraction!( return compute_melt_fraction!(ϕ, phase_ratios, rheology, args) end +function compute_melt_fraction!(ϕ::CuArray, rheology, args) + return compute_melt_fraction!(ϕ, rheology, args) +end + +function compute_melt_fraction!( + ϕ::CuArray, phase_ratios::JustPIC.PhaseRatios, rheology, args +) + return compute_melt_fraction!(ϕ, phase_ratios, rheology, args) +end + # Interpolations function JR2D.temperature2center!(::CUDABackendTrait, thermal::JustRelax.ThermalArrays) return _temperature2center!(thermal) diff --git a/src/ext/CUDA/3D.jl b/src/ext/CUDA/3D.jl index aa5737db..cff47c60 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,7 +238,13 @@ function JR3D.compute_melt_fraction!( return compute_melt_fraction!(ϕ, phase_ratios, rheology, args) end -function JR3D.compute_melt_fraction!(ϕ::CuArray, phase_ratios, rheology, args) +function compute_melt_fraction!(ϕ::CuArray, rheology, args) + return compute_melt_fraction!(ϕ, rheology, args) +end + +function compute_melt_fraction!( + ϕ::CuArray, phase_ratios::JustPIC.PhaseRatios, rheology, args +) return compute_melt_fraction!(ϕ, phase_ratios, rheology, args) end From 38d4f9a138f5ff8b8ef2edcb664c3ab2844c26df Mon Sep 17 00:00:00 2001 From: Pascal Aellig Date: Fri, 4 Oct 2024 11:45:46 +0200 Subject: [PATCH 02/16] indexing error fix? --- src/rheology/Melting.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/rheology/Melting.jl b/src/rheology/Melting.jl index 7a340641..87e47aa6 100644 --- a/src/rheology/Melting.jl +++ b/src/rheology/Melting.jl @@ -13,9 +13,9 @@ end 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!( From 2518c7a2e8912f11aff17f88d4eefcbb60486775 Mon Sep 17 00:00:00 2001 From: Pascal Aellig Date: Fri, 4 Oct 2024 12:00:26 +0200 Subject: [PATCH 03/16] revert ext --- src/ext/CUDA/2D.jl | 10 ---------- 1 file changed, 10 deletions(-) diff --git a/src/ext/CUDA/2D.jl b/src/ext/CUDA/2D.jl index aae5cf04..cd6a1602 100644 --- a/src/ext/CUDA/2D.jl +++ b/src/ext/CUDA/2D.jl @@ -219,16 +219,6 @@ function JR2D.compute_melt_fraction!( return compute_melt_fraction!(ϕ, phase_ratios, rheology, args) end -function compute_melt_fraction!(ϕ::CuArray, rheology, args) - return compute_melt_fraction!(ϕ, rheology, args) -end - -function compute_melt_fraction!( - ϕ::CuArray, phase_ratios::JustPIC.PhaseRatios, rheology, args -) - return compute_melt_fraction!(ϕ, phase_ratios, rheology, args) -end - # Interpolations function JR2D.temperature2center!(::CUDABackendTrait, thermal::JustRelax.ThermalArrays) return _temperature2center!(thermal) From f17e2eddae4ffb82c85e66ce24f69597c2753bd0 Mon Sep 17 00:00:00 2001 From: Pascal Aellig Date: Fri, 4 Oct 2024 12:03:33 +0200 Subject: [PATCH 04/16] completely revert ext changes --- src/ext/AMDGPU/2D.jl | 10 ---------- src/ext/AMDGPU/3D.jl | 10 ---------- src/ext/CUDA/3D.jl | 10 ---------- 3 files changed, 30 deletions(-) diff --git a/src/ext/AMDGPU/2D.jl b/src/ext/AMDGPU/2D.jl index 82c62248..8a99e040 100644 --- a/src/ext/AMDGPU/2D.jl +++ b/src/ext/AMDGPU/2D.jl @@ -231,16 +231,6 @@ function JR2D.compute_melt_fraction!( return compute_melt_fraction!(ϕ, phase_ratios, rheology, args) end -function compute_melt_fraction!(ϕ::ROCArray, rheology, args) - return compute_melt_fraction!(ϕ, rheology, args) -end - -function compute_melt_fraction!( - ϕ::ROCArray, phase_ratios::JustPIC.PhaseRatios, rheology, args -) - return compute_melt_fraction!(ϕ, phase_ratios, rheology, args) -end - # Interpolations function JR2D.temperature2center!(::AMDGPUBackendTrait, thermal::JustRelax.ThermalArrays) return _temperature2center!(thermal) diff --git a/src/ext/AMDGPU/3D.jl b/src/ext/AMDGPU/3D.jl index 651c6f6e..c241ecae 100644 --- a/src/ext/AMDGPU/3D.jl +++ b/src/ext/AMDGPU/3D.jl @@ -236,16 +236,6 @@ function JR3D.compute_melt_fraction!( return compute_melt_fraction!(ϕ, phase_ratios, rheology, args) end -function compute_melt_fraction!(ϕ::ROCArray, rheology, args) - return compute_melt_fraction!(ϕ, rheology, args) -end - -function compute_melt_fraction!( - ϕ::ROCArray, phase_ratios::JustPIC.PhaseRatios, rheology, args -) - return compute_melt_fraction!(ϕ, phase_ratios, rheology, args) -end - # Interpolations function JR3D.temperature2center!(::AMDGPUBackendTrait, thermal::JustRelax.ThermalArrays) return _temperature2center!(thermal) diff --git a/src/ext/CUDA/3D.jl b/src/ext/CUDA/3D.jl index cff47c60..943e79b7 100644 --- a/src/ext/CUDA/3D.jl +++ b/src/ext/CUDA/3D.jl @@ -238,16 +238,6 @@ function JR3D.compute_melt_fraction!( return compute_melt_fraction!(ϕ, phase_ratios, rheology, args) end -function compute_melt_fraction!(ϕ::CuArray, rheology, args) - return compute_melt_fraction!(ϕ, rheology, args) -end - -function compute_melt_fraction!( - ϕ::CuArray, phase_ratios::JustPIC.PhaseRatios, rheology, args -) - return compute_melt_fraction!(ϕ, phase_ratios, rheology, args) -end - # Interpolations function JR3D.temperature2center!(::CUDABackendTrait, thermal::JustRelax.ThermalArrays) return _temperature2center!(thermal) From 79109bb60d5b7f7693fbe265350909eca6e1fe82 Mon Sep 17 00:00:00 2001 From: aelligp Date: Fri, 4 Oct 2024 13:33:00 +0200 Subject: [PATCH 05/16] fix `weno_f`by no inbounds checking --- src/advection/weno5.jl | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/src/advection/weno5.jl b/src/advection/weno5.jl index fb2e9875..c892df8b 100644 --- a/src/advection/weno5.jl +++ b/src/advection/weno5.jl @@ -166,10 +166,12 @@ end end @parallel_indices (i, j) function weno_f!(u, weno, nx, ny) - weno.fL[i, j] = WENO_flux_upwind_x(u, nx, weno, i, j) - weno.fR[i, j] = WENO_flux_downwind_x(u, nx, weno, i, j) - weno.fB[i, j] = WENO_flux_upwind_y(u, ny, weno, i, j) - weno.fT[i, j] = WENO_flux_downwind_y(u, ny, weno, i, j) + @inbounds begin + weno.fL[i, j] = WENO_flux_upwind_x(u, nx, weno, i, j) + weno.fR[i, j] = WENO_flux_downwind_x(u, nx, weno, i, j) + weno.fB[i, j] = WENO_flux_upwind_y(u, ny, weno, i, j) + weno.fT[i, j] = WENO_flux_downwind_y(u, ny, weno, i, j) + end return nothing end From 57739b9f41b60633bc35c6e405fd8e7190d7f1c8 Mon Sep 17 00:00:00 2001 From: aelligp Date: Fri, 4 Oct 2024 15:14:37 +0200 Subject: [PATCH 06/16] udpate tests to disallow scalar, WENO5 w/ backend --- .../stokes2D/Blankenbach2D/Benchmark2D_WENO5.jl | 2 +- .../diffusion/diffusion2D_multiphase.jl | 8 +++++--- miniapps/convection/GlobalConvection2D_WENO5.jl | 2 +- miniapps/convection/WENO5/WENO_convection2D.jl | 2 +- test/test_Blankenbach.jl | 10 +--------- test/test_Interpolations.jl | 7 ++++--- test/test_WENO5.jl | 4 +--- test/test_boundary_conditions2D.jl | 2 -- test/test_boundary_conditions3D.jl | 2 -- test/test_shearheating2D.jl | 2 -- test/test_shearheating3D.jl | 2 -- test/test_sinking_block.jl | 2 -- test/test_stokes_solvi3D.jl | 2 -- test/test_thermalstresses.jl | 6 ++---- 14 files changed, 16 insertions(+), 37 deletions(-) diff --git a/miniapps/benchmarks/stokes2D/Blankenbach2D/Benchmark2D_WENO5.jl b/miniapps/benchmarks/stokes2D/Blankenbach2D/Benchmark2D_WENO5.jl index 3a8fa4ee..82c1aa74 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/benchmarks/thermal_diffusion/diffusion/diffusion2D_multiphase.jl b/miniapps/benchmarks/thermal_diffusion/diffusion/diffusion2D_multiphase.jl index 05afb3cc..75dcfcc7 100644 --- a/miniapps/benchmarks/thermal_diffusion/diffusion/diffusion2D_multiphase.jl +++ b/miniapps/benchmarks/thermal_diffusion/diffusion/diffusion2D_multiphase.jl @@ -6,6 +6,8 @@ using ParallelStencil using GeoParams +using JustPIC, JustPIC._2D +const backend = JustPIC.CPUBackend distance(p1, p2) = mapreduce(x->(x[1]-x[2])^2, +, zip(p1, p2)) |> sqrt @@ -61,7 +63,7 @@ end @parallel_indices (I...) function compute_temperature_source_terms!(H, rheology, phase_ratios, args) args_ij = ntuple_idx(args, I...) - H[I...] = fn_ratio(compute_radioactive_heat, rheology, phase_ratios[I...], args_ij) + H[I...] = fn_ratio(compute_radioactive_heat, rheology, phase_ratios.center[I...], args_ij) return nothing end @@ -73,7 +75,7 @@ function diffusion_2D(; nx=32, ny=32, lx=100e3, ly=100e3, Cp0=1.2e3, K0=3.0) dt = 50 * kyr # physical time step init_mpi = JustRelax.MPI.Initialized() ? false : true - igg = IGG(init_global_grid(nx, ny, 1; init_MPI = init_mpi)...) + # igg = IGG(init_global_grid(nx, ny, 1; init_MPI = init_mpi)...) # Physical domain ni = (nx, ny) @@ -129,7 +131,7 @@ function diffusion_2D(; nx=32, ny=32, lx=100e3, ly=100e3, Cp0=1.2e3, K0=3.0) update_phase_ratios!(phase_ratios, particles, xci, xvi, pPhases) # ---------------------------------------------------- - @parallel (@idx ni) compute_temperature_source_terms!(thermal.H, rheology, phase_ratios.center, args) + @parallel (@idx ni) compute_temperature_source_terms!(thermal.H, rheology, phase_ratios, args) # PT coefficients for thermal diffusion args = (; P=P, T=thermal.Tc) 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 3b804905..83d08ced 100644 --- a/miniapps/convection/WENO5/WENO_convection2D.jl +++ b/miniapps/convection/WENO5/WENO_convection2D.jl @@ -100,7 +100,7 @@ function main2D(igg; ar=8, ny=16, nx=ny*8, figdir="figs2D", do_vtk =false) # ---------------------------------------------------- # Weno model ----------------------------------------- - weno = WENO5(Val(2), (nx,ny).+1) # ni.+1 for Temp + weno = WENO5(backend_JR, Val(2), (nx,ny).+1) # ni.+1 for Temp # ---------------------------------------------------- # Initialize particles ------------------------------- 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 612720e9..8c68603e 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 ecbf1826..342fefde 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 814eaaba..6c5d55e0 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 @@ -317,7 +315,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 @@ -416,7 +414,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 From c0718406b7a2767d71e745d0510c7b8eaf29d6a2 Mon Sep 17 00:00:00 2001 From: aelligp Date: Fri, 4 Oct 2024 17:24:40 +0200 Subject: [PATCH 07/16] fix? --- src/rheology/Melting.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/rheology/Melting.jl b/src/rheology/Melting.jl index 87e47aa6..67662f62 100644 --- a/src/rheology/Melting.jl +++ b/src/rheology/Melting.jl @@ -22,7 +22,7 @@ end ϕ, phase_ratios, rheology, args ) args_ijk = ntuple_idx(args, I...) - ϕ[I...] = compute_melt_frac(rheology, args_ijk, phase_ratios[I...]) + ϕ[I...] = compute_melt_frac(rheology, args_ijk, @cell(phase_ratios[I...])) return nothing end From 05eadfb9281b4a654ff49ea1e640bac7e13a6452 Mon Sep 17 00:00:00 2001 From: Pascal Aellig Date: Tue, 8 Oct 2024 09:27:13 +0200 Subject: [PATCH 08/16] update miniapp --- test/test_thermalstresses.jl | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/test/test_thermalstresses.jl b/test/test_thermalstresses.jl index 6c5d55e0..dffde3ef 100644 --- a/test/test_thermalstresses.jl +++ b/test/test_thermalstresses.jl @@ -220,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) @@ -446,7 +443,14 @@ end @testset "thermal stresses" begin @suppress begin - ϕ, stokes, thermal = main2D(; nx=32, ny=32) + 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 From 197d1ef88582955ad7ad19c800e60f9ce4cabbd6 Mon Sep 17 00:00:00 2001 From: aelligp Date: Tue, 8 Oct 2024 17:41:25 +0200 Subject: [PATCH 09/16] adapt WENO5 for GPU --- Project.toml | 3 +-- src/ext/AMDGPU/2D.jl | 4 ++++ src/ext/AMDGPU/3D.jl | 5 +++++ src/ext/CUDA/2D.jl | 5 +++++ src/ext/CUDA/3D.jl | 4 ++++ src/types/weno.jl | 3 +++ 6 files changed, 22 insertions(+), 2 deletions(-) diff --git a/Project.toml b/Project.toml index 2c5920d2..91a18fff 100644 --- a/Project.toml +++ b/Project.toml @@ -6,9 +6,7 @@ version = "0.3.2" [deps] Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" CellArrays = "d35fcfd7-7af4-4c67-b1aa-d78070614af4" -GMT = "5752ebe1-31b9-557e-87aa-f909b540aa54" GeoParams = "e018b62d-d9de-4a26-8697-af89c310ae38" -GeophysicalModelGenerator = "3700c31b-fa53-48a6-808a-ef22d5a84742" HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f" ImplicitGlobalGrid = "4d7a3746-15be-11ea-1130-334b0c4f5fa0" Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59" @@ -20,6 +18,7 @@ MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221" ParallelStencil = "94395366-693c-11ea-3b26-d9b7aac5d958" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" +SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192" 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 c241ecae..a40c754c 100644 --- a/src/ext/AMDGPU/3D.jl +++ b/src/ext/AMDGPU/3D.jl @@ -381,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 943e79b7..da702989 100644 --- a/src/ext/CUDA/3D.jl +++ b/src/ext/CUDA/3D.jl @@ -381,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/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 From d252f857d53a4df8367fdfd22e9d8d0c35d7eed7 Mon Sep 17 00:00:00 2001 From: Pascal Aellig Date: Tue, 8 Oct 2024 18:57:27 +0200 Subject: [PATCH 10/16] address suggestions --- .../diffusion/diffusion2D_multiphase.jl | 6 +++--- src/advection/weno5.jl | 10 ++++------ 2 files changed, 7 insertions(+), 9 deletions(-) diff --git a/miniapps/benchmarks/thermal_diffusion/diffusion/diffusion2D_multiphase.jl b/miniapps/benchmarks/thermal_diffusion/diffusion/diffusion2D_multiphase.jl index 75dcfcc7..724f7fe8 100644 --- a/miniapps/benchmarks/thermal_diffusion/diffusion/diffusion2D_multiphase.jl +++ b/miniapps/benchmarks/thermal_diffusion/diffusion/diffusion2D_multiphase.jl @@ -63,7 +63,7 @@ end @parallel_indices (I...) function compute_temperature_source_terms!(H, rheology, phase_ratios, args) args_ij = ntuple_idx(args, I...) - H[I...] = fn_ratio(compute_radioactive_heat, rheology, phase_ratios.center[I...], args_ij) + H[I...] = fn_ratio(compute_radioactive_heat, rheology, phase_ratios[I...], args_ij) return nothing end @@ -75,7 +75,7 @@ function diffusion_2D(; nx=32, ny=32, lx=100e3, ly=100e3, Cp0=1.2e3, K0=3.0) dt = 50 * kyr # physical time step init_mpi = JustRelax.MPI.Initialized() ? false : true - # igg = IGG(init_global_grid(nx, ny, 1; init_MPI = init_mpi)...) + igg = IGG(init_global_grid(nx, ny, 1; init_MPI = init_mpi)...) # Physical domain ni = (nx, ny) @@ -131,7 +131,7 @@ function diffusion_2D(; nx=32, ny=32, lx=100e3, ly=100e3, Cp0=1.2e3, K0=3.0) update_phase_ratios!(phase_ratios, particles, xci, xvi, pPhases) # ---------------------------------------------------- - @parallel (@idx ni) compute_temperature_source_terms!(thermal.H, rheology, phase_ratios, args) + @parallel (@idx ni) compute_temperature_source_terms!(thermal.H, rheology, phase_ratios.center, args) # PT coefficients for thermal diffusion args = (; P=P, T=thermal.Tc) diff --git a/src/advection/weno5.jl b/src/advection/weno5.jl index c892df8b..fb2e9875 100644 --- a/src/advection/weno5.jl +++ b/src/advection/weno5.jl @@ -166,12 +166,10 @@ end end @parallel_indices (i, j) function weno_f!(u, weno, nx, ny) - @inbounds begin - weno.fL[i, j] = WENO_flux_upwind_x(u, nx, weno, i, j) - weno.fR[i, j] = WENO_flux_downwind_x(u, nx, weno, i, j) - weno.fB[i, j] = WENO_flux_upwind_y(u, ny, weno, i, j) - weno.fT[i, j] = WENO_flux_downwind_y(u, ny, weno, i, j) - end + weno.fL[i, j] = WENO_flux_upwind_x(u, nx, weno, i, j) + weno.fR[i, j] = WENO_flux_downwind_x(u, nx, weno, i, j) + weno.fB[i, j] = WENO_flux_upwind_y(u, ny, weno, i, j) + weno.fT[i, j] = WENO_flux_downwind_y(u, ny, weno, i, j) return nothing end From 921f05b277c5ce555eaad4b3aa0747e5eab8ec4e Mon Sep 17 00:00:00 2001 From: Pascal Aellig <93140290+aelligp@users.noreply.github.com> Date: Wed, 9 Oct 2024 08:52:57 +0200 Subject: [PATCH 11/16] Update Project.toml Co-authored-by: Albert de Montserrat <58044444+albert-de-montserrat@users.noreply.github.com> --- Project.toml | 1 - 1 file changed, 1 deletion(-) diff --git a/Project.toml b/Project.toml index 91a18fff..076da6b6 100644 --- a/Project.toml +++ b/Project.toml @@ -18,7 +18,6 @@ MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221" ParallelStencil = "94395366-693c-11ea-3b26-d9b7aac5d958" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" -SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192" From 488983e9544b5c61fc200848ce8ac996b9b2e42f Mon Sep 17 00:00:00 2001 From: Pascal Aellig Date: Wed, 9 Oct 2024 11:15:41 +0200 Subject: [PATCH 12/16] rm temp scalar indexing from test --- test/test_thermalstresses.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/test_thermalstresses.jl b/test/test_thermalstresses.jl index dffde3ef..7426e07e 100644 --- a/test/test_thermalstresses.jl +++ b/test/test_thermalstresses.jl @@ -338,8 +338,8 @@ function main2D(igg; 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 From 66632d8b4d0111de82ff7ff005da1ccc0f69b05f Mon Sep 17 00:00:00 2001 From: Albert de Montserrat Date: Wed, 9 Oct 2024 15:53:48 +0200 Subject: [PATCH 13/16] couple of corrections --- src/rheology/Melting.jl | 13 +++---------- 1 file changed, 3 insertions(+), 10 deletions(-) diff --git a/src/rheology/Melting.jl b/src/rheology/Melting.jl index 67662f62..f82ae094 100644 --- a/src/rheology/Melting.jl +++ b/src/rheology/Melting.jl @@ -1,18 +1,15 @@ 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::JustPIC.PhaseRatios, rheology, args) ni = size(ϕ) @parallel (@idx ni) compute_melt_fraction_kernel!(ϕ, phase_ratios.center, rheology, args) @@ -22,10 +19,6 @@ end ϕ, phase_ratios, rheology, args ) args_ijk = ntuple_idx(args, I...) - ϕ[I...] = compute_melt_frac(rheology, args_ijk, @cell(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 From d39678bc3a01b3fbd1e41b63f9069450537b4635 Mon Sep 17 00:00:00 2001 From: aelligp Date: Sat, 12 Oct 2024 12:18:15 +0200 Subject: [PATCH 14/16] fix with specific test? --- test/test_melting.jl | 328 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 328 insertions(+) create mode 100644 test/test_melting.jl diff --git a/test/test_melting.jl b/test/test_melting.jl new file mode 100644 index 00000000..bc20ff52 --- /dev/null +++ b/test/test_melting.jl @@ -0,0 +1,328 @@ +push!(LOAD_PATH, "..") + +@static if ENV["JULIA_JUSTRELAX_BACKEND"] === "AMDGPU" + using AMDGPU +elseif ENV["JULIA_JUSTRELAX_BACKEND"] === "CUDA" + using CUDA +end + +using Test, Suppressor +using GeoParams +using JustRelax, JustRelax.JustRelax2D +using ParallelStencil, ParallelStencil.FiniteDifferences2D + +const backend_JR = @static if ENV["JULIA_JUSTRELAX_BACKEND"] === "AMDGPU" + @init_parallel_stencil(AMDGPU, Float64, 2) + AMDGPUBackend +elseif ENV["JULIA_JUSTRELAX_BACKEND"] === "CUDA" + @init_parallel_stencil(CUDA, Float64, 2) + CUDABackend +else + @init_parallel_stencil(Threads, Float64, 2) + CPUBackend +end + +using JustPIC, JustPIC._2D +# 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 ENV["JULIA_JUSTRELAX_BACKEND"] === "AMDGPU" + JustPIC.AMDGPUBackend +elseif ENV["JULIA_JUSTRELAX_BACKEND"] === "CUDA" + CUDABackend +else + JustPIC.CPUBackend +end + + +# Load script dependencies +using Printf, Statistics, LinearAlgebra, CellArrays, StaticArrays + + +# ----------------------------------------------------------------------------------------- +## SET OF HELPER FUNCTIONS PARTICULAR FOR THIS SCRIPT -------------------------------- +function copyinn_x!(A, B) + @parallel function f_x(A, B) + @all(A) = @inn_x(B) + return nothing + end + + @parallel f_x(A, B) +end + +import ParallelStencil.INDICES +const idx_j = INDICES[2] +macro all_j(A) + return esc(:($A[$idx_j])) +end + +@parallel function init_P!(P, ρg, z) + @all(P) = abs(@all(ρg) * @all_j(z)) * <(@all_j(z), 0.0) + return nothing +end + +function init_phases!(phases, particles, xc_anomaly, yc_anomaly, r_anomaly, sticky_air,top, bottom) + ni = size(phases) + + @parallel_indices (i, j) function init_phases!( + phases, px, py, index, xc_anomaly, yc_anomaly, r_anomaly, sticky_air, top, bottom + ) + @inbounds for ip in cellaxes(phases) + # quick escape + @index(index[ip, i, j]) == 0 && continue + + x = @index px[ip, i, j] + y = -(@index py[ip, i, j]) - sticky_air + if top ≤ y ≤ bottom + @index phases[ip, i, j] = 1.0 # crust + end + + # thermal anomaly - circular + if ((x - xc_anomaly)^2 + (y + yc_anomaly)^2 ≤ r_anomaly^2) + @index phases[ip, i, j] = 2.0 + end + + if y < top + @index phases[ip, i, j] = 3.0 + end + end + return nothing + end + + @parallel (@idx ni) init_phases!( + phases, + particles.coords..., + particles.index, + xc_anomaly, + yc_anomaly, + r_anomaly, + sticky_air, + top, + bottom, + ) +end + +# Initial thermal profile +@parallel_indices (i, j) function init_T!(T, y, sticky_air, top, bottom, dTdz, offset) + depth = -y[j] - sticky_air + + if depth < top + T[i + 1, j] = offset + + elseif top ≤ (depth) < bottom + dTdZ = dTdz + offset = offset + T[i + 1, j] = (depth) * dTdZ + offset + + end + + return nothing +end + +function circular_perturbation!(T, δT, xc_anomaly, yc_anomaly, r_anomaly, xvi, sticky_air) + @parallel_indices (i, j) function _circular_perturbation!( + T, δT, xc_anomaly, yc_anomaly, r_anomaly, x, y, sticky_air + ) + depth = -y[j] - sticky_air + @inbounds if ((x[i] - xc_anomaly)^2 + (depth[j] + yc_anomaly)^2 ≤ r_anomaly^2) + # T[i + 1, j] *= δT / 100 + 1 + T[i + 1, j] = δT + end + return nothing + end + + nx, ny = size(T) + + @parallel (1:(nx - 2), 1:ny) _circular_perturbation!( + T, δT, xc_anomaly, yc_anomaly, r_anomaly, xvi..., sticky_air + ) +end + +function init_rheology(CharDim;) + # plasticity setup + do_DP = true # do_DP=false: Von Mises, do_DP=true: Drucker-Prager (friction angle) + η_reg = 1.0e18Pa * s # regularisation "viscosity" for Drucker-Prager + Coh = 10.0MPa # yield stress. If do_DP=true, τ_y stand for the cohesion: c*cos(ϕ) + ϕ = 30.0 * do_DP # friction angle + G0 = 6.0e11Pa # elastic shear modulus + G_magma = 6.0e11Pa # elastic shear modulus perturbation + + pl = DruckerPrager_regularised(; C=Coh, ϕ=ϕ, η_vp=η_reg, Ψ=0.0) # plasticity + el = SetConstantElasticity(; G=G0, ν=0.5) # elastic spring + el_magma = SetConstantElasticity(; G=G_magma, ν=0.5) # elastic spring + + creep_rock = LinearViscous(; η=1e23 * Pa * s) + creep_magma = LinearViscous(; η=1e18 * Pa * s) + creep_air = LinearViscous(; η=1e18 * Pa * s) + + g = 9.81m/s^2 + rheology = ( + #Name="UpperCrust" + SetMaterialParams(; + Phase = 1, + Density = ConstantDensity(; ρ=2650kg / m^3), + HeatCapacity = ConstantHeatCapacity(; Cp=1050J / kg / K), + Conductivity = ConstantConductivity(; k=3.0Watt / K / m), + LatentHeat = ConstantLatentHeat(; Q_L=350e3J / kg), + ShearHeat = ConstantShearheating(1.0NoUnits), + CompositeRheology = CompositeRheology((creep_rock, el, pl)), + Melting = MeltingParam_Caricchi(), + Gravity = ConstantGravity(; g=g), + Elasticity = el, + CharDim = CharDim, + ), + + #Name="Magma" + SetMaterialParams(; + Phase = 1, + Density = ConstantDensity(; ρ=2650kg / m^3), + HeatCapacity = ConstantHeatCapacity(; Cp=1050J / kg / K), + Conductivity = ConstantConductivity(; k=1.5Watt / K / m), + LatentHeat = ConstantLatentHeat(; Q_L=350e3J / kg), + ShearHeat = ConstantShearheating(0.0NoUnits), + CompositeRheology = CompositeRheology((creep_magma, el_magma)), + Melting = MeltingParam_Caricchi(), + Gravity = ConstantGravity(; g=g), + Elasticity = el_magma, + CharDim = CharDim, + ), + + #Name="Sticky Air" + SetMaterialParams(; + Phase = 1, + Density = ConstantDensity(ρ=1kg/m^3,), + HeatCapacity = ConstantHeatCapacity(; Cp=1000J / kg / K), + Conductivity = ConstantConductivity(; k=15Watt / K / m), + LatentHeat = ConstantLatentHeat(; Q_L=0.0J / kg), + ShearHeat = ConstantShearheating(0.0NoUnits), + CompositeRheology = CompositeRheology((creep_air,)), + Gravity = ConstantGravity(; g=g), + CharDim = CharDim, + ), + ) + +end + + +function main2D(igg; nx=32, ny=32) + + # Characteristic lengths + CharDim = GEO_units(;length=12.5km, viscosity=1e21, temperature = 1e3C) + + #-------JustRelax parameters------------------------------------------------------------- + # Domain setup for JustRelax + sticky_air = nondimensionalize(1.5km, CharDim) # thickness of the sticky air layer + ly = nondimensionalize(12.5km,CharDim) + sticky_air # domain length in y-direction + lx = nondimensionalize(15.5km, CharDim) # domain length in x-direction + li = lx, ly # domain length in x- and y-direction + ni = nx, ny # number of grid points in x- and y-direction + di = @. li / ni # grid step in x- and y-direction + origin = nondimensionalize(0.0km,CharDim), -ly # origin coordinates of the domain + grid = Geometry(ni, li; origin=origin) + εbg = nondimensionalize(0.0 / s,CharDim) # background strain rate + (; xci, xvi) = grid # nodes at the center and vertices of the cells + #--------------------------------------------------------------------------------------- + + # Physical Parameters + rheology = init_rheology(CharDim; ) + cutoff_visc = nondimensionalize((1e16Pa*s, 1e24Pa*s),CharDim) + κ = (4 / (rheology[1].HeatCapacity[1].Cp * rheology[1].Density[1].ρ)) + dt = dt_diff = (0.5 * min(di...)^2 / κ / 2.01) # diffusive CFL timestep limiter + + # Initialize particles ------------------------------- + nxcell, max_xcell, min_xcell = 20, 40, 15 + particles = init_particles(backend, nxcell, max_xcell, min_xcell, xvi...) + subgrid_arrays = SubgridDiffusionCellArrays(particles) + # velocity grids + grid_vx, grid_vy = velocity_grids(xci, xvi, di) + # temperature + pT, pPhases = init_cell_arrays(particles, Val(2)) + particle_args = (pT, pPhases) + + # Circular temperature anomaly ----------------------- + x_anomaly = lx * 0.5 + y_anomaly = nondimensionalize(-5km,CharDim) # origin of the small thermal anomaly + r_anomaly = nondimensionalize(1.5km, CharDim) # radius of perturbation + anomaly = nondimensionalize((750 + 273)K, CharDim) # thermal perturbation (in K) + init_phases!(pPhases, particles, x_anomaly, y_anomaly, r_anomaly, sticky_air, nondimensionalize(0.0km,CharDim), nondimensionalize(20km,CharDim)) + phase_ratios = PhaseRatios(backend, length(rheology), ni) + update_phase_ratios!(phase_ratios, particles, xci, xvi, pPhases) + + # Initialisation of thermal profile + thermal = ThermalArrays(backend_JR, ni) # initialise thermal arrays and boundary conditions + thermal_bc = TemperatureBoundaryConditions(; + no_flux = (left=true, right=true, top=false, bot=false), + ) + @parallel (@idx ni .+ 1) init_T!( + thermal.T, xvi[2], + sticky_air, + nondimensionalize(0e0km,CharDim), + nondimensionalize(15km,CharDim), + nondimensionalize((723 - 273)K,CharDim) / nondimensionalize(15km,CharDim), + nondimensionalize(273K,CharDim) + ) + circular_perturbation!( + thermal.T, anomaly, x_anomaly, y_anomaly, r_anomaly, xvi, sticky_air + ) + thermal_bcs!(thermal, thermal_bc) + temperature2center!(thermal) + + # STOKES --------------------------------------------- + # Allocate arrays needed for every Stokes problem + stokes = StokesArrays(backend_JR, ni) # initialise stokes arrays with the defined regime + pt_stokes = PTStokesCoeffs(li, di; ϵ = 1e-4, CFL = 1 / √2.1) + # ---------------------------------------------------- + + args = (; T=thermal.Tc, P=stokes.P, dt=dt) + pt_thermal = PTThermalCoeffs( + backend_JR, rheology, phase_ratios, args, dt, ni, di, li; ϵ=1e-5, CFL=0.8 / √2.1 + ) + + # Pure shear far-field boundary conditions + stokes.V.Vx .= PTArray(backend_JR)([ + εbg * (x - lx * 0.5) / (lx / 2) / 2 for x in xvi[1], _ in 1:(ny + 2) + ]) + stokes.V.Vy .= PTArray(backend_JR)([ + (abs(y) - sticky_air) * εbg * (abs(y) > sticky_air) for _ in 1:(nx + 2), y in xvi[2] + ]) + + flow_bcs = VelocityBoundaryConditions(; + free_slip = (left=true, right=true, top=true, bot=true), + free_surface = true, + ) + flow_bcs!(stokes, flow_bcs) + + compute_viscosity!(stokes, phase_ratios, args, rheology, cutoff_visc) + + ϕ = @zeros(ni...) + + # 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 + for _ in 1:5 + compute_ρg!(ρg[2], phase_ratios, rheology, (T=thermal.Tc, P=stokes.P)) + @parallel init_P!(stokes.P, ρg[2], xci[2]) + compute_melt_fraction!( + ϕ, phase_ratios, rheology, (T=thermal.Tc, P=stokes.P) + ) + end + + return ϕ, stokes, thermal +end + +@testset "thermal stresses" begin + # @suppress begin + 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.8035 rtol = 1e-2 + @test Array(ϕ)[nx_T >>> 1 + 1, ny_T >>> 1 + 1] ≈ 0.10152 rtol = 1e-1 + + # end +end From 49dbaa05da1ce944704ef3552f89528dcb99cb98 Mon Sep 17 00:00:00 2001 From: aelligp Date: Sat, 12 Oct 2024 13:20:48 +0200 Subject: [PATCH 15/16] delete test again after local pass --- test/test_melting.jl | 328 ------------------------------------------- 1 file changed, 328 deletions(-) delete mode 100644 test/test_melting.jl diff --git a/test/test_melting.jl b/test/test_melting.jl deleted file mode 100644 index bc20ff52..00000000 --- a/test/test_melting.jl +++ /dev/null @@ -1,328 +0,0 @@ -push!(LOAD_PATH, "..") - -@static if ENV["JULIA_JUSTRELAX_BACKEND"] === "AMDGPU" - using AMDGPU -elseif ENV["JULIA_JUSTRELAX_BACKEND"] === "CUDA" - using CUDA -end - -using Test, Suppressor -using GeoParams -using JustRelax, JustRelax.JustRelax2D -using ParallelStencil, ParallelStencil.FiniteDifferences2D - -const backend_JR = @static if ENV["JULIA_JUSTRELAX_BACKEND"] === "AMDGPU" - @init_parallel_stencil(AMDGPU, Float64, 2) - AMDGPUBackend -elseif ENV["JULIA_JUSTRELAX_BACKEND"] === "CUDA" - @init_parallel_stencil(CUDA, Float64, 2) - CUDABackend -else - @init_parallel_stencil(Threads, Float64, 2) - CPUBackend -end - -using JustPIC, JustPIC._2D -# 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 ENV["JULIA_JUSTRELAX_BACKEND"] === "AMDGPU" - JustPIC.AMDGPUBackend -elseif ENV["JULIA_JUSTRELAX_BACKEND"] === "CUDA" - CUDABackend -else - JustPIC.CPUBackend -end - - -# Load script dependencies -using Printf, Statistics, LinearAlgebra, CellArrays, StaticArrays - - -# ----------------------------------------------------------------------------------------- -## SET OF HELPER FUNCTIONS PARTICULAR FOR THIS SCRIPT -------------------------------- -function copyinn_x!(A, B) - @parallel function f_x(A, B) - @all(A) = @inn_x(B) - return nothing - end - - @parallel f_x(A, B) -end - -import ParallelStencil.INDICES -const idx_j = INDICES[2] -macro all_j(A) - return esc(:($A[$idx_j])) -end - -@parallel function init_P!(P, ρg, z) - @all(P) = abs(@all(ρg) * @all_j(z)) * <(@all_j(z), 0.0) - return nothing -end - -function init_phases!(phases, particles, xc_anomaly, yc_anomaly, r_anomaly, sticky_air,top, bottom) - ni = size(phases) - - @parallel_indices (i, j) function init_phases!( - phases, px, py, index, xc_anomaly, yc_anomaly, r_anomaly, sticky_air, top, bottom - ) - @inbounds for ip in cellaxes(phases) - # quick escape - @index(index[ip, i, j]) == 0 && continue - - x = @index px[ip, i, j] - y = -(@index py[ip, i, j]) - sticky_air - if top ≤ y ≤ bottom - @index phases[ip, i, j] = 1.0 # crust - end - - # thermal anomaly - circular - if ((x - xc_anomaly)^2 + (y + yc_anomaly)^2 ≤ r_anomaly^2) - @index phases[ip, i, j] = 2.0 - end - - if y < top - @index phases[ip, i, j] = 3.0 - end - end - return nothing - end - - @parallel (@idx ni) init_phases!( - phases, - particles.coords..., - particles.index, - xc_anomaly, - yc_anomaly, - r_anomaly, - sticky_air, - top, - bottom, - ) -end - -# Initial thermal profile -@parallel_indices (i, j) function init_T!(T, y, sticky_air, top, bottom, dTdz, offset) - depth = -y[j] - sticky_air - - if depth < top - T[i + 1, j] = offset - - elseif top ≤ (depth) < bottom - dTdZ = dTdz - offset = offset - T[i + 1, j] = (depth) * dTdZ + offset - - end - - return nothing -end - -function circular_perturbation!(T, δT, xc_anomaly, yc_anomaly, r_anomaly, xvi, sticky_air) - @parallel_indices (i, j) function _circular_perturbation!( - T, δT, xc_anomaly, yc_anomaly, r_anomaly, x, y, sticky_air - ) - depth = -y[j] - sticky_air - @inbounds if ((x[i] - xc_anomaly)^2 + (depth[j] + yc_anomaly)^2 ≤ r_anomaly^2) - # T[i + 1, j] *= δT / 100 + 1 - T[i + 1, j] = δT - end - return nothing - end - - nx, ny = size(T) - - @parallel (1:(nx - 2), 1:ny) _circular_perturbation!( - T, δT, xc_anomaly, yc_anomaly, r_anomaly, xvi..., sticky_air - ) -end - -function init_rheology(CharDim;) - # plasticity setup - do_DP = true # do_DP=false: Von Mises, do_DP=true: Drucker-Prager (friction angle) - η_reg = 1.0e18Pa * s # regularisation "viscosity" for Drucker-Prager - Coh = 10.0MPa # yield stress. If do_DP=true, τ_y stand for the cohesion: c*cos(ϕ) - ϕ = 30.0 * do_DP # friction angle - G0 = 6.0e11Pa # elastic shear modulus - G_magma = 6.0e11Pa # elastic shear modulus perturbation - - pl = DruckerPrager_regularised(; C=Coh, ϕ=ϕ, η_vp=η_reg, Ψ=0.0) # plasticity - el = SetConstantElasticity(; G=G0, ν=0.5) # elastic spring - el_magma = SetConstantElasticity(; G=G_magma, ν=0.5) # elastic spring - - creep_rock = LinearViscous(; η=1e23 * Pa * s) - creep_magma = LinearViscous(; η=1e18 * Pa * s) - creep_air = LinearViscous(; η=1e18 * Pa * s) - - g = 9.81m/s^2 - rheology = ( - #Name="UpperCrust" - SetMaterialParams(; - Phase = 1, - Density = ConstantDensity(; ρ=2650kg / m^3), - HeatCapacity = ConstantHeatCapacity(; Cp=1050J / kg / K), - Conductivity = ConstantConductivity(; k=3.0Watt / K / m), - LatentHeat = ConstantLatentHeat(; Q_L=350e3J / kg), - ShearHeat = ConstantShearheating(1.0NoUnits), - CompositeRheology = CompositeRheology((creep_rock, el, pl)), - Melting = MeltingParam_Caricchi(), - Gravity = ConstantGravity(; g=g), - Elasticity = el, - CharDim = CharDim, - ), - - #Name="Magma" - SetMaterialParams(; - Phase = 1, - Density = ConstantDensity(; ρ=2650kg / m^3), - HeatCapacity = ConstantHeatCapacity(; Cp=1050J / kg / K), - Conductivity = ConstantConductivity(; k=1.5Watt / K / m), - LatentHeat = ConstantLatentHeat(; Q_L=350e3J / kg), - ShearHeat = ConstantShearheating(0.0NoUnits), - CompositeRheology = CompositeRheology((creep_magma, el_magma)), - Melting = MeltingParam_Caricchi(), - Gravity = ConstantGravity(; g=g), - Elasticity = el_magma, - CharDim = CharDim, - ), - - #Name="Sticky Air" - SetMaterialParams(; - Phase = 1, - Density = ConstantDensity(ρ=1kg/m^3,), - HeatCapacity = ConstantHeatCapacity(; Cp=1000J / kg / K), - Conductivity = ConstantConductivity(; k=15Watt / K / m), - LatentHeat = ConstantLatentHeat(; Q_L=0.0J / kg), - ShearHeat = ConstantShearheating(0.0NoUnits), - CompositeRheology = CompositeRheology((creep_air,)), - Gravity = ConstantGravity(; g=g), - CharDim = CharDim, - ), - ) - -end - - -function main2D(igg; nx=32, ny=32) - - # Characteristic lengths - CharDim = GEO_units(;length=12.5km, viscosity=1e21, temperature = 1e3C) - - #-------JustRelax parameters------------------------------------------------------------- - # Domain setup for JustRelax - sticky_air = nondimensionalize(1.5km, CharDim) # thickness of the sticky air layer - ly = nondimensionalize(12.5km,CharDim) + sticky_air # domain length in y-direction - lx = nondimensionalize(15.5km, CharDim) # domain length in x-direction - li = lx, ly # domain length in x- and y-direction - ni = nx, ny # number of grid points in x- and y-direction - di = @. li / ni # grid step in x- and y-direction - origin = nondimensionalize(0.0km,CharDim), -ly # origin coordinates of the domain - grid = Geometry(ni, li; origin=origin) - εbg = nondimensionalize(0.0 / s,CharDim) # background strain rate - (; xci, xvi) = grid # nodes at the center and vertices of the cells - #--------------------------------------------------------------------------------------- - - # Physical Parameters - rheology = init_rheology(CharDim; ) - cutoff_visc = nondimensionalize((1e16Pa*s, 1e24Pa*s),CharDim) - κ = (4 / (rheology[1].HeatCapacity[1].Cp * rheology[1].Density[1].ρ)) - dt = dt_diff = (0.5 * min(di...)^2 / κ / 2.01) # diffusive CFL timestep limiter - - # Initialize particles ------------------------------- - nxcell, max_xcell, min_xcell = 20, 40, 15 - particles = init_particles(backend, nxcell, max_xcell, min_xcell, xvi...) - subgrid_arrays = SubgridDiffusionCellArrays(particles) - # velocity grids - grid_vx, grid_vy = velocity_grids(xci, xvi, di) - # temperature - pT, pPhases = init_cell_arrays(particles, Val(2)) - particle_args = (pT, pPhases) - - # Circular temperature anomaly ----------------------- - x_anomaly = lx * 0.5 - y_anomaly = nondimensionalize(-5km,CharDim) # origin of the small thermal anomaly - r_anomaly = nondimensionalize(1.5km, CharDim) # radius of perturbation - anomaly = nondimensionalize((750 + 273)K, CharDim) # thermal perturbation (in K) - init_phases!(pPhases, particles, x_anomaly, y_anomaly, r_anomaly, sticky_air, nondimensionalize(0.0km,CharDim), nondimensionalize(20km,CharDim)) - phase_ratios = PhaseRatios(backend, length(rheology), ni) - update_phase_ratios!(phase_ratios, particles, xci, xvi, pPhases) - - # Initialisation of thermal profile - thermal = ThermalArrays(backend_JR, ni) # initialise thermal arrays and boundary conditions - thermal_bc = TemperatureBoundaryConditions(; - no_flux = (left=true, right=true, top=false, bot=false), - ) - @parallel (@idx ni .+ 1) init_T!( - thermal.T, xvi[2], - sticky_air, - nondimensionalize(0e0km,CharDim), - nondimensionalize(15km,CharDim), - nondimensionalize((723 - 273)K,CharDim) / nondimensionalize(15km,CharDim), - nondimensionalize(273K,CharDim) - ) - circular_perturbation!( - thermal.T, anomaly, x_anomaly, y_anomaly, r_anomaly, xvi, sticky_air - ) - thermal_bcs!(thermal, thermal_bc) - temperature2center!(thermal) - - # STOKES --------------------------------------------- - # Allocate arrays needed for every Stokes problem - stokes = StokesArrays(backend_JR, ni) # initialise stokes arrays with the defined regime - pt_stokes = PTStokesCoeffs(li, di; ϵ = 1e-4, CFL = 1 / √2.1) - # ---------------------------------------------------- - - args = (; T=thermal.Tc, P=stokes.P, dt=dt) - pt_thermal = PTThermalCoeffs( - backend_JR, rheology, phase_ratios, args, dt, ni, di, li; ϵ=1e-5, CFL=0.8 / √2.1 - ) - - # Pure shear far-field boundary conditions - stokes.V.Vx .= PTArray(backend_JR)([ - εbg * (x - lx * 0.5) / (lx / 2) / 2 for x in xvi[1], _ in 1:(ny + 2) - ]) - stokes.V.Vy .= PTArray(backend_JR)([ - (abs(y) - sticky_air) * εbg * (abs(y) > sticky_air) for _ in 1:(nx + 2), y in xvi[2] - ]) - - flow_bcs = VelocityBoundaryConditions(; - free_slip = (left=true, right=true, top=true, bot=true), - free_surface = true, - ) - flow_bcs!(stokes, flow_bcs) - - compute_viscosity!(stokes, phase_ratios, args, rheology, cutoff_visc) - - ϕ = @zeros(ni...) - - # 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 - for _ in 1:5 - compute_ρg!(ρg[2], phase_ratios, rheology, (T=thermal.Tc, P=stokes.P)) - @parallel init_P!(stokes.P, ρg[2], xci[2]) - compute_melt_fraction!( - ϕ, phase_ratios, rheology, (T=thermal.Tc, P=stokes.P) - ) - end - - return ϕ, stokes, thermal -end - -@testset "thermal stresses" begin - # @suppress begin - 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.8035 rtol = 1e-2 - @test Array(ϕ)[nx_T >>> 1 + 1, ny_T >>> 1 + 1] ≈ 0.10152 rtol = 1e-1 - - # end -end From bfaad2028ecee89c77ddb9c667a824f725dc4690 Mon Sep 17 00:00:00 2001 From: aelligp Date: Sat, 12 Oct 2024 13:33:14 +0200 Subject: [PATCH 16/16] format; add stop condition for CI on new push --- .github/workflows/ci.yml | 6 ++++++ src/rheology/Melting.jl | 4 +++- 2 files changed, 9 insertions(+), 1 deletion(-) 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/src/rheology/Melting.jl b/src/rheology/Melting.jl index f82ae094..04ed5c82 100644 --- a/src/rheology/Melting.jl +++ b/src/rheology/Melting.jl @@ -12,7 +12,9 @@ end function compute_melt_fraction!(ϕ, phase_ratios::JustPIC.PhaseRatios, rheology, args) ni = size(ϕ) - @parallel (@idx ni) compute_melt_fraction_kernel!(ϕ, phase_ratios.center, rheology, args) + @parallel (@idx ni) compute_melt_fraction_kernel!( + ϕ, phase_ratios.center, rheology, args + ) end @parallel_indices (I...) function compute_melt_fraction_kernel!(