From bfda8d7a1368c56b2d6bc98cbc57545ef0cf87dc Mon Sep 17 00:00:00 2001 From: Albert de Montserrat <58044444+albert-de-montserrat@users.noreply.github.com> Date: Sat, 9 Nov 2024 20:12:18 +0100 Subject: [PATCH] Bump to JustPIC v0.5.3 (#272) * update stress kernels to new phase ratios from JP * bump JP * format * update 3D shear band miniapp * typos * fix 3D shear heating test * fix indexing bug in thermal BCs * fix indexing bug * make test pass again --- Project.toml | 6 +- .../stokes3D/shear_band/ShearBand3D.jl | 55 ++++++---- src/boundaryconditions/BoundaryConditions.jl | 6 +- src/stokes/Stokes2D.jl | 3 + src/stokes/Stokes3D.jl | 3 + src/stokes/StressKernels.jl | 102 ++++++------------ test/test_shearheating3D.jl | 30 +++--- test/test_thermalstresses.jl | 2 +- 8 files changed, 100 insertions(+), 107 deletions(-) diff --git a/Project.toml b/Project.toml index 8555574d6..e9f32bbd3 100644 --- a/Project.toml +++ b/Project.toml @@ -6,7 +6,9 @@ version = "0.4.0" [deps] Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" CellArrays = "d35fcfd7-7af4-4c67-b1aa-d78070614af4" +ExactFieldSolutions = "2a6b1ac7-2de5-4354-bffb-bf009f6a4bef" GeoParams = "e018b62d-d9de-4a26-8697-af89c310ae38" +GeophysicalModelGenerator = "3700c31b-fa53-48a6-808a-ef22d5a84742" HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f" ImplicitGlobalGrid = "4d7a3746-15be-11ea-1130-334b0c4f5fa0" JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" @@ -34,11 +36,13 @@ AMDGPU = "0.8, 0.9, 1" Adapt = "4" CUDA = "5" CellArrays = "0.2" +ExactFieldSolutions = "0.1.5" GeoParams = "0.6.4" +GeophysicalModelGenerator = "0.7.7" HDF5 = "0.17.1" ImplicitGlobalGrid = "0.15" JLD2 = "0.4, 0.5" -JustPIC = "0.5" +JustPIC = "0.5.3" MPI = "0.20" MuladdMacro = "0.2" ParallelStencil = "0.13.6" diff --git a/miniapps/benchmarks/stokes3D/shear_band/ShearBand3D.jl b/miniapps/benchmarks/stokes3D/shear_band/ShearBand3D.jl index ad5d511bf..b2dfbece5 100644 --- a/miniapps/benchmarks/stokes3D/shear_band/ShearBand3D.jl +++ b/miniapps/benchmarks/stokes3D/shear_band/ShearBand3D.jl @@ -1,39 +1,42 @@ +using CUDA using JustRelax, JustRelax.JustRelax3D, JustRelax.DataIO -const backend_JR = CPUBackend - +const backend_JR = CUDABackend using Printf, GeoParams, GLMakie, CellArrays using ParallelStencil -@init_parallel_stencil(Threads, Float64, 3) +@init_parallel_stencil(CUDA, Float64, 3) using JustPIC, JustPIC._3D -const backend = JustPIC.CPUBackend +const backend = CUDABackend # HELPER FUNCTIONS --------------------------------------------------------------- solution(ε, t, G, η) = 2 * ε * η * (1 - exp(-G * t / η)) # Initialize phases on the particles -function init_phases!(phase_ratios, xci, xvi, radius) - ni = size(phase_ratios.center) +function init_phases!(phases, particles, radius) + ni = size(phases) origin = 0.5, 0.5, 0.5 - @parallel_indices (i, j, k) function init_phases!(phases, xc, yc, zc, o_x, o_y, o_z) - x, y, z = xc[i], yc[j], zc[k] - if ((x-o_x)^2 + (y-o_y)^2 + (z-o_z)^2) > radius - @index phases[1, i, j, k] = 1.0 - @index phases[2, i, j, k] = 0.0 + @parallel_indices (I...) function init_phases!(phases, index, xc, yc, zc, o_x, o_y, o_z) + + for ip in cellaxes(xc) + (@index index[ip, I...]) || continue - else - @index phases[1, i, j, k] = 0.0 - @index phases[2, i, j, k] = 1.0 + x = @index xc[ip, I...] + y = @index yc[ip, I...] + z = @index zc[ip, I...] + if ((x-o_x)^2 + (y-o_y)^2 + (z-o_z)^2) > radius^2 + @index phases[ip, I...] = 1.0 + else + @index phases[ip, I...] = 2.0 + end end return nothing end - @parallel (@idx ni) init_phases!(phase_ratios.center, xci..., origin...) - @parallel (@idx ni .+ 1) init_phases!(phase_ratios.vertex, xvi..., origin...) + @parallel (@idx ni) init_phases!(phases, particles.index, particles.coords...,origin...) return nothing end @@ -91,9 +94,15 @@ function main(igg; nx=64, ny=64, nz=64, figdir="model_figs") ) # Initialize phase ratios ------------------------------- - radius = 0.1 + nxcell, max_xcell, min_xcell = 125, 150, 75 + particles = init_particles(backend, nxcell, max_xcell, min_xcell, xvi, di, ni) + radius = 0.1 + phase_ratios = PhaseRatios(backend, length(rheology), ni) + pPhases, = init_cell_arrays(particles, Val(1)) + # Assign particles phases anomaly + init_phases!(pPhases, particles, radius) phase_ratios = PhaseRatios(backend, length(rheology), ni) - init_phases!(phase_ratios, xci, xvi, radius) + update_phase_ratios!(phase_ratios, particles, xci, xvi, pPhases) # STOKES --------------------------------------------- # Allocate arrays needed for every Stokes problem @@ -167,9 +176,9 @@ function main(igg; nx=64, ny=64, nz=64, figdir="model_figs") velocity2vertex!(Vx_v, Vy_v, Vz_v, @velocity(stokes)...) data_v = (; - τII = Array(stokes.τ.II), - εII = Array(stokes.ε.II), - εII = Array(stokes.ε_pl.II), + τII = Array(stokes.τ.II), + εII = Array(stokes.ε.II), + εII_pl = Array(stokes.ε_pl.II), ) data_c = (; P = Array(stokes.P), @@ -214,7 +223,7 @@ function main(igg; nx=64, ny=64, nz=64, figdir="model_figs") return nothing end -n = 128 +n = 64 nx = ny = nz = n figdir = "ShearBand3D" igg = if !(JustRelax.MPI.Initialized()) @@ -222,4 +231,4 @@ igg = if !(JustRelax.MPI.Initialized()) else igg end -main(igg; figdir = figdir, nx = nx, ny = ny, nz = nz); +main(igg; figdir = figdir, nx = nx, ny = ny, nz = nz); \ No newline at end of file diff --git a/src/boundaryconditions/BoundaryConditions.jl b/src/boundaryconditions/BoundaryConditions.jl index c5006f7b2..766d84182 100644 --- a/src/boundaryconditions/BoundaryConditions.jl +++ b/src/boundaryconditions/BoundaryConditions.jl @@ -4,9 +4,11 @@ include("free_surface.jl") include("no_slip.jl") include("pure_shear.jl") -@inline bc_index(x::NTuple{2,T}) where {T} = mapreduce(xi -> max(size(xi)...), max, x) -@inline bc_index(x::T) where {T<:AbstractArray} = max(size(x)...) +@inline bc_index(x::T) where {T<:AbstractArray{_T,2} where {_T}} = max(size(x)...) +@inline bc_index(x::T) where {T<:AbstractArray{_T,3} where {_T}} = + max(size(x)...), max(size(x)...) +@inline bc_index(x::NTuple{2,T}) where {T} = mapreduce(xi -> max(size(xi)...), max, x) @inline function bc_index(x::NTuple{3,T}) where {T} n = mapreduce(xi -> max(size(xi)...), max, x) return n, n diff --git a/src/stokes/Stokes2D.jl b/src/stokes/Stokes2D.jl index ad2eb38ae..8faf8c31f 100644 --- a/src/stokes/Stokes2D.jl +++ b/src/stokes/Stokes2D.jl @@ -598,6 +598,9 @@ function _solve!( rheology, phase_ratios.center, phase_ratios.vertex, + phase_ratios.xy, + phase_ratios.yz, + phase_ratios.xz, ) update_halo!(stokes.τ.xy) diff --git a/src/stokes/Stokes3D.jl b/src/stokes/Stokes3D.jl index d4617339e..3f603bbe9 100644 --- a/src/stokes/Stokes3D.jl +++ b/src/stokes/Stokes3D.jl @@ -491,6 +491,9 @@ function _solve!( rheology, phase_ratios.center, phase_ratios.vertex, + phase_ratios.xy, + phase_ratios.yz, + phase_ratios.xz, ) update_halo!(stokes.τ.yz) update_halo!(stokes.τ.xz) diff --git a/src/stokes/StressKernels.jl b/src/stokes/StressKernels.jl index 63d6c5019..bf3c1bad1 100644 --- a/src/stokes/StressKernels.jl +++ b/src/stokes/StressKernels.jl @@ -587,6 +587,9 @@ end rheology, phase_center, phase_vertex, + phase_xy, + phase_yz, + phase_xz, ) τyzv, τxzv, τxyv = τshear_v τyzv_old, τxzv_old, τxyv_old = τshear_ov @@ -622,7 +625,7 @@ end τxyv_old_ij = av_clamped_yz_z(τxyv_old, Ic...) # vertex parameters - phase = @inbounds phase_vertex[I...] + phase = @inbounds phase_yz[I...] is_pl, Cv, sinϕv, cosϕv, sinψv, η_regv = plastic_params_phase( rheology, EIIv_ij, phase ) @@ -632,24 +635,12 @@ end dτ_rv = inv(θ_dτ + ηv_ij * _Gvdt + 1.0) # stress increments @ vertex - dτxxv = - (-(τxxv_ij - τxxv_old_ij) * ηv_ij * _Gvdt - τxxv_ij + 2.0 * ηv_ij * εxxv_ij) * - dτ_rv - dτyyv = - (-(τyyv_ij - τyyv_old_ij) * ηv_ij * _Gvdt - τyyv_ij + 2.0 * ηv_ij * εyyv_ij) * - dτ_rv - dτzzv = - (-(τzzv_ij - τzzv_old_ij) * ηv_ij * _Gvdt - τzzv_ij + 2.0 * ηv_ij * εzzv_ij) * - dτ_rv - dτyzv = - (-(τyzv_ij - τyzv_old_ij) * ηv_ij * _Gvdt - τyzv_ij + 2.0 * ηv_ij * εyzv_ij) * - dτ_rv - dτxzv = - (-(τxzv_ij - τxzv_old_ij) * ηv_ij * _Gvdt - τxzv_ij + 2.0 * ηv_ij * εxzv_ij) * - dτ_rv - dτxyv = - (-(τxyv_ij - τxyv_old_ij) * ηv_ij * _Gvdt - τxyv_ij + 2.0 * ηv_ij * εxyv_ij) * - dτ_rv + dτxxv = compute_stress_increment(τxxv_ij, τxxv_old_ij, ηv_ij, εxxv_ij, _Gvdt, dτ_rv) + dτyyv = compute_stress_increment(τyyv_ij, τyyv_old_ij, ηv_ij, εyyv_ij, _Gvdt, dτ_rv) + dτzzv = compute_stress_increment(τzzv_ij, τzzv_old_ij, ηv_ij, εzzv_ij, _Gvdt, dτ_rv) + dτyzv = compute_stress_increment(τyzv_ij, τyzv_old_ij, ηv_ij, εyzv_ij, _Gvdt, dτ_rv) + dτxzv = compute_stress_increment(τxzv_ij, τxzv_old_ij, ηv_ij, εxzv_ij, _Gvdt, dτ_rv) + dτxyv = compute_stress_increment(τxyv_ij, τxyv_old_ij, ηv_ij, εxyv_ij, _Gvdt, dτ_rv) dτijv = dτxxv, dτyyv, dτzzv, dτyzv, dτxzv, dτxyv τijv = τxxv_ij, τyyv_ij, τzzv_ij, τyzv_ij, τxzv_ij, τxyv_ij @@ -697,7 +688,7 @@ end τxyv_old_ij = av_clamped_xz_z(τxyv_old, Ic...) # vertex parameters - phase = @inbounds phase_vertex[I...] + phase = @inbounds phase_xz[I...] is_pl, Cv, sinϕv, cosϕv, sinψv, η_regv = plastic_params_phase( rheology, EIIv_ij, phase ) @@ -707,24 +698,12 @@ end dτ_rv = inv(θ_dτ + ηv_ij * _Gvdt + 1.0) # stress increments @ vertex - dτxxv = - (-(τxxv_ij - τxxv_old_ij) * ηv_ij * _Gvdt - τxxv_ij + 2.0 * ηv_ij * εxxv_ij) * - dτ_rv - dτyyv = - (-(τyyv_ij - τyyv_old_ij) * ηv_ij * _Gvdt - τyyv_ij + 2.0 * ηv_ij * εyyv_ij) * - dτ_rv - dτzzv = - (-(τzzv_ij - τzzv_old_ij) * ηv_ij * _Gvdt - τzzv_ij + 2.0 * ηv_ij * εzzv_ij) * - dτ_rv - dτyzv = - (-(τyzv_ij - τyzv_old_ij) * ηv_ij * _Gvdt - τyzv_ij + 2.0 * ηv_ij * εyzv_ij) * - dτ_rv - dτxzv = - (-(τxzv_ij - τxzv_old_ij) * ηv_ij * _Gvdt - τxzv_ij + 2.0 * ηv_ij * εxzv_ij) * - dτ_rv - dτxyv = - (-(τxyv_ij - τxyv_old_ij) * ηv_ij * _Gvdt - τxyv_ij + 2.0 * ηv_ij * εxyv_ij) * - dτ_rv + dτxxv = compute_stress_increment(τxxv_ij, τxxv_old_ij, ηv_ij, εxxv_ij, _Gvdt, dτ_rv) + dτyyv = compute_stress_increment(τyyv_ij, τyyv_old_ij, ηv_ij, εyyv_ij, _Gvdt, dτ_rv) + dτzzv = compute_stress_increment(τzzv_ij, τzzv_old_ij, ηv_ij, εzzv_ij, _Gvdt, dτ_rv) + dτyzv = compute_stress_increment(τyzv_ij, τyzv_old_ij, ηv_ij, εyzv_ij, _Gvdt, dτ_rv) + dτxzv = compute_stress_increment(τxzv_ij, τxzv_old_ij, ηv_ij, εxzv_ij, _Gvdt, dτ_rv) + dτxyv = compute_stress_increment(τxyv_ij, τxyv_old_ij, ηv_ij, εxyv_ij, _Gvdt, dτ_rv) dτijv = dτxxv, dτyyv, dτzzv, dτyzv, dτxzv, dτxyv τijv = τxxv_ij, τyyv_ij, τzzv_ij, τyzv_ij, τxzv_ij, τxyv_ij @@ -774,7 +753,7 @@ end τxyv_old_ij = τxyv_old[I...] # vertex parameters - phase = @inbounds phase_vertex[I...] + phase = @inbounds phase_xy[I...] is_pl, Cv, sinϕv, cosϕv, sinψv, η_regv = plastic_params_phase( rheology, EIIv_ij, phase ) @@ -784,25 +763,12 @@ end dτ_rv = inv(θ_dτ + ηv_ij * _Gvdt + 1.0) # stress increments @ vertex - dτxxv = - (-(τxxv_ij - τxxv_old_ij) * ηv_ij * _Gvdt - τxxv_ij + 2.0 * ηv_ij * εxxv_ij) * - dτ_rv - dτyyv = - (-(τyyv_ij - τyyv_old_ij) * ηv_ij * _Gvdt - τyyv_ij + 2.0 * ηv_ij * εyyv_ij) * - dτ_rv - dτzzv = - (-(τzzv_ij - τzzv_old_ij) * ηv_ij * _Gvdt - τzzv_ij + 2.0 * ηv_ij * εzzv_ij) * - dτ_rv - dτyzv = - (-(τyzv_ij - τyzv_old_ij) * ηv_ij * _Gvdt - τyzv_ij + 2.0 * ηv_ij * εyzv_ij) * - dτ_rv - dτxzv = - (-(τxzv_ij - τxzv_old_ij) * ηv_ij * _Gvdt - τxzv_ij + 2.0 * ηv_ij * εxzv_ij) * - dτ_rv - dτxyv = - (-(τxyv_ij - τxyv_old_ij) * ηv_ij * _Gvdt - τxyv_ij + 2.0 * ηv_ij * εxyv_ij) * - dτ_rv - + dτxxv = compute_stress_increment(τxxv_ij, τxxv_old_ij, ηv_ij, εxxv_ij, _Gvdt, dτ_rv) + dτyyv = compute_stress_increment(τyyv_ij, τyyv_old_ij, ηv_ij, εyyv_ij, _Gvdt, dτ_rv) + dτzzv = compute_stress_increment(τzzv_ij, τzzv_old_ij, ηv_ij, εzzv_ij, _Gvdt, dτ_rv) + dτyzv = compute_stress_increment(τyzv_ij, τyzv_old_ij, ηv_ij, εyzv_ij, _Gvdt, dτ_rv) + dτxzv = compute_stress_increment(τxzv_ij, τxzv_old_ij, ηv_ij, εxzv_ij, _Gvdt, dτ_rv) + dτxyv = compute_stress_increment(τxyv_ij, τxyv_old_ij, ηv_ij, εxyv_ij, _Gvdt, dτ_rv) dτijv = dτxxv, dτyyv, dτzzv, dτyzv, dτxzv, dτxyv τijv = τxxv_ij, τyyv_ij, τzzv_ij, τyzv_ij, τxzv_ij, τxyv_ij τIIv_ij = second_invariant(τijv .+ dτijv) @@ -895,6 +861,9 @@ end rheology, phase_center, phase_vertex, + phase_xy, + phase_yz, + phase_xz, ) τxyv = τshear_v[1] τxyv_old = τshear_ov[1] @@ -921,15 +890,11 @@ end dτ_rv = inv(θ_dτ + ηv_ij * _Gvdt + 1.0) # stress increments @ vertex - dτxxv = - (-(τxxv_ij - τxxv_old_ij) * ηv_ij * _Gvdt - τxxv_ij + 2.0 * ηv_ij * εxxv_ij) * dτ_rv - dτyyv = - (-(τyyv_ij - τyyv_old_ij) * ηv_ij * _Gvdt - τyyv_ij + 2.0 * ηv_ij * εyyv_ij) * dτ_rv - dτxyv = - ( - -(τxyv[I...] - τxyv_old[I...]) * ηv_ij * _Gvdt - τxyv[I...] + - 2.0 * ηv_ij * ε[3][I...] - ) * dτ_rv + dτxxv = compute_stress_increment(τxxv_ij, τxxv_old_ij, ηv_ij, εxxv_ij, _Gvdt, dτ_rv) + dτyyv = compute_stress_increment(τyyv_ij, τyyv_old_ij, ηv_ij, εyyv_ij, _Gvdt, dτ_rv) + dτxyv = compute_stress_increment( + τxyv[I...], τxyv_old[I...], ηv_ij, ε[3][I...], _Gvdt, dτ_rv + ) τIIv_ij = √(0.5 * ((τxxv_ij + dτxxv)^2 + (τyyv_ij + dτyyv)^2) + (τxyv[I...] + dτxyv)^2) # yield function @ center @@ -964,7 +929,8 @@ end εij_ve = @. εij + 0.5 * τij_o * _Gdt εII_ve = GeoParams.second_invariant(εij_ve) # stress increments @ center - dτij = @. (-(τij - τij_o) * ηij * _Gdt - τij .+ 2.0 * ηij * εij) * dτ_r + # dτij = @. (-(τij - τij_o) * ηij * _Gdt - τij .+ 2.0 * ηij * εij) * dτ_r + dτij = compute_stress_increment(τij, τij_o, ηij, εij, _Gdt, dτ_r) τII_ij = GeoParams.second_invariant(dτij .+ τij) # yield function @ center F = τII_ij - C - Pr[I...] * sinϕ diff --git a/test/test_shearheating3D.jl b/test/test_shearheating3D.jl index 68b29d583..939322e63 100644 --- a/test/test_shearheating3D.jl +++ b/test/test_shearheating3D.jl @@ -82,6 +82,7 @@ function Shearheating3D(igg; nx=16, ny=16, nz=16) # Initialize particles ------------------------------- nxcell, max_xcell, min_xcell = 20, 40, 10 particles = init_particles(backend, nxcell, max_xcell, min_xcell, xvi...) + subgrid_arrays = SubgridDiffusionCellArrays(particles) # velocity grids grid_vx, grid_vy, grid_vz = velocity_grids(xci, xvi, di) # temperature @@ -107,7 +108,8 @@ function Shearheating3D(igg; nx=16, ny=16, nz=16) # TEMPERATURE PROFILE -------------------------------- thermal = ThermalArrays(backend_JR, ni) thermal_bc = TemperatureBoundaryConditions(; - no_flux = (left = true , right = true , top = false, bot = false, front = true , back = true), + no_flux = (left = true , right = true , top = true, bot = true, front = true , back = true), + # no_flux = (left = true , right = true , top = false, bot = false, front = true , back = true), ) # Initialize constant temperature @@ -136,7 +138,7 @@ function Shearheating3D(igg; nx=16, ny=16, nz=16) no_slip = (left = false, right = false, top = false, bot = false, front = false, back = false), ) ## Compression and not extension - fix this - εbg = 5e-14 + εbg = 5e-14 stokes.V.Vx .= PTArray(backend_JR)([ -(x - lx/2) * εbg for x in xvi[1], _ in 1:ny+2, _ in 1:nz+2]) stokes.V.Vy .= PTArray(backend_JR)([ -(y - ly/2) * εbg for _ in 1:nx+2, y in xvi[2], _ in 1:nz+2]) stokes.V.Vz .= PTArray(backend_JR)([ (lz - abs(z)) * εbg for _ in 1:nx+2, _ in 1:ny+2, z in xvi[3]]) @@ -144,12 +146,17 @@ function Shearheating3D(igg; nx=16, ny=16, nz=16) update_halo!(@velocity(stokes)...) grid2particle!(pT, xvi, thermal.T, particles) + dt₀ = similar(stokes.P) # Time loop t, it = 0.0, 0 local iters while it < 5 + # interpolate fields from particle to grid vertices + particle2grid!(thermal.T, pT, xvi, particles) + temperature2center!(thermal) + # Stokes solver ---------------- iters = solve!( stokes, @@ -165,18 +172,14 @@ function Shearheating3D(igg; nx=16, ny=16, nz=16) kwargs = ( iterMax = 100e3, nout=1e3, - viscosity_cutoff=(-Inf, Inf), + viscosity_cutoff=(1e18, 1e22), verbose=false, ) ) tensor_invariant!(stokes.ε) - dt = compute_dt(stokes, di, dt_diff) + dt = compute_dt(stokes, di, dt_diff) * 0.1 # ------------------------------ - # interpolate fields from particle to grid vertices - particle2grid!(thermal.T, pT, xvi, particles) - temperature2center!(thermal) - compute_shear_heating!( thermal, stokes, @@ -202,6 +205,13 @@ function Shearheating3D(igg; nx=16, ny=16, nz=16) verbose = false, ) ) + 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 -------------------- @@ -211,8 +221,6 @@ function Shearheating3D(igg; nx=16, ny=16, nz=16) ) # advect particles in memory move_particles!(particles, xvi, particle_args) - # interpolate fields from grid vertices to particles - grid2particle_flip!(pT, xvi, thermal.T, thermal.Told, particles) # check if we need to inject particles inject_particles_phase!(particles, pPhases, (pT,), (thermal.T,), xvi) # update phase ratios @@ -223,8 +231,6 @@ function Shearheating3D(igg; nx=16, ny=16, nz=16) # ------------------------------ end - finalize_global_grid(; finalize_MPI=true) - return iters, thermal end diff --git a/test/test_thermalstresses.jl b/test/test_thermalstresses.jl index 9dc806e4c..56a9099b1 100644 --- a/test/test_thermalstresses.jl +++ b/test/test_thermalstresses.jl @@ -430,7 +430,7 @@ function main2D(igg; nx=32, ny=32) thermal_bcs!(thermal, thermal_bc) temperature2center!(thermal) thermal.ΔT .= thermal.T .- thermal.Told - vertex2center!(thermal.ΔTc, thermal.ΔT) + vertex2center!(thermal.ΔTc, thermal.ΔT[2:end-1, :]) @show it += 1 t += dt