Lecture 4

Agenda
📚 Multi-physics and going 2D
🚧 Exercises:

• Implement advection-diffusion solver in 2D

• Use Makie.jl for visualisation

---

Content

👉 get started with exercises

---

Coupled multi-physics in 2D

The goal of this lecture 4 is to:

• Learn the difference between different time integration schemes in the PT method

• Solve partial differential equations in 2D

• Better understand the coupling between physical processes

• Plot advanced graphics with Makie.jl

Explicit vs semi-implicit fluxes

In Lecture 3, you learned how to solve elliptic PDEs using the pseudo-transient (PT) method by augmenting the definition of the diffusion flux with a pseudo-time derivative:

\[ \rho\frac{\partial q}{\partial t} + \frac{q}{D} = -\frac{\partial C}{\partial x}~. \]

When discretising the first term on the left-hand side, we use first-order finite differences to approximate the pseudo-time derivative:

\[ \rho\frac{\partial q}{\partial t} \approx \rho\frac{q^{n+1} - q^n}{\Delta\tau}~. \]

For the second term in Eq. (1), \(q/D\), there are two possible choices:

1. Use the flux from the current pseudo-time layer \(n\);

2. Use the flux from the next pseudo-time layer \(n+1\).

In the first case, we use an explicit flux discretization; in the second case, we use a semi-implicit discretization.

In Lecture 3, we used the explicit discretization of fluxes. Let’s now implement the semi-implicit version!

1. 👉 Use your script for the 1D steady diffusion problem from the previous lecture, or start from this script;

2. Create a new file called steady_diffusion_implicit_flux_1d.jl for this exercise;

3. Think about how to compute the flux when using \(q^{n+1}/D\) in the flux update rule, and implement the semi-implicit scheme.

This script should produce the same final result as the explicit version. So why bother with another scheme?

👉 Modify the definition of the pseudo-time step. Replace this line:

dτ      = dx / sqrt(1 / ρ) / 1.1

with

dτ      = dx / sqrt(1 / ρ)

Observe how the PT iterations converge in the semi-implicit case with the theoretically maximal pseudo-time step, while the explicit flux discretization diverges.

Going 2D

Converting the 1D code to higher dimensions is remarkably easy thanks to the explicit time integration scheme. First, we define the domain size and the number of grid points in the y-direction:

# physics
lx, ly  = 20.0, 20.0
# ...
# numerics
nx, ny  = 100, 101

Next, we compute the grid spacing, the coordinates of grid cell centers, and update the pseudo-time step to satisfy the 2D stability criterion:

# derived numerics
dx, dy  = lx / nx, ly / ny
xc, yc  = LinRange(dx / 2, lx - dx / 2, nx), LinRange(dy / 2, ly - dy / 2, ny)
dτ      = min(dx, dy) / sqrt(1 / ρ) / sqrt(2)

We now allocate 2D arrays for the concentration field and the fluxes:

# array initialisation
C       = @. 1.0 + exp(-(xc - lx / 4)^2 - (yc' - ly / 4)^2) - xc / lx
qx, qy  = zeros(nx-1, ny), zeros(nx, ny-1)

Finally, we add the update rules for the second dimension:

while err >= ϵtol && iter <= maxiter
    #qx                 .-= ...
    #qy                 .-= ...
    #C[2:end-1,2:end-1] .-= ...
    # ...
end

Last thing to fix is the visualisation, as now we want the top-down view of the computational domain:

p1 = heatmap(xc, yc, C'; xlims=(0, lx), ylims=(0, ly), clims=(0, 1), aspect_ratio=1,
             xlabel="lx", ylabel="ly", title="iter/nx=$(round(iter / nx, sigdigits=3))")

Let's run the simulation:

Coupled systems of PDEs

You’ve learned how to solve second-order partial differential equations with a single independent variable — nice work!

Now, let’s take a step further and look at a slightly more advanced example with two independent variables.

This is a small but important move toward real-world applications.

Here's the system of equations:

In this system, the two main variables are \(P\) and \(T\), which we interpret as the pressure and temperature of a fluid:

1. The flux \(q\) is defined so that the fluid flows from regions of higher pressure to lower pressure.

2. The flow is enhanced by temperature perturbations — hotter fluid is more buoyant.

3. The temperature changes due to thermal conduction, but heat is also transported by the moving fluid.

👉 If we ignore coupling terms, what PDE types are the equations for \(P\) and \(T\)?

We’ll use an operator splitting approach, similar to Exercise 2 from Lecture 3.

• Step 1: Solve the elliptic equation for pressure \(P\), assuming the temperature \(T\) remains fixed.

• Step 2: Update \(T\) using the flux \(q\) computed in Step 1.

👉 Start from making a copy of your own 1D steady diffusion script or use this one

Rename the file to double_diffusion_1D.jl, and the main function to double_diffusion_1D() accordingly.

First, rename variables C and qx to P and qDx, respectively. Rename the diffusion coefficient dc to k.

Add new physical parameters:

# physics
lx      = 20.0
λ       = 0.001
k       = 1.0
α       = 1.0

Next, we will streamline a bit the PT parameters (it will be helpful in the next lecture).

Rename and replace the PT parameters:

qx         .-= dτ ./ (ρ * dc .+ dτ) .* (qx .+ dc .* diff(C) ./ dx)
C[2:end-1] .-= dτ .* diff(qx) ./ dx

with

qDx        .-= (qDx .+ k .* diff(P) ./ dx) ./ (θ_dτ_D + 1.0)
P[2:end-1] .-= (diff(qDx) ./ dx) ./ β_dτ_D

Use the following definitions for the new parameters in the # derived numerics section:

cfl     = 0.99
re_D    = 2π
θ_dτ_D  = lx / re_D / (cfl * dx)
β_dτ_D  = k * re_D / (cfl * dx * lx)

Implement the physical time stepping. Add the number of time steps and visualisation frequency into # numerics section of the script:

...
ncheck  = ceil(Int, 0.25nx)
nt      = 50
nvis    = 5

Remove the arrays storing the error evolution history, and wrap the iterative PT loop with the physical time loop:

for it in 1:nt
    @printf("it = %d\n", it)
    iter = 1; err = 2ϵtol
    while err >= ϵtol && iter <= maxiter
        #qDx        .-= ...
        #P[2:end-1] .-= ...
        if iter % ncheck == 0
            #err = ...
            @printf("  iter = %.1f × N, err = %1.3e\n", iter / nx, err)
        end
        iter += 1
    end
    # TODO
end

Add temperature arrays; keep pressure and fluid flux zero:

# temperature
T   = @. exp(-(xc - lx/4)^2)
T_i = copy(T)
# pressure
P   = zeros(nx)
qDx = zeros(Float64, nx - 1)

After the iterative loop for the pressure:

• Add the computation of the stable time step

• Implement diffusion and advection of temperature as two separate substeps:

dt  = min(dta, dtd)
# temperature
#T[2:end-1] .+= ...
#T[2:end-1] .-= ...
if it % nvis == 0
    # visualisation
    p1 = plot(xc, [T_i, T]; xlims=(0, lx), ylabel="Temperature", title="iter/nx=$(round(iter/nx,sigdigits=3))")
    p2 = plot(xc, P       ; xlims=(0, lx), xlabel="lx", ylabel="Pressure")
    display(plot(p1, p2; layout=(2, 1)))
end

Finally, add the coupling between the fluid flux qDx and the temperature T in the form of the term \(\alpha T\). Run the script.

Can you explain what you see?

⤴ back to Content

---

Visualisation with Makie.jl

Until now, we used Plots.jl for visualisation.

The cool thing about Plots.jl is how simple it is to plot something.

For more control and advanced graphics there is Makie.jl.

The docs are here.

Like Plots.jl, Makie.jl supports multiple backends, e.g.:

• GLMakie.jl for hardware-accelerated graphics, supports most features, but requires dedicated GPU;

• CairoMakie.jl for headless systems and publication ready vector graphics;

• WGLMakie.jl for interactive graphics on the web.

Basic usage

For basic usage, Makie.jl offers a simple high-level API:

using CairoMakie

plot(1:3)  # a scatter plot, for a line use `line`
A = rand(50, 50);
#try heatmap

Advanced usage

For more control, you can create and manipulate figures, axes, plots etc. as separate objects.

f = Figure()
scatter(f[1, 1], rand(100, 2))
lines(f[1, 2], cumsum(randn(100)))

ax = Axis(f[2, 1]; xlabel="x", ylabel="y", title="subplot")
lines!(ax, cumsum(randn(20)); label="line", linewidth=3, color=:red)
scatter!(ax, cumsum(randn(20)); label="scatter", marker=:cross, markersize=rand(5:20, 20))
axislegend(ax)

Animations

To create simple videos, you can use the record function

f = Figure()
ax = Axis(f[1, 1][1, 1]; aspect=DataAspect())
hm = heatmap!(ax, rand(10, 10))
cb = Colorbar(f[1, 1][1, 2], hm)

record(f, "anim.mp4"; fps=30) do io
    for i = 1:100
        hm[1] = rand(10,10) # simple way to update the plot in-place
        recordframe!(io)
    end
end

Live "animation"

To simply plot display figure in a computational loop like we did with Plots.jl, use display function:

for i = 1:10
    hm[1] = rand(10,10) # simple way to update the plot in-place
    display(f)
end

⤴ back to Content

Exercises - lecture 4

Exercise 1 — Advection-diffusion in 2D

👉 See Logistics for submission details.

The goal of this exercise is to

• Extend the advection-diffusion solver with implicit diffusion step from 1D to 2D

• Implement the upwind advection scheme in 2D

• Modify the problem configuration

Create a code implicit_advection_diffusion_2D.jl for this exercise and add it to the lecture4 folder in your private GitHub repo. Report the results of this exercise within a new section in the README.

Getting started

1. Duplicate the the implicit_diffusion_1D.jl code you created in Exercise 1 from Homework 3 and name it implicit_advection_diffusion_2D.jl.

2. Extend the 1D calculations to 2D

3. Add advection as in Exercise 2 from Homework 3 and name it implicit_advection_diffusion_2D.jl.

Modify the initial conditions to include following parameters in the # physics section:

# physics
lx, ly  = 10.0, 10.0
dc      = 1.0
vx      = 10.0
vy      = -10.0

and implement following initial condition

# array initialisation
C       = @. exp(-(xc - lx / 4)^2 - (yc' - 3ly / 4)^2)
C_old   = copy(C)

Choose the time step according to the following (stability) criterion:

dt = min(dx / abs(vx), dy / abs(vy)) / 2

Also make sure to use the following numerical parameters (in number of grid points nx,ny)

# numerics
nx, ny  = 200, 201
ϵtol    = 1e-8
maxiter = 10nx
ncheck  = ceil(Int, 0.02nx)
nt      = 100

Note that the iterative pseudo-timestep limitation should be updated to

dτ = min(dx, dy) / sqrt(1 / ρ) / sqrt(2)

for 2D configurations.

Task 1

Repeat the steps from the Exercise 1 in Homework 3 to create the implicit time-dependent diffusion solver but in 2D. Do not include advection yet. Pay attention to add information relative to the second dimension whenever it's needed.

Make a short animation showing the time evolution of the concentration field C during nt = 100 physical time steps. The figure should contain 2 subplots, the first displaying the heatmap of the C field and the second the evolution of the by nx normalised iteration count:

# visualisation
p1 = heatmap(xc, yc, C'; xlims=(0, lx), ylims=(0, ly), clims=(0, 1), aspect_ratio=1,
             xlabel="lx", ylabel="ly", title="iter/nx=$(round(iter/nx,sigdigits=3))")
p2 = plot(iter_evo, err_evo; xlabel="iter/nx", ylabel="err",
          yscale=:log10, grid=true, markershape=:circle, markersize=10)
display(plot(p1, p2; layout=(2, 1)))

Task 2

Include now the advection step in a similar way as in the 1D case from the previous exercise, i.e., adding them after the iteration loop within the time loop. Use advection velocities and parameters listed above, taking care in implementing the "upwind" strategy discussed in Lecture 2.

Make a short animation showing the time evolution of the concentration field C during nt = 50 physical time steps using the same figure layout as for Task 1.

⤴ back to Content

---

Exercise 2 — Using Makie.jl for visualisation

👉 See Logistics for submission details.

The goals of this exercise are to:

• learn visualisation techniques in scientific computing;

• use Makie.jl to visualise the 2D simulation.

Create a script implicit_advection_diffusion_makie_2D.jl for this exercise and add it to the lecture4 folder in your private GitHub repository. Summarise your results in a new section of the README.

Getting started

1. Duplicate the implicit_advection_diffusion_2D.jl code you created in Exercise 1 and rename it to implicit_advection_diffusion_makie_2D.jl.

2. Remove all visualisation that uses Plots.jl.

3. Recreate the animation from Exercise 1 with Makie.jl, and add an arrow (quiver) plot to visualise the flux vector field.

Create the figure, two axes, a heatmap, a quiver plot, and a line plot with markers before the time loop:

# time loop
fig = Figure(size=(400, 650))
ax1 = Axis(...)
ax2 = Axis(...)
hm  = heatmap!(ax1, ...)
cb  = Colorbar(...)
ar  = arrows2d!(ax1, ...)
plt = scatterlines!(ax2, Float64[], Float64[])
record(fig, "heatmap_arrows.mp4"; fps=20) do io
    for it = 1:nt
        # ...
        # visualisation
        # update plots here
        recordframe!(io)
    end
end

In an arrow plot, showing one arrow per grid point would result in too much visual noise. Show an arrow to every 10th cell in x and y directions.

Hints

• Use the colormap and colorrange attributes to configure the heatmap.

• To update plot data, index into the plot’s arguments. For example, for a heatmap created with hm = heatmap!(ax, x, y, z), update the field with hm[3] = z_new. Here, hm[1] and hm[2] correspond to the x and y data, respectively.

• To get every nth value of a 2D array, use syntax Q[1:n:end, 1:n:end]. Same applies to 1D and 3D arrays.

Task 1

Make a short animation showing the time evolution of the concentration field C during nt = 100 physical time steps. The numerical algorithm should be the same as in Exercise 1.

The figure should contain 2 subplots, the first displaying the heatmap of the C field with the "roma" colormap, and the quiver plot showing magnitude and direction of the flux vector q. The second subplot should show the evolution of the by nx normalised iteration count. In the time loop, only update the existing plot, don't create new figure every time step.

You should get an animation like this:

⤴ back to Content