Skip to content

Commit

Permalink
up miniapps
Browse files Browse the repository at this point in the history
  • Loading branch information
albert-de-montserrat committed Apr 23, 2024
1 parent cc92215 commit 7afa625
Show file tree
Hide file tree
Showing 6 changed files with 593 additions and 41 deletions.
Original file line number Diff line number Diff line change
Expand Up @@ -204,7 +204,7 @@ function main2D(igg; ar=8, ny=16, nx=ny*8, figdir="figs2D", do_vtk =false)
ax1 = Axis(fig[1,1], aspect = 2/3, title = "T")
ax2 = Axis(fig[1,2], aspect = 2/3, title = "log10(η)")
scatter!(ax1, Array(thermal.T[2:end-1,:][:]), Yv)
scatter!(ax2, Array(log10.(A[:])), Y)
scatter!(ax2, Array(log10.(η[:])), Y)
ylims!(ax1, minimum(xvi[2]), 0)
ylims!(ax2, minimum(xvi[2]), 0)
hideydecorations!(ax2)
Expand Down Expand Up @@ -395,4 +395,4 @@ else
end

# run main script
main2D(igg; figdir = figdir, ar = ar, nx = nx, ny = ny, do_vtk = do_vtk);
# main2D(igg; figdir = figdir, ar = ar, nx = nx, ny = ny, do_vtk = do_vtk);
16 changes: 13 additions & 3 deletions src/stokes/Stokes2D.jl
Original file line number Diff line number Diff line change
Expand Up @@ -556,6 +556,7 @@ function JustRelax.solve!(

wtime0 += @elapsed begin
compute_maxloc!(ητ, η; window=(1, 1))
# compute_maxloc!(ητ, η_vep; window=(1, 1))
update_halo!(ητ)

@parallel (@idx ni) compute_∇V!(stokes.∇V, @velocity(stokes)..., _di...)
Expand Down Expand Up @@ -628,16 +629,25 @@ function JustRelax.solve!(
)

@hide_communication b_width begin # communication/computation overlap
# @parallel compute_V!(
# @velocity(stokes)...,
# Vx_on_Vy,
# θ,
# @stress(stokes)...,
# pt_stokes.ηdτ,
# ρg...,
# ητ,
# _di...,
# dt * free_surface,
# )
@parallel compute_V!(
@velocity(stokes)...,
Vx_on_Vy,
θ,
stokes.P,
@stress(stokes)...,
pt_stokes.ηdτ,
ρg...,
ητ,
_di...,
dt * free_surface,
)
# apply boundary conditions
free_surface_bcs!(stokes, flow_bcs, η, rheology, phase_ratios, dt, di)
Expand Down
56 changes: 28 additions & 28 deletions subduction/LAMEM2D.jl
Original file line number Diff line number Diff line change
@@ -1,21 +1,21 @@
using CUDA
# using CUDA
using JustRelax, JustRelax.DataIO
import JustRelax.@cell
using ParallelStencil
# @init_parallel_stencil(Threads, Float64, 2)
@init_parallel_stencil(CUDA, Float64, 2)
@init_parallel_stencil(Threads, Float64, 2)
# @init_parallel_stencil(CUDA, 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 = CPUBackend # Options: CPUBackend, CUDABackend, AMDGPUBackend
const backend = CUDABackend # Options: CPUBackend, CUDABackend, AMDGPUBackend
const backend = CPUBackend # Options: CPUBackend, CUDABackend, AMDGPUBackend
# const backend = CUDABackend # Options: CPUBackend, CUDABackend, AMDGPUBackend

# setup ParallelStencil.jl environment
# model = PS_Setup(:cpu, Float64, 2) # or (:CUDA, Float64, 3) or (:AMDGPU, Float64, 3)
model = PS_Setup(:CUDA, Float64, 2) # or (:CUDA, Float64, 3) or (:AMDGPU, Float64, 3)
model = PS_Setup(:cpu, Float64, 2) # or (:CUDA, Float64, 3) or (:AMDGPU, Float64, 3)
# model = PS_Setup(:CUDA, Float64, 2) # or (:CUDA, Float64, 3) or (:AMDGPU, Float64, 3)
environment!(model)

# Load script dependencies
Expand Down Expand Up @@ -93,8 +93,8 @@ function main(li, origin, phases_GMG, igg; nx=16, ny=16, figdir="figs2D", do_vtk
# ----------------------------------------------------

# TEMPERATURE PROFILE --------------------------------
Ttop = 20
Tbot = 1565.0
Ttop = 20 + 273
Tbot = 1565.0 + 273
thermal = ThermalArrays(ni)
@views thermal.T[2:end-1, :] .= PTArray(T_GMG)
thermal_bc = TemperatureBoundaryConditions(;
Expand Down Expand Up @@ -134,22 +134,22 @@ function main(li, origin, phases_GMG, igg; nx=16, ny=16, figdir="figs2D", do_vtk
flow_bcs!(stokes, flow_bcs) # apply boundary conditions
update_halo!(@velocity(stokes)...)

# # Plot initial T and η profiles
# let
# Yv = [y for x in xvi[1], y in xvi[2]][:]
# Y = [y for x in xci[1], y in xci[2]][:]
# fig = Figure(size = (1200, 900))
# ax1 = Axis(fig[1,1], aspect = 2/3, title = "T")
# ax2 = Axis(fig[1,2], aspect = 2/3, title = "log10(η)")
# scatter!(ax1, Array(thermal.T[2:end-1,:][:]), Yv./1e3)
# scatter!(ax2, Array(log10.(η[:])), Y./1e3)
# # scatter!(ax2, Array(ρg[2][:]), Y./1e3)
# ylims!(ax1, minimum(xvi[2])./1e3, 0)
# ylims!(ax2, minimum(xvi[2])./1e3, 0)
# hideydecorations!(ax2)
# # save(joinpath(figdir, "initial_profile.png"), fig)
# fig
# end
# Plot initial T and η profiles
let
Yv = [y for x in xvi[1], y in xvi[2]][:]
Y = [y for x in xci[1], y in xci[2]][:]
fig = Figure(size = (1200, 900))
ax1 = Axis(fig[1,1], aspect = 2/3, title = "T")
ax2 = Axis(fig[1,2], aspect = 2/3, title = "log10(η)")
scatter!(ax1, Array(thermal.T[2:end-1,:][:]), Yv./1e3)
scatter!(ax2, Array(log10.(η[:])), Y./1e3)
# scatter!(ax2, Array(ρg[2][:]), Y./1e3)
ylims!(ax1, minimum(xvi[2])./1e3, 0)
ylims!(ax2, minimum(xvi[2])./1e3, 0)
hideydecorations!(ax2)
# save(joinpath(figdir, "initial_profile.png"), fig)
fig
end

# IO -------------------------------------------------
# if it does not exist, make folder where figures are stored
Expand Down Expand Up @@ -213,8 +213,8 @@ function main(li, origin, phases_GMG, igg; nx=16, ny=16, figdir="figs2D", do_vtk
args,
dt,
igg;
iterMax = 150e3,
nout = 2e3,
iterMax = 100e3,
nout = 1e3,
viscosity_cutoff = (1e18, 1e24),
free_surface = false,
# viscosity_relaxation = 1e-5
Expand Down Expand Up @@ -354,5 +354,5 @@ else
igg
end

main(li, origin, phases_GMG, igg; figdir = figdir, nx = nx, ny = ny, do_vtk = do_vtk);
# main(li, origin, phases_GMG, igg; figdir = figdir, nx = nx, ny = ny, do_vtk = do_vtk);

Loading

0 comments on commit 7afa625

Please sign in to comment.