From cdc2f779b7d02de9ed4709d53dabce9ab80fa871 Mon Sep 17 00:00:00 2001 From: aelligp Date: Thu, 4 Apr 2024 11:38:16 +0200 Subject: [PATCH] add Pressure plot, add arrows --- .../Thermal_Stress_Magma_Chamber.jl | 69 +++++++++++++------ .../Thermal_Stress_Magma_Chamber_nondim.jl | 62 ++++++++++++----- 2 files changed, 94 insertions(+), 37 deletions(-) diff --git a/miniapps/benchmarks/thermal_stress/Thermal_Stress_Magma_Chamber.jl b/miniapps/benchmarks/thermal_stress/Thermal_Stress_Magma_Chamber.jl index c9bcb03a..b539689b 100644 --- a/miniapps/benchmarks/thermal_stress/Thermal_Stress_Magma_Chamber.jl +++ b/miniapps/benchmarks/thermal_stress/Thermal_Stress_Magma_Chamber.jl @@ -1,18 +1,17 @@ -using CUDA using JustRelax, JustRelax.DataIO import JustRelax.@cell using ParallelStencil -@init_parallel_stencil(CUDA, Float64, 2) #or (CUDA, Float64, 2) or (AMDGPU, Float64, 2) +@init_parallel_stencil(Threads, Float64, 2) #or (CUDA, Float64, 2) or (AMDGPU, Float64, 2) using JustPIC using 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 = CUDABackend # Options: CPUBackend, CUDABackend, AMDGPUBackend +const backend = CPUBackend # Options: CPUBackend, CUDABackend, AMDGPUBackend # setup ParallelStencil.jl environment -model = PS_Setup(:CUDA, Float64, 2) #or (:CUDA, Float64, 2) or (:AMDGPU, Float64, 2) +model = PS_Setup(:Threads, Float64, 2) #or (:CUDA, Float64, 2) or (:AMDGPU, Float64, 2) environment!(model) using Printf, Statistics, LinearAlgebra, GeoParams, CairoMakie, CellArrays @@ -566,7 +565,8 @@ function main2D(igg; figdir=figdir, nx=nx, ny=ny, do_vtk=false) ax2 = Axis( fig[2, 2][1, 1]; aspect=ar, - title=L"ΔP [MPa]", + title=L"Viscosity [\mathrm{Pa s}]", + xlabel="Width [km]", titlesize=40, yticklabelsize=25, xticklabelsize=25, @@ -575,8 +575,7 @@ function main2D(igg; figdir=figdir, nx=nx, ny=ny, do_vtk=false) ax3 = Axis( fig[3, 1][1, 1]; aspect=ar, - title=L"Viscosity [\mathrm{Pa s}]", - xlabel="Width [km]", + title=L"ΔP [MPa]", titlesize=40, yticklabelsize=25, xticklabelsize=25, @@ -585,8 +584,7 @@ function main2D(igg; figdir=figdir, nx=nx, ny=ny, do_vtk=false) ax4 = Axis( fig[3, 2][1, 1]; aspect=ar, - title=L"\log_{10}(\dot{\varepsilon}_{\textrm{II}}) [\mathrm{s}^{-1}]", - xlabel="Width [km]", + title=L"P [MPa]", titlesize=40, yticklabelsize=25, xticklabelsize=25, @@ -595,6 +593,16 @@ function main2D(igg; figdir=figdir, nx=nx, ny=ny, do_vtk=false) ax5 = Axis( fig[4, 1][1, 1]; aspect=ar, + title=L"\log_{10}(\dot{\varepsilon}_{\textrm{II}}) [\mathrm{s}^{-1}]", + xlabel="Width [km]", + titlesize=40, + yticklabelsize=25, + xticklabelsize=25, + xlabelsize=25, + ) + ax6 = Axis( + fig[4, 2][1, 1]; + aspect=ar, title=L"\tau_{\textrm{II}} [MPa]", xlabel="Width [km]", titlesize=40, @@ -610,37 +618,53 @@ function main2D(igg; figdir=figdir, nx=nx, ny=ny, do_vtk=false) Array(thermal.T[2:(end - 1), :] .- 273); colormap=:batlow, ) - # Plot Pressure difference + # Plot effective viscosity p2 = heatmap!( ax2, xci[1] .* 1e-3, xci[2] .* 1e-3, - Array((stokes.P .- P_init) ./ 1e6); - colormap=:roma, + Array(log10.(η_vep)); + colormap=:glasgow, + colorrange=(log10(1e16), log10(1e24)), ) - # Plot effective viscosity + arrows!( + ax2, + xvi[1][1:5:end-1]./1e3, xvi[2][1:5:end-1]./1e3, + Array.((Vx_v[1:5:end-1, 1:5:end-1], + Vy_v[1:5:end-1, 1:5:end-1]))..., + lengthscale = 2 / max(maximum(Vx_v), maximum(Vy_v)), + color = :darkblue, + ) + # Plot Pressure difference p3 = heatmap!( ax3, xci[1] .* 1e-3, xci[2] .* 1e-3, - Array(log10.(η_vep)); - colormap=:glasgow, - colorrange=(log10(1e16), log10(1e24)), + Array((stokes.P .- P_init) ./ 1e6); + colormap=:roma, ) - # Plot 2nd invariant of strain rate + # Plot Pressure difference p4 = heatmap!( ax4, xci[1] .* 1e-3, xci[2] .* 1e-3, - Array(log10.(stokes.ε.II)); + Array((stokes.P) ./ 1e6); colormap=:roma, ) - # Plot 2nd invariant of stress + # Plot 2nd invariant of strain rate p5 = heatmap!( ax5, xci[1] .* 1e-3, xci[2] .* 1e-3, - Array((stokes.τ.II)); + Array(log10.(stokes.ε.II)); + colormap=:roma, + ) + # Plot 2nd invariant of stress + p6 = heatmap!( + ax6, + xci[1] .* 1e-3, + xci[2] .* 1e-3, + Array((stokes.τ.II)./1e6); colormap=:batlow, ) hidexdecorations!(ax1) @@ -661,6 +685,9 @@ function main2D(igg; figdir=figdir, nx=nx, ny=ny, do_vtk=false) Colorbar( fig[4, 1][1, 2], p5; height=Relative(0.7), ticklabelsize=25, ticksize=15 ) + Colorbar( + fig[4, 2][1, 2], p6; height=Relative(0.7), ticklabelsize=25, ticksize=15 + ) rowgap!(fig.layout, 1) colgap!(fig.layout, 1) colgap!(fig.layout, 1) @@ -679,7 +706,7 @@ function main2D(igg; figdir=figdir, nx=nx, ny=ny, do_vtk=false) scatter!(ax1, Array(thermal.T[2:(end - 1), :][:] .- 273.0), Yv) lines!(ax2, Array(stokes.P[:]), Y) - scatter!(a3, Array((stokes.τ.II[:])), Y) + scatter!(a3, Array((stokes.τ.II[:])./1e6), Y) hideydecorations!(ax2) save(joinpath(figdir, "pressure_profile_$it.png"), fig) 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 a6439381..3f9ab99f 100644 --- a/miniapps/benchmarks/thermal_stress/Thermal_Stress_Magma_Chamber_nondim.jl +++ b/miniapps/benchmarks/thermal_stress/Thermal_Stress_Magma_Chamber_nondim.jl @@ -604,11 +604,12 @@ function main2D(igg; figdir=figdir, nx=nx, ny=ny, do_vtk=false) yticklabelsize=25, xticklabelsize=25, xlabelsize=25, - ) + ) ax2 = Axis( fig[2, 2][1, 1]; aspect=ar, - title=L"ΔP [MPa]", + title=L"Viscosity [\mathrm{Pa s}]", + xlabel="Width [km]", titlesize=40, yticklabelsize=25, xticklabelsize=25, @@ -617,8 +618,7 @@ function main2D(igg; figdir=figdir, nx=nx, ny=ny, do_vtk=false) ax3 = Axis( fig[3, 1][1, 1]; aspect=ar, - title=L"Viscosity [\mathrm{Pa s}]", - xlabel="Width [km]", + title=L"ΔP [MPa]", titlesize=40, yticklabelsize=25, xticklabelsize=25, @@ -627,8 +627,7 @@ function main2D(igg; figdir=figdir, nx=nx, ny=ny, do_vtk=false) ax4 = Axis( fig[3, 2][1, 1]; aspect=ar, - title=L"\log_{10}(\dot{\varepsilon}_{\textrm{II}}) [\mathrm{s}^{-1}]", - xlabel="Width [km]", + title=L"P [MPa]", titlesize=40, yticklabelsize=25, xticklabelsize=25, @@ -637,6 +636,16 @@ function main2D(igg; figdir=figdir, nx=nx, ny=ny, do_vtk=false) ax5 = Axis( fig[4, 1][1, 1]; aspect=ar, + title=L"\log_{10}(\dot{\varepsilon}_{\textrm{II}}) [\mathrm{s}^{-1}]", + xlabel="Width [km]", + titlesize=40, + yticklabelsize=25, + xticklabelsize=25, + xlabelsize=25, + ) + ax6 = Axis( + fig[4, 2][1, 1]; + aspect=ar, title=L"\tau_{\textrm{II}} [MPa]", xlabel="Width [km]", titlesize=40, @@ -652,36 +661,54 @@ function main2D(igg; figdir=figdir, nx=nx, ny=ny, do_vtk=false) ustrip.(dimensionalize((Array(thermal.T[2:(end - 1), :])),C,CharDim)); colormap=:batlow, ) - # Plot Pressure difference + # Plot effective viscosity p2 = heatmap!( ax2, ustrip.(dimensionalize(xci[1],km,CharDim)), ustrip.(dimensionalize(xci[2],km,CharDim)), - ustrip.(dimensionalize((Array((stokes.P .- P_init))),MPa,CharDim)); - colormap=:roma, + ustrip.(dimensionalize((Array(log10.(η_vep))),Pa*s,CharDim)); + colormap=:glasgow, + colorrange=(log10(1e16), log10(1e24)), ) - # Plot effective viscosity + arrows!( + ax2, + ustrip.(dimensionalize(xvi[1],km,CharDim))[1:5:end-1], + ustrip.(dimensionalize(xvi[2],km,CharDim))[1:5:end-1], + Array.((ustrip.(dimensionalize(Vx_v, cm/yr, CharDim))[1:5:end-1, 1:5:end-1], + ustrip.(dimensionalize(Vy_v, cm/yr, CharDim))[1:5:end-1, 1:5:end-1]))..., + lengthscale = 1 / max(maximum(ustrip.(dimensionalize(Vx_v, cm/yr, CharDim))), + maximum(ustrip.(dimensionalize(Vy_v, cm/yr, CharDim)))), + color = :red, + ) + # Plot Pressure difference p3 = heatmap!( ax3, ustrip.(dimensionalize(xci[1],km,CharDim)), ustrip.(dimensionalize(xci[2],km,CharDim)), - ustrip.(dimensionalize((Array(log10.(η_vep))),Pa*s,CharDim)); - colormap=:glasgow, - colorrange=(log10(1e16), log10(1e24)), + ustrip.(dimensionalize((Array((stokes.P .- P_init))),MPa,CharDim)); + colormap=:roma, ) - # Plot 2nd invariant of strain rate + # Plot Pressure difference p4 = heatmap!( ax4, ustrip.(dimensionalize(xci[1],km,CharDim)), ustrip.(dimensionalize(xci[2],km,CharDim)), - log10.(ustrip.(dimensionalize(Array((stokes.ε.II)),s^-1,CharDim))); + ustrip.(dimensionalize((Array((stokes.P))),MPa,CharDim)); colormap=:roma, ) - # Plot 2nd invariant of stress + # Plot 2nd invariant of strain rate p5 = heatmap!( ax5, ustrip.(dimensionalize(xci[1],km,CharDim)), ustrip.(dimensionalize(xci[2],km,CharDim)), + log10.(ustrip.(dimensionalize(Array((stokes.ε.II)),s^-1,CharDim))); + colormap=:roma, + ) + # Plot 2nd invariant of stress + p6 = heatmap!( + ax6, + ustrip.(dimensionalize(xci[1],km,CharDim)), + ustrip.(dimensionalize(xci[2],km,CharDim)), ustrip.(dimensionalize(Array((stokes.τ.II)),MPa,CharDim)); colormap=:batlow, ) @@ -703,6 +730,9 @@ function main2D(igg; figdir=figdir, nx=nx, ny=ny, do_vtk=false) Colorbar( fig[4, 1][1, 2], p5; height=Relative(0.7), ticklabelsize=25, ticksize=15 ) + Colorbar( + fig[4, 2][1, 2], p6; height=Relative(0.7), ticklabelsize=25, ticksize=15 + ) rowgap!(fig.layout, 1) colgap!(fig.layout, 1) colgap!(fig.layout, 1)