Porous Convection

This Porous Convection model is based on the equations, constants and mesh structure given by ETH Zurich's course on Solving Partial Differential Equations in Parallel on GPUs Lecture 4.

Dependencies

using CairoMakie
using Catlab
using CombinatorialSpaces
using ComponentArrays
using Decapodes
using DiagrammaticEquations
using Distributions
using GeometryBasics: Point2, Point3
using LinearAlgebra
using MLStyle
using OrdinaryDiffEq
using SparseArrays
using StaticArrays
using StatsBase

import CombinatorialSpaces.DiscreteExteriorCalculus: eval_constant_primal_form

The overall physics the model wishes to capture here are those of a fluid at different temperatures. Differences in the temperature of the fluid led to differences in the densities of that fluid which leads to convection of the fluid.

To drive the convection, the model has a cooling element at the top and a heating element at the bottom of our mesh. To avoid the need to determine the density of the fluid, the Boussinesq approximation is used.

The model generates a divergence-free pressure field by solving Poisson's equation and the Darcy flux is a combination of the forces caused by pressure differences and buoyant forces. This Darcy flux leads to advection of the fluid and along with some heat diffusion, the overall change in temperature is captured.

Finer details of the physics can be found at the source listed above.

Porous_Convection = @decapode begin
  (λ_ρ₀Cp, αρ₀, k_ηf, ϕ)::Constant
  (P, T, Adv, bound_T, bound_Ṫ)::Form0
  (g, qD)::Form1

  bound_T == adiabatic(T)

  # Darcy flux
  ρ == g ∧ (αρ₀ * bound_T)
  P == Δ⁻¹(δ(ρ))
  qD == -k_ηf * (d(P) - ρ)

  Adv == ⋆(interpolate(∧ᵈᵖ₁₁(⋆(d(bound_T)), qD)))
  Ṫ == -1/ϕ * Adv + λ_ρ₀Cp * Δ(bound_T)

  bound_Ṫ == tb_bc(Ṫ)

  ∂ₜ(T) == bound_Ṫ
end
infer_types!(Porous_Convection)
resolve_overloads!(Porous_Convection)
to_graphviz(Porous_Convection)

The Mesh

Our mesh will be a triangulated grid mesh of width 40 units and height 20. We ues the circumcenter subdivision method as the triangulated grid's triangles are well-behaved and thus we can enjoy faster solve times.

lx, ly = 40.0, 20.0
dx = dy = 0.4
s = triangulated_grid(lx, ly, dx, dy, Point3{Float64});
sd = EmbeddedDeltaDualComplex2D{Bool, Float64, Point2{Float64}}(s);
subdivide_duals!(sd, Circumcenter());

Operators and Boundary Conditions

We set up our custom operators first. Since our system is fairly small, we can directly factorize our Laplacian to make solve the Poisson Equation fast and accurate.

As you may have noticed, all of our data is located on the primal mesh, including our temperature. While this means we can keep the Darcy flux on the primal elements, the Lie derivative needs to compute the advection is slightly complicated now. However, we can employ interpolation to faithfully map data from primal to dual triangle elements to sidestep this problem. For fine enough meshes, the error in doing so in negligible.

Along with our operators, we implement adiabatic conditions on the side of our mesh as well as no-change conditions on the top/bottom of our mesh, where our cooling/heating elements reside. The adiabatic condition works by giving the boundary the values of the appropriate horizontal neighbor, which easy to do on the triangulated grid as its structure allows for easy lookup of this neighbor.

Δ0 = Δ(0,sd);
fΔ0 = LinearAlgebra.factorize(Δ0);

mat = p2_d2_interpolation(sd)

left_wall_idxs = findall(p -> p[1] < dx, s[:point]);
right_wall_idxs = findall(p -> p[1] > lx - dx, s[:point]);

# For adiabatic conditions:
next_left_wall_idxs = left_wall_idxs .+ 1;
next_right_wall_idxs = right_wall_idxs .- 1;

# For no-change conditions
bottom_wall_idxs= findall(p -> p[2] == 0, s[:point]);
top_wall_idxs = findall(p -> p[2] == ly, s[:point]);

apply_tb_bc(x) = begin x[bottom_wall_idxs] .= 0; x[top_wall_idxs] .= 0; return x; end

function generate(sd, my_symbol; hodge=GeometricHodge())
  op = @match my_symbol begin
    :Δ⁻¹ => x -> begin
      y = fΔ0 \ x
      # Constant changes in solution are valid
      y .-= minimum(y)
    end
    :adiabatic => x -> begin
      x[left_wall_idxs] .= x[next_left_wall_idxs]
      x[right_wall_idxs] .= x[next_right_wall_idxs]
      return x
    end
    :tb_bc => apply_tb_bc
    :interpolate => x -> mat * x
    _ => error("No operator $my_symbol found.")
  end
  return op
end

sim = eval(gensim(Porous_Convection))
f = sim(sd, generate, DiagonalHodge())

Initial Conditions

We set up a Gaussian heat disturbance at the center of our mesh to produce interesting convective behavior. We also set the cooling/heating elements.

ΔT = 200.0

T_dist = MvNormal([lx/2.0, ly/2.0], [1/sqrt(2), 1/sqrt(2)])
T = [2 * ΔT * pdf(T_dist, [p[1], p[2]]) for p in sd[:point]]
T[top_wall_idxs] .= -ΔT/2
T[bottom_wall_idxs] .= ΔT/2

Since we depend on gravity to drive the convection, as this is why lower densities rise against higher densities, we define gravity as a force moving in the downward y-direction and use eval_constant_primal_form to generate the appropriate 1-Form to represent this force.

Afterwards we establish a variety of physical constants.

# Gravity
accl_g = 9.81
grav = SVector{3}([0.0, -accl_g, 0.0])
g = eval_constant_primal_form(sd, grav)

# Physical constants
Ra = 750
k_ηf = 1.0
αρ₀ = (1.0/accl_g)
ϕ = 0.1
λ_ρ₀Cp = 1/Ra*(accl_g*αρ₀*k_ηf*ΔT*ly/ϕ)

u₀ = ComponentArray(T=T, g=g)
constants = (k_ηf = k_ηf, αρ₀ = αρ₀, ϕ = ϕ, λ_ρ₀Cp = λ_ρ₀Cp)

(k_ηf = 1.0, αρ₀ = 0.1019367991845056, ϕ = 0.1, λ_ρ₀Cp = 53.33333333333333)

Solving the Porous Convection Equations

We simply plug in our initial conditions and generate simulation code and solve it using Tsit5(). Below is an output of the full simulation.

tₑ = 0.7
prob = ODEProblem(f, u₀, (0, tₑ), constants)
soln = solve(prob, Tsit5(); saveat = 0.005)
soln.retcode

ReturnCode.Success = 1

save_dynamics("Porous_Convection.mp4", 120)

"Porous_Convection.mp4"

[ Info: Page built in 67 seconds.
[ Info: This page was last built at 2025-03-15T01:43:41.681.

Porous Convection

Dependencies