11) Variable coefficients

11) Variable coefficients#

Last time:

Conservation Laws
Intro to Finite Volume Methods
Nonhomogeneous equations
Cell-centered grids and vertex-centered grids

Today:

Validation and Verification
Non constant coefficients
Manufactured solutions for variable coefficients
Discretizing in conservative form
Sparsity
5.1 Upwind advection-diffusion solver using COO

1. Verification and Validation (V&V)#

In Numerical Analysis, the concepts of Verification and Validation can be oversimplified in a succinct manner by saying that “verification is checking that you are solving the equations right” (verifying the numerics) and “validation is checking that you are solving the right equations” (for physical or engineering problems this involves verifying the physics - often done by comparing model results with actual data given by observations).

In essence:

Verification: solving the problem right#

It answers:

✅ Is our code correct?

It does not answer:

❌ Does the implemented model have predictive value/is it fit for the purpose?

Validation: solving the right problem#

It answers:

✅ Are the equations a suitable model for the relevant processes?

It does not answer:

❌ Are we getting closer to the analytical solution?

Similarly, in Software Engineering, we can say that verification is asking the question: “are we building the product right (i.e., correctly)?” and validation: “are we building the right product (i.e., does it match the original plan/design)?”
From a modeling and simulation point of view in scientific computing, we can have best practices like checking the convergence or order of accuracy of a numerical scheme and this would be part of verification. While, if we are modeling a physical system, using conservation laws, then checking that the quantities of interest are really conserved is part of validation (in fact, we would check that our assumptions are not violated).

This might remind you of the differences between forward/backward error.

2. Variable coefficients#

Problems that can be described with variable coefficients:

Heat conduction: steel, brick, wood, foam
Electrical conductivity: copper, rubber, air
Elasticity: steel, rubber, concrete, rock
Linearization of nonlinear materials
- ketchup, glacier ice, rocks (mantle/lithosphere)

using Plots
using LinearAlgebra
using LaTeXStrings
default(linewidth=3)
default(legendfontsize=12)

# Recall our utility/convenience `vander` and  `fdstencil` from Lecture 7

function vander(x, k=nothing)
    if k === nothing
        k = length(x)
    end
    V = ones(length(x), k)
    for j = 2:k
        V[:, j] = V[:, j-1] .* x
    end
    V
end

function fdstencil(source, target, k)
    "kth derivative stencil from source to target"
    x = source .- target
    V = vander(x)
    rhs = zero(x)'
    rhs[k+1] = factorial(k)
    rhs / V
end

function my_spy(A)
    cmax = norm(vec(A), Inf)
    s = max(1, ceil(120 / size(A, 1)))
    spy(A, marker=(:square, s), c=:diverging_rainbow_bgymr_45_85_c67_n256, clims=(-cmax, cmax))
end

my_spy (generic function with 1 method)

kappa_step(x) = .1 + .9 * (x > 0)
kappa_smooth(x) = .55 + .45 * sin(pi*x/2)
plot([kappa_step, kappa_smooth], xlims=(-1, 1), ylims=(0, 1), label=["κ step" "κ smooth"])

(16)#\[\begin{align} -\big(\kappa(x) u_x\big)_x &= 0 & u(-1) &= 0 & \kappa u_x(1) &= 1 \end{align}\]

What physical scenario could this represent?
Qualitatively, what would a solution look like?

A naive finite difference solver#

Conservative (divergence) form	Non-divergence form
\(-(\kappa u_x)_x = 0\)	\(-\kappa u_{xx} - \kappa_x u_x = 0\)

function poisson_nondivergence(x, spoints, kappa, forcing; leftbc=(0, zero), rightbc=(0, zero))
    n = length(x)
    L = zeros(n, n)
    rhs = forcing.(x)
    kappax = kappa.(x)
    for i in 2:n-1
        jleft = min(max(1, i-spoints÷2), n-spoints+1)
        js = jleft : jleft + spoints - 1
        kappa_x = fdstencil(x[js], x[i], 1) * kappax[js]
        L[i, js] = -fdstencil(x[js], x[i], 2) .* kappax[i] - fdstencil(x[js], x[i], 1) * kappa_x
    end
    L[1,1:spoints] = fdstencil(x[1:spoints], x[1], leftbc[1])
    if leftbc[1] == 1
        L[1, :] *= kappax[1]
    end
    L[n,n-spoints+1:n] = fdstencil(x[n-spoints+1:n], x[n], rightbc[1])
    if rightbc[1] == 1
        L[n, :] *= kappax[n]
    end
    rhs[1] = leftbc[2](x[1])
    rhs[n] = rightbc[2](x[n])
    L, rhs
end

poisson_nondivergence (generic function with 1 method)

Let’s try it.

x = LinRange(-1, 1, 60)
kappa = kappa_step
L, rhs = poisson_nondivergence(x, 3, kappa, zero, rightbc=(1, one))
u = L \ rhs
plot(x, u, label="numerical solution")
plot!(x -> 5*kappa(x), label="κ step")

3. Manufactured solutions for variable coefficients#

manufactured(x) = tanh(2x)
d_manufactured(x) = 2*cosh(2x)^-2
d_kappa_smooth(x) = .45*pi/2 * cos(pi*x/2)
flux_manufactured_kappa_smooth(x) = kappa_smooth(x) * d_manufactured(x)
function forcing_manufactured_kappa_smooth(x)
    8 * tanh(2x) / cosh(2x)^2 * kappa_smooth(x) -
     d_kappa_smooth(x) * d_manufactured(x)
end
x = LinRange(-1, 1, 40)
L, rhs = poisson_nondivergence(x, 3, kappa_smooth,
    forcing_manufactured_kappa_smooth,
    leftbc=(0, manufactured), rightbc=(1, flux_manufactured_kappa_smooth))
u = L \ rhs;

plot(x, u, marker=:circle, legend=:bottomright, label= L"u_h")
plot!([manufactured flux_manufactured_kappa_smooth forcing_manufactured_kappa_smooth kappa_smooth],
label=["u_manufactured" "flux" "forcing" "κ smooth"], ylim=(-2.5, 2))

Convergence#

Let’s try first with the smooth \(\kappa\) function for the non constant coefficients.

function poisson_error(n, spoints=3)
    x = LinRange(-1, 1, n)
    L, rhs = poisson_nondivergence(x, spoints, kappa_smooth,
        forcing_manufactured_kappa_smooth,
        leftbc=(0, manufactured), rightbc=(1, flux_manufactured_kappa_smooth))
    u = L \ rhs
    norm(u - manufactured.(x), Inf)
end
ns = 2 .^ (3:10)
plot(ns, poisson_error.(ns, 5), marker=:auto, xscale=:log10, yscale=:log10, label = "error")
plot!([n -> n^-2, n -> n^-4], label=["\$n^{-2}\$" "\$n^{-4}\$"])

✅ Verified!#

What would happen with the discontinuous coefficients step function#

x = LinRange(-1, 1, 20)
L, rhs = poisson_nondivergence(x, 3, kappa_step,
    zero,
    leftbc=(0, zero), rightbc=(0, one))
u = L \ rhs
plot(x, u, marker=:auto, legend=:bottomright, label = L"u_h")

4. Discretizing in conservative form#

Conservative (divergence) form	Non-divergence form
\(-(\kappa u_x)_x = 0\)	\(-\kappa u_{xx} - \kappa_x u_x = 0\)

function poisson_conservative(n, kappa, forcing; leftbc=(0, zero), rightbc=(0, zero))
    x = LinRange(-1, 1, n)
    xstag = (x[1:end-1] + x[2:end]) / 2
    L = zeros(n, n)
    rhs = forcing.(x)
    kappa_stag = kappa.(xstag)
    for i in 2:n-1
        flux_L = kappa_stag[i-1] * fdstencil(x[i-1:i], xstag[i-1], 1)
        flux_R = kappa_stag[i] * fdstencil(x[i:i+1], xstag[i], 1)
        js = i-1:i+1
        weights = -fdstencil(xstag[i-1:i], x[i], 1)
        L[i, i-1:i+1] = weights[1] * [flux_L..., 0] + weights[2] * [0, flux_R...]
    end
    if leftbc[1] == 0
        L[1, 1] = 1
        rhs[1] = leftbc[2](x[1])
        rhs[2:end] -= L[2:end, 1] * rhs[1]
        L[2:end, 1] .= 0
    end
    if rightbc[1] == 0
        L[end,end] = 1
        rhs[end] = rightbc[2](x[end])
        rhs[1:end-1] -= L[1:end-1,end] * rhs[end]
        L[1:end-1,end] .= 0
    end
    x, L, rhs
end

poisson_conservative (generic function with 1 method)

Compare conservative vs non-divergence forms#

forcing = zero # try one
x, L, rhs = poisson_conservative(20, kappa_step,
    forcing, leftbc=(0, zero), rightbc=(0, one))
u = L \ rhs
plot(x, u, marker=:circle, legend=:bottomright, label = L"u_h", title = "Conservative form")

x = LinRange(-1, 1, 20)
L, rhs = poisson_nondivergence(x, 3, kappa_step,
    forcing, leftbc=(0, zero), rightbc=(0, one))
u = L \ rhs
plot(x, u, marker=:circle, legend=:bottomright, label = L"u_h", title = "Non-divergence form" )

Continuity of flux#

forcing = zero

L, rhs = poisson_nondivergence(x, 3, kappa_step,
    forcing, leftbc=(0, zero), rightbc=(0, one))
u = L \ rhs
plot(x, u, marker=:circle, legend=:bottomright, label = L"u_h", title = "Non-divergence form" )

xstag = (x[1:end-1] + x[2:end]) ./ 2 # a staggered grid
du = (u[1:end-1] - u[2:end]) ./ diff(x)
plot(xstag, [du .* kappa_step.(xstag)], marker=:circle, ylims=[-1, 1], label = L"D(u) \kappa \textrm{-step}")

Manufactured solutions with discontinuous coefficients#

We need to be able to evaluate derivatives of the flux \(-\kappa u_x\).
A physically-realizable solution would have continuous flux, but we we’d have to be making a physical solution to have that in verification.
Idea: replace the discontinuous function with a continuous one with a rapid transition (something that can “mollify” the discontinuity).

kappa_tanh(x, epsilon=.1) = .55 + .45 * tanh(x / epsilon)
d_kappa_tanh(x, epsilon=.1) = .45/epsilon * cosh(x/epsilon)^-2
plot([kappa_tanh], label = L"\kappa = \tanh")

Solving with the smoothed step function \(\kappa\):

kappa_tanh(x, epsilon=.02)= .55 + .45 * tanh(x / epsilon)
d_kappa_tanh(x, epsilon=.02) = .45/epsilon * cosh(x/epsilon)^-2
flux_manufactured_kappa_tanh(x) = kappa_tanh(x) * d_manufactured(x)
function forcing_manufactured_kappa_tanh(x)
    8 * tanh(2x) / cosh(2x)^2 * kappa_tanh(x) -
    d_kappa_tanh(x) * d_manufactured(x)
end
x, L, rhs = poisson_conservative(200, kappa_tanh,
    forcing_manufactured_kappa_tanh,
    leftbc=(0, manufactured), rightbc=(0, manufactured))
u = L \ rhs
plot(x, u, marker=:circle, legend=:bottomright, title=L"\kappa = tanh", label = L"u_h")
plot!([manufactured flux_manufactured_kappa_tanh forcing_manufactured_kappa_tanh kappa_tanh],
label=["u_manufactured" "flux" "forcing" L"\kappa"], ylim=(-5, 5))

Error  = (norm(u - manufactured.(x), Inf))

0.0005721837233710336

Convergence#

function poisson_error(n, spoints=3)
    x, L, rhs = poisson_conservative(n, kappa_tanh,
        forcing_manufactured_kappa_tanh,
        leftbc=(0, manufactured), rightbc=(0, manufactured))
    u = L \ rhs
    norm(u - manufactured.(x), Inf)
end
ns = 2 .^ (3:10)
plot(ns, poisson_error.(ns, 3), marker=:auto, xscale=:log10, yscale=:log10, label = "error")
plot!(n -> n^-2, label="\$1/n^2\$")

Outlook#

Manufactured solutions can be misleading for verification in the presence of rapid coefficient variation.
- Forcing terms may be poorly resolved on the grid
Non-conservative methods are scary
- Can look great during verification on smooth problems, then break catastrophically
Conservative methods are generally robust, but now we have choices to make
- What if \(\kappa(x)\) is gridded data? How do we evaluate at staggered points?
- What does this mean for boundary conditions, \(\kappa u_x = g\)
- What do staggered points mean on unstructured/adaptive meshes?
- Which terms should we evaluate where?

5. Sparsity#

When we use a direct system solver like

u = L \ rhs

how expensive is it?

For instance, we have seen the advection-diffusion problem:

function advdiff_upwind(n, kappa, wind, forcing)
    x = LinRange(-1, 1, n) # cell values
    xstag = (x[1:end-1] + x[2:end]) / 2 # edge/face values
    L = zeros(n, n)
    rhs = forcing.(x)
    kappa_stag = kappa.(xstag)
    for i in 2:n-1 # for the interior nodes
        # first-order upwinding for the advective flux: use the value from the upstream side of the face
        flux_L = -kappa_stag[i-1] * fdstencil(x[i-1:i], xstag[i-1], 1) +
            wind * (wind > 0 ? [1 0] : [0 1]) # for w>0, wu ≈ w u_{i-1}, for w<0, wu ≈ w u_{i+1} 
        flux_R = -kappa_stag[i] * fdstencil(x[i:i+1], xstag[i], 1) +
            wind * (wind > 0 ? [1 0] : [0 1]) # for w>0, wu ≈ w u_{i+1}, for w<0, wu ≈ w u_{i-1} 
        js = i-1:i+1
        weights = fdstencil(xstag[i-1:i], x[i], 1)
        L[i, i-1:i+1] = weights[1] * [flux_L..., 0] + weights[2] * [0, flux_R...] # hence the flux is computed using F_{i-1/2} and F_{i+1/2}
    end
    L[1, 1] = 1
    rhs[1] = 0
    L[2:end, 1] .= 0
    L[end,end] = 1
    rhs[end] = 0
    L[1:end-1,end] .= 0
    x, L, rhs
end

advdiff_upwind (generic function with 1 method)

n = 30; h = 2/n
kappa = 1
wind = -1000
x, L, rhs = advdiff_upwind(2000, x -> kappa, wind, one)
@timev L \ rhs;

  0.072349 seconds (4 allocations: 30.548 MiB, 11.43% gc time)
elapsed time (ns):  72349023
gc time (ns):       8271299
bytes allocated:    32032304
pool allocs:        1
non-pool GC allocs: 2
malloc() calls:     1
free() calls:       0
minor collections:  1
full collections:   0

function elapsed_solve(n)
    x, L, rhs = advdiff_upwind(n, x -> kappa, wind, one)
    @elapsed L \ rhs
end
ns = (1:6) * 1000
times = elapsed_solve.(ns)

6-element Vector{Float64}:
030589824
250805679
197368013
392390769
594081925
2253659

scatter(ns, times, xscale=:log10, yscale=:log10)
plot!(ns, [(ns/3000).^2, (ns/3000).^3], label=["n^2" "n^3"], xrange=(1e3, 1e4), xlabel = "n", ylabel = "time")

Sparsity#

The matrices derived via Finite Difference Schemes are very sparse. They are banded matrices with only a few non-zero entries along the diagonals.

x, L, rhs = advdiff_upwind(50, one, 1, one)
my_spy(L)

using Pkg
Pkg.add("SparseArrays")
using SparseArrays
sparse(L)

   Resolving package versions...

  No Changes to `~/.julia/environments/v1.10/Project.toml`
  No Changes to `~/.julia/environments/v1.10/Manifest.toml`

┌ Warning: Circular dependency detected.
│ Precompilation will be skipped for dependencies in this cycle:
│  ┌ Symbolics → SymbolicsForwardDiffExt
│  └─ Symbolics → SymbolicsPreallocationToolsExt
└ @ Pkg.API.Precompilation ~/julia-1.10.10/share/julia/stdlib/v1.10/Pkg/src/precompilation.jl:583

50×50 SparseMatrixCSC{Float64, Int64} with 144 stored entries:
⎡⠱⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎤
⎢⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⎥
⎣⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠑⎦

How does solving with the sparse matrix scale?#

function elapsed_solve(n)
    x, L, rhs = advdiff_upwind(n, x -> kappa, wind, one)
    L = sparse(L)
    @elapsed L \ x
end
ns = (1:10) * 1000
times = elapsed_solve.(ns)
scatter(ns, times, xscale=:log10, yscale=:log10)
plot!([n -> 1e-2*n/1000, n -> 1e-2*(n/1000)^2, n -> 1e-2*(n/1000)^3],
  label=["n" "n^2" "n^3"])

SparseArrays Compressed Sparse Column (CSC)#

# starting with a dense matrix
A = sparse([1 0 0 -2; 0 3.14 0 .4; 5.1 0 0 .06])

3×4 SparseMatrixCSC{Float64, Int64} with 6 stored entries:
 1.0   ⋅     ⋅   -2.0
  ⋅   3.14   ⋅    0.4
 5.1   ⋅     ⋅    0.06

@show A.colptr
@show A.rowval
A.nzval

A.colptr = [1, 3, 4, 4, 7]
A.rowval = [1, 3, 2, 1, 2, 3]

6-element Vector{Float64}:
0
1
14
 -2.0
4
06

CSC (or, more commonly outside Matlab/Julia, row-based CSR) are good for matrix operations (multiplication, solve), but sometimes inconvenient to assemble.

Matrix assembly using COO format#

3×4 SparseMatrixCSC{Float64, Int64} with 6 stored entries:
 1.0   ⋅     ⋅   -2.0
  ⋅   3.14   ⋅    0.4
 5.1   ⋅     ⋅    0.06

sparse([1, 1, 2, 2, 3, 3, 1],
       [1, 4, 2, 4, 1, 4, 1],
       [1, -2, 3.14, .4, 5.1, .06, 10])

3×4 SparseMatrixCSC{Float64, Int64} with 6 stored entries:
 11.0   ⋅     ⋅   -2.0
   ⋅   3.14   ⋅    0.4
  5.1   ⋅     ⋅    0.06

5.1 Upwind advection-diffusion solver using COO#

function advdiff_sparse(n, kappa, wind, forcing; rightbc=0)
    x = LinRange(-1, 1, n)
    xstag = (x[1:end-1] + x[2:end]) / 2
    rhs = forcing.(x)
    kappa_stag = kappa.(xstag)
    rows = [1, n]
    cols = [1, n]
    vals = [1., 1.] # diagonals entries (float)
    rhs[[1,n]] = [0, rightbc] # Dirichlet boundary condition
    for i in 2:n-1
        flux_L = -kappa_stag[i-1] * fdstencil(x[i-1:i], xstag[i-1], 1) +
            wind * (wind > 0 ? [1 0] : [0 1])
        flux_R = -kappa_stag[i] * fdstencil(x[i:i+1], xstag[i], 1) +
            wind * (wind > 0 ? [1 0] : [0 1])
        weights = fdstencil(xstag[i-1:i], x[i], 1)
        append!(rows, [i,i,i])
        append!(cols, i-1:i+1)
        append!(vals, weights[1] *  [flux_L..., 0] + weights[2] * [0, flux_R...])
    end
    L = sparse(rows, cols, vals)
    x, L, rhs
end

advdiff_sparse (generic function with 1 method)

Assembly cost:

n = 4000; h = 2/n
kappa = 1
wind = 100
x, L, rhs = advdiff_sparse(n, x -> kappa, wind, one)
@show minimum(diff(x))
plot(x, L \ rhs, legend=:none, title="Pe_h $(wind*h/kappa)")

minimum(diff(x)) = 0.0005001250312576255

n = 10000
@time advdiff_upwind(n, one, 1, one);
@time advdiff_sparse(n, one, 1, one);

  0.349460 seconds (598.37 k allocations: 794.344 MiB, 12.84% gc time)
  0.031009 seconds (588.92 k allocations: 35.243 MiB)

# It's also possible to dynamically insert.
# (But performance this way is generally poor.)
A = spzeros(5, 5)
A[1,1] = 3
A

5×5 SparseMatrixCSC{Float64, Int64} with 1 stored entry:
 3.0   ⋅    ⋅    ⋅    ⋅ 
  ⋅    ⋅    ⋅    ⋅    ⋅ 
  ⋅    ⋅    ⋅    ⋅    ⋅ 
  ⋅    ⋅    ⋅    ⋅    ⋅ 
  ⋅    ⋅    ⋅    ⋅    ⋅ 

11) Variable coefficients

Contents

11) Variable coefficients#

1. Verification and Validation (V&V)#

Verification: solving the problem right#

Validation: solving the right problem#

2. Variable coefficients#

A naive finite difference solver#

3. Manufactured solutions for variable coefficients#

Convergence#

✅ Verified!#

What would happen with the discontinuous coefficients step function#

4. Discretizing in conservative form#

Compare conservative vs non-divergence forms#

Continuity of flux#

Manufactured solutions with discontinuous coefficients#

Convergence#

Outlook#

5. Sparsity#

Sparsity#

How does solving with the sparse matrix scale?#

SparseArrays Compressed Sparse Column (CSC)#

Matrix assembly using COO format#

5.1 Upwind advection-diffusion solver using COO#