diff --git a/Cavity2D.jl b/Cavity2D.jl new file mode 100644 index 00000000..d7e1d57a --- /dev/null +++ b/Cavity2D.jl @@ -0,0 +1,396 @@ +using Printf +using JustRelax, JustRelax.JustRelax2D +import JustRelax.JustRelax2D as JR +const backend_JR = CPUBackend + +using JustPIC, JustPIC._2D +const backend = JustPIC.CPUBackend + +using ParallelStencil, ParallelStencil.FiniteDifferences2D +@init_parallel_stencil(Threads, Float64, 2) + +# Load script dependencies +using LinearAlgebra, GeoParams, GLMakie + +include("mask.jl") +include("MiniKernels.jl") + +# Velocity helper grids for the particle advection +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) + esc(:($A[$idx_j])) +end + +# Initial pressure profile - not accurate +@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) + ni = size(phases) + + @parallel_indices (i, j) function init_phases!(phases, px, py, index) + r = 5 + + @inbounds for ip in cellaxes(phases) + # quick escape + @index(index[ip, i, j]) == 0 && continue + + x = @index px[ip, i, j] + depth = abs(@index py[ip, i, j]) + @index phases[ip, i, j] = 2.0 + + if 0e0 ≤ depth ≤ 10e0 || ((x - 25)^2 + (depth - 55)^2 ≤ r^2) + @index phases[ip, i, j] = 1.0 + end + + end + return nothing + end + + @parallel (@idx ni) init_phases!(phases, particles.coords..., particles.index) +end +## END OF HELPER FUNCTION ------------------------------------------------------------ + +# (Path)/folder where output data and figures are stored +n = 101 +nx = n +ny = n +igg = if !(JustRelax.MPI.Initialized()) # initialize (or not) MPI grid + IGG(init_global_grid(nx, ny, 1; init_MPI= true)...) +else + igg +end + +## BEGIN OF MAIN SCRIPT -------------------------------------------------------------- +function main(igg, nx, ny) + + # Physical domain ------------------------------------ + thick_air = 5e0 # thickness of sticky air layer + ly = 50e0 + thick_air # domain length in y + lx = 50e0 # domain length in x + ni = nx, ny # number of cells + li = lx, ly # domain length in x- and y- + di = @. li / ni # grid step in x- and -y + origin = 0.0, -ly # origin coordinates (15km f sticky air layer) + grid = Geometry(ni, li; origin = origin) + (; xci, xvi) = grid # nodes at the center and vertices of the cells + # ---------------------------------------------------- + + # Physical properties using GeoParams ---------------- + rheology = rheology = ( + # Name = "Air", + SetMaterialParams(; + Phase = 1, + Density = ConstantDensity(; ρ=2.7e3), + CompositeRheology = CompositeRheology((LinearViscous(; η=1e21),)), + Gravity = ConstantGravity(; g=9.81), + ), + # Name = "Rock", + SetMaterialParams(; + Phase = 2, + Density = ConstantDensity(; ρ=2.7e3), + CompositeRheology = CompositeRheology((LinearViscous(; η=1e21),)), + Gravity = ConstantGravity(; g=9.81), + ), + ) + # ---------------------------------------------------- + + # Initialize particles ------------------------------- + nxcell, max_xcell, min_xcell = 120, 120, 100 + particles = init_particles( + backend, nxcell, max_xcell, min_xcell, xvi, di, ni + ) + # velocity grids + grid_vx, grid_vy = velocity_grids(xci, xvi, di) + # temperature + pT, pPhases = init_cell_arrays(particles, Val(2)) + particle_args = (pT, pPhases) + + # Elliptical temperature anomaly + init_phases!(pPhases, particles) + phase_ratios = PhaseRatios(backend, length(rheology), ni) + update_phase_ratios!(phase_ratios, particles, xci, xvi, pPhases) + # ---------------------------------------------------- + + # RockRatios + air_phase = 1 + ϕ = RockRatio(ni...) + update_rock_ratio!(ϕ, phase_ratios, air_phase) + + # STOKES --------------------------------------------- + # Allocate arrays needed for every Stokes problem + stokes = StokesArrays(backend_JR, ni) + pt_stokes = PTStokesCoeffs(li, di; ϵ=1e-4, Re=3e0*π, r=0.7, CFL = 0.95 / √2.1) + # ---------------------------------------------------- + + # TEMPERATURE PROFILE -------------------------------- + thermal = ThermalArrays(backend_JR, ni) + # ---------------------------------------------------- + + # Buoyancy forces & rheology + ρg = @zeros(ni...), @zeros(ni...) + args = (; T = thermal.Tc, P = stokes.P, dt = Inf) + compute_ρg!(ρg[2], phase_ratios, rheology, (T=thermal.Tc, P=stokes.P)) + # @parallel init_P!(stokes.P, ρg[2], xci[2]) + compute_viscosity!(stokes, phase_ratios, args, rheology, (-Inf, Inf)) + + # Boundary conditions + flow_bcs = VelocityBoundaryConditions(; + free_slip = (left = true, right = true, top = true, bot = false), + no_slip = (left = false, right = false, top = false, bot = true), + ) + + Vx_v = @zeros(ni.+1...) + Vy_v = @zeros(ni.+1...) + + figdir = "FreeSurfacePlume" + take(figdir) + + # Time loop + t, it = 0.0, 0 + dt = 1e3 * (3600 * 24 * 365.25) + + _di = inv.(di) + (; ϵ, r, θ_dτ, ηdτ) = pt_stokes + (; η, η_vep) = stokes.viscosity + ni = size(stokes.P) + iterMax = 5e3 + nout = 1e3 + viscosity_cutoff = (-Inf, Inf) + free_surface = false + ητ = @zeros(ni...) + while it < 1 + + ## variational solver + + # errors + err = 2 * ϵ + iter = 0 + err_evo1 = Float64[] + err_evo2 = Float64[] + norm_Rx = Float64[] + norm_Ry = Float64[] + norm_∇V = Float64[] + sizehint!(norm_Rx, Int(iterMax)) + sizehint!(norm_Ry, Int(iterMax)) + sizehint!(norm_∇V, Int(iterMax)) + sizehint!(err_evo1, Int(iterMax)) + sizehint!(err_evo2, Int(iterMax)) + + # solver loop + @copy stokes.P0 stokes.P + wtime0 = 0.0 + relλ = 0.2 + θ = deepcopy(stokes.P) + λ = @zeros(ni...) + λv = @zeros(ni .+ 1...) + η0 = deepcopy(η) + do_visc = true + + for Aij in @tensor_center(stokes.ε_pl) + Aij .= 0.0 + end + # Vx_on_Vy = @zeros(size(stokes.V.Vy)) + + # compute buoyancy forces and viscosity + compute_ρg!(ρg[end], phase_ratios, rheology, args) + compute_viscosity!(stokes, phase_ratios, args, rheology, viscosity_cutoff) + + # for a in 1:100 + while iter ≤ iterMax + err < ϵ && break + # for _ in 1:100 + JR.compute_maxloc!(ητ, η; window=(1, 1)) + # update_halo!(ητ) + + @parallel (@idx ni) compute_∇V!(stokes.∇V, @velocity(stokes), ϕ, _di) + + @parallel (@idx ni) compute_P!( + θ, + stokes.P0, + stokes.R.RP, + stokes.∇V, + ητ, + rheology, + phase_ratios.center, + ϕ, + dt, + pt_stokes.r, + pt_stokes.θ_dτ + ) + + JR.update_ρg!(ρg[2], phase_ratios, rheology, args) + + @parallel (@idx ni .+ 1) JR.compute_strain_rate!( + @strain(stokes)..., stokes.∇V, @velocity(stokes)..., _di... + ) + + # if rem(iter, nout) == 0 + # @copy η0 η + # end + # if do_visc + # update_viscosity!( + # stokes, + # phase_ratios, + # args, + # rheology, + # viscosity_cutoff; + # relaxation=viscosity_relaxation, + # ) + # end + + @parallel (@idx ni .+ 1) update_stresses_center_vertex_ps!( + @strain(stokes), + @tensor_center(stokes.ε_pl), + stokes.EII_pl, + @tensor_center(stokes.τ), + (stokes.τ.xy,), + @tensor_center(stokes.τ_o), + (stokes.τ_o.xy,), + θ, + stokes.P, + stokes.viscosity.η, + λ, + λv, + stokes.τ.II, + stokes.viscosity.η_vep, + relλ, + dt, + θ_dτ, + rheology, + phase_ratios.center, + phase_ratios.vertex, + ϕ, + ) + # update_halo!(stokes.τ.xy) + + # @parallel (1:(size(stokes.V.Vy, 1) - 2), 1:size(stokes.V.Vy, 2)) JR.interp_Vx∂ρ∂x_on_Vy!( + # Vx_on_Vy, stokes.V.Vx, ρg[2], _di[1] + # ) + + # @hide_communication b_width begin # communication/computation overlap + @parallel (@idx ni.+1) compute_V!( + @velocity(stokes)..., + stokes.R.Rx, + stokes.R.Ry, + stokes.P, + @stress(stokes)..., + ηdτ, + ρg..., + ητ, + ϕ, + # ϕ.Vx, + # ϕ.Vy, + _di..., + ) + # apply boundary conditions + # velocity2displacement!(stokes, dt) + # JR.free_surface_bcs!(stokes, flow_bcs, η, rheology, phase_ratios, dt, di) + flow_bcs!(stokes, flow_bcs) + # end + # f,ax,h=heatmap(stokes.V.Vy) + # # f,ax,h=heatmap(stokes.V.Vx) + # Colorbar(f[1,2], h, label="Vy"); f + # update_halo!(@velocity(stokes)...) + # end + + iter += 1 + + if iter % nout == 0 && iter > 1 + # er_η = norm_mpi(@.(log10(η) - log10(η0))) + # er_η < 1e-3 && (do_visc = false) + # errs = maximum_mpi.((abs.(stokes.R.Rx), abs.(stokes.R.Ry), abs.(stokes.R.RP))) + errs = ( + norm(@views stokes.R.Rx[2:(end - 1), 2:(end - 1)]) / length(stokes.R.Rx), + norm(@views stokes.R.Ry[2:(end - 1), 2:(end - 1)]) / length(stokes.R.Ry), + norm(stokes.R.RP) / length(stokes.R.RP), + ) + push!(norm_Rx, errs[1]) + push!(norm_Ry, errs[2]) + push!(norm_∇V, errs[3]) + err = maximum(errs) + push!(err_evo1, err) + push!(err_evo2, iter) + + # if igg.me == 0 #&& ((verbose && err > ϵ) || iter == iterMax) + @printf( + "Total steps = %d, err = %1.3e [norm_Rx=%1.3e, norm_Ry=%1.3e, norm_∇V=%1.3e] \n", + iter, + err, + norm_Rx[end], + norm_Ry[end], + norm_∇V[end] + ) + # end + isnan(err) && error("NaN(s)") + isinf(err) && error("Inf(s)") + end + end + heatmap(stokes.V.Vy) + + dt = compute_dt(stokes, di) / 2 + # ------------------------------ + + # Advection -------------------- + # advect particles in space + advection!(particles, RungeKutta2(), @velocity(stokes), (grid_vx, grid_vy), dt) + # advect particles in memory + move_particles!(particles, xvi, particle_args) + # check if we need to inject particles + inject_particles_phase!(particles, pPhases, (), (), xvi) + # update phase ratios + update_phase_ratios!(phase_ratios, particles, xci, xvi, pPhases) + update_rock_ratio!(ϕ, phase_ratios, air_phase) + + @show it += 1 + t += dt + + if it == 1 || rem(it, 1) == 0 + velocity2vertex!(Vx_v, Vy_v, @velocity(stokes)...) + nt = 5 + fig = Figure(size = (900, 900), title = "t = $t") + ax = Axis(fig[1,1], aspect = 1, title = " t=$(round.(t/(1e3 * 3600 * 24 *365.25); digits=3)) Kyrs") + # heatmap!(ax, xci[1].*1e-3, xci[2].*1e-3, Array([argmax(p) for p in phase_ratios.vertex]), colormap = :grayC) + # arrows!( + # ax, + # xvi[1][1:nt:end-1]./1e3, xvi[2][1:nt:end-1]./1e3, Array.((Vx_v[1:nt:end-1, 1:nt:end-1], Vy_v[1:nt:end-1, 1:nt:end-1]))..., + # lengthscale = 25 / max(maximum(Vx_v), maximum(Vy_v)), + # color = :red, + # ) + # heatmap!(ax, stokes.V.Vy) + # heatmap!(ax, stokes.τ.xx) + + # ind = iszero.(stokes.V.Vy) + # stokes.V.Vy[ind] .= NaN + heatmap!(ax, stokes.V.Vy) + display(fig) + # save(joinpath(figdir, "$(it).png"), fig) + + end + end + return stokes, ϕ +end +# ## END OF MAIN SCRIPT ---------------------------------------------------------------- +stoke, ϕ = main(igg, nx, ny); + +# # # (Path)/folder where output data and figures are stored +# # n = 100 +# # nx = n +# # ny = n +# # igg = if !(JustRelax.MPI.Initialized()) # initialize (or not) MPI grid +# # IGG(init_global_grid(nx, ny, 1; init_MPI= true)...) +# # else +# # igg +# # end + diff --git a/Duretz2016.jl b/Duretz2016.jl index 9f7fab9a..674a7dcb 100644 --- a/Duretz2016.jl +++ b/Duretz2016.jl @@ -43,6 +43,7 @@ function init_phases!(phases, particles) @parallel_indices (i, j) function init_phases!(phases, px, py, index) r=0.5 + rinc=0.1 @inbounds for ip in cellaxes(phases) # quick escape @@ -51,6 +52,7 @@ function init_phases!(phases, particles) x = @index px[ip, i, j] depth = (@index py[ip, i, j]) + # h = √(r^2 - x^2) - 0.4330127018922193*2 h = -√(r^2 - x^2) if depth ≤ h @index phases[ip, i, j] = 2.0 @@ -58,6 +60,15 @@ function init_phases!(phases, particles) @index phases[ip, i, j] = 1.0 end + # h = √(r^2 - x^2) - 0.55*2 + # h = -√(r^2 - x^2) - 0.2 + # if depth ≤ h + # @index phases[ip, i, j] = 1.0 + # end + + # if (x^2 + (depth + 0.75)^2 ≤ rinc^2) + # @index phases[ip, i, j] = 1.0 + # end end return nothing @@ -65,6 +76,14 @@ function init_phases!(phases, particles) @parallel (@idx ni) init_phases!(phases, particles.coords..., particles.index) end + +# init_phases!(pPhases, particles) +# phase_ratios = PhaseRatios(backend, length(rheology), ni) +# update_phase_ratios!(phase_ratios, particles, xci, xvi, pPhases) +# p=[argmax(p) for p in phase_ratios.center] +# heatmap(p) + + ## END OF HELPER FUNCTION ------------------------------------------------------------ # (Path)/folder where output data and figures are stored @@ -78,7 +97,7 @@ else end ## BEGIN OF MAIN SCRIPT -------------------------------------------------------------- -# function main(igg, nx, ny) +function main(igg, nx, ny) # Physical domain ------------------------------------ thick_air = 0 # thickness of sticky air layer @@ -115,7 +134,7 @@ end # ---------------------------------------------------- # Initialize particles ------------------------------- - nxcell, max_xcell, min_xcell = 30, 40, 15 + nxcell, max_xcell, min_xcell = 40, 70, 15 particles = init_particles( backend, nxcell, max_xcell, min_xcell, xvi, di, ni ) @@ -173,7 +192,7 @@ end (; ϵ, r, θ_dτ, ηdτ) = pt_stokes (; η, η_vep) = stokes.viscosity ni = size(stokes.P) - iterMax = 15e3 + iterMax = 5e3 nout = 1e3 viscosity_cutoff = (-Inf, Inf) free_surface = false @@ -347,50 +366,73 @@ end dt = compute_dt(stokes, di) / 2 # ------------------------------ - - # Advection -------------------- - # advect particles in space - advection!(particles, RungeKutta2(), @velocity(stokes), (grid_vx, grid_vy), dt) - # advect particles in memory - move_particles!(particles, xvi, particle_args) - # check if we need to inject particles - inject_particles_phase!(particles, pPhases, (), (), xvi) - # update phase ratios - update_phase_ratios!(phase_ratios, particles, xci, xvi, pPhases) - update_rock_ratio!(ϕ, phase_ratios, air_phase) - - @show it += 1 - t += dt - - if it == 1 || rem(it, 1) == 0 - velocity2vertex!(Vx_v, Vy_v, @velocity(stokes)...) - nt = 5 - fig = Figure(size = (900, 900), title = "t = $t") - ax = Axis(fig[1,1], aspect = 1, title = " t=$(round.(t/(1e3 * 3600 * 24 *365.25); digits=3)) Kyrs") - heatmap!(ax, xci[1].*1e-3, xci[2].*1e-3, Array([argmax(p) for p in phase_ratios.vertex]), colormap = :grayC) - # arrows!( - # ax, - # xvi[1][1:nt:end-1]./1e3, xvi[2][1:nt:end-1]./1e3, Array.((Vx_v[1:nt:end-1, 1:nt:end-1], Vy_v[1:nt:end-1, 1:nt:end-1]))..., - # lengthscale = 25 / max(maximum(Vx_v), maximum(Vy_v)), - # color = :red, - # ) - fig - save(joinpath(figdir, "$(it).png"), fig) - - end - # end - return nothing -# end + heatmap(stokes.τ.xy) + + # # Advection -------------------- + # # advect particles in space + # advection!(particles, RungeKutta2(), @velocity(stokes), (grid_vx, grid_vy), dt) + # # advect particles in memory + # move_particles!(particles, xvi, particle_args) + # # check if we need to inject particles + # inject_particles_phase!(particles, pPhases, (), (), xvi) + # # update phase ratios + # update_phase_ratios!(phase_ratios, particles, xci, xvi, pPhases) + # update_rock_ratio!(ϕ, phase_ratios, air_phase) + + # @show it += 1 + # t += dt + + # if it == 1 || rem(it, 1) == 0 + # velocity2vertex!(Vx_v, Vy_v, @velocity(stokes)...) + # nt = 5 + # fig = Figure(size = (900, 900), title = "t = $t") + # ax = Axis(fig[1,1], aspect = 1, title = " t=$(round.(t/(1e3 * 3600 * 24 *365.25); digits=3)) Kyrs") + # heatmap!(ax, xci[1].*1e-3, xci[2].*1e-3, Array([argmax(p) for p in phase_ratios.vertex]), colormap = :grayC) + # # arrows!( + # # ax, + # # xvi[1][1:nt:end-1]./1e3, xvi[2][1:nt:end-1]./1e3, Array.((Vx_v[1:nt:end-1, 1:nt:end-1], Vy_v[1:nt:end-1, 1:nt:end-1]))..., + # # lengthscale = 25 / max(maximum(Vx_v), maximum(Vy_v)), + # # color = :red, + # # ) + # fig + # save(joinpath(figdir, "$(it).png"), fig) + + # end + # # end + # return nothing + # heatmap(stokes.τ.xy) + ind = iszero.(stokes.P) + stokes.P[ind] .= NaN + f,ax,=heatmap(xci...,stokes.P) + ind = iszero.(stokes.V.Vy) + stokes.V.Vy[ind] .= NaN + f,ax,=heatmap(xci[1],xvi[2], stokes.V.Vy[2:end-1,:]) + + px, py = particles.coords + scatter!( + px.data[:], py.data[:],color=pPhases.data[:], + markersize=3) + h = @. -√(0.5^2 - xci[1]^2) + lines!(ax,xci[1], h, color=:red) + xlims!(ax, -0.25, -0.15) + ylims!(ax, -0.5, -0.4) + f + # ind = iszero.(stokes.V.Vx) + # stokes.V.Vx[ind] .= NaN + # heatmap(stokes.V.Vx) + # stokes +end # ## END OF MAIN SCRIPT ---------------------------------------------------------------- main(igg, nx, ny) -# # # (Path)/folder where output data and figures are stored -# # n = 100 -# # nx = n -# # ny = n -# # igg = if !(JustRelax.MPI.Initialized()) # initialize (or not) MPI grid -# # IGG(init_global_grid(nx, ny, 1; init_MPI= true)...) -# # else -# # igg -# # end +# heatmap(stokes.τ.xy) +# # ind = iszero.(stokes.P) +# # stokes.P[ind] .= NaN +# heatmap(stokes.P) +# ind = iszero.(stokes.V.Vy) +# stokes.V.Vy[ind] .= NaN +# heatmap(stokes.V.Vy) +# # ind = iszero.(stokes.V.Vx) +# # stokes.V.Vx[ind] .= NaN +# # heatmap(stokes.V.Vx) \ No newline at end of file diff --git a/VariationalSolver.jl b/VariationalSolver.jl index f6c70ad1..b8441c29 100644 --- a/VariationalSolver.jl +++ b/VariationalSolver.jl @@ -270,17 +270,29 @@ Colorbar(f[1,2], h, label="Vy") lines(stokes.P[1,:]) lines(stokes.V.Vy[1,:]) +xc = [x for x in xci[1], y in xci[2]] +yc = [y for x in xci[1], y in xci[2]] + v = [isvalid_c(ϕ, i, j) for i in 1:nx, j in 1:ny] heatmap(xci..., v) +heatmap(xci..., ϕ.center) +h = @. -√(0.5^2 - xci[1]^2) +lines!(xci[1], h, color=:red) + v = [isvalid_v(ϕ, i, j) for i in 1:nx, j in 1:ny] heatmap(xvi..., v) +heatmap(xvi..., ϕ.vertex) +h = @. -√(0.5^2 - xvi[1]^2) +lines!(xvi[1], h, color=:red) -vx = [isvalid_vx(ϕ, i, j) for i in 1:nx+1, j in 1:ny] -heatmap(vx) +vx = [isvalid_vx(ϕ, i+1, j) for i in 1:nx-1, j in 1:ny] +heatmap(xvi[1], xci[2], vx) +h = @. -√(0.5^2 - xci[1]^2) +lines!(xci[1], h, color=:red) -vy = [isvalid_vy(ϕ, i, j) for i in 1:nx, j in 1:ny+1] +vy = [isvalid_vy(ϕ, i, j+1) for i in 1:nx, j in 1:ny-1] heatmap(vy) heatmap( diff --git a/mask.jl b/mask.jl index cf7d1a3e..edc040ac 100644 --- a/mask.jl +++ b/mask.jl @@ -122,55 +122,48 @@ end function isvalid_c(ϕ::RockRatio, i, j) vx = (ϕ.Vx[i, j] > 0) * (ϕ.Vx[i + 1, j] > 0) vy = (ϕ.Vy[i, j] > 0) * (ϕ.Vy[i, j + 1] > 0) - div = vx * vy - return div * (ϕ.center[i, j] > 0) + v = vx * vy + # return v * (ϕ.center[i, j] > 0) + return true end function isvalid_v(ϕ::RockRatio, i, j) nx, ny = size(ϕ.Vx) - il = clamp(i - 1, 1, nx) - i0 = clamp(i, 1, nx) - j0 = clamp(j, 1, ny) - vx = (ϕ.Vx[i0, j0] > 0) * (ϕ.Vx[il, j0] > 0) + j_bot = max(j - 1, 1) + j0 = min(j, ny) + vx = (ϕ.Vx[i, j0] > 0) * (ϕ.Vx[i, j_bot] > 0) - nx, ny = size(ϕ.Vy) - jl = max(j - 1, 1) - i0 = clamp(i, 1, nx) - j0 = clamp(j, 1, ny) - vy = (ϕ.Vy[i0, j0] > 0) * (ϕ.Vy[i0, jl] > 0) - div = vx * vy - return div * (ϕ.vertex[i, j] > 0) + nx, ny = size(ϕ.Vy) + i_left = max(i - 1, 1) + i0 = min(i, nx) + vy = (ϕ.Vy[i0, j] > 0) * (ϕ.Vy[i_left, j] > 0) + v = vx * vy + # return v * (ϕ.vertex[i, j] > 0) + return true end function isvalid_vx(ϕ::RockRatio, i, j) - nx, ny = size(ϕ.center) - il = clamp(i - 1, 1, nx) - i0 = clamp(i, 1, nx) - j0 = clamp(j, 1, ny) - c = (ϕ.center[i0, j0] > 0) * (ϕ.center[il, j0] > 0) - - _, ny = size(ϕ.vertex) - jt = clamp(j - 1, 1, ny) - v = (ϕ.vertex[i, j] > 0) * (ϕ.vertex[i, jt] > 0) - div = c * v - return div * (ϕ.Vx[i, j] > 0) + c = (ϕ.center[i, j] > 0) * (ϕ.center[i-1, j] > 0) + v = (ϕ.vertex[i, j] > 0) * (ϕ.vertex[i, j+1] > 0) + cv = c * v + return cv * (ϕ.Vx[i, j] > 0) + # c = (ϕ.center[i, j] > 0) || (ϕ.center[i-1, j] > 0) + # v = (ϕ.vertex[i, j] > 0) || (ϕ.vertex[i, j+1] > 0) + # cv = c || v + # return cv || (ϕ.Vx[i, j] > 0) end function isvalid_vy(ϕ::RockRatio, i, j) - nx, ny = size(ϕ.center) - jb = clamp(j - 1, 1, ny) - i0 = clamp(i, 1, nx) - j0 = clamp(j, 1, ny) - c = (ϕ.center[i0, j0] > 0) * (ϕ.center[i0, jb] > 0) - - nx, ny = size(ϕ.vertex) - ir = clamp(i - 1, 1, ny) - v = (ϕ.vertex[i, j] > 0) * (ϕ.vertex[ir, j] > 0) - div = c * v - return div * (ϕ.Vy[i, j] > 0) + c = (ϕ.center[i, j] > 0) * (ϕ.center[i, j - 1] > 0) + v = (ϕ.vertex[i, j] > 0) * (ϕ.vertex[i + 1, j] > 0) + cv = c * v + return cv * (ϕ.Vy[i, j] > 0) + # c = (ϕ.center[i, j] > 0) || (ϕ.center[i, j - 1] > 0) + # v = (ϕ.vertex[i, j] > 0) || (ϕ.vertex[i + 1, j] > 0) + # cv = c || v + # return cv || (ϕ.Vy[i, j] > 0) end - # function isvalid(A::T, I::Vararg{Integer, N}) where {N, T<:AbstractArray} # v = true # Base.@nexprs 2 j -> begin @@ -374,7 +367,7 @@ end harm_ya(A) = _av_ya(A, i, j) if all((i, j) .< size(Vx) .- 1) - if isvalid_vx(ϕ, i, j) + if isvalid_vx(ϕ, i + 1, j) Rx[i, j]= R_Vx = (-d_xa(P, ϕ.center) + d_xa(τxx, ϕ.center) + d_yi(τxy, ϕ.vertex) - av_xa(ρgx)) Vx[i + 1, j + 1] += R_Vx * ηdτ / av_xa(ητ) else @@ -384,7 +377,7 @@ end end if all((i, j) .< size(Vy) .- 1) - if isvalid_vy(ϕ, i, j) + if isvalid_vy(ϕ, i, j + 1) Ry[i, j] = R_Vy = -d_ya(P, ϕ.center) + d_ya(τyy, ϕ.center) + d_xi(τxy, ϕ.vertex) - av_ya(ρgy) Vy[i + 1, j + 1] += R_Vy * ηdτ / av_ya(ητ) else @@ -396,40 +389,40 @@ end return nothing end -# @parallel_indices (i, j) function compute_V!( -# Vx::AbstractArray{T,2}, Vy, Rx, Ry, P, τxx, τyy, τxy, ηdτ, ρgx, ρgy, ητ, ϕ_Vx, ϕ_Vy, _dx, _dy -# ) where {T} -# d_xi(A, ϕ) = _d_xi(A, ϕ, _dx, i, j) -# d_xa(A, ϕ) = _d_xa(A, ϕ, _dx, i, j) -# d_yi(A, ϕ) = _d_yi(A, ϕ, _dy, i, j) -# d_ya(A, ϕ) = _d_ya(A, ϕ, _dy, i, j) -# av_xa(A) = _av_xa(A, i, j) -# av_ya(A) = _av_ya(A, i, j) -# harm_xa(A) = _av_xa(A, i, j) -# harm_ya(A) = _av_ya(A, i, j) +@parallel_indices (i, j) function compute_V!( + Vx::AbstractArray{T,2}, Vy, Rx, Ry, P, τxx, τyy, τxy, ηdτ, ρgx, ρgy, ητ, ϕ_Vx, ϕ_Vy, _dx, _dy +) where {T} + d_xi(A, ϕ) = _d_xi(A, ϕ, _dx, i, j) + d_xa(A, ϕ) = _d_xa(A, ϕ, _dx, i, j) + d_yi(A, ϕ) = _d_yi(A, ϕ, _dy, i, j) + d_ya(A, ϕ) = _d_ya(A, ϕ, _dy, i, j) + av_xa(A) = _av_xa(A, i, j) + av_ya(A) = _av_ya(A, i, j) + harm_xa(A) = _av_xa(A, i, j) + harm_ya(A) = _av_ya(A, i, j) -# if all((i, j) .< size(Vx) .- 1) -# if iszero(ϕ_Vx[i + 1, j]) -# Rx[i, j] = zero(T) -# Vx[i + 1, j + 1] = zero(T) -# else -# Rx[i, j]= R_Vx = (-d_xa(P, ϕ.center) + d_xa(τxx, ϕ.center) + d_yi(τxy, ϕ.vertex) - av_xa(ρgx)) -# Vx[i + 1, j + 1] += R_Vx * ηdτ / av_xa(ητ) -# end -# end + if all((i, j) .< size(Vx) .- 1) + if iszero(ϕ_Vx[i + 1, j]) + Rx[i, j] = zero(T) + Vx[i + 1, j + 1] = zero(T) + else + Rx[i, j]= R_Vx = (-d_xa(P, ϕ.center) + d_xa(τxx, ϕ.center) + d_yi(τxy, ϕ.vertex) - av_xa(ρgx)) + Vx[i + 1, j + 1] += R_Vx * ηdτ / av_xa(ητ) + end + end -# if all((i, j) .< size(Vy) .- 1) -# if iszero(ϕ_Vy[i, j + 1]) -# Ry[i, j] = zero(T) -# Vy[i + 1, j + 1] = zero(T) -# else -# Ry[i, j] = R_Vy = -d_ya(P, ϕ.center) + d_ya(τyy, ϕ.center) + d_xi(τxy, ϕ.vertex) - av_ya(ρgy) -# Vy[i + 1, j + 1] += R_Vy * ηdτ / av_ya(ητ) -# end -# end + if all((i, j) .< size(Vy) .- 1) + if iszero(ϕ_Vy[i, j + 1]) + Ry[i, j] = zero(T) + Vy[i + 1, j + 1] = zero(T) + else + Ry[i, j] = R_Vy = -d_ya(P, ϕ.center) + d_ya(τyy, ϕ.center) + d_xi(τxy, ϕ.vertex) - av_ya(ρgy) + Vy[i + 1, j + 1] += R_Vy * ηdτ / av_ya(ητ) + end + end -# return nothing -# end + return nothing +end # @parallel_indices (i, j) function compute_V!( # Vx::AbstractArray{T,2}, Vy, Vx_on_Vy, P, τxx, τyy, τxy, ηdτ, ρgx, ρgy, ητ, ϕ_Vx, ϕ_Vy, _dx, _dy, dt