Unit 1: Introduction

Published

26/06/2026

TipWhat this unit is for

This unit is motivation, not curriculum. We solve a real-looking coastal-oceanography problem end-to-end up front so that everything that follows has a destination: the PDE theory in Unit 6, the deep-learning toolkit in Unit 2, the SciML stack in Units 3–4, the PINN machinery in Units 5 and 7, and the AIMS thermistor-chain capstone in Units 8–10. You are not expected to write or understand this code yet — the goal is to watch the workflow and build a mental picture of where the course lands.

This unit opens the course with a concrete, hands-on problem — a Brisbane River surge propagating into Moreton Bay — that threads together every theme you will revisit over the next nine units: PDEs, Julia, computation, real measurements, the inverse problem, ocean / water modelling, and PINNs. Then it previews the larger AIMS capstone in §1.3.

1.1 An example scientific question

Before we touch any code, let’s set up a concrete problem and the gap between what we can measure in it and what we actually want to know. The next two subsections sketch (a) the physical setting — Brisbane’s outflow into Moreton Bay — and (b) the asymmetry between cheap downstream observations and the hard-to-measure upstream forcing that drives them. That asymmetry is what the rest of the unit will try to close.

Where the Brisbane River meets Moreton Bay

Most of Brisbane’s storm runoff and a meaningful slice of South-East Queensland’s wastewater eventually drains down the Brisbane River and exits at Pile Light into Moreton Bay — a shallow, semi-enclosed body of water bounded by the mainland to the west, Moreton Island to the east, and North Stradbroke Island to the south-east. The bay is a Ramsar-listed estuarine nursery, a major shipping channel for the Port of Brisbane, and the recreational heart of South-East Queensland; what comes down the river matters to ecology, public health, fisheries, and shipping.

After heavy rain in the upper catchment, the river runs hard for a day or two and pushes a freshwater surge — a low-density plume carrying sediment, nutrients, and (occasionally) sewage overflow — out into the bay. The bay responds: tide gauges across the bay record a small, transient surface-elevation anomaly riding on top of the usual tide, even though the wind is calm. The Healthy Land & Water Moreton Bay catchment report cards have been tracking the downstream impact for years; the bay’s seagrass beds in the western flats still bear scars from the 2011 and 2022 floods.

What we actually measure, and what we don’t

Here is the operational catch: the bay has plenty of tide gauges (Mud Island, Wellington Point, Tangalooma anchorage, Russell / Macleay channel — Brisbane Bar / Pile Light too, but it sits so close to the river mouth that it doesn’t constrain the source much more than a direct measurement would), but the river outflow itself is hard to instrument in flood conditions — current meters get washed out, banks become inaccessible, suspended-sediment loads kill optical probes within hours. We rarely have a reliable time-history of what came out of the river during an event, only of how the bay responded over the days that followed.

A real tide gauge in Moreton Bay records the total water level

\eta_{\text{gauge}}(t) \;=\; \underbrace{\eta_{\text{tide}}(t)}_{\sim\, 1.5\,\text{m}} \;+\; \underbrace{\eta_{\text{surge}}(t)}_{\sim\, 0.45\,\text{m}} \;+\; \text{noise.}

The tidal part is huge compared to the surge we care about — roughly 1.5 m of astronomical tide vs. the ~45 cm freshwater surge — but it’s also predictable. From a long calibration record at each gauge, oceanographers fit a harmonic tidal model (a sum of \sim 30 sinusoids at known astronomical frequencies: M_2 at 12.42 h, S_2 at 12.00 h, K_1 at 23.93 h, O_1 at 25.82 h, and so on); \eta_{\text{tide}}(t) can then be reconstructed and subtracted from the live record. What remains is the de-tided residual, \eta_{\text{gauge}}(t) - \eta_{\text{tide}}(t) \approx \eta_{\text{surge}}(t) + \text{noise} — which is what we treat as the synthetic mooring data through this unit (with \sigma = 1.5\,\text{cm} Gaussian noise to imitate the harmonic-fit residual after de-tiding).

Two practical points worth carrying forward:

  • De-tiding is itself an inverse problem of sorts. Errors in the harmonic constants (especially the shallow-water tidal components M_4, M_6 near the river mouth) leak directly into the surge residual. Any source-recovery inversion built on tide-gauge data inherits that calibration uncertainty.
  • Surge timescale ≪ tidal timescale, but not by much. A 3-hour freshwater surge has substantial energy at periods near the diurnal and semi-diurnal tidal bands. Aggressive de-tiding (notch filters around the tidal lines) risks smearing surge content; sloppy de-tiding leaves tide leakage in the residual. The Brisbane River Hydraulic Model the Department of Environment uses for operational forecasts hits this trade-off head-on.

That sets up the scientific question we will solve in this unit:

TipThe Unit 1 question

Given sparse, noisy tide-gauge timeseries scattered across Moreton Bay, can we reconstruct the time-history of the surge at the river mouth that produced them?

If we can, the same machinery — a forward physical model plus a learned source — lets us infer contaminant load, freshwater discharge, and bay-flushing timescales from measurements that are already routinely collected. The river that you cannot instrument in a flood becomes legible through the gauges that survived in the bay.

The same shape of question recurs throughout coastal oceanography: given downstream observations, infer the upstream forcing. It is the classical inverse problem of geophysics — and it is notoriously ill-posed, because many different upstream scenarios can produce essentially the same downstream signal once the physics has smoothed, delayed, and noisified it. We will solve it twice in this unit so the ill-posedness is visible before we spend Units 5–9 learning the modern (PINN-based) tools that disambiguate it.

Build your own miniature de-tiding pipeline. Generate a synthetic 5-day gauge record at 3-minute sampling: a “tide” made of two sinusoids (M_2 at 12.42 h, amplitude 0.8 m; S_2 at 12.00 h, amplitude 0.3 m), plus the Unit 1 surge pulse 0.45\exp[-((t - 50\,\text{h})/3\,\text{h})^2], plus \sigma = 1.5\,\text{cm} Gaussian noise. Then:

  1. Fit the harmonic constants (amplitude and phase of both constituents) by linear least squares using only the first 36 hours — a “calm” calibration window before the surge.
  2. Reconstruct the tide over the full record, subtract it, and compare the residual to the true surge pulse.
  3. Repeat the fit using the whole record (surge included) as the calibration window. How much surge energy leaks into the harmonic fit, and what does that do to the recovered pulse?

💡 Hint

With the constituent frequencies known, the fit is linear least squares. Build the design matrix design(t) = hcat(ones(length(t)), sin.(ω[1].*t), cos.(ω[1].*t), sin.(ω[2].*t), cos.(ω[2].*t)) with ω = 2π ./ [12.42, 12.0], and solve with Julia’s backslash (A \ y). Step 1: select the calm window with a boolean mask calm = t .< 36.0, fit the coefficients on t[calm], then reconstruct the tide over the whole record from those coefficients and subtract it. Step 3: refit on the entire record and overlay both residuals on the true surge — watch the full-record fit pull surge energy near the tidal band into its constituents. The amplitude/phase form A\sin(\omega t + \phi) is equivalent to the sin+cos form, which is what keeps the problem linear.

Go to solution →

1.2 The worked example

This section walks the full pipeline end-to-end: we set up a tractable model of the bay, write down the physics, run a forward simulation to generate synthetic measurements, then try to invert those measurements back to the river-mouth source two different ways. Each subsubsection is one stage of that pipeline. All code lives under units/unit_01/scripts/ — the runtime + role table at the end of §1.2 gives the full inventory.

The geography, made tractable

Site map for the Unit 1 worked example. (a) The SWE-modelling domain (blue rectangle) covers 50 km east–west by 95 km north–south, takes in the mainland coast, Moreton Island, North & South Stradbroke, the Russell / Macleay / Karragarra cluster in the southern bay, and the open ocean to the east of Moreton. The four tide-gauge sites (gold) carry the names of real Moreton Bay sites, but their positions are approximate (see the note below). The red triangle near the western coast is the Brisbane River mouth, where the unknown surge ψ(t) is imposed. (b) Inset zoom on the river entrance / Pile Light area; the nearest informative gauge (G1, Mud Island) sits ~18 km north of the mouth in mid-bay, deliberately not on top of the source — see the Why not put G1 at Brisbane Bar? discussion.

We choose a domain that is small enough to evaluate in seconds on a laptop CPU but big enough to capture the geography that matters:

  • A 50 \times 95\,\text{km} rectangle at 500\,\text{m} resolution — so the SWE mesh is 100 \times 190 = 19\,000 cells, of which about half are water.
  • A hand-built bathymetry that echoes the real bay: a \sim 6\,\text{m} inner-bay flat, an eastern deep channel (the NE channel) rising to \sim 25\,\text{m} between Moreton Island and the mainland, the dredged \sim 14\,\text{m} Brisbane River shipping channel, and a 30–50\,\text{m} shelf east of Moreton Island.
  • A land mask for the mainland coast, Moreton Island, North Stradbroke Island, the South Stradbroke spit, the Russell / Macleay / Karragarra island cluster in the southern bay, and a Peel / Mud cluster mid-bay. The southern bay genuinely is half-blocked off by islands and mangroves — modelling it as open water (as we did in an earlier draft of this unit) made the surge propagation visibly wrong.

Bathymetry used by the model — depths in metres, land in grey. Note the sloped basin, the deeper eastern channel (NE channel), the dredged Brisbane River channel running east from the source, the shoals around the river mouth, the eastern shelf outside Moreton Island, and the way the southern bay is constricted by the Russell / Macleay / Karragarra island cluster and the South Stradbroke spit. This is what the SWE solver actually sees. The map is drawn in landscape with north to the left (← N) and east up, so the bay’s long 95 km north–south axis runs across the page.
NoteA geographical apology

This bay is a caricature, not a chart. The bathymetry is hand-built — plausible depths and the right overall shape, but not survey data — and the four gauges carry the names of real Moreton Bay sites only for flavour. Their positions are approximate, and a couple are frankly off: they were nudged onto the idealised grid to give a clean spread of informative observers (mid-bay, across from the river, in the island lee, down south) and then snapped to the nearest water cell, so several sit a few kilometres from their namesakes — “Tangalooma Roads” and “Russell Channel” most of all. Please don’t navigate by this map; it exists to make the inverse problem legible, not to locate anything. The physics it carries — wave speeds, travel times, the way the bay smooths the source — is faithful; the cartography is impressionistic.

NoteWhat this isn’t

This is not a bathymetric chart. For real work, students should use AusBathyTopo (Geoscience Australia) or the GEBCO global grids — both are open and downloadable. The hand-built model here is the smallest one that gets the shape of the answer right; it takes the place a real bathymetry file would in production.

The bay model is produced by scripts/build_bay.jl and stored as plain CSVs in units/unit_01/data/:

using DelimitedFiles, CSV, DataFrames
H_grid = readdlm("units/unit_01/data/bay_bathymetry.csv", ',', Float64)
mask   = readdlm("units/unit_01/data/bay_mask.csv",       ',', Int)
gauges = CSV.read("units/unit_01/data/bay_gauges.csv",   DataFrame)
river  = CSV.read("units/unit_01/data/river_source.csv", DataFrame)

The CSVs are committed in the repo, so you have three options:

  1. Clone the whole course (recommended for hands-on work): git clone https://github.com/open-AIMS/Julia_PINN_training_2026
  2. Download individual files straight from GitHub — click any of the linked filenames above and use Raw → Save as.
  3. Regenerate them with julia --project=. units/unit_01/scripts/build_bay.jl from the repo root. The script is fast (a few seconds) and the only inputs are the constants at the top of the file, so you can change the grid, coastline, or gauge layout and rebuild.

Throughout the course, paths shown in backticks (e.g. scripts/build_bay.jl) are repo-relative; the same path also works as a GitHub link when rendered.

The physics in one paragraph

The water motion is described by the linearised shallow-water equations (de Wolff et al. 2021, eq. 6–7; reproduced in detail in §7.2). With H(x,y) the local depth, g gravity, and b a small linear drag, the surface elevation \eta(x, y, t) and depth-averaged horizontal velocity \mathbf{u}(x, y, t) = (u, v) satisfy

\underbrace{\frac{\partial \eta}{\partial t} + \nabla\!\cdot\!\bigl(H\,\mathbf{u}\bigr) = 0}_{\textbf{continuity (mass)}} \tag{1}

\underbrace{\frac{\partial \mathbf{u}}{\partial t} + g\,\nabla \eta = -b\,\mathbf{u}}_{\textbf{momentum}} \tag{2}

Those two equations combine into a single, more familiar one. If you differentiate the mass equation in time and substitute the momentum equation (the algebra is in the callout), and drop the small drag, the surface height \eta obeys a 2-D wave equation:

\frac{\partial^2 \eta}{\partial t^2} \;=\; \nabla\!\cdot\!\bigl(g\,H\,\nabla \eta\bigr). \tag{3}

In words: a bump in the water surface spreads outward as a gravity wave — gravity pulls the raised water back down, it overshoots, and so the disturbance travels. The one number that matters is its speed,

c(x, y) \;=\; \sqrt{g\,H(x, y)},

which depends only on the local depth H: waves travel faster where the water is deeper. In Moreton Bay that is roughly 8\,\text{m/s} over the shallow \sim 6\,\text{m} inner flats and about 15\,\text{m/s} along the deeper eastern trough — so the surge does not reach every gauge at the same time, and that difference in timing is a big part of what later lets us reconstruct the source. This wave equation is the physics the inverse solver below is told to obey.

Take \partial/\partial t of the continuity equation Equation 1 and substitute \partial \mathbf{u}/\partial t from the momentum equation Equation 2:

\frac{\partial^2 \eta}{\partial t^2} \;=\; -\nabla\!\cdot\!\Bigl(H\,\tfrac{\partial \mathbf{u}}{\partial t}\Bigr) \;=\; \nabla\!\cdot\!\bigl(g\,H\,\nabla \eta\bigr) \;+\; \nabla\!\cdot\!\bigl(H\,b\,\mathbf{u}\bigr).

The last term is the drag; it is small over the few-hour surge window, and dropping it leaves Equation 3. The restoring force here is gravity, not the fluid’s compressibility — that is what separates these surface gravity waves from acoustic (“pressure”) waves.

NoteWhy “linearised”?

We have dropped the Coriolis, advective (\mathbf{u}\!\cdot\!\nabla\mathbf{u}), and Reynolds-stress terms from the full SWE. For surge events of a few hours over a region \lesssim 100\,\text{km}, with no strong wind, the dominant balance is between pressure gradient and free-surface acceleration — exactly the same balance that gives a tsunami or tidal wave its speed. We will revisit when each dropped term matters in §7.2.

We close the system with a Dirichlet boundary condition at the river-mouth cell — a “Dirichlet” BC just means we pin the value of \eta at that point, as opposed to a “Neumann” BC which would pin its derivative:

\eta(x_r, y_r, t) \;=\; \psi(t),

where \psi(t) is the surface-elevation anomaly imposed by the river inflow. In the forward problem we prescribe \psi(t); in the inverse problem we will recover it from gauge data. Solid coasts get the no-flux condition \mathbf{u}\cdot\hat{\mathbf{n}}=0 (water doesn’t cross the coast).

The bay is not closed, though: it exchanges water with the open ocean through its wide northern entrance (the Northwest Channel), the strip of sea east of Moreton Island, and the southern Broadwater outlet. Modelling those as walls would bounce the surge straight back into the gauges. We don’t have a true radiation boundary condition here — instead we approximate one with a thin Rayleigh-damped sponge: a strip of cells along each of those three open rims where we add an artificial -\sigma\eta relaxation term, ramping up toward the edge, so outgoing surge waves are bled to zero before they reach the domain rim and reflect. (The western edge is solid mainland, so it stays a reflecting coast — which is physically correct.) A genuine open boundary is one of the things a fuller model would improve; here the sponge is enough to keep spurious reflections out of the synthetic gauge records.

A synthetic surge, propagated by finite differences

To generate the data we will work with, we run the linearised SWE forward on a staggered Arakawa-C grid — a layout where \eta lives at cell centres and the velocity components u, v live at the corresponding cell faces, which avoids the spurious checkerboard modes a co-located grid would suffer from — with explicit time stepping at \Delta t = 12\,\text{s} for 10 hours. The full driver is scripts/generate_surge_data.jl; Unit 6 covers staggered-grid FD schemes properly. The “ground-truth” surge is a two-pulse profile

\psi_{\text{truth}}(t) \;=\; 0.45 \exp\!\left[-\left(\tfrac{t - 2.0\,\text{h}}{0.55\,\text{h}}\right)^2\right] \;+\; 0.18 \exp\!\left[-\left(\tfrac{t - 4.3\,\text{h}}{0.55\,\text{h}}\right)^2\right]

— a \sim 45\,\text{cm} leading freshwater pulse from heavy catchment rain, followed by a \sim 18\,\text{cm} follow-on bulge about two hours later. These amplitudes are realistic for the de-tided residual of a moderate Brisbane River flood.

The inner FD step is just a discretisation of Equation 1 and Equation 2 — the listing below is what the code actually executes:

# Inner FD step (linearised SWE on an Arakawa-C grid).
# η lives at cell centres; u, v on the vertical/horizontal faces.
for step in 1:NT
    t = step * DT

    # 1. Mass continuity:  ∂η/∂t  =  -∂(Hu·u)/∂x  -  ∂(Hv·v)/∂y
    @inbounds for j in 2:NY-1, i in 2:NX-1
        mask[j, i] == 1 || continue
        flux_x_e = Hu[j, i  ] * u[j, i  ]
        flux_x_w = Hu[j, i-1] * u[j, i-1]
        flux_y_n = Hv[j,   i] * v[j,   i]
        flux_y_s = Hv[j-1, i] * v[j-1, i]
        η[j, i] -= DT * ((flux_x_e - flux_x_w) / DX +
                         (flux_y_n - flux_y_s) / DY)
    end

    # 2. Sponge layers on the open rims (north / east / south) absorb
    #    outgoing waves; `sponge` is pre-built with the per-cell damping rate.
    @inbounds for j in 1:NY, i in 1:NX
        sponge[j, i] > 0 && (η[j, i] -= DT * sponge[j, i] * η[j, i])
    end

    # 3. River-mouth source (Dirichlet on η at the source cell)
    η[JR, IR] = ψ_truth(t)

    # 4. Momentum (using the new η):  ∂u/∂t = -g ∂η/∂x - b·u
    @inbounds for j in 1:NY, i in 1:NX-1
        if u_open[j, i]
            dηdx = (η[j, i+1] - η[j, i]) / DX
            u[j, i] += DT * (-G_GRAV * dηdx - B_DRAG * u[j, i])
        else
            u[j, i] = 0.0
        end
    end
    # ... v update mirrors u, omitted for brevity ...
end

The full forward solve takes \sim 20\,\text{s} on a laptop at 500\,\text{m} resolution.

Watching the surge cross the bay

The movie below plays through the FD solution: 121 snapshots, 5 minutes of real time per step, total 10 hours. It loops automatically (click the image to pause or resume); watch the surge rise from rest, through the two peaks, and then the slow drain-out as the piled-up water leaks back to the ocean through the bay’s open rims, until by the 10-hour mark the bay has nearly returned to rest. The red triangle is the Brisbane River source, the gold discs are the four tide gauges, and grey is land.

Brisbane River surge propagating across Moreton Bay
t = 0.00 h · ψ(t) = +0.000 m
A slow movie of the surge — it loops automatically; click the image to pause or resume. The red triangle marks the unknown source ψ(t).

A few things to notice as it plays:

  1. Wave speed depends on depth. c = \sqrt{g\,H}. In the inner bay (H \approx 6\,\text{m}) that is about 8\,\text{m/s}; along the eastern deep trough you can see the wavefront noticeably faster.
  2. Geometry matters. The surge wraps around the Peel / Mud cluster and reflects off Moreton Island; the southern bay is slow to fill because of the Russell / Macleay / Karragarra constriction.
  3. Different gauges see different signals. G1 (Mud Island, ~18 km NNE of the mouth) responds first and strongest. G2 (Wellington Point, ~26 km ESE) sees a clean propagated wavefront across the central bay. G3 (Tangalooma, ~35 km NE) and G4 (Russell Channel, ~50 km south) see attenuated, delayed, reflection-laden signals. It’s the contrast between gauges, not any one of them, that constrains ψ(t).

The measurements we get

The four gauges produce timeseries \eta_g(t_k) at one-minute sampling. We add \sigma = 1.5\,\text{cm} Gaussian observation noise to imitate real de-tided tide-gauge residuals, then thin to the kind of cadence (\sim 3 min) that operational systems actually use.

Gauge Location Lat, Lon Cell Depth
G1 Mud Island -27.210°S, 153.215°E (39, 143) 6.0 m
G2 Wellington Point -27.480°S, 153.236°E (43, 84) 6.0 m
G3 Tangalooma Roads -27.205°S, 153.320°E (60, 144) 18.2 m
G4 Russell Channel -27.660°S, 153.290°E (60, 48) 6.0 m

These four timeseries, the bathymetry, and the physics are all the information the inverse solver gets. The river-mouth source \psi(t) is hidden.

NoteWhy not put G1 at Brisbane Bar?

The obvious choice for the closest gauge would be the Brisbane Bar / Pile Light tide gauge — it’s a real BoM station, right at the river entrance. We deliberately don’t use it: at ~1.5 km from the model’s river-mouth cell, Brisbane Bar essentially measures ψ(t) directly (travel time ~3 min, no smoothing), so the inversion from that gauge alone reduces to a one-to-one read-off and the broader four-gauge machinery looks like overkill.

Put another way: recovering \psi(t) from the four bay gauges is exactly the problem you would face if the Brisbane Bar gauge were broken or offline — reconstructing the near-source signal it would have measured from the other, more distant gauges. That sensor-outage scenario is the realistic one, and it is what makes the bay’s smoothing (and the whole four-gauge inversion) worth the trouble.

Mud Island (~18 km NNE of the mouth, real mid-bay tide-gauge site) is the closest gauge where the Green’s-function smoothing, geometric spreading, and reflection paths actually do something to the signal. With G1 at Mud Island, all four gauges are meaningful observers and the inverse-problem story is honest.

If you want to try the Brisbane-Bar version anyway, change the G1 row in scripts/build_bay.jl back to (-27.367, 153.166) and ./build.sh execute 1 to regenerate.

The inverse problem, two ways

We try two complementary approaches that bracket the spectrum covered by the rest of the course. You are not expected to read either listing in detail yet — they are sketches of what the inverse problem looks like once you have Units 2–9 in your toolbox.

What is a PINN, in one breath? A physics-informed neural network is just a small neural network used as a stand-in for the unknown field — here, the water height \eta(x, y, t) everywhere in the bay and at every moment. Rather than march the physics forward in time the way the simulator does, we train the network so that it does two things at once: (1) it matches the handful of gauge readings we actually have, and (2) it obeys the wave equation Equation 3 at thousands of random check-points scattered through space and time. That second part is the “physics-informed” bit — the equation is folded straight into the training objective, so the network is only allowed to settle on a water surface that is physically possible, not just any curve that threads the gauges. Once trained, the network has effectively filled in the whole bay, and we read off the river source simply by asking it for the height at the river-mouth cell, \psi_\theta(t) = \eta_\theta(x_r, y_r, t). Units 5 and 7 build this properly; for now it is enough to picture a neural network that is forced to respect the physics.

NoteWhat’s a collocation point?

Those “random check-points” have a name: collocation points. A collocation point is just a location in space and time (x, y, t) where we ask the network to obey the equation — nothing more. Unlike a gauge reading, there is no measured value there: we evaluate the wave-equation residual on the network’s own output and push it toward zero. Unit 1’s PINN scatters a few hundred of them at random through the water part of the bay (land excluded, and a small circle around the river mouth left out, since the source lives there), drawing a fresh batch every training step. They are essentially free — you can place as many as you like, anywhere — because checking the physics needs only the equation we already know, never data and never the true source. That is the mesh-free heart of a PINN: sample the physics wherever you want.

Method A — Naive PINN. Train a smooth MLP (multi-layer perceptron, the basic feed-forward neural network — covered in Unit 2) \eta_\theta(x, y, t) on inputs scaled to [-1, 1]^3 to match the gauge observations, satisfy the wave equation (Equation 3) at random water collocation points (the random space-time points where we evaluate the PDE residual as a soft constraint on the network — covered in Unit 5), and respect the zero initial condition. Treat the recovered source as \psi_\theta(t) := \eta_\theta(x_r, y_r, t). The only data the PINN ever sees are those same four noisy gauge time-series (g_{\text{obs}}, sub-sampled every few minutes) that the adjoint method uses — the wave equation, the zero initial condition and the smoothness preference are constraints, not data. Like the adjoint, it never sees the true source; the recovered \psi_\theta is just the trained network read off at the river-mouth cell. This is what de Wolff et al. (2021) classified as the “vanilla PINN” baseline and what §7.2 of this course re-derives in detail. The full training driver is scripts/train_inverse_pinn.jl.

# Sketch of the naive PINN loss — see scripts/train_inverse_pinn.jl.
function total_loss(θ, Xphys, Kx, Ky, Xz, Xt, Xsm)
    # Data: gauge readings ↔ η_θ at the gauge cells
    L_data = mean((η_batch(GAUGE_INPUT, θ) .- GAUGE_TARGET).^2)

    # Physics: ∂_ττ η  =  K_x ∂_ξξ η + K_y ∂_ζζ η   in normalised coords,
    #   sampled at random water-cell collocation points (excluding a small
    #   radius around the Dirichlet source so we don't fight the BC).
    η_p = reshape(η_batch(Xphys, θ), 7, :)
    η_ξξ = (η_p[2, :] .- 2 .* η_p[1, :] .+ η_p[3, :]) ./ H_NX^2
    η_ζζ = (η_p[4, :] .- 2 .* η_p[1, :] .+ η_p[5, :]) ./ H_NY^2
    η_ττ = (η_p[6, :] .- 2 .* η_p[1, :] .+ η_p[7, :]) ./ H_NT^2
    L_phys = mean((η_ττ .- Kx .* η_ξξ .- Ky .* η_ζζ) .^ 2)

    # IC + smoothness terms — see script
    return LAM_DATA*L_data + LAM_PHYS*L_phys + LAM_IC*L_ic + ...
end

Method B — the simulator as a forward map, run backwards. Here is the plain-English version before the equations. We already have a simulator of the bay: hand it a source signal \psi(t) at the river mouth and it predicts what each tide gauge will read. Call that the forward map — source in, gauge readings out. The inverse problem is to run it the other way: given the gauge readings, find the source that must have produced them.

Two facts make this clean. First, the bay physics is linear: double the source and every gauge reading doubles, and the response to two pushes is just the sum of the two separate responses. Second — and because of that — the entire simulator is captured by a single experiment: its response to one sharp ping at the source (a brief pulse). That ping-response is the bay’s Green’s function. Once you have it, the gauge reading for any source is just that ping-response slid along in time and added up — one copy for each moment of the source. So “predict the gauges” becomes a matrix multiply, and “find the source” becomes solving a linear system.

NoteThe trick: none of this used the unknown source

Notice what building G did not need: the true \psi. We poked the bay with a probe pulse we chose — a short, known blip — and recorded how the physics responded. That single experiment converts “I don’t know the source” into “I’ve measured how any source maps to the gauges.” The bay’s response to the real, unknown \psi is then simply G\psi — which is exactly the thing we now invert.

The only snag is that solving it exactly is unstable, and the cure — regularisation — is the idea worth dwelling on. The precise version:

Adjoint inverse, precisely. (The background machinery here — impulse response, Green’s function, convolution and the Toeplitz matrix that follow from the bay being a linear time-invariant system — is built up from scratch in Unit 6 §6.8; the L^{(2)} Tikhonov cost and how \lambda is chosen live in the regularisation subsection.) Because the linearised SWE is linear in \psi, the forward map G: \psi(\cdot) \mapsto \eta_g(\cdot) is a linear operator. Its columns are the Green’s function of the bay — the response at every gauge to a unit-impulse source at the river mouth. Compute it once (a single FD pass with a narrow Gaussian pulse), assemble G as a Toeplitz convolution matrix (the standard “discrete convolution as a matrix” trick: each row of G is the impulse response shifted in time by one sample), and solve the Tikhonov-regularised linear least-squares problem

\hat\psi \;=\; \arg\min_\psi \; \tfrac{1}{2}\,\|G\psi - g_{\text{obs}}\|_2^2 \;+\; \tfrac{\lambda^2}{2}\,\|\mathrm{L}^{(2)} \psi\|_2^2,

Despite the \arg\min notation, solving this is not like training the naive PINN: it has a one-step closed-form answer — a single linear solve (the A \ b in the code sketch below), with no neural network and no gradient descent at all. (The naive PINN of Method A minimises a similar-looking objective, but by gradient descent through a network; here we solve the same kind of problem exactly, in one shot.)

The first term says “the predicted gauges should match the observations.” The second term is the regularisation, and it is worth dwelling on, because every inverse problem in this course needs something like it. On its own, the first term is unstable: because the bay smooths and delays the signal, a wildly wobbly source — full of fast wiggles that nearly cancel out by the time they reach the gauges — can fit the noisy data just as well as the true, smooth one. Minimising the misfit alone therefore hands back a jagged, nonsensical \psi. The regularisation term cures this by adding a small penalty \lambda^2 \|\mathrm{L}^{(2)}\psi\|^2 on the curvature of \psi (the operator \mathrm{L}^{(2)} is just the discrete second derivative), which says: of all the sources that fit the gauges, prefer the smoothest one. The knob \lambda sets how hard we lean on that preference — too small and the wiggles return, too large and we flatten the real peaks.

This raises the obvious question: how do we set \lambda when we can’t peek at the true \psi to check which value works best? The trick is to lean on the one thing we do know — how noisy the gauges are (here about 1.5 cm). The reasoning is just common sense: a good reconstruction should match the gauges about as well as the noise allows — no better, no worse. Fitting them better than the noise means we are bending \psi to chase meaningless measurement jitter; fitting them worse means we have smoothed so hard that we have wiped out real signal. So the code tries a range of \lambda values and keeps the one whose leftover mismatch sits right at the noise level. (This rule has a name — the discrepancy principle — but the idea is entirely that common-sense balance, and it never once looks at the true source.)

This idea — fit the data, but prefer a simple answer — is exactly the “prior” the naive PINN instead supplies implicitly, through the smoothness of its network. No neural network here, but it leans on the same physics the PINN tries to learn. The full driver is scripts/train_inverse_adjoint.jl.

# Sketch of the adjoint inverse — see scripts/train_inverse_adjoint.jl.
h_g = fd_forward(narrow_gaussian_pulse)            # impulse response per gauge
G   = vcat([toeplitz_from_impulse(h_g[:, k]) for k in 1:NG]...)
A   = vcat(G, λ_smooth * L2, λ_anchor * e1')       # stacked LHS
b   = vcat(g_obs, zeros(size(L2, 1)), 0.0)
ψ̂   = A \ b                                        # closed-form Tikhonov
TipWhere — and only where — the true source is used

It is easy to suspect we are quietly peeking at the answer somewhere. We are not. Here is what each step actually feeds on:

Step What it uses Uses the true \psi?
Build G (the bay’s response) the simulator + a probe pulse we chose no
Solve for the source G and the observed gauges, plus “prefer a smooth answer” no
Choose \lambda the leftover mismatch vs how noisy the gauges are no
Score the result the recovered \hat\psi and the true \psi yes — only here

The true \psi shows up in exactly one place: the final accuracy number we print to grade ourselves, because this is a teaching example and we happen to have the ground truth to compare against. Delete that one line and everything above still runs — unchanged — on real gauge data where the source is genuinely unknown. That is the whole point: we reconstruct a source nobody measured, from its ripples at four gauges, the physics, and an honest sense of how noisy the measurements are.

What we recover

Top panel: the true ψ(t) (red) vs the adjoint inverse (blue, using the SWE solver as forward map) vs the naive PINN (orange dashed, smooth-MLP η-network with ψ read at the source cell). Bottom 2×2: predicted gauge timeseries at all four bay gauges — G1 (Mud Island) and G2 (Wellington Point) along the top, G3 (Tangalooma) and G4 (Russell Channel) along the bottom — overlaid on the noisy observations both methods see. Both methods fit the gauges, only the adjoint inverse recovers the source.

The plot is the centrepiece of Unit 1. Two observations matter pedagogically:

  1. Both methods fit the gauge observations. Look at any of the four gauge panels: the adjoint prediction (blue) and the naive-PINN prediction (dashed orange) both lie near or inside the noisy cloud of observations at G1–G4. By the metric the algorithms are optimising, both succeed.
  2. They disagree wildly on what produced those gauges. The adjoint inverse recovers both surge peaks with the right timing and two-pulse shape, but overshoots their amplitude (\sim 0.53\,\text{m} and \sim 0.22\,\text{m}, vs truth 0.45 / 0.18). The naive PINN gets the timing of the main pulse but flattens its amplitude to \sim 0.05\,\text{m} — almost an order of magnitude under.

This is the classic signature of an ill-posed inverse problem: many qualitatively different sources can produce the same downstream signal once it has been spread by the physics, smoothed by the gauges, and contaminated by noise. The remedy is to add a prior — the adjoint inverse adds it explicitly through Tikhonov smoothness; the naive PINN implicitly priors-in the smooth-MLP ansatz, which turns out to bias the answer in the wrong direction here.

A fair question: if the naive PINN does worse, why is this a PINN course? Because this problem is the classical method’s home turf. The bay here is linear, and for a linear inverse the Green’s-function + Tikhonov solve is essentially optimal — and almost free. Note that both methods saw the same four noisy gauges and the same physics; the adjoint wins not by seeing more, but because linearity lets it solve the problem exactly, in one shot, while the naive PINN takes on the much harder job of learning the entire field \eta(x,y,t) by gradient descent — with none of the tricks the rest of the course adds. We show that honestly, rather than rigging an easy PINN victory, so you can trust the comparison: on a problem this simple, reach for the classical tool. The interesting question is what happens when the problem stops being this simple.

de Wolff et al. (2021) documented exactly this failure mode on the same linearised SWE problem. §7.2–§7.3 walk through the modern fixes — Fourier feature embeddings, causal training, hard boundary-condition enforcement, adaptive loss weighting — that close the gap, so that by Unit 9 the PINN can do what the adjoint does here without the luxury of precomputing G.

TipSo why PINNs at all? Where the classical method can’t follow

We made this bay linear on purpose, so you could watch both methods on a level field. That linearity is exactly what hands the win to the Green’s-function solve — and exactly what real problems rarely give you. A PINN keeps working in the places the classical recipe breaks down:

  • Nonlinear surges. Once the water runs shallow (wetting/drying, currents that advect the wave), superposition fails — there is no impulse response to build G from. The PINN just changes the residual it minimises.
  • Unknown physics. Don’t know the bathymetry, or a friction coefficient? Make it a trainable parameter and recover it alongside the source from the same gauges — something the Green’s-function solve, which assumes the whole model, simply can’t do. (Unit 7 does exactly this.)
  • No simulator required. The PINN learns the field straight from data + the equation; it never needs a working forward solver in order to invert. The adjoint method is built from one.
  • Anywhere, any resolution. Being mesh-free, it returns \eta(x,y,t) at any continuous point, and copes with irregular or moving sensors.

Unit 1’s bay sits in the one corner where you don’t need any of that — which is precisely why it’s the honest place to start. The rest of the course lives in the corners where only the PINN can go.

The same solve, on a GPU

The forward solve above runs in a second or two at the bay’s native 100\times190 resolution. But the moment you want finer structure — sharper features around the islands, or a higher-resolution surge front — the cell count and the number of time steps both climb (the CFL condition ties the time step to the grid spacing), and the CPU solve slows to minutes.

This is the canonical case for a GPU. The shallow-water update is a stencil: every water cell’s new surface height depends only on its immediate neighbours, so thousands of cells can update at once. The script below is the same linearised SWE model, rewritten so every step is a whole-array broadcast instead of a scalar loop — which means the identical code runs on the CPU (Array) or the GPU (CuArray). It refines the real Moreton Bay bathymetry by an integer factor and races the two:

units/unit_01/scripts/surge_gpu.jl
#!/usr/bin/env julia
# ===========================================================================
# Unit 1 — the Moreton Bay surge solver, on the GPU.
#
# This is the SAME linearised shallow-water model as scripts/generate_surge_data.jl
# (Arakawa C-grid, forward-backward stepping, land mask, eastern sponge, a
# Dirichlet river-mouth source), but rewritten in a VECTORISED, backend-agnostic
# style: every update is a whole-array broadcast, no scalar loops. The exact same
# code therefore runs on the CPU (`Array`) or the GPU (`CuArray`) — the only
# difference is where the arrays live.
#
# Why bother? The reference solver runs a 100x190 grid in a few seconds on a CPU.
# Real bathymetric studies want metre-scale resolution over the whole bay, which
# is 10-100x more cells and 10-100x more time steps (CFL ties dt to dx). That is
# exactly the regime where a GPU's thousands of cores win: this script refines the
# bay by an integer factor and shows the CPU vs GPU wall-clock crossover.
#
# Run on the GPU hub (the @pinn env has CUDA.jl):
#   julia --project=@pinn units/unit_01/scripts/surge_gpu.jl
#
# It prints a benchmark table and writes figures/surge_gpu_field.png.
# Nothing here executes during `quarto render` — the .qmd shows it `eval: false`
# and includes the captured output.
# ===========================================================================

using Printf, Statistics, JSON3

# CUDA is optional: on the GPU hub it drives the real benchmark; on a CPU-only
# machine (e.g. a laptop regenerating just the figure) we skip it gracefully and
# the identical vectorised solver runs on `Array`. The field is the same either way.
const HAVE_CUDA = try
    @eval using CUDA
    true
catch
    @info "CUDA.jl not available — running CPU-only (benchmark GPU column blank; figure still rendered)"
    false
end
const HAVE_GPU = HAVE_CUDA && CUDA.functional()

# --- Bay geometry ----------------------------------------------------------
# Prefer the real Unit 1 bathymetry/mask if we can find them (on the hub the
# course repo is at /home/efs/_shared/course-materials); otherwise synthesise an
# idealised bay with the same character: a sloping basin, a shallower western
# shelf near the river mouth, a barrier-island gap to the east, land to the west.

function load_or_make_bay()
    candidates = String[]
    haskey(ENV, "GPU_DATA_DIR") && push!(candidates, ENV["GPU_DATA_DIR"])
    push!(candidates, joinpath(@__DIR__, "..", "data"))
    push!(candidates, "/home/efs/_shared/course-materials/units/unit_01/data")
    for d in candidates
        bf = joinpath(d, "bay_bathymetry.csv")
        mf = joinpath(d, "bay_mask.csv")
        if isfile(bf) && isfile(mf)
            H = _readcsv_float(bf)
            M = round.(Int, _readcsv_float(mf))
            H[isnan.(H)] .= 5.0
            src0 = _read_source(joinpath(d, "river_source.csv"), M)
            return (Float32.(H), M, src0, "real Moreton Bay bathymetry ($(size(H,1))x$(size(H,2)))")
        end
    end
    # --- idealised fallback (NY x NX = 100 x 190, same aspect as the real bay)
    NY, NX = 100, 190
    H = fill(5.0f0, NY, NX)
    M = ones(Int, NY, NX)
    @inbounds for j in 1:NY, i in 1:NX
        # depth deepens to the east; shallow shelf in the west near the source
        depth = 3.0f0 + 42.0f0 * (i / NX)^1.3f0
        H[j, i] = depth
        # western land wedge (mainland coast) + a barrier island band near i≈0.8NX
        west_land = i < (6 + 10 * sin(3.0 * j / NY))
        barrier   = (abs(i - round(Int, 0.82NX)) ≤ 1) && !(38 ≤ j ≤ 52)  # gap = tidal inlet
        if west_land || barrier
            M[j, i] = 0
            H[j, i] = 1.0f0
        end
    end
    jr0 = clamp(round(Int, 0.33NY), 2, NY-1)
    return (H, M, (jr0, _first_water_col(M, jr0)), "idealised bay (100x190)")
end

function _readcsv_float(path)
    rows = Vector{Vector{Float64}}()
    for ln in eachline(path)
        isempty(strip(ln)) && continue
        push!(rows, [(t == "NaN" || t == "nan" || isempty(t)) ? NaN : parse(Float64, t)
                     for t in split(ln, ',')])
    end
    return reduce(vcat, [permutedims(r) for r in rows])
end

_first_water_col(M, row) = (for i in 1:size(M, 2); M[row, i] == 1 && return max(i, 3); end; 3)

# The REAL Brisbane River-mouth cell (row = iy = north, col = ix = east) read from
# river_source.csv, so the GPU surge starts where the actual source is — not at a
# synthetic 1/3-up-the-coast guess. Falls back to that heuristic if the file is absent.
function _read_source(path, M)
    if isfile(path)
        lines = readlines(path)
        if length(lines) >= 2
            cols = split(strip(lines[2]), ',')
            if length(cols) >= 6
                ix = tryparse(Int, cols[5]); iy = tryparse(Int, cols[6])
                ix !== nothing && iy !== nothing && return (iy, ix)
            end
        end
    end
    jr = clamp(round(Int, 0.33size(M, 1)), 2, size(M, 1) - 1)
    return (jr, _first_water_col(M, jr))
end

# Nearest-neighbour integer refinement of a field (each cell -> r x r block).
refine_field(A, r) = r == 1 ? A : repeat(A, inner = (r, r))

# --- Build all the static arrays a solve needs, on a given backend ----------
# `dev` is `identity` for CPU or `CuArray` for GPU.
function build_state(Hc, Mc, src0, refine, dev)
    H  = Float32.(refine_field(Hc, refine))
    M  = refine_field(Mc, refine)
    NY, NX = size(H)
    maskf = Float32.(M)

    # face depths and open-face masks (water-water faces only)
    Hu = 0.5f0 .* (H[:, 1:NX-1] .+ H[:, 2:NX])
    Hv = 0.5f0 .* (H[1:NY-1, :] .+ H[2:NY, :])
    uopen = Float32.((M[:, 1:NX-1] .== 1) .& (M[:, 2:NX] .== 1))
    vopen = Float32.((M[1:NY-1, :] .== 1) .& (M[2:NY, :] .== 1))

    # absorbing sponge on every open rim — northern entrance, eastern ocean
    # strip, southern Broadwater outlet (west is solid mainland → no sponge).
    sponge = zeros(Float32, NY, NX)
    SPONGE_W = 5 * refine
    @inbounds for j in 1:NY, i in 1:NX
        de = i - (NX - SPONGE_W)
        dn = j - (NY - SPONGE_W)
        ds = (SPONGE_W + 1) - j
        s = 0.0f0
        de > 0 && (s = max(s, (1.0f0 / 60) * (de / SPONGE_W)^2))
        dn > 0 && (s = max(s, (1.0f0 / 60) * (dn / SPONGE_W)^2))
        ds > 0 && (s = max(s, (1.0f0 / 60) * (ds / SPONGE_W)^2))
        sponge[j, i] = s
    end
    spu = 0.5f0 .* (sponge[:, 1:NX-1] .+ sponge[:, 2:NX])
    spv = 0.5f0 .* (sponge[1:NY-1, :] .+ sponge[2:NY, :])

    # river-mouth source cell — the REAL mouth (src0 = base-grid row/col), mapped
    # into the refined grid and nudged onto water if the centre lands on land.
    jr0, ir0 = src0
    half = cld(refine, 2)
    jr = clamp((jr0 - 1) * refine + half, 2, NY - 1)
    ir = clamp((ir0 - 1) * refine + half, 2, NX - 1)
    if M[jr, ir] == 0
        best = (jr, ir); bestd = typemax(Int)
        @inbounds for jj in 1:NY, ii in 1:NX
            M[jj, ii] == 1 || continue
            d = (jj - jr)^2 + (ii - ir)^2
            d < bestd && (bestd = d; best = (jj, ii))
        end
        jr, ir = best
    end
    src = zeros(Float32, NY, NX); src[jr, ir] = 1.0f0

    return (; H, maskf, Hu, Hv, uopen, vopen, sponge, spu, spv,
            src = dev(src), NY, NX,
            maskf_d = dev(maskf), Hu_d = dev(Hu), Hv_d = dev(Hv),
            uopen_d = dev(uopen), vopen_d = dev(vopen),
            sponge_d = dev(sponge), spu_d = dev(spu), spv_d = dev(spv))
end

# Surge profile imposed at the source (two flood pulses), in metres.
psi(t) = 0.45f0 * exp(-((t - 2.0f0*3600) / (0.55f0*3600))^2) +
         0.18f0 * exp(-((t - 4.3f0*3600) / (0.55f0*3600))^2)

# --- One fully-vectorised time step (works for Array OR CuArray) ------------
const G_GRAV = 9.81f0
const B_DRAG = 5.0f-5

function swe_solve(Hc, Mc, src0, refine, dev; dx0 = 500.0f0, t_end = 3*3600.0f0,
                   nwarm = 5, capture::Int = 0)
    s = build_state(Hc, Mc, src0, refine, dev)
    NY, NX = s.NY, s.NX
    dx = dx0 / refine; dy = dx
    cmax = sqrt(G_GRAV * maximum(s.H))
    dt = 0.45f0 / (cmax * sqrt(1/dx^2 + 1/dy^2))        # CFL-safe
    nt = Int(floor(t_end / dt))
    frames = Matrix{Float32}[]; ftimes = Float64[]      # for the movie (capture>0)
    cap_every = capture > 0 ? max(1, nt ÷ capture) : typemax(Int)

    η = dev(zeros(Float32, NY, NX))
    u = dev(zeros(Float32, NY, NX-1))
    v = dev(zeros(Float32, NY-1, NX))
    divx = similar(η); divy = similar(η)
    fill!(divx, 0); fill!(divy, 0)

    function step!(tt)
        Fx = s.Hu_d .* u
        Fy = s.Hv_d .* v
        @views divx[:, 2:NX-1] .= (Fx[:, 2:NX-1] .- Fx[:, 1:NX-2]) ./ dx
        @views divy[2:NY-1, :] .= (Fy[2:NY-1, :] .- Fy[1:NY-2, :]) ./ dy
        η .-= dt .* (divx .+ divy) .* s.maskf_d          # continuity (water only)
        η .*= (1f0 .- dt .* s.sponge_d)                  # sponge on η
        η .= η .* (1f0 .- s.src) .+ (psi(tt) .* s.src)   # river-mouth Dirichlet
        dηdx = (η[:, 2:NX] .- η[:, 1:NX-1]) ./ dx
        dηdy = (η[2:NY, :] .- η[1:NY-1, :]) ./ dy
        u .= (u .+ dt .* (.-G_GRAV .* dηdx .- B_DRAG .* u)) .* s.uopen_d
        v .= (v .+ dt .* (.-G_GRAV .* dηdy .- B_DRAG .* v)) .* s.vopen_d
        u .*= (1f0 .- dt .* s.spu_d)                     # sponge on velocities
        v .*= (1f0 .- dt .* s.spv_d)
        return nothing
    end

    for n in 1:nwarm; step!(n*dt); end                   # warm up / compile
    fill!(η, 0); fill!(u, 0); fill!(v, 0)
    dev === identity || CUDA.synchronize()
    t0 = time()
    for n in 1:nt
        step!(n*dt)
        if capture > 0 && n % cap_every == 0
            push!(frames, Array(η)); push!(ftimes, n*dt)
        end
    end
    dev === identity || CUDA.synchronize()
    elapsed = time() - t0

    env = maximum(abs, Array(η))
    return (; field = Array(η), elapsed, nt, dt, NY, NX, ncells = NY*NX, env,
            frames, ftimes)
end

# ---------------------------------------------------------------------------
println("="^64)
println("Unit 1 — Moreton Bay shallow-water surge on CPU vs GPU")
println("="^64)
Hc, Mc, src0, src_desc = load_or_make_bay()
@printf("geometry: %s\n", src_desc)
@printf("GPU available: %s%s\n", HAVE_GPU,
        HAVE_GPU ? "  ($(CUDA.name(CUDA.device())), $(round(CUDA.totalmem(CUDA.device())/2^30; digits=1)) GiB)" : "")
println()

@printf("%-8s %-12s %-9s %-7s %10s %10s %9s\n",
        "refine", "grid", "cells", "steps", "CPU (s)", "GPU (s)", "speedup")
println("-"^72)

refines = [1, 2, 3, 4]
fine_field = nothing
for r in refines
    cpu = swe_solve(Hc, Mc, src0, r, identity)
    gpu = HAVE_GPU ? swe_solve(Hc, Mc, src0, r, CuArray) : nothing
    spd = gpu === nothing ? NaN : cpu.elapsed / gpu.elapsed
    @printf("%-8d %-12s %-9d %-7d %10.2f %10.2f %8.1fx\n",
            r, "$(cpu.NY)x$(cpu.NX)", cpu.ncells, cpu.nt,
            cpu.elapsed, gpu === nothing ? NaN : gpu.elapsed, spd)
    if gpu !== nothing
        @assert maximum(abs, cpu.field .- gpu.field) < 1f-2 "CPU/GPU fields diverged"
    end
    global fine_field = (gpu === nothing ? cpu : gpu).field
end

println()
println("CPU and GPU fields agree to < 1e-2 m at every resolution (same code, same physics).")

# --- Figure + movie: the surge on the finest grid --------------------------
# Rotated landscape, North ← left, land greyed, km axes — same look as the rest
# of Unit 1's bay maps (units/unit_01/scripts/_mapfig.jl). We re-run the finest
# grid once more (NOT timed) capturing frames, then render a static field PNG and
# a sequence of movie frames that unit_01.qmd plays back with a JS widget.
try
    include(joinpath(@__DIR__, "_mapfig.jl"))
    rfac      = last(refines)
    fine_mask = refine_field(Mc, rfac)
    dxkm      = 0.5 / rfac                       # 500 m base grid, refined ×rfac

    run  = swe_solve(Hc, Mc, src0, rfac, HAVE_GPU ? CuArray : identity;
                     t_end = 10*3600.0f0, capture = 60)  # movie runs to t = 10 h (matches the CPU walkthrough)
    NYf, NXf = size(run.field)
    vlim = max(0.05f0, 0.6f0 * maximum(maximum(abs, F) for F in run.frames))

    figdir = get(ENV, "GPU_FIG_DIR", joinpath(@__DIR__, "..", "figures"))
    isdir(figdir) || mkpath(figdir)

    # static field: the peak-surge frame — the final 10 h state has drained back
    # to rest, so the last frame would be a near-empty bay.
    pk = argmax([maximum(abs, F) for F in run.frames])
    p = bay_map(run.frames[pk], fine_mask, dxkm;
        clims = (-vlim, vlim), cmap = :balance, clabel = "η  (m)",
        title = @sprintf("Surge η near peak (t = %.1f h) — refined %d×%d grid",
                         run.ftimes[pk] / 3600, NYf, NXf))
    savefig(p, joinpath(figdir, "surge_gpu_field.png"))
    println("wrote figures/surge_gpu_field.png")

    # movie frames + metadata (played back by the widget in unit_01.qmd)
    moviedir = joinpath(figdir, "surge_gpu_frames")
    isdir(moviedir) || mkpath(moviedir)
    for f in readdir(moviedir; join = true); endswith(f, ".png") && rm(f); end
    recs = Dict{String, Any}[]
    for (k, (Fk, tk)) in enumerate(zip(run.frames, run.ftimes))
        pp = bay_map(Fk, fine_mask, dxkm;
            clims = (-vlim, vlim), cmap = :balance, clabel = "η  (m)",
            title = @sprintf("GPU surge (refined %d×%d) — t = %.2f h", NYf, NXf, tk/3600))
        fn = @sprintf("frame_%03d.png", k - 1)
        savefig(pp, joinpath(moviedir, fn))
        push!(recs, Dict("idx" => k - 1, "file" => fn, "t_hr" => tk/3600))
    end
    open(joinpath(moviedir, "frames_meta.json"), "w") do io
        JSON3.pretty(io, Dict("nframes" => length(recs),
                              "grid" => "$(NYf)x$(NXf)", "frames" => recs))
    end
    println("wrote ", length(recs), " GPU movie frames into figures/surge_gpu_frames/")
catch e
    println("(figure/movie skipped: ", e, ")")
end

Captured on the workshop GPU hub (NVIDIA A10G):

================================================================
Unit 1 — Moreton Bay shallow-water surge on CPU vs GPU
================================================================
geometry: real Moreton Bay bathymetry (190x100)
GPU available: true  (NVIDIA A10G, 22.0 GiB)

refine   grid         cells     steps      CPU (s)    GPU (s)   speedup
------------------------------------------------------------------------
1        190x100      19000     1462          1.77       0.37       4.7x
2        380x200      76000     2925          6.86       0.63      10.9x
3        570x300      171000    4388         26.99       1.39      19.5x
4        760x400      304000    5851         69.04       2.90      23.8x

CPU and GPU fields agree to < 1e-2 m at every resolution (same code, same physics).
wrote figures/surge_gpu_field.png

At the bay’s native grid the GPU barely helps — the problem is too small to fill the card. But each refinement multiplies the work, and by a 760\times400 grid the GPU is more than 20× faster while the CPU solve has passed a minute. The two never disagree by more than a centimetre: it is genuinely the same physics, just spread across thousands of cores. The surge at the finest resolution:

Surface elevation \eta near the peak of the synthetic surge (~2 h), on the refined 760\times400 Moreton Bay grid, computed on the GPU. Same landscape view as the other bay maps — north to the left (← N), east up — with land (islands, mainland coast) in grey; the surge has spread from the river mouth across the western bay.

And because it is the same solver stepping through time, we can watch it run. Here is the refined-grid surge as a movie — the same Moreton Bay as the CPU walkthrough above, just resolved four times finer in each direction:

GPU surge on the refined Moreton Bay grid
t = 0.00 h
The same shallow-water solver on the refined 760x400 grid, run on the GPU — it loops automatically; click to pause or resume.

What you have just seen

In one worked example we have touched, at an introductory depth, every theme of the course:

Theme Where it appeared
PDEs Linearised SWE, derived as a 2-D wave equation
Julia The whole forward + inverse stack, ~\!400 LOC
Computation Staggered-grid FD solve in 12\,\text{s}
Measurements Four noisy tide-gauge timeseries
Inverse problem Recover \psi(t) from gauge data
Ocean / water modelling Real geography (Moreton Bay), realistic surge amplitudes
PINNs Naive PINN baseline, its failure mode, the path to fix it

The scripts for the actual SWE solver, animation, and both inversions all live under units/unit_01/scripts/. Click any filename to jump to the source on GitHub:

Script Role Runtime
build_bay.jl bathymetry + land mask + gauge layout + bathymetry figure ~\!1 s
generate_site_map.py OSM/CartoDB site-map figure ~\!10 s
generate_surge_data.jl linearised-SWE forward solve, gauge timeseries, snapshots ~\!12 s
generate_surge_frames.jl per-snapshot PNG frames for the slider ~\!90 s
generate_surge_animation.jl optional GIF assembled from the surge frames ~\!30 s
train_inverse_pinn.jl naive PINN inverse (Lux.jl for the network, Zygote.jl for autodiff — both covered in Unit 2 / Unit 5) ~\!3 min
train_inverse_adjoint.jl Green’s-function / Tikhonov inverse ~\!15 s
render_recovery_plot.jl the ψ + 2×2 gauge comparison plot ~\!2 s

The whole pipeline is offline-reproducible on a single laptop CPU. Everything is written so it ports cleanly to a GPU once the rest of the course has motivated why we’d want one.

The shallow-water wave speed is c = \sqrt{g\,H}. Using the gauge table above (distances from the river mouth: G1 ≈ 18 km, G2 ≈ 26 km, G3 ≈ 35 km, G4 ≈ 50 km) and a representative depth along each path, estimate the arrival time of the leading surge peak at each gauge. Then watch the §1.2 movie and read off when each gauge first responds. Two follow-ups:

  • Which gauge’s predicted arrival is most sensitive to your choice of “representative depth”, and why does the bathymetry figure make that obvious?
  • The surge peaks at the mouth at t = 2 h. Roughly when should the G4 (Russell Channel) peak arrive, and why is the constant-depth estimate likely to be late or early there?

💡 Hint

Everything follows from c = \sqrt{gH} with g = 9.81: a 6 m inner-bay flat gives c ≈ 7.7 m/s, the ~25 m NE channel gives c ≈ 15 m/s. Arrival = distance / c (e.g. G1: 18 km / 7.7 ≈ 39 min), then add the t = 2 h source peak to get each gauge’s peak time. For the G3 path, run it twice — all-shallow vs all-channel — and quote the bracket rather than a single number (that spread is why G3 is the most depth-sensitive). For G4, ask what a straight-line constant-depth path leaves out: the detour around the constriction and the shallow southern flats both delay the real arrival, so the naive estimate lands early.

Go to solution →

1.3 The capstone preview: an AIMS thermistor column

The Moreton Bay surge is a 2-D shallow-water problem. The capstone across Units 8–10 keeps the inverse-problem flavour but moves the physics: a 1-D vertical column of ocean at three AIMS-monitored reef sites along the central Great Barrier Reef, with a thermistor chain sampling temperature at five depths every hour.

After a storm passes through, the deeper sensors record an unexpected cooling pattern. The competing hypotheses:

  1. The storm intensified upwelling — cold water pumped up from depth.
  2. The storm enhanced vertical mixing — heat drawn down from the surface faster than usual.
  3. Surface heating reduced — cloud cover and evaporation cut the net flux into the ocean.

Each leaves a distinct fingerprint in the depth–time record. The full model specification — Task A (single-site, CPU) and Task B (three-site joint, GPU) — lives in Unit 9; the worked solution that recovers which mechanism dominated is in Unit 10.

Same intellectual shape as today’s worked example — physics + sparse data → unknown driver — but with a different geometry, different physics, and a richer inverse-problem story.

Before reading any of the capstone units, commit to a prediction. For each of the three competing mechanisms (intensified upwelling, enhanced vertical mixing, reduced surface heating), sketch — on paper, axes z vs t — the pattern of temperature change you’d expect a five-sensor thermistor chain to record over the days after the storm. For each mechanism note: which sensors cool, which (if any) warm, and which respond first. Keep the sketch; you’ll grade it against Unit 10’s mechanism-partition plots.

💡 Hint

For each mechanism, ask which term of the column equation it perturbs: w\,\partial_z T (upwelling), \partial_z(\kappa\,\partial_z T) (mixing), or the surface flux/source (heating) — then where in the column that term bites hardest. The key discriminator: mixing only redistributes heat (it cools above the thermocline but warms below it, with column heat roughly conserved), while upwelling and reduced heating both remove heat, so no sensor warms. A warming deep sensor is the giveaway for mixing. For the other two, the tell is which end of the chain responds first in time.

Go to solution →

1.4 Course roadmap

Each remaining unit picks up one ingredient that today’s worked example assumed without explanation. In order:

  • Unit 2 — Machine-learning + deep-learning basics. The function-approximation toolkit (linear models, MLPs, autodiff, optimisers) you’ll need to even write a PINN.
  • Unit 3 — Scientific-ML landscape. A survey of physics-informed ML approaches, how PINNs fit in, and what other tools exist (Neural ODEs, operator learning, SINDy, …).
  • Unit 4 — ODEs via universal approximation. Neural ODEs and Universal Differential Equations (UDEs) — the gateway between classical dynamical systems and learned components.
  • Unit 5 — Your first PINN. Basic ODE and PDE models, the loss-as-physics-constraint pattern, the collocation-points machinery used above.
  • Unit 6 — Formal PDE theory + classical numerics. What a wave equation is, where finite differences come from, why the staggered Arakawa-C grid we used isn’t arbitrary.
  • Unit 7 — Modern PINN techniques for harder PDEs. Fourier features, causal training, hard BC enforcement, adaptive loss weighting — the fixes for the exact failure mode the naive PINN exhibited above.
  • Unit 8 — PDE modelling in key AIMS domains. The reef-scale ocean physics the capstone column rests on: the column idealisation, advection vs diffusion in fluids, the surface energy budget, stratification and mixing, dimensionless regimes.
  • Unit 9 — Capstone specification. Both versions of the project: Task A (single-site, CPU-friendly) and Task B (three-site joint inverse, GPU-class). Every equation, parameter, deliverable, and success criterion.
  • Unit 10 — Capstone solution. Worked solutions for both tasks, behind a triple-click reveal so you can attempt the project first. Closes the loop by recovering the unknown driver from sparse mooring data — exactly the workflow you saw in §1.2.

75 % Julia, 25 % Python

Julia is the primary language: Lux.jl, NeuralPDE.jl, MethodOfLines.jl, the broader SciML stack. Python parallels appear in Unit 2 (scikit-learn, PyTorch) and Unit 10 (DeepXDE) so participants leave with a sense of the cross-ecosystem landscape.