L2 error of semidiscrete heat finite elements

import PDETextbook.Base

open MeasureTheory
open scoped BigOperators ENNReal

noncomputable section

namespace PDETextbook.LarssonThomee.Thm10_1

open Filter Topology

local notation "E2" => Euclidean 2

/-- Bounded convex domain in R^2 with smooth (C^infty) boundary Gamma, via a smooth defining function. -/
def IsBddConvexSmoothDomain (Ω : Set E2) : Prop :=
  IsOpen Ω ∧ Ω.Nonempty ∧ Bornology.IsBounded Ω ∧ Convex ℝ Ω ∧
  ∃ ρ : E2 → ℝ, ContDiff ℝ (⊤ : ℕ∞) ρ ∧ Ω = {x | ρ x < 0} ∧
    frontier Ω = {x | ρ x = 0} ∧ ∀ x, ρ x = 0 → fderiv ℝ ρ x ≠ 0

/-- A triangle given by its three vertices; K is their convex hull. -/
abbrev Tri := Fin 3 → E2

def triSet (K : Tri) : Set E2 := convexHull ℝ (Set.range K)

/-- Omega_h, the union of the triangles. -/
def unionSet (T : Finset Tri) : Set E2 := ⋃ K ∈ T, triSet K

/-- Sect. 5.2 triangulation of a smooth convex Omega: nondegenerate triangles in closure Omega,
conforming (two triangles meet in the hull of their common vertices: empty, vertex or edge),
boundary vertices of Omega_h lying on Gamma. -/
def IsTriangulation (Ω : Set E2) (T : Finset Tri) : Prop :=
  T.Nonempty ∧
  (∀ K ∈ T, AffineIndependent ℝ K ∧ triSet K ⊆ closure Ω) ∧
  (∀ K ∈ T, ∀ K' ∈ T, K ≠ K' →
    triSet K ∩ triSet K' = convexHull ℝ (Set.range K ∩ Set.range K')) ∧
  (∀ K ∈ T, ∀ i, K i ∈ frontier (unionSet T) → K i ∈ frontier Ω)

/-- Family {T_h}_{0<h<1}: h = max_K diam K, all angles bounded below independently of h. -/
def IsAdmissibleFamily (Ω : Set E2) (T : ℝ → Finset Tri) : Prop :=
  (∀ h ∈ Set.Ioo (0:ℝ) 1, IsTriangulation Ω (T h) ∧
    (∀ K ∈ T h, Metric.diam (triSet K) ≤ h) ∧ ∃ K ∈ T h, Metric.diam (triSet K) = h) ∧
  ∃ θ₀ : ℝ, 0 < θ₀ ∧ ∀ h ∈ Set.Ioo (0:ℝ) 1, ∀ K ∈ T h, ∀ i j k : Fin 3,
    i ≠ j → j ≠ k → i ≠ k → θ₀ ≤ EuclideanGeometry.angle (K i) (K j) (K k)

/-- S_h: continuous on closure Omega, linear (affine) on each K, vanishing on Gamma and on Omega \ Omega_h. -/
def InSh (Ω : Set E2) (T : Finset Tri) (χ : E2 → ℝ) : Prop :=
  ContinuousOn χ (closure Ω) ∧
  (∀ K ∈ T, ∃ g : E2 →ᵃ[ℝ] ℝ, Set.EqOn χ g (triSet K)) ∧
  ∀ x ∈ closure Ω \ interior (unionSet T), χ x = 0

/-- H^2 norm ||w||_2 from weak derivatives up to order 2; equals ⊤ if w ∉ H^2. -/
def h2Norm (Ω : Set E2) (w : E2 → ℝ) : ℝ≥0∞ :=
  ⨅ (G : Fin 2 → E2 → ℝ) (H : Fin 2 → Fin 2 → E2 → ℝ)
    (_ : HasH1Gradient Ω w G ∧ ∀ i, HasH1Gradient Ω (G i) (H i)),
    (eLpNorm w 2 (domainMeasure Ω) ^ 2 + ∑ i, eLpNorm (G i) 2 (domainMeasure Ω) ^ 2 +
      ∑ i, ∑ j, eLpNorm (H i j) 2 (domainMeasure Ω) ^ 2) ^ (1/2 : ℝ)

/-- w ∈ H^2 with weak Laplacian equal to L a.e. on Omega. -/
def HasWeakLaplacian (Ω : Set E2) (w L : E2 → ℝ) : Prop :=
  ∃ (G : Fin 2 → E2 → ℝ) (H : Fin 2 → Fin 2 → E2 → ℝ),
    HasH1Gradient Ω w G ∧ (∀ i, HasH1Gradient Ω (G i) (H i)) ∧
    ∀ᵐ x ∂(domainMeasure Ω), L x = ∑ i, H i i x

/-- Time derivative in L^2(Omega) at t. -/
def HasL2TimeDeriv (Ω : Set E2) (u : ℝ → E2 → ℝ) (ut : E2 → ℝ) (t : ℝ) : Prop :=
  Tendsto (fun s => eLpNorm (fun x => (s - t)⁻¹ * (u s x - u t x) - ut x) 2 (domainMeasure Ω))
    (𝓝[≠] t) (𝓝 0)

def L2RightContAtZero (Ω : Set E2) (u : ℝ → E2 → ℝ) : Prop :=
  Tendsto (fun s => eLpNorm (fun x => u s x - u 0 x) 2 (domainMeasure Ω)) (𝓝[>] 0) (𝓝 0)

/-- (10.1): u_t - Δu = f in Omega x R_+, u = 0 on Gamma (H^1_0), u(.,0) = v; ut is u_t. -/
def IsSolution10_1 (Ω : Set E2) (f : ℝ → E2 → ℝ) (v : E2 → ℝ) (u ut : ℝ → E2 → ℝ) : Prop :=
  (∀ t > 0, HasL2TimeDeriv Ω u (ut t) t ∧ InH10Closure Ω (u t) ∧
    HasWeakLaplacian Ω (u t) (fun x => ut t x - f t x)) ∧
  (∀ x ∈ Ω, u 0 x = v x) ∧ L2RightContAtZero Ω u

/-- (10.3): (u_{h,t},χ) + a(u_h,χ) = (f,χ) for all χ ∈ S_h, t > 0; u_h(0) = v_h;
a(v,w) = ∫ ∇v·∇w with weak gradients. -/
def IsSemidiscrete10_3 (Ω : Set E2) (T : Finset Tri) (f : ℝ → E2 → ℝ) (vh : E2 → ℝ)
    (uh uht : ℝ → E2 → ℝ) : Prop :=
  (∀ t ≥ 0, InSh Ω T (uh t)) ∧ (∀ x ∈ Ω, uh 0 x = vh x) ∧ L2RightContAtZero Ω uh ∧
  ∀ t > 0, HasL2TimeDeriv Ω uh (uht t) t ∧
    ∀ χ, InSh Ω T χ → ∃ (G Gχ : Fin 2 → E2 → ℝ),
      HasH1Gradient Ω (uh t) G ∧ HasH1Gradient Ω χ Gχ ∧
      (∫ x in Ω, uht t x * χ x) + (∫ x in Ω, ∑ i, G i x * Gχ i x) = ∫ x in Ω, f t x * χ x

/-- Theorem 10.1 (Larsson–Thomée). C is independent of h (may depend on Omega and the family).
Norms are ENNReal-valued so that a missing H^2 regularity gives ⊤ rather than a Lean junk value. -/
theorem thm_10_1 (Ω : Set E2) (hΩ : IsBddConvexSmoothDomain Ω)
    (T : ℝ → Finset Tri) (hT : IsAdmissibleFamily Ω T) :
    ∃ C : ℝ, ∀ h ∈ Set.Ioo (0:ℝ) 1, ∀ (f : ℝ → E2 → ℝ) (v vh : E2 → ℝ)
      (u ut uh uht : ℝ → E2 → ℝ),
      InSh Ω (T h) vh → IsSolution10_1 Ω f v u ut → IsSemidiscrete10_3 Ω (T h) f vh uh uht →
      ∀ t ≥ 0, eLpNorm (fun x => uh t x - u t x) 2 (domainMeasure Ω) ≤
        eLpNorm (fun x => vh x - v x) 2 (domainMeasure Ω) +
        ENNReal.ofReal (C * h ^ 2) * (h2Norm Ω v + ∫⁻ s in Set.Ioc 0 t, h2Norm Ω (ut s)) := by
  sorry

end PDETextbook.LarssonThomee.Thm10_1