Skip to content

Commit

Permalink
add Pressure plot, add arrows
Browse files Browse the repository at this point in the history
  • Loading branch information
@aelligp
aelligp committed Apr 4, 2024
1 parent a8c1be9 commit cdc2f77
Show file tree
Hide file tree
Showing 2 changed files with 94 additions and 37 deletions.
69 changes: 48 additions & 21 deletions miniapps/benchmarks/thermal_stress/Thermal_Stress_Magma_Chamber.jl
View file
Original file line number Diff line number Diff line change
@@ -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
Expand Down Expand Up @@ -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,
Expand All @@ -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,
Expand All @@ -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,
Expand All @@ -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,
Expand All @@ -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)
Expand All @@ -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)
Expand All @@ -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)
Expand Down
62 changes: 46 additions & 16 deletions miniapps/benchmarks/thermal_stress/Thermal_Stress_Magma_Chamber_nondim.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand All @@ -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,
Expand All @@ -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,
Expand All @@ -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,
Expand All @@ -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,
)
Expand All @@ -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)
Expand Down

0 comments on commit cdc2f77

Please sign in to comment.