From 1d3b1bd446d2111def415be2dd8728d2b8275372 Mon Sep 17 00:00:00 2001 From: Albert de Montserrat Date: Sun, 5 May 2024 23:34:08 +0200 Subject: [PATCH] first commit --- Subduction2D.jl | 396 +++++++++++++++++++ Subduction2D_GMG.jl | 119 ++++++ Subduction2D_rheology.jl | 182 +++++++++ src/boundaryconditions/BoundaryConditions.jl | 2 +- src/stokes/Stokes2D.jl | 151 +------ 5 files changed, 700 insertions(+), 150 deletions(-) create mode 100644 Subduction2D.jl create mode 100644 Subduction2D_GMG.jl create mode 100644 Subduction2D_rheology.jl diff --git a/Subduction2D.jl b/Subduction2D.jl new file mode 100644 index 00000000..ba320400 --- /dev/null +++ b/Subduction2D.jl @@ -0,0 +1,396 @@ +const isCUDA = false + +@static if isCUDA + using CUDA +end + +using JustRelax, JustRelax.JustRelax2D, JustRelax.DataIO +import JustRelax.@cell +using ParallelStencil, ParallelStencil.FiniteDifferences2D +@static if isCUDA + @init_parallel_stencil(CUDA, Float64, 2) +else + @init_parallel_stencil(Threads, Float64, 2) +end + +const backend = @static if isCUDA + CUDABackend # Options: CPUBackend, CUDABackend, AMDGPUBackend +else + CPUBackend # Options: CPUBackend, CUDABackend, AMDGPUBackend +end + +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_JP = @static if isCUDA + CUDABackend # Options: CPUBackend, CUDABackend, AMDGPUBackend +else + JustPIC.CPUBackend # Options: CPUBackend, CUDABackend, AMDGPUBackend +end + +# Load script dependencies +using GeoParams, GLMakie, CellArrays + +# Load file with all the rheology configurations +include("Subduction2D_rheology.jl") +include("Subduction2D_GMG.jl") + +## SET OF HELPER FUNCTIONS PARTICULAR FOR THIS SCRIPT -------------------------------- + +import ParallelStencil.INDICES +const idx_k = INDICES[2] +macro all_k(A) + esc(:($A[$idx_k])) +end + +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 + +# Initial pressure profile - not accurate +@parallel function init_P!(P, ρg, z) + @all(P) = abs(@all(ρg) * @all_k(z)) * <(@all_k(z), 0.0) + return nothing +end +## END OF HELPER FUNCTION ------------------------------------------------------------ + +## BEGIN OF MAIN SCRIPT -------------------------------------------------------------- +function main(li, origin, phases_GMG, igg; nx=16, ny=16, figdir="figs2D", do_vtk =false) + + # Physical domain ------------------------------------ + ni = nx, ny # number of cells + di = @. li / ni # grid steps + grid = Geometry(ni, li; origin = origin) + (; xci, xvi) = grid # nodes at the center and vertices of the cells + # ---------------------------------------------------- + + # Physical properties using GeoParams ---------------- + ρbg = 2.7e3 + rheology_dev = init_rheologies(; ρbg = ρbg) + rheology = init_rheologies(; ρbg = 0e0) + dt = 50e3 * 3600 * 24 * 365 # diffusive CFL timestep limiter + # ---------------------------------------------------- + + # Initialize particles ------------------------------- + nxcell = 40 + max_xcell = 60 + min_xcell = 20 + particles = init_particles( + backend_JP, nxcell, max_xcell, min_xcell, xvi, di, ni + ) + subgrid_arrays = SubgridDiffusionCellArrays(particles) + # velocity grids + grid_vxi = velocity_grids(xci, xvi, di) + # temperature + pPhases, pT = init_cell_arrays(particles, Val(2)) + # τxx_p, τyy_p, τxy_p = init_cell_arrays(particles, Val(3)) + # vorticity_p, = init_cell_arrays(particles, Val(1)) + # particle_args = (pT, τxx_p, τyy_p, τxy_p, vorticity_p, pPhases) + particle_args = (pT, pPhases) + + # Assign particles phases anomaly + phases_device = PTArray(backend)(phases_GMG) + init_phases!(pPhases, phases_device, particles, xvi) + phase_ratios = PhaseRatio(backend, ni, length(rheology)) + phase_ratios_center(phase_ratios, particles, grid, pPhases) + # ---------------------------------------------------- + + # STOKES --------------------------------------------- + # Allocate arrays needed for every Stokes problem + stokes = StokesArrays(backend, ni) + pt_stokes = PTStokesCoeffs(li, di; ϵ=1e-4, Re=3π, r=1e0, CFL = 1 / √2.1) # Re=3π, r=0.7 + # ---------------------------------------------------- + + # TEMPERATURE PROFILE -------------------------------- + thermal = ThermalArrays(backend, ni) + @views thermal.T[2:end-1, :] .= PTArray(backend)(T_GMG) + Ttop, Tbot = extrema(T_GMG) + thermal_bc = TemperatureBoundaryConditions(; + no_flux = (left = true, right = true, top = false, bot = false), + ) + thermal_bcs!(thermal, thermal_bc) + @views thermal.T[:, end] .= Ttop + @views thermal.T[:, 1] .= Tbot + temperature2center!(thermal) + # ---------------------------------------------------- + + # Buoyancy forces + args = (; T = thermal.Tc, P = zeros(ni...), Plitho = zeros(ni...), dt = Inf) + ρg = ntuple(_ -> @zeros(ni...), Val(2)) + # for _ in 1:2 + # compute_ρg!(ρg[2], phase_ratios, rheology, args) + # # JustRelax.Stokes2D.init_P!(stokes.P, ρg[2], xci[2]) + # end + compute_ρg!(ρg[2], phase_ratios, rheology_dev, args) + ρg_bg = ρbg * 9.81 + args.P .= reverse(cumsum(reverse((ρg[2] .+ ρg_bg).* di[2], dims=2), dims=2), dims=2) + + # Rheology + viscosity_cutoff = (1e18, 1e24) + compute_viscosity!(stokes, phase_ratios, args, rheology, viscosity_cutoff) + + # PT coefficients for thermal diffusion + pt_thermal = PTThermalCoeffs( + backend, rheology, phase_ratios, args, dt, ni, di, li; ϵ=1e-5, CFL=1e-3 / √3 + ) + + # Boundary conditions + flow_bcs = FlowBoundaryConditions(; + free_slip = (left = true , right = true , top = true , bot = true), + free_surface = false, + ) + 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(stokes.P[:]), Y./1e3) + # # scatter!(ax2, Array(Plitho[:]), 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 + if do_vtk + vtk_dir = joinpath(figdir, "vtk") + take(vtk_dir) + end + take(figdir) + # ---------------------------------------------------- + + local Vx_v, Vy_v + if do_vtk + Vx_v = @zeros(ni.+1...) + Vy_v = @zeros(ni.+1...) + end + + T_buffer = @zeros(ni.+1) + Told_buffer = similar(T_buffer) + dt₀ = similar(stokes.P) + for (dst, src) in zip((T_buffer, Told_buffer), (thermal.T, thermal.Told)) + copyinn_x!(dst, src) + end + grid2particle!(pT, xvi, T_buffer, particles) + + # Time loop + t, it = 0.0, 0 + + # Vertice arrays of normal components of the stress tensor + # τxx_vertex, τyy_vertex = @zeros(ni.+1...), @zeros(ni.+1...) + # τxx_o_vertex, τyy_o_vertex = @zeros(ni.+1...), @zeros(ni.+1...) + + while it < 1000 # run only for 5 Myrs + # while (t/(1e6 * 3600 * 24 *365.25)) < 5 # run only for 5 Myrs + + # interpolate fields from particle to grid vertices + particle2grid!(T_buffer, pT, xvi, particles) + @views T_buffer[:, end] .= Ttop + @views T_buffer[:, 1] .= Tbot + @views thermal.T[2:end-1, :] .= T_buffer + thermal_bcs!(thermal, thermal_bc) + temperature2center!(thermal) + + # Update buoyancy and viscosity - + args.Plitho .= reverse(cumsum(reverse((ρg[2] .+ ρg_bg).* di[2], dims=2), dims=2), dims=2) + args.P .= stokes.P .- minimum(stokes.P) .+ args.Plitho + + # args = (; T = thermal.Tc, P = Plitho, dt=Inf) + # args = (; T = thermal.Tc, P = stokes.P, dt=Inf) + compute_viscosity!(stokes, phase_ratios, args, rheology, viscosity_cutoff) + compute_ρg!(ρg[2], phase_ratios, rheology_dev, args) + + # Stokes solver ---------------- + t_stokes = @elapsed begin + out = solve!( + stokes, + pt_stokes, + di, + flow_bcs, + ρg, + phase_ratios, + rheology_dev, + args, + dt, + igg; + kwargs = ( + iterMax = 50e3, + nout = 1e3, + viscosity_cutoff = (1e18, 1e24), + free_surface = false, + viscosity_relaxation = 1e-2 + ) + ); + end + println("Stokes solver time ") + println(" Total time: $t_stokes s") + println(" Time/iteration: $(t_stokes / out.iter) s") + tensor_invariant!(stokes.ε) + dt = compute_dt(stokes, di) + # ------------------------------ + + # Thermal solver --------------- + heatdiffusion_PT!( + thermal, + pt_thermal, + thermal_bc, + rheology, + args, + dt, + di; + kwargs = ( + igg = igg, + phase = phase_ratios, + iterMax = 150e3, + nout = 2e3, + verbose = true + ), + ) + subgrid_characteristic_time!( + subgrid_arrays, particles, dt₀, phase_ratios, rheology, thermal, stokes, xci, di + ) + centroid2particle!(subgrid_arrays.dt₀, xci, dt₀, particles) + subgrid_diffusion!( + pT, thermal.T, thermal.ΔT, subgrid_arrays, particles, xvi, di, dt + ) + # ------------------------------ + + # Advection -------------------- + # advect particles in space + advection!(particles, RungeKutta2(), @velocity(stokes), grid_vxi, dt) + # advect particles in memory + move_particles!(particles, xvi, particle_args) + # JustRelax.Stokes2D.rotate_stress_particles!( + # stokes, + # τxx_vertex, + # τyy_vertex, + # τxx_o_vertex, + # τyy_o_vertex, + # τxx_p, + # τyy_p, + # τxy_p, + # vorticity_p, + # particles, + # grid, + # dt + # ) + + # check if we need to inject particles + # inject_particles_phase!(particles, pPhases, (pT, τxx_p, τyy_p, τxy_p, vorticity_p), (T_buffer, τxx_vertex, τyy_vertex, stokes.τ.xy, stokes.ω.xy_v), xvi) + inject_particles_phase!(particles, pPhases, (pT, ), (T_buffer,), xvi) + # update phase ratios + phase_ratios_center(phase_ratios, particles, grid, pPhases) + + @show it += 1 + t += dt + + # Data I/O and plotting --------------------- + if it == 1 || rem(it, 5) == 0 + # checkpointing(figdir, stokes, thermal.T, η, t) + (; η, η_vep) = stokes.viscosity + # if do_vtk + # velocity2vertex!(Vx_v, Vy_v, @velocity(stokes)...) + # data_v = (; + # T = Array(T_buffer), + # τxy = Array(stokes.τ.xy), + # εxy = Array(stokes.ε.xy), + # Vx = Array(Vx_v), + # Vy = Array(Vy_v), + # ) + # data_c = (; + # P = Array(stokes.P), + # τxx = Array(stokes.τ.xx), + # τyy = Array(stokes.τ.yy), + # εxx = Array(stokes.ε.xx), + # εyy = Array(stokes.ε.yy), + # η = Array(η_vep), + # ) + # velocity_v = ( + # Array(Vx_v), + # Array(Vy_v), + # ) + # save_vtk( + # joinpath(vtk_dir, "vtk_" * lpad("$it", 6, "0")), + # xvi, + # xci, + # data_v, + # data_c, + # velocity_v + # ) + # end + + # Make particles plottable + p = particles.coords + ppx, ppy = p + pxv = ppx.data[:]./1e3 + pyv = ppy.data[:]./1e3 + clr = pPhases.data[:] + # clr = pT.data[:] + idxv = particles.index.data[:]; + + # Make Makie figure + ar = 3 + fig = Figure(size = (1200, 900), title = "t = $t") + ax1 = Axis(fig[1,1], aspect = ar, title = "T [K] (t=$(t/(1e6 * 3600 * 24 *365.25)) Myrs)") + ax2 = Axis(fig[2,1], aspect = ar, title = "Vy [m/s]") + ax3 = Axis(fig[1,3], aspect = ar, title = "log10(εII)") + ax4 = Axis(fig[2,3], aspect = ar, title = "log10(η)") + # Plot temperature + h1 = heatmap!(ax1, xvi[1].*1e-3, xvi[2].*1e-3, Array(thermal.T[2:end-1,:]) , colormap=:batlow) + # Plot particles phase + h2 = scatter!(ax2, Array(pxv[idxv]), Array(pyv[idxv]), color=Array(clr[idxv]), markersize = 1) + # Plot 2nd invariant of strain rate + h3 = heatmap!(ax3, xci[1].*1e-3, xci[2].*1e-3, Array(log10.(stokes.ε.II)) , colormap=:batlow) + # Plot effective viscosity + h4 = heatmap!(ax4, xci[1].*1e-3, xci[2].*1e-3, Array(log10.(η_vep)) , colormap=:batlow) + hidexdecorations!(ax1) + hidexdecorations!(ax2) + hidexdecorations!(ax3) + Colorbar(fig[1,2], h1) + Colorbar(fig[2,2], h2) + Colorbar(fig[1,4], h3) + Colorbar(fig[2,4], h4) + linkaxes!(ax1, ax2, ax3, ax4) + fig + save(joinpath(figdir, "$(it).png"), fig) + end + # ------------------------------ + end + + return nothing +end + +## END OF MAIN SCRIPT ---------------------------------------------------------------- +do_vtk = true # set to true to generate VTK files for ParaView +figdir = "Subduction_LAMEM_2D" +# nx, ny = 512, 256 +# nx, ny = 512, 128 +n = 64 +nx, ny = n*6, n +# nx, ny = (512, 128) +nx, ny = 256, 64 +li, origin, phases_GMG, T_GMG = GMG_subduction_2D(nx+1, ny+1) +igg = if !(JustRelax.MPI.Initialized()) # initialize (or not) MPI grid + IGG(init_global_grid(nx, ny, 1; init_MPI= true)...) +else + igg +end + +# main(li, origin, phases_GMG, igg; figdir = figdir, nx = nx, ny = ny, do_vtk = do_vtk); diff --git a/Subduction2D_GMG.jl b/Subduction2D_GMG.jl new file mode 100644 index 00000000..8260d0be --- /dev/null +++ b/Subduction2D_GMG.jl @@ -0,0 +1,119 @@ +using GeophysicalModelGenerator + +function GMG_subduction_2D(nx, ny) + model_depth = 660 + # Our starting basis is the example above with ridge and overriding slab + nx, nz = nx, ny + x = range(-2000, 2000, nx); + z = range(-model_depth, 20, nz); + Grid2D = CartData(xyz_grid(x,0,z)) + Phases = zeros(Int64, nx, 1, nz); + Temp = fill(1280.0, nx, 1, nz); + air_thickness = 20.0 + lith = LithosphericPhases(Layers=[20 80], Phases=[1 2 0]) + + # Add left oceanic plate + add_box!( + Phases, + Temp, + Grid2D; + xlim =(-2000, 0), + zlim =(-model_depth, 0.0), + Origin = nothing, StrikeAngle=0, DipAngle=0, + phase = lith, + T = SpreadingRateTemp( + Tsurface = 20, + Tmantle = 1280.0, + MORside = "left", + SpreadingVel= 0.5, + AgeRidge = 0.01; + maxAge = 80.0 + ) + ) + + # Add right oceanic plate + add_box!( + Phases, + Temp, + Grid2D; + xlim =(1500, 2000), + zlim =(-model_depth, 0.0), + Origin = nothing, StrikeAngle=0, DipAngle=0, + phase = lith, + T = SpreadingRateTemp( + Tsurface = 20, + Tmantle = 1280.0, + MORside = "right", + SpreadingVel= 0.5, + AgeRidge = 0.01; + maxAge = 80.0 + ) + ) + + # Add overriding plate margin + add_box!( + Phases, + Temp, + Grid2D; + xlim =(0, 400), + zlim =(-model_depth, 0.0), + Origin = nothing, StrikeAngle=0, DipAngle=0, + phase = LithosphericPhases(Layers=[25 90], Phases=[3 4 0] ), + T = HalfspaceCoolingTemp( + Tsurface = 20, + Tmantle = 1280.0, + Age = 80.0 + ) + ) + + # Add overriding plate craton + add_box!( + Phases, + Temp, + Grid2D; + xlim =(400, 1500), + zlim =(-model_depth, 0.0), + Origin = nothing, StrikeAngle=0, DipAngle=0, + phase = LithosphericPhases(Layers=[35 100], Phases=[3 4 0] ), + T = HalfspaceCoolingTemp( + Tsurface = 20, + Tmantle = 1280.0, + Age = 120.0 + ) + ) + # Add slab + add_box!( + Phases, + Temp, + Grid2D; + xlim =(0, 300), + zlim =(-model_depth, 0.0), + Origin = nothing, StrikeAngle=0, DipAngle=30, + phase = LithosphericPhases(Layers=[30 80], Phases=[1 2 0], Tlab=1250 ), + T = HalfspaceCoolingTemp( + Tsurface = 20, + Tmantle = 1280.0, + Age = 120.0 + ) + ) + + Adiabat = 0.4 + @. Temp = Temp - Grid2D.z.val .* Adiabat + + # Lithosphere-asthenosphere boundary: + # ind = findall(Temp .> 1250 .&& (Phases.==2 .|| Phases.==5)); + # Phases[ind] .= 0; + + surf = Grid2D.z.val .> 0.0 + Temp[surf] .= 20.0 + Phases[surf] .= 5 + + Grid2D = addfield(Grid2D,(;Phases, Temp)) + + li = (abs(last(x)-first(x)), abs(last(z)-first(z))).* 1e3 + origin = (x[1], z[1]) .* 1e3 + + ph = Phases[:,1,:] .+ 1 + + return li, origin, ph, Temp[:,1,:].+273 +end diff --git a/Subduction2D_rheology.jl b/Subduction2D_rheology.jl new file mode 100644 index 00000000..3d23e260 --- /dev/null +++ b/Subduction2D_rheology.jl @@ -0,0 +1,182 @@ +using GeoParams.Dislocation +using GeoParams.Diffusion + +function init_rheologies(; ρbg = 0e0) + disl_dry_olivine = SetDislocationCreep(Dislocation.dry_olivine_Hirth_2003; V = 14.5e-6) + disl_oceanic_crust = SetDislocationCreep(Dislocation.plagioclase_An75_Ji_1993) + # disl_oceanic_litho = SetDislocationCreep(Dislocation.plagioclase_An75_Ji_1993) + disl_cont_crust = SetDislocationCreep(Dislocation.wet_quartzite_Kirby_1983) + + diff_dry_olivine = SetDiffusionCreep(Diffusion.dry_olivine_Hirth_2003; V = 14.5e-6) + + ϕ_dry_olivine = sind(20) + C_dry_olivine = 30e6 + + ϕ_oceanic_crust = sind(0) + C_oceanic_crust = 5e6 + + ϕ_oceanic_litho = sind(10) + C_oceanic_litho = 5e6 + + ϕ_cont_crust = sind(20) + C_cont_crust = 30e6 + + soft_C = LinearSoftening((C_oceanic_litho*0.05, C_oceanic_litho), (0.1, 0.5)) + + elasticity = ConstantElasticity(; G=5e10, ν=0.4) + # common physical properties + α = 3e-5 # 1 / K + Cp = 1000 # J / kg K + # C = 3e6 # Pa + η_reg = 1e20 + # ρbg = 2700 # kg / m^3 + + # Define rheolgy struct + rheology = ( + # Name = "dry olivine - Hirth_Kohlstedt_2003", + SetMaterialParams(; + Phase = 1, + Density = PT_Density(; ρ0=3.3e3-ρbg, α = α, β = 0e0, T0 = 273), + HeatCapacity = ConstantHeatCapacity(; Cp=Cp), + Conductivity = ConstantConductivity(; k = 3), + CompositeRheology = CompositeRheology( + ( + disl_dry_olivine, + diff_dry_olivine, + elasticity, + DruckerPrager_regularised(; C = C_dry_olivine, ϕ=ϕ_dry_olivine, η_vp=η_reg, Ψ=0.0, softening_C = soft_C) # non-regularized plasticity + ) + ), + RadioactiveHeat = ConstantRadioactiveHeat(6.6667e-12), + Elasticity = elasticity, + Gravity = ConstantGravity(; g=9.81), + ), + # Name = "oceanic crust", + SetMaterialParams(; + Phase = 2, + Density = PT_Density(; ρ0=3.3e3-ρbg, α = α, β = 0e0, T0 = 273), + HeatCapacity = ConstantHeatCapacity(; Cp=Cp), + Conductivity = ConstantConductivity(; k =3 ), + CompositeRheology = CompositeRheology( + ( + disl_oceanic_crust, + elasticity, + DruckerPrager_regularised(; C = C_oceanic_crust, ϕ = ϕ_oceanic_crust, η_vp=η_reg, Ψ=0.0, softening_C = soft_C) # non-regularized plasticity + ) + ), + RadioactiveHeat = ConstantRadioactiveHeat(2.333e-10), + ), + # Name = "oceanic lithosphere", + SetMaterialParams(; + Phase = 3, + Density = PT_Density(; ρ0=3.3e3-ρbg, α = α, β = 0e0, T0 = 273), + HeatCapacity = ConstantHeatCapacity(; Cp=Cp), + Conductivity = ConstantConductivity(; k = 3), + CompositeRheology = CompositeRheology( + ( + disl_dry_olivine, + diff_dry_olivine, + elasticity, + DruckerPrager_regularised(; C = C_dry_olivine, ϕ=ϕ_dry_olivine, η_vp=η_reg, Ψ=0.0, softening_C = soft_C) # non-regularized plasticity + ) + ), + RadioactiveHeat = ConstantRadioactiveHeat(6.6667e-12), + Elasticity = elasticity, + ), + # Name = "continental crust", + SetMaterialParams(; + Phase = 4, + Density = PT_Density(; ρ0=2.7e3-ρbg, α = α, β = 0e0, T0 = 273), + RadioactiveHeat = ConstantRadioactiveHeat(5.3571e-10), + HeatCapacity = ConstantHeatCapacity(; Cp=Cp), + Conductivity = ConstantConductivity(; k =3 ), + CompositeRheology = CompositeRheology( + ( + disl_cont_crust, + elasticity, + DruckerPrager_regularised(; C = C_cont_crust, ϕ = ϕ_cont_crust, η_vp=η_reg, Ψ=0.0, softening_C = soft_C) # non-regularized plasticity + ) + ), + Elasticity = elasticity, + ), + # Name = "continental lithosphere", + SetMaterialParams(; + Phase = 5, + Density = PT_Density(; ρ0=3.3e3-ρbg, α = α, β = 0e0, T0 = 273), + HeatCapacity = ConstantHeatCapacity(; Cp=Cp), + Conductivity = ConstantConductivity(; k = 3), + CompositeRheology = CompositeRheology( + ( + disl_dry_olivine, + diff_dry_olivine, + elasticity, + DruckerPrager_regularised(; C = C_dry_olivine, ϕ=ϕ_dry_olivine, η_vp=η_reg, Ψ=0.0, softening_C = soft_C) # non-regularized plasticity + ) + ), + RadioactiveHeat = ConstantRadioactiveHeat(6.6667e-12), + Elasticity = elasticity, + ), + # Name = "StickyAir", + SetMaterialParams(; + Phase = 6, + Density = ConstantDensity(; ρ=1-ρbg), # water density + HeatCapacity = ConstantHeatCapacity(; Cp=3e3), + RadioactiveHeat = ConstantRadioactiveHeat(0.0), + Conductivity = ConstantConductivity(; k=3.0), + CompositeRheology = CompositeRheology( + ( + LinearViscous(; η=1e21), + # elasticity + ) + ), + # Elasticity = elasticity, + ), + ) +end + +function init_phases!(phases, phase_grid, particles, xvi) + ni = size(phases) + @parallel (@idx ni) _init_phases!(phases, phase_grid, particles.coords, particles.index, xvi) +end + +@parallel_indices (I...) function _init_phases!(phases, phase_grid, pcoords::NTuple{N, T}, index, xvi) where {N,T} + + ni = size(phases) + + for ip in JustRelax.cellaxes(phases) + # quick escape + @cell(index[ip, I...]) == 0 && continue + + pᵢ = ntuple(Val(N)) do i + @cell pcoords[i][ip, I...] + end + + d = Inf # distance to the nearest particle + particle_phase = -1 + for offi in 0:1, offj in 0:1 + ii = I[1] + offi + jj = I[2] + offj + + !(ii ≤ ni[1]) && continue + !(jj ≤ ni[2]) && continue + + xvᵢ = ( + xvi[1][ii], + xvi[2][jj], + ) + d_ijk = √(sum((pᵢ[i] - xvᵢ[i])^2 for i in 1:N)) + if d_ijk < d + d = d_ijk + particle_phase = phase_grid[ii, jj] + end + end + @cell phases[ip, I...] = Float64(particle_phase) + + if pᵢ[2] ≥ 0.0 + @cell phases[ip, I...] = 6.0 + end + + end + + return nothing +end \ No newline at end of file diff --git a/src/boundaryconditions/BoundaryConditions.jl b/src/boundaryconditions/BoundaryConditions.jl index bdc9efeb..cdca024d 100644 --- a/src/boundaryconditions/BoundaryConditions.jl +++ b/src/boundaryconditions/BoundaryConditions.jl @@ -20,7 +20,7 @@ Apply the prescribed heat boundary conditions `bc` on the `T` """ thermal_bcs!(thermal, bcs) = thermal_bcs!(backend(thermal), thermal, bcs) function thermal_bcs!( - ::CPUBackendTrait, thermal::JustRelax.ThermalArrays, bcs::FlowBoundaryConditions + ::CPUBackendTrait, thermal::JustRelax.ThermalArrays, bcs::TemperatureBoundaryConditions ) return thermal_bcs!(thermal.T, bcs) end diff --git a/src/stokes/Stokes2D.jl b/src/stokes/Stokes2D.jl index b792b520..bdc80cc9 100644 --- a/src/stokes/Stokes2D.jl +++ b/src/stokes/Stokes2D.jl @@ -558,7 +558,7 @@ function _solve!( @strain(stokes), @tensor_center(stokes.ε_pl), stokes.EII_pl, - stokes.P, + args.P, θ, η, η_vep, @@ -568,7 +568,6 @@ function _solve!( dt, θ_dτ, ) - # free_surface_bcs!(stokes, flow_bcs) center2vertex!(stokes.τ.xy, stokes.τ.xy_c) update_halo!(stokes.τ.xy) @@ -643,7 +642,7 @@ function _solve!( end end - stokes.P .= θ + # stokes.P .= θ @parallel (@idx ni .+ 1) multi_copy!(@tensor(stokes.τ_o), @tensor(stokes.τ)) @parallel (@idx ni) multi_copy!(@tensor_center(stokes.τ_o), @tensor_center(stokes.τ)) @@ -660,149 +659,3 @@ function _solve!( norm_∇V=norm_∇V, ) end - -function _solve!( - stokes::JustRelax.StokesArrays, - thermal::JustRelax.ThermalArrays, - pt_stokes, - di::NTuple{2,T}, - flow_bcs, - ϕ, - ρg, - phase_v, - phase_c, - args_η, - rheology::NTuple{N,MaterialParams}, - dt, - igg::IGG; - iterMax=10e3, - nout=500, - b_width=(4, 4, 1), - verbose=true, - kwargs..., -) where {N,T} - - # unpack - - _di = inv.(di) - (; ϵ, r, θ_dτ, ηdτ) = pt_stokes - (; η, η_vep) = stokes.viscosity - ni = size(stokes.P) - # ~preconditioner - ητ = deepcopy(η) - # @hide_communication b_width begin # communication/computation overlap - compute_maxloc!(ητ, η; window=(1, 1)) - update_halo!(ητ) - # end - - # errors - err = 2 * ϵ - iter = 0 - err_evo1 = Float64[] - err_evo2 = Float64[] - norm_Rx = Float64[] - norm_Ry = Float64[] - norm_∇V = Float64[] - - # solver loop - wtime0 = 0.0 - while iter < 2 || (err > ϵ && iter ≤ iterMax) - wtime0 += @elapsed begin - @parallel (@idx ni) compute_∇V!(stokes.∇V, @velocity(stokes)..., _di...) - @parallel (@idx ni) compute_P!( - stokes.P, - stokes.P0, - stokes.R.RP, - stokes.∇V, - η, - rheology, - phase_c, - dt, - r, - θ_dτ, - ) - @parallel (@idx ni .+ 1) compute_strain_rate!( - @strain(stokes)..., stokes.∇V, @velocity(stokes)..., _di... - ) - @parallel (@idx ni) compute_ρg!( - ρg[end], ϕ, rheology, (T=thermal.Tc, P=stokes.P) - ) - @parallel (@idx ni) compute_τ_gp!( - @tensor_center(stokes.τ), - stokes.τ.II, - @tensor(stokes.τ_o), - @strain(stokes), - η, - η_vep, - thermal.T, - phase_v, - phase_c, - args_η, - rheology, # needs to be a tuple - dt, - θ_dτ, - ) - @parallel center2vertex!(stokes.τ.xy, stokes.τ.xy_c) - @hide_communication b_width begin # communication/computation overlap - @parallel compute_V!( - @velocity(stokes)..., - stokes.P, - @stress(stokes)..., - pt_stokes.ηdτ, - ρg..., - ητ, - _di..., - dt, - ) - # apply boundary conditions boundary conditions - flow_bcs!(stokes, flow_bcs) - update_halo!(stokes.V.Vx, stokes.V.Vy) - end - end - - iter += 1 - if iter % nout == 0 && iter > 1 - @parallel (@idx ni) compute_Res!( - stokes.R.Rx, - stokes.R.Ry, - @velocity(stokes)..., - stokes.P, - @stress(stokes)..., - ρg..., - _di..., - dt, - ) - errs = maximum_mpi.((abs.(stokes.R.Rx), abs.(stokes.R.Ry), abs.(stokes.R.RP))) - push!(norm_Rx, errs[1]) - push!(norm_Ry, errs[2]) - push!(norm_∇V, errs[3]) - err = maximum_mpi(errs) - push!(err_evo1, err) - push!(err_evo2, iter) - if igg.me == 0 && (verbose || 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)") - end - - if igg.me == 0 && err ≤ ϵ - println("Pseudo-transient iterations converged in $iter iterations") - end - end - - return ( - iter=iter, - err_evo1=err_evo1, - err_evo2=err_evo2, - norm_Rx=norm_Rx, - norm_Ry=norm_Ry, - norm_∇V=norm_∇V, - ) -end