2021-10-08 Finite Difference in 2D

Last time

  • Accuracy vs cost tradeoffs for ODE integrators

  • Stiff and non-stiff problems, discussion of costs

Today

  • FD methods in 2D

  • Cost profile

  • The need for fast algebraic solvers

using Plots
using LinearAlgebra
using SparseArrays
using Zygote

default(linewidth=4)

function newton(residual, jacobian, u0; maxits=20)
    u = u0
    uhist = [copy(u)]
    normhist = []
    for k in 1:maxits
        f = residual(u)
        push!(normhist, norm(f))
        J = jacobian(u)
        delta_u = - J \ f
        u += delta_u
        push!(uhist, copy(u))
    end
    uhist, normhist
end

function plot_stability(Rz, title; xlims=(-2, 2), ylims=(-2, 2))
    x = LinRange(xlims[1], xlims[2], 100)
    y = LinRange(ylims[1], ylims[2], 100)
    heatmap(x, y, (x, y) -> abs(Rz(x + 1im*y)), c=:bwr, clims=(0, 2), aspect_ratio=:equal, title=title)
end
plot_stability (generic function with 1 method)

Extending advection-diffusion to 2D

1 dimension

(19)\[\begin{align} u_t + (- \kappa u_x + w u)_x &= f(x) & \text{ on } \Omega &= (a,b) \\ u(a) &= g_D(a) & u'(b) &= g_N(b) \end{align}\]
  • Cell Peclet number \(\mathrm{Pe}_h = \frac{\lvert w \rvert h}{\kappa}\)

    • \(\mathrm{Pe}_h \lesssim 1\) avoids oscillations

    • \(\mathrm{Pe}_h \gtrsim 1\) is non-stiff for time-dependent model

  • Centered versus upwind for advection

  • Need uniformly bounded \(\kappa \ge \epsilon > 0\)

  • “Strong form” not defined at discontinuities in \(\kappa\)

    • Works okay using divergence form and fluxes at staggered points

2 dimensions

(20)\[\begin{align} u_t + \nabla\cdot\big(- \kappa \nabla u + \mathbf{w} u \big) &= f(x,y) & \text{ on } \Omega & \subset \mathbb R^2 \\ u|_{\Gamma_D} &= g_D(x,y) \\ (-\kappa \nabla u + \mathbf w u) \cdot \hat n|_{\Gamma_N} &= g_N(x,y) \end{align}\]
  • \(\Omega\) is some well-connected open set (we will assume simply connected) and the Dirichlet boundary \(\Gamma_D \subset \partial \Omega\) is nonempty.

  • Finite difference methods don’t have an elegant/flexible way to specify boundaries

  • We’ll choose \(\Omega = (-1, 1) \times (-1, 1)\)

On finite difference grids

  • Non-uniform grids can mesh “special” domains

    • Rare in 3D; overset grids, immersed boundary methods

  • Concept of staggering is complicated/ambiguous

Wesseling 11.4: A boundary-fitted grid around an airfoil.

Time-independent advection-diffusion

Advection

\[ \nabla\cdot(\mathbf w u) = \mathbf w \cdot \nabla u + (\nabla\cdot \mathbf w) u\]

If we choose divergence-free \(\mathbf w\), we can use the stencil

\[\begin{split} \mathbf w \cdot \nabla u \approx \frac 1 h \begin{bmatrix} & w_2 & \\ -w_1 & & w_1 \\ & -w_2 & \end{bmatrix} \!:\! \begin{bmatrix} u_{i-1, j+1} & u_{i, j+1} & u_{i+1,j+1} \\ u_{i-1, j} & u_{i, j} & u_{i+1,j} \\ u_{i-1, j-1} & u_{i, j-1} & u_{i+1,j-1} \\ \end{bmatrix} \end{split}\]

Diffusion

\[ -\nabla\cdot(\kappa \nabla u) = -\kappa \nabla\cdot \nabla u - \nabla\kappa\cdot \nabla u\]
  • When would you trust this decomposition?

  • If we have constant \(\kappa\), we can write

    \[\begin{split} -\kappa \nabla\cdot \nabla u \approx \frac{\kappa}{h^2} \begin{bmatrix} & -1 & \\ -1 & 4 & -1 \\ & -1 & \end{bmatrix} \!:\! \begin{bmatrix} u_{i-1, j+1} & u_{i, j+1} & u_{i+1,j+1} \\ u_{i-1, j} & u_{i, j} & u_{i+1,j} \\ u_{i-1, j-1} & u_{i, j-1} & u_{i+1,j-1} \\ \end{bmatrix} \end{split}\]

Advection-diffusion in code

function advdiff_matrix(n; kappa=1, wind=[1, 1]/sqrt(2))
    h = 2 / n
    rows = Vector{Int64}()
    cols = Vector{Int64}()
    vals = Vector{Float64}()
    idx(i, j) = (i-1)*n + j
    stencil_advect = [-wind[1], -wind[2], 0, wind[1], wind[2]] / h
    stencil_diffuse = [-1, -1, 4, -1, -1] * kappa / h^2
    for i in 1:n
        for j in 1:n
            if i in [1, n] || j in [1, n]
                push!(rows, idx(i, j))
                push!(cols, idx(i, j))
                push!(vals, 1.)
            else
                append!(rows, repeat([idx(i,j)], 5))
                append!(cols, [idx(i-1,j), idx(i,j-1), idx(i,j), idx(i+1,j), idx(i,j+1)])
                append!(vals, stencil_advect + stencil_diffuse)
            end
        end
    end
    sparse(rows, cols, vals)
end
advdiff_matrix (generic function with 1 method)

Spy the matrix

A = advdiff_matrix(10, wind=[1, 3], kappa=.1)
clim = norm(vec(A), Inf)
spy(A, marker=(:square, 3), c=:bwr, clims=(-clim, clim))
../_images/2021-10-08-fd-2d_14_0.svg
A = advdiff_matrix(20, wind=[0, 0], kappa=.1)
ev = eigvals(Matrix(A))
scatter(real(ev), imag(ev))
../_images/2021-10-08-fd-2d_15_0.svg

Plot a solution

n = 100
x = LinRange(-1, 1, n)
y = x
f = cos.(pi*x/2) * cos.(pi*y/2)'
heatmap(x, y, f, aspect_ratio=:equal)
../_images/2021-10-08-fd-2d_17_0.svg
A = advdiff_matrix(n, wind=[3,1], kappa=.01)
u = A \ vec(f)
heatmap(x, y, reshape(u, n, n), aspect_ratio=:equal)
../_images/2021-10-08-fd-2d_18_0.svg
  • What happens when advection dominates?

  • As you refine the grid?

Cost breadown and optimization

using ProfileSVG
function assemble_and_solve(n)
    A = advdiff_matrix(n)
    x = LinRange(-1, 1, n)
    f = cos.(pi*x/2) * cos.(pi*x/2)'
    u = A \ vec(f)
end

@profview assemble_and_solve(600)
Profile results in :-1 #15 in task.jl:411 eventloop in eventloop.jl:8 invokelatest in essentials.jl:706 #invokelatest#2 in essentials.jl:708 execute_request in execute_request.jl:67 softscope_include_string in SoftGlobalScope.jl:65 include_string in loading.jl:1116 eval in boot.jl:360 assemble_and_solve in In[133]:3 advdiff_matrix in In[131]:2 #advdiff_matrix#21 in In[131]:11 vect in array.jl:108 _array_for in array.jl:670 _array_for in array.jl:672 similar in abstractarray.jl:784 similar in abstractarray.jl:785 Array in boot.jl:457 Array in boot.jl:448 umfpack_free_symbolic in umfpack.jl:591 umfpack_free_numeric in umfpack.jl:606 umfpack_dl_free_numeric in umfpack.jl:602 setindex! in array.jl:843 in in operators.jl:1136 iterate in array.jl:781 < in int.jl:448 < in int.jl:441 #advdiff_matrix#21 in In[131]:16 repeat in abstractarraymath.jl:224 repeat##kw in abstractarraymath.jl:267 #repeat#1 in abstractarraymath.jl:269 repeat_inner_outer in abstractarraymath.jl:327 repeat_outer in abstractarraymath.jl:347 similar in array.jl:357 Array in boot.jl:448 repeat_outer in abstractarraymath.jl:348 repeat_outer in abstractarraymath.jl:350 setindex! in array.jl:870 unsafe_copyto! in array.jl:266 pointer in abstractarray.jl:1119 unsafe_convert in pointer.jl:65 unsafe_copyto! in array.jl:271 repeat_outer in abstractarraymath.jl:352 append! in array.jl:980 _growend! in array.jl:888 append! in array.jl:981 copyto! in array.jl:304 vect in array.jl:108 _array_for in array.jl:670 _array_for in array.jl:672 similar in abstractarray.jl:784 similar in abstractarray.jl:785 Array in boot.jl:457 Array in boot.jl:448 #advdiff_matrix#21 in In[131]:17 append! in array.jl:980 _growend! in array.jl:888 append! in array.jl:981 copyto! in array.jl:304 + in int.jl:87 vect in array.jl:108 _array_for in array.jl:670 _array_for in array.jl:672 similar in abstractarray.jl:784 similar in abstractarray.jl:785 Array in boot.jl:457 Array in boot.jl:448 setindex! in array.jl:843 #advdiff_matrix#21 in In[131]:18 append! in array.jl:980 _growend! in array.jl:888 append! in array.jl:981 copyto! in array.jl:304 _copyto_impl! in array.jl:313 unsafe_copyto! in array.jl:271 + in arraymath.jl:47 broadcast_preserving_zero_d in broadcast.jl:872 materialize in broadcast.jl:883 copy in broadcast.jl:908 copyto! in broadcast.jl:936 copyto! in broadcast.jl:983 macro expansion in simdloop.jl:77 macro expansion in broadcast.jl:984 setindex! in array.jl:843 similar in broadcast.jl:196 similar in broadcast.jl:197 similar in abstractarray.jl:784 similar in abstractarray.jl:785 Array in boot.jl:465 Array in boot.jl:457 Array in boot.jl:448 #advdiff_matrix#21 in In[131]:22 sparse in sparsematrix.jl:952 dimlub in sparsematrix.jl:948 maximum in reducedim.jl:873 #maximum#687 in reducedim.jl:873 _maximum in reducedim.jl:877 #_maximum#689 in reducedim.jl:877 _maximum in reducedim.jl:878 #_maximum#690 in reducedim.jl:878 mapreduce in reducedim.jl:310 #mapreduce#672 in reducedim.jl:310 _mapreduce_dim in reducedim.jl:318 _mapreduce in reduce.jl:415 sparse in sparsematrix.jl:956 sparse in sparsematrix.jl:756 Array in boot.jl:448 sparse in sparsematrix.jl:767 sparse! in sparsematrix.jl:841 iterate in range.jl:674 == in promotion.jl:410 sparse! in sparsematrix.jl:849 setindex! in array.jl:843 sparse! in sparsematrix.jl:857 > in operators.jl:305 < in int.jl:83 sparse! in sparsematrix.jl:861 >= in operators.jl:352 <= in int.jl:442 sparse! in sparsematrix.jl:862 setindex! in array.jl:843 sparse! in sparsematrix.jl:864 getindex in array.jl:805 iterate in range.jl:674 == in promotion.jl:410 sparse! in sparsematrix.jl:875 fill! in array.jl:335 setindex! in array.jl:843 sparse! in sparsematrix.jl:885 != in operators.jl:204 == in promotion.jl:410 sparse! in sparsematrix.jl:889 + in int.jl:87 sparse! in sparsematrix.jl:920 - in int.jl:86 sparse! in sparsematrix.jl:922 getindex in array.jl:805 sparse! in sparsematrix.jl:926 setindex! in array.jl:843 iterate in range.jl:674 == in promotion.jl:410 repeat_outer in abstractarraymath.jl:345 repeat##kw in abstractarraymath.jl:266 assemble_and_solve in In[133]:6 \ in linalg.jl:1564 ishermitian in sparsematrix.jl:3428 is_hermsym in sparsematrix.jl:3437 copy in array.jl:349 \ in linalg.jl:1567 \ in factorization.jl:104 ldiv! in umfpack.jl:497 ldiv! in umfpack.jl:510 _Aq_ldiv_B! in umfpack.jl:527 _AqldivB_kernel! in umfpack.jl:532 solve! in umfpack.jl:366 lu in umfpack.jl:155 #lu#1 in umfpack.jl:156 copy in array.jl:349 decrement in SuiteSparse.jl:15 copy in array.jl:349 decrement! in SuiteSparse.jl:12 setindex! in array.jl:843 iterate in range.jl:674 == in promotion.jl:410 #lu#1 in umfpack.jl:161 umfpack_numeric! in umfpack.jl:324 #umfpack_numeric!#13 in umfpack.jl:325 umfpack_symbolic! in umfpack.jl:304 #umfpack_numeric!#13 in umfpack.jl:327

What’s left?

  • Symmetric Dirichlet boundary conditions

  • Symmetric Neumann boundary conditions

  • Verification with method of manufactured solutions

  • Non-uniform grids

  • Upwinding for advection-dominated problems

  • Variable coefficients

  • Time-dependent problems

  • Fast algebraic solvers