add a DG SW model exemple by yolhan83 · Pull Request #1283 · Ferrite-FEM/Ferrite.jl

yolhan83 · 2026-02-11T13:55:46Z

This PR adds a shallow water RKDG exemple to showcase how this can be done in Ferrite.
I won't go further on this description since everything is in the md, please suggest anykind of clarification.

yolhan83 · 2026-02-11T14:25:47Z

Manifest should be in the gitignore ?

yolhan83 · 2026-02-11T15:07:08Z

Ok for review, sorry some tutos didnt work on my pc couldn't test doc locally, only could preview. Seems to work now

termi-official · 2026-02-11T16:28:36Z

Really awesome, thanks for the addition!

Manifest should be in the gitignore ?

We have it in the CI to have a fully reproducible state for the docs, as we rely on a lot of packages in the docs. I will review later.

codecov · 2026-02-11T16:35:54Z

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 94.25%. Comparing base (f245839) to head (87cfa1e).

Additional details and impacted files

@@           Coverage Diff           @@
##           master    #1283   +/-   ##
=======================================
  Coverage   94.25%   94.25%           
=======================================
  Files          40       40           
  Lines        6750     6750           
=======================================
  Hits         6362     6362           
  Misses        388      388

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

termi-official

First feedback from my side.

[...] some tutos didnt work on my pc [...]

What failure do you get?

docs/src/literate-tutorials/dg_shallow_water.jl

yolhan83 · 2026-02-12T06:55:09Z

First feedback from my side.

[...] some tutos didnt work on my pc [...]

What failure do you get?

I get :

┌ Error: failed to run `@example` block in docs\src\tutorials\stokes-flow.md:459-461
│ ```@example stokes-flow
│ main()
│ ```
│   exception =
│ 
│    Stacktrace:
│      [1] error(s::String)
│        @ Base .\error.jl:44
│      [2] open
│        @ C:\Users\yolha\.julia\artifacts\649584fb27316a4395dbb47fde904c83626feae5\lib\gmsh.jl:113 [inlined]
│      [3] togrid(filename::String; domain::String)
│        @ FerriteGmsh C:\Users\yolha\.julia\packages\FerriteGmsh\MH6r7\src\FerriteGmsh.jl:185
│      [4] togrid
│        @ C:\Users\yolha\.julia\packages\FerriteGmsh\MH6r7\src\FerriteGmsh.jl:176 [inlined]
│      [5] (::Main.var"__atexample__named__stokes-flow".var"#setup_grid##0#setup_grid##1")(dir::String)
│        @ Main.var"__atexample__named__stokes-flow" .\stokes-flow.md:213
│      [6] mktempdir(fn::Main.var"__atexample__named__stokes-flow".var"#setup_grid##0#setup_grid##1", parent::String; prefix::String)
│        @ Base.Filesystem .\file.jl:936
│      [7] mktempdir (repeats 2 times)
│        @ .\file.jl:932 [inlined]
│      [8] setup_grid(h::Float64)
│        @ Main.var"__atexample__named__stokes-flow" .\stokes-flow.md:210
│      [9] main()
│        @ Main.var"__atexample__named__stokes-flow" .\stokes-flow.md:462
│     [10] top-level scope
│        @ stokes-flow.md:460
│     [11] eval(m::Module, e::Any)
│        @ Core .\boot.jl:489
│     [12] #61
│        @ C:\Users\yolha\.julia\packages\Documenter\HdXI4\src\expander_pipeline.jl:856 [inlined]
│     [13] cd(f::Documenter.var"#61#62"{Module, Expr}, dir::String)
│        @ Base.Filesystem .\file.jl:101
│     [14] (::Documenter.var"#59#60"{Documenter.Page, Module, Expr})()
│        @ Documenter C:\Users\yolha\.julia\packages\Documenter\HdXI4\src\expander_pipeline.jl:855
│     [15] (::IOCapture.var"#12#13"{Type{InterruptException}, Documenter.var"#59#60"{Documenter.Page, Module, Expr}, IOContext{Base.PipeEndpoint}, IOContext{Base.PipeEndpoint}, Base.TTY, Base.TTY})()
│        @ IOCapture C:\Users\yolha\.julia\packages\IOCapture\Y5rEA\src\IOCapture.jl:170
│     [16] with_logstate(f::IOCapture.var"#12#13"{Type{InterruptException}, Documenter.var"#59#60"{Documenter.Page, Module, Expr}, IOContext{Base.PipeEndpoint}, IOContext{Base.PipeEndpoint}, Base.TTY, Base.TTY}, logstate::Base.CoreLogging.LogState)  
│        @ Base.CoreLogging .\logging\logging.jl:542
│     [17] with_logger(f::Function, logger::Base.CoreLogging.ConsoleLogger)
│        @ Base.CoreLogging .\logging\logging.jl:653
│     [18] capture(f::Documenter.var"#59#60"{Documenter.Page, Module, Expr}; rethrow::Type, color::Bool, passthrough::Bool, capture_buffer::IOBuffer, io_context::Vector{Any})
│        @ IOCapture C:\Users\yolha\.julia\packages\IOCapture\Y5rEA\src\IOCapture.jl:167
│     [19] runner(::Type{Documenter.Expanders.ExampleBlocks}, node::MarkdownAST.Node{Nothing}, page::Documenter.Page, doc::Documenter.Document)
│        @ Documenter C:\Users\yolha\.julia\packages\Documenter\HdXI4\src\expander_pipeline.jl:854
└ @ Documenter C:\Users\yolha\.julia\packages\Documenter\HdXI4\src\utilities\utilities.jl:47

termi-official · 2026-02-12T11:23:54Z

docs/src/literate-tutorials/dg_shallow_water.jl

+    Ui = @view U[idofs]
+    for qp in 1:getnquadpoints(iv)
+        n = getnormal(iv, qp; here = true)
+        dF = getdetJdV(iv, qp)
+        hL = function_value(iv, qp, Ui, dr_h, dr_h .+ ndpc; here = true)
+        qxL = function_value(iv, qp, Ui, dr_qx, dr_qx .+ ndpc; here = true)
+        qyL = function_value(iv, qp, Ui, dr_qy, dr_qy .+ ndpc; here = true)
+        hR = function_value(iv, qp, Ui, dr_h, dr_h .+ ndpc; here = false)
+        qxR = function_value(iv, qp, Ui, dr_qx, dr_qx .+ ndpc; here = false)
+        qyR = function_value(iv, qp, Ui, dr_qy, dr_qy .+ ndpc; here = false)
+        UL = Vec{3}((hL, qxL, qyL))
+        UR = Vec{3}((hR, qxR, qyR))
+        fhat = LF(UL, UR, n)
+        for i in 1:nb
+            ΦᵢL = shape_value(iv, qp, i; here = true)
+            ΦᵢR = shape_value(iv, qp, i + nb; here = false)
+            Ri[dr_h[i]] += (-ΦᵢL * fhat[1]) * dF
+            Ri[dr_qx[i]] += (-ΦᵢL * fhat[2]) * dF
+            Ri[dr_qy[i]] += (-ΦᵢL * fhat[3]) * dF
+            Ri[dr_h[i] + ndpc] += (ΦᵢR * fhat[1]) * dF
+            Ri[dr_qx[i] + ndpc] += (ΦᵢR * fhat[2]) * dF
+            Ri[dr_qy[i] + ndpc] += (ΦᵢR * fhat[3]) * dF
+        end
+    end


@AbdAlazezAhmed can you take a look at this part here? Any chance that we miss some interface for this, so users do not need to manually index this stuff? Also, aren't we missing some 1/h term?

I found that weird using it too especially, I thought it should be either

hL = function_value(iv, qp, Ui, dr_h; here=true) hR = function_value(iv, qp, Ui, dr_h; here=false)

or

hL = function_value(iv, qp, Ui, dr_h) hR = function_value(iv, qp, Ui, dr_h .+ ndpc)

but having to do both seems weird. Maybe differencing Ui_here and Ui_there since the beggening would be better ?

docs/src/literate-tutorials/dg_shallow_water.jl

termi-official · 2026-02-12T11:34:18Z

docs/src/literate-tutorials/dg_shallow_water.jl

+            Rb[dr_h[i]] -= (Φᵢ * fhat[1]) * dF
+            Rb[dr_qx[i]] -= (Φᵢ * fhat[2]) * dF
+            Rb[dr_qy[i]] -= (Φᵢ * fhat[3]) * dF


Correct me if I am wrong, but don't we need some element size term here to enforce the Dirichlet condition weakly?

I'm using a ghost element method for the week forcing at the boundary, so the code is exactly similar to the inner interface beside that we dont compute the dofs for the ghost element and we preomcute what we want or wish at the outer trace. Maybe you're thinking of a strong one with a forcing term of order 1/size but that can lead to non-physical solution in this framework. It is ok when DG is used on parabolic / eliptic pdes though.

Interesting, thanks for sharing this. Do you know some paper analyzing this? I have not read the theory part at the top yet, bit if it is not explained there yet, could you add the reference there?

The clearest I know is this one :
The Runge–Kutta discontinuous Galerkin method for conservation laws V: multidimensional systems](https://www.sciencedirect.com/science/article/pii/S0021999198958922
by B Cockburn and CW Shu its quite old 1998 but its just too well explained.

Of course can't do that without citing a friend who did his phd at the same time as I did,
Poussel, Camille, Mehmet Ersoy, and Frederic Golay. "Wetting and drying treatments with mesh adaptation for shallow water equations using a Runge–Kutta discontinuous Galerkin method." International Journal for Numerical Methods in Fluids 97.5 (2025): 692-712.

won't cite myself since I'm more on the modelling side but I used this method for every papers I did (with a twist for advection-diffusion problem)

also if you want the theoric part to be clearer I can do that but I thought I would go really fast to the code since the point stay Ferrite

I always struggle with a good balance for how much theory should be stated in the docs. @AbdAlazezAhmed what would you put in the docs for this example?

We also might want to put this tutorial before the Navier-Stokes tutorial, because it shows OrdinaryDiffEq withouth the funky parts to strongly enforce the Dirichlet boundary condition and discontinuity handling.

docs/src/tutorials/index.md

yolhan83 · 2026-02-12T16:11:14Z

comment for myself : if docs pass, remove the vtk gen froom CI

docs/src/literate-tutorials/dg_shallow_water.jl

termi-official · 2026-02-12T18:27:34Z

docs/src/literate-tutorials/dg_shallow_water.jl

+            Rb[dr_h[i]] -= (Φᵢ * fhat[1]) * dF
+            Rb[dr_qx[i]] -= (Φᵢ * fhat[2]) * dF
+            Rb[dr_qy[i]] -= (Φᵢ * fhat[3]) * dF


I always struggle with a good balance for how much theory should be stated in the docs. @AbdAlazezAhmed what would you put in the docs for this example?

docs/src/literate-tutorials/dg_shallow_water.jl

termi-official · 2026-02-12T18:32:10Z

docs/src/literate-tutorials/dg_shallow_water.jl

+            Rb[dr_h[i]] -= (Φᵢ * fhat[1]) * dF
+            Rb[dr_qx[i]] -= (Φᵢ * fhat[2]) * dF
+            Rb[dr_qy[i]] -= (Φᵢ * fhat[3]) * dF


We also might want to put this tutorial before the Navier-Stokes tutorial, because it shows OrdinaryDiffEq withouth the funky parts to strongly enforce the Dirichlet boundary condition and discontinuity handling.

docs/src/tutorials/index.md

yolhan83 · 2026-02-12T20:58:59Z

I think its ok for review again, besides waiting for some answer on the amount of math.
little note too : the first lines aren't hide in preview yet this is how litterate should handle this ?

lijas

Nice work! I dont know anything about the shallow water equations, but from reading the introduction it seems like the fields are grouped in to a vector U = [h, qx, qy] ? In that case, I think it would be better to treat this as a vector problem instead of three scalar equations. The code will become cleaner and follow the equations much closer. Below I will give suggestion on how to update parts of the code, but you have to update other parts of the code accordingly... This requires some refactoring but I think it can be worth it.

lijas · 2026-02-15T11:54:28Z

docs/src/literate-tutorials/dg_shallow_water.jl

+droptol!(M, 1.0e-12)
+M = Diagonal(M)
+# ### Volume integral on each element
+function element_volume_integral!(R, cell, U, dr, param)


Element rotine becomes:

function element_volume_integral!(R, cell, U, param) Re = param.Re fill!(Re, 0.0) cv = param.cv nb = getnbasefunctions(cv) Ferrite.reinit!(cv, cell) dofs = celldofs(cell) Ue = @view U[dofs] # Note: For performance reasons, we don’t recommend using this outside the DG framework; here, the DOFs are guaranteed to be stored contiguously in memory. coords = getcoordinates(cell) for qp in 1:getnquadpoints(cv) dE = getdetJdV(cv, qp) Uval = function_value(cv, qp, Ue) F = F_function(Uval) X = spatial_coordinate(cv, qp, coords) SS = S(Uval, X) for i in 1:nb Φᵢ = shape_value(cv, qp, i) ∇Φᵢ = shape_gradient(cv, qp, i) Re[i] += (F ⊡ ∇Φᵢ+ SS ⋅ Φᵢ) * dE end end Ferrite.assemble!(R, dofs, Re) return nothing end

with

# Define double contraction operator for StaticArrays matrices ⊡(A::StaticMatrix, B::StaticMatrix) = sum(A .* B)

lijas · 2026-02-15T11:54:31Z

docs/src/literate-tutorials/dg_shallow_water.jl

+
+#-
+# ### Mass matrix (block-diagonal per element)
+function assemble_mass_matrix!(M, dh, cv)


Mass matrix routine simplifies to:

function assemble_mass_matrix!(M, dh, cv) nb = getnbasefunctions(cv) Me = zeros(nb, nb) asM = start_assemble(M) for cell in CellIterator(dh) Ferrite.reinit!(cv, cell) dofs = celldofs(cell) fill!(Me, 0.0) for qp in 1:getnquadpoints(cv) dE = getdetJdV(cv, qp) for i in 1:nb Φᵢ = shape_value(cv, qp, i) for j in 1:nb Φⱼ = shape_value(cv, qp, j) m = (Φᵢ ⋅ Φⱼ) * dE Me[i, j] += m end end end assemble!(asM, dofs, inv(Me)) end return M end

lijas · 2026-02-15T11:55:35Z

docs/src/literate-tutorials/dg_shallow_water.jl

+
+# ### Interface integral on each interface
+function interface_integral!(R, ic, U, dr, param)
+    dr_h, dr_qx, dr_qy = dr


function interface_integral!(R, ic, U, param) Ri = param.Ri fill!(Ri, 0.0) iv = param.iv nb = getnbasefunctions(param.cv) Ferrite.reinit!(iv, ic) idofs = interfacedofs(ic) Ui = @view U[idofs] for qp in 1:getnquadpoints(iv) n = getnormal(iv, qp; here = true) dF = getdetJdV(iv, qp) UL = function_value(iv, qp, Ui; here = true) UR = function_value(iv, qp, Ui; here = false) #Ujmp = function_value_jump(iv, qp, U) #c = max(max_speed(UL, n), max_speed(UR, n)) #fhat = 0.5(F(UL)⋅n + F(UR)⋅n) - c*Ujmp fhat = LF(UL, UR, n) for i in 1:nb Φᵢ = 2*shape_value_jump(iv, qp, i) Ri[i] += (Φᵢ ⋅ fhat) * dF end end Ferrite.assemble!(R, idofs, Ri) return nothing end

lijas · 2026-02-15T11:56:34Z

docs/src/literate-tutorials/dg_shallow_water.jl

+batpx(X) = -2 * (X[1] - 5) * 0.5 * exp(-(X[1] - 5)^2 - (X[2] - 5)^2)
+batpy(X) = -2 * (X[2] - 5) * 0.5 * exp(-(X[1] - 5)^2 - (X[2] - 5)^2)
+
+function F(U)


We need F to ouput and 3x2 matrix, which also matches your equations

function F(U) h, qx, qy = U ux = qx / h uy = qy / h p = g * h * h / 2 return @SMatrix [qx qy; qx*ux+p qy*ux; qx*uy qy*uy+p] end

Really nice thank you I will update afaik, I remember now why I didn't do that, I was a bit scared to mix SMatrix with Tensors but since you think it's all right 👍

If we will merge Ferrite-FEM/Tensors.jl#236 in the future then mixing the SMatrix in here also won't be necessary anymore. :)

yolhan83 added 6 commits February 11, 2026 14:52

add a DG SW model exemple

7ed3d63

format runic

544df97

format runic

fc49777

format runic

d2a688d

format pre

51d6d03

fix usings

f3a1b27

yolhan83 added 3 commits February 11, 2026 15:31

forgot a +

6543875

dont print the param

0370c9a

dont export vtk here and misplaced parenthesis

6087d24

termi-official self-requested a review February 11, 2026 16:28

termi-official reviewed Feb 11, 2026

View reviewed changes

yolhan83 added 6 commits February 12, 2026 08:11

apply suggestions

e9cd5ab

apply suggestions

8c6b716

Runic

224441c

Runic

71aaba8

Runic

a6e5440

Add titles back

46552a1

termi-official reviewed Feb 12, 2026

View reviewed changes

adapt to comments

b6299b6

KristofferC reviewed Feb 12, 2026

View reviewed changes

docs/src/tutorials/index.md Outdated Show resolved Hide resolved

adapt to comments

2513dfc

rm from CI

9292d2b

termi-official reviewed Feb 12, 2026

View reviewed changes

apply comments

5c15819

yolhan83 added 3 commits February 12, 2026 21:50

fix inline weirdness

b0ef7c6

fix runic

11e6b2d

fix runic

87cfa1e

lijas reviewed Feb 15, 2026

View reviewed changes

Conversation

yolhan83 commented Feb 11, 2026

Uh oh!

yolhan83 commented Feb 11, 2026

Uh oh!

yolhan83 commented Feb 11, 2026

Uh oh!

termi-official commented Feb 11, 2026

Uh oh!

codecov bot commented Feb 11, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Codecov Report

Uh oh!

termi-official left a comment • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

yolhan83 commented Feb 12, 2026

Uh oh!

Choose a reason for hiding this comment

Uh oh!

yolhan83 Feb 12, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

yolhan83 Feb 12, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

termi-official Feb 12, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

yolhan83 commented Feb 12, 2026

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

yolhan83 commented Feb 12, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

lijas left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

codecov bot commented Feb 11, 2026 •

edited

Loading

termi-official left a comment •

edited

Loading

yolhan83 Feb 12, 2026 •

edited

Loading

yolhan83 Feb 12, 2026 •

edited

Loading

termi-official Feb 12, 2026 •

edited

Loading

yolhan83 commented Feb 12, 2026 •

edited

Loading