Square Bells
This example runs a small electron flow problem in the "square bells" geometry from (arXiv:2603.11175). We timestep the equations from a trivial initial condition until it reaches a steady state. Then we show how to save and visualize the solution.
To run this example, from the repo root run the following command which should take around 10 minutes on a laptop:
julia --project=. docs/src/tutorials/square_bells.jlIf needed, first add the script packages to the repo-root environment:
import Pkg
Pkg.add(["CairoMakie", "Interpolations", "OrdinaryDiffEqSSPRK"])First, load required packages.
using FermiSea
using CairoMakie
using HDF5
using Interpolations
using Logging
using OrdinaryDiffEqSSPRK
using StaticArrays
using TrixiNext, load the mesh. Trixi can load an external 2D Abaqus-style mesh with P4estMesh{2}. The boundary names in the mesh file are the symbols we use later in the boundary condition tuple.
mesh = P4estMesh{2}(
joinpath(@__DIR__, "..", "..", "..", "assets", "square_bells","square_bells.inp");
polydeg=1,
boundary_symbols=[:contact_bottom, :contact_top, :walls]
);Next, we choose which equations we want to solve. IsotropicFermiHarmonics2D keeps angular harmonics in the order a0, a1, b1, a2, b2, .... The first mode is the density-like mode, while a1 and b1 are proportional to the current. Let's keep 20 modes for now, which is plenty to get a good approximation in the hydrodynamic limit, though we want up to 200 if we have weaker collisions.
equations = IsotropicFermiHarmonics2D(20; v_fermi=1.0);Next let's add the collision term. A collision operator is just a Trixi source term. If you only need collisions, pass it directly as source_terms. If you want to combine collisions with a magnetic field source, wrap them with SourceTerms(collision, magnetic_source).
collisions = LinearCollisionMatrix(equations; gamma_mr=0.0, gamma_mc=100.0);as an initial condition, let's just start with zero everywhere.
initial_condition(x, t, equations) =
zero(SVector{nvariables(equations), typeof(equations.v_fermi)});Next, we specify our boundary conditions. CurrentContactBC(I) enforces a total integrated current I through the whole named contact by solving for one shared contact potential. MaxwellWallBC(1) gives fully diffuse walls.
boundary_conditions = (;
contact_bottom = CurrentContactBC(-0.1),
contact_top = CurrentContactBC(0.1),
walls = MaxwellWallBC(1.0),
);Next, we set up Trixi's solver, which is discontinuous Galerkin spectral element method (DGSEM), with a Lax-Friedrichs surface flux. The polydeg argument sets the polynomial degree for the DG method.
solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs);then build the descretized system:
semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver;
source_terms=collisions,
boundary_conditions=boundary_conditions);which then just reduces to an ODE problem that we can solve with any ODE solver in Julia, like this.
ode = semidiscretize(semi, (0.0, 100.0));It is often the case that we want to do things during the solve, like printing log messages, calcualting diagnostics, saving the solution, etc. These operations, just functions that we want to call at various points during the solve, are called "callbacks". Here we set up a few, with the most important being the SteadyStateCallback, which will stop the solve once the solution has reached a steady state with some tolerance.
callbacks = CallbackSet(
SummaryCallback(),
SteadyStateCallback(abstol=1e-2, reltol=0.0),
AnalysisCallback(semi; interval=2000, analysis_errors=Symbol[]),
StepsizeCallback(; cfl=0.5),
);finally, we just solve!
sol = solve(ode, SSPRK43(); adaptive=false, dt=1.0e-3, callback=callbacks,
save_everystep=false);
████████╗██████╗ ██╗██╗ ██╗██╗
╚══██╔══╝██╔══██╗██║╚██╗██╔╝██║
██║ ██████╔╝██║ ╚███╔╝ ██║
██║ ██╔══██╗██║ ██╔██╗ ██║
██║ ██║ ██║██║██╔╝ ██╗██║
╚═╝ ╚═╝ ╚═╝╚═╝╚═╝ ╚═╝╚═╝
┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ SemidiscretizationHyperbolic │
│ ════════════════════════════ │
│ #spatial dimensions: ……………………………………… 2 │
│ mesh: ……………………………………………………………………………… P4estMesh{2, 2, Float64} │
│ equations: ………………………………………………………………… IsotropicFermiHarmonics2D │
│ initial condition: …………………………………………… initial_condition │
│ boundary conditions: ……………………………………… 3 │
│ │ contact_bottom: ……………………………………………… FermiSea.CurrentContactBC{Float64} │
│ │ contact_top: ……………………………………………………… FermiSea.CurrentContactBC{Float64} │
│ │ walls: ……………………………………………………………………… FermiSea.MaxwellWallBC{Float64} │
│ source terms: ………………………………………………………… FermiSea.LinearCollisionMatri…0 0.0 0.0 0.0 0.0 0.0 100.0]) │
│ solver: ………………………………………………………………………… DG │
│ total #DOFs per field: ………………………………… 7648 │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘
┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ P4estMesh{2, 2, Float64} │
│ ════════════════════════ │
│ #trees: ………………………………………………………………………… 478 │
│ current #cells: …………………………………………………… 478 │
│ polydeg: ……………………………………………………………………… 1 │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘
┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ IsotropicFermiHarmonics2D │
│ ═════════════════════════ │
│ #variables: ……………………………………………………………… 41 │
│ │ variable 1: ………………………………………………………… a0 │
│ │ variable 2: ………………………………………………………… a1 │
│ │ variable 3: ………………………………………………………… b1 │
│ │ variable 4: ………………………………………………………… a2 │
│ │ variable 5: ………………………………………………………… b2 │
│ │ variable 6: ………………………………………………………… a3 │
│ │ variable 7: ………………………………………………………… b3 │
│ │ variable 8: ………………………………………………………… a4 │
│ │ variable 9: ………………………………………………………… b4 │
│ │ variable 10: ……………………………………………………… a5 │
│ │ variable 11: ……………………………………………………… b5 │
│ │ variable 12: ……………………………………………………… a6 │
│ │ variable 13: ……………………………………………………… b6 │
│ │ variable 14: ……………………………………………………… a7 │
│ │ variable 15: ……………………………………………………… b7 │
│ │ variable 16: ……………………………………………………… a8 │
│ │ variable 17: ……………………………………………………… b8 │
│ │ variable 18: ……………………………………………………… a9 │
│ │ variable 19: ……………………………………………………… b9 │
│ │ variable 20: ……………………………………………………… a10 │
│ │ variable 21: ……………………………………………………… b10 │
│ │ variable 22: ……………………………………………………… a11 │
│ │ variable 23: ……………………………………………………… b11 │
│ │ variable 24: ……………………………………………………… a12 │
│ │ variable 25: ……………………………………………………… b12 │
│ │ variable 26: ……………………………………………………… a13 │
│ │ variable 27: ……………………………………………………… b13 │
│ │ variable 28: ……………………………………………………… a14 │
│ │ variable 29: ……………………………………………………… b14 │
│ │ variable 30: ……………………………………………………… a15 │
│ │ variable 31: ……………………………………………………… b15 │
│ │ variable 32: ……………………………………………………… a16 │
│ │ variable 33: ……………………………………………………… b16 │
│ │ variable 34: ……………………………………………………… a17 │
│ │ variable 35: ……………………………………………………… b17 │
│ │ variable 36: ……………………………………………………… a18 │
│ │ variable 37: ……………………………………………………… b18 │
│ │ variable 38: ……………………………………………………… a19 │
│ │ variable 39: ……………………………………………………… b19 │
│ │ variable 40: ……………………………………………………… a20 │
│ │ variable 41: ……………………………………………………… b20 │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘
┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ DG{Float64} │
│ ═══════════ │
│ basis: …………………………………………………………………………… LobattoLegendreBasis{Float64}(polydeg=3) │
│ mortar: ………………………………………………………………………… LobattoLegendreMortarL2{Float64}(polydeg=3) │
│ surface integral: ……………………………………………… SurfaceIntegralWeakForm │
│ │ surface flux: …………………………………………………… FluxLaxFriedrichs(max_abs_speed) │
│ volume integral: ………………………………………………… VolumeIntegralWeakForm │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘
┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ SteadyStateCallback │
│ ═══════════════════ │
│ absolute tolerance: ………………………………………… 0.01 │
│ relative tolerance: ………………………………………… 0.0 │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘
┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ AnalysisCallback │
│ ════════════════ │
│ interval: …………………………………………………………………… 2000 │
│ analyzer: …………………………………………………………………… LobattoLegendreAnalyzer{Float64}(polydeg=6) │
│ │ integral 1: ………………………………………………………… FermiSea.SteadyStateResidual() │
│ save analysis to file: ………………………………… no │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘
┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ StepsizeCallback │
│ ════════════════ │
│ CFL Advective: ……………………………………………………… Returns{Float64}(0.5) │
│ CFL Diffusive: ……………………………………………………… Returns{Float64}(0.0) │
│ Interval: …………………………………………………………………… 1 │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘
┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ Time integration │
│ ════════════════ │
│ Start time: ……………………………………………………………… 0.0 │
│ Final time: ……………………………………………………………… 100.0 │
│ time integrator: ………………………………………………… SSPRK43 │
│ adaptive: …………………………………………………………………… false │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘
┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ Environment information │
│ ═══════════════════════ │
│ #threads: …………………………………………………………………… 1 │
│ threading backend: …………………………………………… polyester │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘
────────────────────────────────────────────────────────────────────────────────────────────────────
Simulation running 'IsotropicFermiHarmonics2D' with DGSEM(polydeg=3)
────────────────────────────────────────────────────────────────────────────────────────────────────
#timesteps: 0 run time: 6.79200000e-06 s
Δt: 1.00000000e-03 └── GC time: 0.00000000e+00 s (0.000%)
sim. time: 0.00000000e+00 (0.000%) time/DOF/rhs!: NaN s
PID: Inf s
#DOFs per field: 7648 alloc'd memory: 617.269 MiB
#elements: 478
steady state: 3.95031195e+02
────────────────────────────────────────────────────────────────────────────────────────────────────
────────────────────────────────────────────────────────────────────────────────────────────────────
Simulation running 'IsotropicFermiHarmonics2D' with DGSEM(polydeg=3)
────────────────────────────────────────────────────────────────────────────────────────────────────
#timesteps: 2000 run time: 8.51244396e+01 s
Δt: 9.30743933e-04 └── GC time: 1.08920667e+01 s (12.795%)
sim. time: 1.86148787e+00 (1.861%) time/DOF/rhs!: 1.35787102e-06 s
PID: 1.36758066e-06 s
#DOFs per field: 7648 alloc'd memory: 172.548 MiB
#elements: 478
steady state: 1.83668766e+00
────────────────────────────────────────────────────────────────────────────────────────────────────
────────────────────────────────────────────────────────────────────────────────────────────────────
Simulation running 'IsotropicFermiHarmonics2D' with DGSEM(polydeg=3)
────────────────────────────────────────────────────────────────────────────────────────────────────
#timesteps: 4000 run time: 1.67981876e+02 s
Δt: 9.30743933e-04 └── GC time: 2.13986610e+01 s (12.739%)
sim. time: 3.72297573e+00 (3.723%) time/DOF/rhs!: 1.32658067e-06 s
PID: 1.35404522e-06 s
#DOFs per field: 7648 alloc'd memory: 187.061 MiB
#elements: 478
steady state: 1.10197483e+00
────────────────────────────────────────────────────────────────────────────────────────────────────
────────────────────────────────────────────────────────────────────────────────────────────────────
Simulation running 'IsotropicFermiHarmonics2D' with DGSEM(polydeg=3)
────────────────────────────────────────────────────────────────────────────────────────────────────
#timesteps: 6000 run time: 2.51128380e+02 s
Δt: 9.30743933e-04 └── GC time: 3.22621243e+01 s (12.847%)
sim. time: 5.58446360e+00 (5.584%) time/DOF/rhs!: 1.33148898e-06 s
PID: 1.35878285e-06 s
#DOFs per field: 7648 alloc'd memory: 170.369 MiB
#elements: 478
steady state: 4.20215417e-01
────────────────────────────────────────────────────────────────────────────────────────────────────
────────────────────────────────────────────────────────────────────────────────────────────────────
Simulation running 'IsotropicFermiHarmonics2D' with DGSEM(polydeg=3)
────────────────────────────────────────────────────────────────────────────────────────────────────
#timesteps: 8000 run time: 3.34734739e+02 s
Δt: 9.30743933e-04 └── GC time: 4.30750828e+01 s (12.868%)
sim. time: 7.44595147e+00 (7.446%) time/DOF/rhs!: 1.33916190e-06 s
PID: 1.36629301e-06 s
#DOFs per field: 7648 alloc'd memory: 206.290 MiB
#elements: 478
steady state: 1.83975763e-01
────────────────────────────────────────────────────────────────────────────────────────────────────
────────────────────────────────────────────────────────────────────────────────────────────────────
Simulation running 'IsotropicFermiHarmonics2D' with DGSEM(polydeg=3)
────────────────────────────────────────────────────────────────────────────────────────────────────
#timesteps: 10000 run time: 4.18532220e+02 s
Δt: 9.30743933e-04 └── GC time: 5.40895639e+01 s (12.924%)
sim. time: 9.30743933e+00 (9.307%) time/DOF/rhs!: 1.34220617e-06 s
PID: 1.36940582e-06 s
#DOFs per field: 7648 alloc'd memory: 201.157 MiB
#elements: 478
steady state: 1.82935443e-01
────────────────────────────────────────────────────────────────────────────────────────────────────
────────────────────────────────────────────────────────────────────────────────────────────────────
Simulation running 'IsotropicFermiHarmonics2D' with DGSEM(polydeg=3)
────────────────────────────────────────────────────────────────────────────────────────────────────
#timesteps: 12000 run time: 5.02519386e+02 s
Δt: 9.30743933e-04 └── GC time: 6.51064065e+01 s (12.956%)
sim. time: 1.11689272e+01 (11.169%) time/DOF/rhs!: 1.34514697e-06 s
PID: 1.37251517e-06 s
#DOFs per field: 7648 alloc'd memory: 202.972 MiB
#elements: 478
steady state: 1.06104563e-01
────────────────────────────────────────────────────────────────────────────────────────────────────
────────────────────────────────────────────────────────────────────────────────────────────────────
Simulation running 'IsotropicFermiHarmonics2D' with DGSEM(polydeg=3)
────────────────────────────────────────────────────────────────────────────────────────────────────
#timesteps: 14000 run time: 5.86250552e+02 s
Δt: 9.30743933e-04 └── GC time: 7.61357352e+01 s (12.987%)
sim. time: 1.30304151e+01 (13.030%) time/DOF/rhs!: 1.34102201e-06 s
PID: 1.36832090e-06 s
#DOFs per field: 7648 alloc'd memory: 190.450 MiB
#elements: 478
steady state: 4.01763516e-02
────────────────────────────────────────────────────────────────────────────────────────────────────
────────────────────────────────────────────────────────────────────────────────────────────────────
Simulation running 'IsotropicFermiHarmonics2D' with DGSEM(polydeg=3)
────────────────────────────────────────────────────────────────────────────────────────────────────
#timesteps: 16000 run time: 6.69441244e+02 s
Δt: 9.30743933e-04 └── GC time: 8.72602285e+01 s (13.035%)
sim. time: 1.48919029e+01 (14.892%) time/DOF/rhs!: 1.33233520e-06 s
PID: 1.35949907e-06 s
#DOFs per field: 7648 alloc'd memory: 186.580 MiB
#elements: 478
steady state: 2.53491105e-02
────────────────────────────────────────────────────────────────────────────────────────────────────
────────────────────────────────────────────────────────────────────────────────────────────────────
Simulation running 'IsotropicFermiHarmonics2D' with DGSEM(polydeg=3)
────────────────────────────────────────────────────────────────────────────────────────────────────
#timesteps: 18000 run time: 7.52662810e+02 s
Δt: 9.30743933e-04 └── GC time: 9.84505173e+01 s (13.080%)
sim. time: 1.67533908e+01 (16.753%) time/DOF/rhs!: 1.33274042e-06 s
PID: 1.36000817e-06 s
#DOFs per field: 7648 alloc'd memory: 172.857 MiB
#elements: 478
steady state: 2.26718041e-02
────────────────────────────────────────────────────────────────────────────────────────────────────
┌ Info: Steady state tolerance reached
│ steady_state_callback = Trixi.SteadyStateCallback{Float64}(0.01, 0.0)
└ t = 17.746494573641414
────────────────────────────────────────────────────────────────────────────────────────────────────
Simulation running 'IsotropicFermiHarmonics2D' with DGSEM(polydeg=3)
────────────────────────────────────────────────────────────────────────────────────────────────────
#timesteps: 19067 run time: 7.97700896e+02 s
Δt: 9.30743933e-04 └── GC time: 1.04361815e+02 s (13.083%)
sim. time: 1.77464946e+01 (17.746%) time/DOF/rhs!: 1.33378146e-06 s
PID: 1.37942730e-06 s
#DOFs per field: 7648 alloc'd memory: 212.449 MiB
#elements: 478
steady state: 9.99792111e-03
────────────────────────────────────────────────────────────────────────────────────────────────────
─────────────────────────────────────────────────────────────────────────────────
Trixi.jl Time Allocations
─────────────────────── ────────────────────────
Tot / % measured: 798s / 97.9% 2.93TiB / 98.5%
Section ncalls time %tot avg alloc %tot avg
─────────────────────────────────────────────────────────────────────────────────
rhs! 76.3k 779s 99.7% 10.2ms 2.88TiB 100.0% 39.6MiB
boundary flux 76.3k 357s 45.6% 4.68ms 1.94TiB 67.2% 26.6MiB
source terms 76.3k 228s 29.2% 2.99ms 435GiB 14.7% 5.83MiB
volume integral 76.3k 86.0s 11.0% 1.13ms 365GiB 12.4% 4.90MiB
interface flux 76.3k 66.4s 8.5% 871μs 168GiB 5.7% 2.26MiB
prolong2interfaces 76.3k 19.8s 2.5% 259μs 0.00B 0.0% 0.00B
surface integral 76.3k 11.0s 1.4% 145μs 0.00B 0.0% 0.00B
Jacobian 76.3k 4.88s 0.6% 64.0μs 0.00B 0.0% 0.00B
reset ∂u/∂t 76.3k 4.46s 0.6% 58.4μs 0.00B 0.0% 0.00B
prolong2boundaries 76.3k 1.78s 0.2% 23.3μs 0.00B 0.0% 0.00B
~rhs!~ 76.3k 204ms 0.0% 2.68μs 4.78KiB 0.0% 0.06B
prolong2mortars 76.3k 14.8ms 0.0% 194ns 0.00B 0.0% 0.00B
mortar flux 76.3k 12.1ms 0.0% 158ns 0.00B 0.0% 0.00B
analyze solution 11 1.71s 0.2% 156ms 661MiB 0.0% 60.1MiB
calculate dt 19.1k 260ms 0.0% 13.6μs 0.00B 0.0% 0.00B
─────────────────────────────────────────────────────────────────────────────────
For visualization it is convenient to sample the DG solution on a Cartesian grid. The package can write that grid to HDF5; this is also a handy format for post-processing in Python, MATLAB, ParaView-adjacent scripts, or another Julia session.
cartesian_file = joinpath(@__DIR__, "square_bells_cartesian.h5");
with_logger(NullLogger()) do
save_cartesian(sol, semi, cartesian_file; nvisnodes=200)
end;then we read it as follows and visualize in Makie
x, y, a1, b1, mask = h5open(cartesian_file, "r") do file
read(file["x"]), read(file["y"]), read(file["a1"]), read(file["b1"]),
read(file["mask"])
end;
current_magnitude = sqrt.(a1.^2 .+ b1.^2);
current_magnitude[.!mask] .= NaN;
stream_x = a1;
stream_y = b1;
stream_x[.!mask] .= 0;
stream_y[.!mask] .= 0;
stream_x_itp = linear_interpolation((x, y), stream_x; extrapolation_bc=0.0);
stream_y_itp = linear_interpolation((x, y), stream_y; extrapolation_bc=0.0);
current_vector(px, py) =
Point2f(stream_x_itp(px, py), stream_y_itp(px, py));
fig = Figure();
ax = Axis(fig[1, 1]; aspect=DataAspect(), xlabel="x", ylabel="y",
title="Current streamlines");
hm = heatmap!(ax, x, y, current_magnitude; colormap=:viridis);
x_interval = CairoMakie.Makie.IntervalSets.ClosedInterval(minimum(x), maximum(x));
y_interval = CairoMakie.Makie.IntervalSets.ClosedInterval(minimum(y), maximum(y));
streamplot!(ax, current_vector, x_interval, y_interval;
color=Returns(:white), linewidth=1.2, density=0.8,
gridsize=(28, 28), arrow_size=8);
Colorbar(fig[1, 2], hm; label="|j|");
fig
This page was generated using Literate.jl.